Corso di Calcolo Numerico II profAl df.Alessandra d’Al id’Alessio …unina.stidue.net/Calcolo...

23/01/2009

1

Corso di Calcolo Numerico IIf Al d d’Al iprof.Alessandra d’Alessio

Lezioni di Calcolo Numerico e MatlabEd. Liguori- II edizione

Dispense OnLine

CALCOLO NUMERICO A.D'ALESSIO 1

Dispense OnLine

disciplinadisciplina attaatta aa risolvererisolvere problemiproblemi realirealicolcol calcolatorecalcolatore

obiettivoobiettivofornirefornire strumentistrumenti softwaresoftware

perper lala risoluzionerisoluzione effettivaeffettiva didi problemiproblemidelladella scienzascienza ee dalladalla tecnicatecnica


23/01/2009

2

svilupposviluppo didi metodimetodi computazionalicomputazionali ee didi algoritmialgoritmiperper lala risoluzionerisoluzione didi problemiproblemi

neinei piùpiù svariatisvariati campicampi applicativiapplicativi


Fisica Chimica

Medicina

Ingegneria Economia

Scienze della terra

Calcolo Calcolo NumericoNumerico

Computer Computer ScienceScience

CALCOLO SCIENTIFICOCALCOLO SCIENTIFICO


pStrumenti softwaresoftware per la risoluzione effettivala risoluzione effettiva

di problemi realireali

23/01/2009

3

La tecnologia avanzata è tecnologia matematica

senza il calcolo scientifico moltefunzionalità utilizzate quotidianamente

non sarebbero disponibili


Modelli di previsioni del tempoModelli epidemiologiciModelli economiciTAC, RNMTAC, RNMComputer GraphicFotografia digitaleMusica digitale:CD,mp3,ipodFilm e tv digitaliTelefonia mobileVolo automaticoCrittografia:bancomat


Crittografia:bancomat,..Indagini statistiche,exit pollInternetModelli finanziari, di codeAnalisi di rischio:assicurazioni

23/01/2009

4

Matematica e mondo reale


Studio del vento lungo un tratto di costa


23/01/2009

5

Calcolo della temperatura e del degassamento di CO2nel cono del Vesuvio


frattali


23/01/2009

6

calcolo PageRank Google


RisolvereRisolvere unun problemaproblema concon unun CalcolatoreCalcolatore


23/01/2009

7

Problema reale Problema reale PP

Modello Matematico Modello Matematico M(P)M(P)Modello Matematico Modello Matematico M(P)M(P)

Modello Numerico Modello Numerico MN(P)MN(P)


AlgoritmoAlgoritmo SoftwareSoftware

IndividuareEntità note (dati)Entità incognite (soluzioni)gLeggi di dipendenza tra esse


Attenzione a che il modello sia adeguato

23/01/2009

8

di tipo discreto e finitodiscreto e finito

i dati sono numeri realirealile relazioni matematiche sono


m mun numero finitofinito di operazioni aritmetiche

RisoluzioneRisoluzione didi MN(P)MN(P) concon unun

Sequenza Sequenza finitafinita di operazioni che calcolano a partire dai dati di input i dati di output

E d ll’ l i


Esecutore dell’algoritmosistema aritmetico a precisione

finita di un elaboratore

23/01/2009

9

Implementazione dell’algoritmo in

uno specifico linguaggio di programmazione

uno specifico ambiente di elaborazione

p g

influenza di


sistema aritmeticocompilatoresistema operativogestione memoria

Soluzione del problema P

Risultato dell’esecuzione del softwarein uno

specifico ambiente di elaborazione


23/01/2009

10

Problema reale Problema reale P :P :Previsioni Metereologiche

Previsioni del tempo per le successive 72hpartendo dai rilevamenti attuali

(migliaia di osservazioni metereologiche)

V i bili d l l i t t t Variabili da calcolare:pressione,vento,temperatura, umidità, nuvolosità,….


Modello Matematico(semplificato) :sistema di PDE

Modello Numerico: discretizzazione


23/01/2009

11

discretizzazione del dominio: griglie


Sistema lineare di grandi dimensioni e sparsoSistema lineare di grandi dimensioni e sparso

Griglia di 2.2Km : Arco Alpino: 10080000 punti

Sviluppo di software accurato ed efficientegrandi risorse di calcolo ad alte prestazioni

Earth Simulator NEC (Giappone)35 Tflops (1TFlop=1012 flops)


RoadRunner IBM-Los Alamos1 Pflops (1mil. di miliardi al sec.)

35 Tflops (1TFlop=1012 flops)previsioni metereologiche

23/01/2009

12

Esempi di risultati : vento, umidità, precipitazioni


Problema reale Problema reale PP

Modello Matematico Modello Matematico M(P)M(P)

Fonti di erroreFonti di errore

Errori di semplificazioneErrori di semplificazione

Modello Matematico Modello Matematico M(P)M(P)

Modello Numerico Modello Numerico MN(P)MN(P)

Errore di troncamentoErrore di troncamento

E di dE di d ffff


AlgoritmoAlgoritmo SoftwareSoftware

Errore di roundErrore di round--offoff

23/01/2009

13

ERRORE ERRORE DIDI TRONCAMENTO (DISCRETIZZAZIONE)TRONCAMENTO (DISCRETIZZAZIONE)

Differenza tra il risultato vero e quello prodotto dal modello numerico (in R)modello numerico (in R)

Troncare una serieTerminare un processo iterativo prima della convergenza…………..


ERRORE di ROUND OFFERRORE di ROUND OFF

Differenza tra il risultato prodotto da un particolare particolare algoritmo algoritmo in R e quello prodotto dallo stesso algoritmo nelsistema aritmetico a precisione finita di un elaboratore

dipende da criterio di rappresentazione in memoria dei dati gli algoritmi per le operazioni aritmetiche

SISTEMA ARITMETICO DEL CALCOLATORESISTEMA ARITMETICO DEL CALCOLATORE



23/01/2009

14

Risoluzione di un problema col calcolatoreErrori che non si possono evitare

Non si ha la soluzione esatta

Errore computazionale = Errore computazionale = EEt.a.t.a.+errore round+errore round--offoff

VALUTARE L’VALUTARE L’ACCURATEZZA ACCURATEZZA DEL CALCOLODEL CALCOLO


VALUTARE L’VALUTARE L’ACCURATEZZA ACCURATEZZA DEL CALCOLODEL CALCOLO

ERRORE ASSOLUTOEa = |x-y|

ERRORE RELATIVOEr = |x-y|/|x|

y è un’approssimazione di x (valore vero)

se y ed x hanno pp cifre se y ed x(#0) hanno pp cifrese y ed x hanno pp cifredecimalidecimali ugualiEa= |x-y| < 10--pp

y ( ) pp fsignificativesignificative uguali

Er= |x-y|/|x| < 10--p+p+11

L’errore relativo tiene conto dell’ordine di grandezzaInformazione più significativa

Il l è i it


Il valore vero è sempre incognitoL’errore si può solo stimare

23/01/2009

15

Modello matematico

∑∞

=

++++==0

32

...!3!2

1!n

nx xxx

nxe

Modello NumericoModello Numerico

!...

!3!21

!)(

0

32

Nxxxx

nxxS

NN

n

n

N +++++== ∑=


EEtata(N)(N) == |e|exx-- SSNN||

Modello matematico

Modello Numerico

hxfhxf

dxxdf

h

)()(lim)(0

−+=

→

Modello Numerico

hxfhxf

dxxdf )()()( −+

≅


EEtata(h)(h) == |f’(x)|f’(x)-- [f([f(x+hx+h))--f(x)]/f(x)]/h|h|

23/01/2009

16

per N→∞ Eta(N) = ex- SN→0

per h→0 Eta(h) = |f’(x)- [f(x+h)-f(x)]/h| →0

il modello numerico convergeconverge

È interessante misurare lala velocitàvelocità di convergenza


valutazione dell’efficienzadell’efficienza

f(x+h)=f(x)+hf’(x)+(h2/2)f”(ξ) formula di Taylor

f'( ) [f( h) f( )]/h (h/2)f”(ξ)f'(x)=[f(x+h)-f(x)]/h-(h/2)f”(ξ)

Eta(h) = (h/2)f”(ξ) =O(h)

per h→0 Et (h) = →0 come h


per h→0 Eta(h) = →0 come h

Ordine di convergenza 1 LINEARE

23/01/2009

17

f continua fino alla derivata terza formula di Taylor

f(x+h)=f(x)+hf’(x)+h2f’’(x)+h3f”’(ξ)2 6

f(x-h)=f(x)-hf’(x)+h2f’’(x)-h3f”’(η)2 62 6

diverso modello numericof'(x)≅[f(x+h)-f(x-h)]/2h

Eta(h) = (h2/12)|f”’(ξ)+f”’(η)| =O(h2)


ta( ) ( )| (ξ) (η)| ( )

Ordine di convergenza 2 PIU’ VELOCE

modello 1 O(h)dimezzando h l’errore diminuisce della metà

modello 2 O(h2)dimezzando h l’errore diminuisc di 1/4modello 2 O(h ) diminuisce di 1/4

il modello 2 converge più velocementeossia

a parità di h è più accurato


23/01/2009

18

h mod 1 mod 2 er mod1 er mod2

derivata di ex per x=0 f’(0)=e0 =1dimezzando h con i due modelli

h mod.1 mod.2 er.mod1 er.mod25.0000e-01 1.2974e+00 1.0422e+00 2.9744e-01 4.2191e-022.5000e-01 1.1361e+00 1.0104e+00 1.3610e-01 1.0449e-021.2500e-01 1.0652e+00 1.0026e+00 6.5188e-02 2.6062e-036.2500e-02 1.0319e+00 1.0007e+00 3.1911e-02 6.5117e-043.1250e-02 1.0158e+00 1.0002e+00 1.5789e-02 1.6277e-041.5625e-02 1.0079e+00 1.0000e+00 7.8533e-03 4.0691e-057.8125e-03 1.0039e+00 1.0000e+00 3.9164e-03 1.0173e-053.9063e-03 1.0020e+00 1.0000e+00 1.9557e-03 2.5431e-061 9531e-03 1 0010e+00 1 0000e+00 9 7720e-04 6 3578e-07


1.9531e 03 1.0010e+00 1.0000e+00 9.7720e 04 6.3578e 079.7656e-04 1.0005e+00 1.0000e+00 4.8844e-04 1.5895e-074.8828e-04 1.0002e+00 1.0000e+00 2.4418e-04 3.9736e-082.4414e-04 1.0001e+00 1.0000e+00 1.2208e-04 9.9340e-09

rapporto tra gli errori dei due modelli

4.5757e-001 2.4767e-0014.7896e-001 2.4941e-001

4.8953e-001 2.4985e-0014.9478e-001 2.4996e-0014.9739e-001 2.4999e-0014.9870e-001 2.5000e-0014.9935e-001 2.5000e-0014.9967e-001 2.5000e-0014.9984e-001 2.5000e-0014.9992e-001 2.5000e-0014.9996e-001 2.5000e-001

0

0.05

0.1

0.15

0.2

0.25 err. mod 1err. mod 2


1 2 3 4 5 6 7 8 9 10 11 120

Grafico degli errori dei 2 modelli

23/01/2009

19

ERRORE di ROUND OFFERRORE di ROUND OFF


Dipende dal

Criterio di rappresentazioneCriterio di rappresentazioneFLOATING POINT NORMALIZZATO

A PRECISIONE FINITA

Per ogni x∈ℝ, x ≠ 0esponente

x= mm ×β×βee


ββbase

mantissa|m|

23/01/2009

20

x≠0, x∈F x= + . d1d2…dt×βe

L < e < U0

23/01/2009

21

x=0.d1d2…dtdt+1…×βe ,L

23/01/2009

22

F(β,t L,U) (con arr.) si dimostra

|fl|fl(x)(x) ––x)|x)|

23/01/2009

23

Esempio errore di round-off

>> format hex>> 0.3ans =

3fd3333333333333>> 0.1+0.2ans =

3fd3333333333334>> 0.3/0.1 %risultato esatto 3ans =

4007ffffffffffff

CALCOLO NUMERICO-A.d'Alessio 45

4007ffffffffffff>> 3ans =

4008000000000000

Per ogni x∈F esiste un insieme di numeri macchinache è “zero” per x

Per ogni x∈F si definisceεx = ε|x| = epsilon relativo ad xx p

ε εx

0 1 x

CALCOLO NUMERICO A.d'ALESSIO 46

23/01/2009

24

32 bit

s.m. Espon. Mantissa1bit 1bit 8bit8bit 23 bit23 bit

p

Bit implicito : il primo bit della mantissanon è memorizzato è sempre 1

Esponente in memoria: intero ∈[0,255]t l i i 127 (bi )


vero esponente = valore in memoria – 127 (bias)e ∈[-127,128]

Fornisce un meccanismo per proseguire il calcolo in presenza di situazioni di errore

op. non corretta (∞∗0, 0/0,…) (not a number) NaNNaN

Divisione per zero , overflow ±∞ ++infinfpc/+0 =+inf , c/+inf =0

underflow graduale denormalizzazione mantissasubnormal numbers

minr = (1.00..0)2×2-126≈1.2×10-38


maxr = (1.11..1)2×2127≈3.4×1038

23/01/2009

25

Gli esponenti minimo e massimo sono usatiper gestire le situazioni eccezionali

v=NaN e=255 m#0v= +inf e=255 m=0v = + 0 e=0 m=0 v = + 0 e=0 m=0 v=subnormal numbers e=0 m#0

V=numero macchina 0

23/01/2009

26

aritmetica IEEE singola precisione

ε = 1 2-23+1=2-23 ≅ 1.9209×10--772

aritmetica IEEE doppia precisionearitmetica IEEE doppia precisione

ε = 1 2-52+1=2-52 ≅ 2.2204×10--16162

In matlab:eps(‘single’), eps , eps(x)


criterio di arresto di un processo iterativo convergentex1,x2,… → xxn+1=axn+cn formula ricorrente

Quando arrestare il processo iterativo?

| x x | |x | l ti x Criterio di

la successione “converge” numericamente

Quando xn+1 e xn hanno tutte le t cifre significative uguali( uguali a quelle di x )


| xn+1-xn| < ε|xn|= ε relativo a xnCriterio diarresto naturale

23/01/2009

27

calcolo del limite della successione an=an-1+1/2n a0=0lim an = 1

| an+1-an| < ε|an|= ε relativo ad an Criterio di arresto naturale

eseguendo i calcoli in matlab eps ≈ 10-16

n an1 0.000000000000000

2 0.500000000000000

3 0.750000000000000

4 0.875000000000000

n an8 0.992187500000000

9 0.996093750000000

10 0.998046875000000

………


4 0.875000000000000

5 0.937500000000000

6 0.968750000000000

7 0.984375000000000

50 0.999999999999998

51 0.999999999999999

52 1.000000000000000

%esempio criterio di arrestoa(1)=0;a(2)=1/2;n=2;while(abs(a(n)-a(n-1))>eps(a(n))&&n

23/01/2009

28

Criterio di arresto naturale

richiedereil massimo numero

di cifre significativeesatte del sistema

richiesta meno impegnativap gk cifre significative esatte

xn+1e xn con k cifre significative uguali k

23/01/2009

29

TOL è una quantità “piccola”

se TOLX= TOL|xn| va in underflow

|xn+1-xn|< 0 impossibile!

calcolo di TOLXche evita l’underflow

if ⏐xn⏐>minr/TOL thenTOLX =⏐xn⏐×TOL

else


che evita l underflow elseTOLX = minr

endif

Inevitabile introduzione diErrore di round off nei dati di input può

essere amplificato dal problemaCONDIZIONAMENTOCONDIZIONAMENTO

E di d ff d ll i i ò ss Errore di round off delle operazioni può essere amplificato dall’algoritmo

STABILITA’STABILITA’

tali errori possono propagarsi ed amplificarsi

stimare e controllare la crescita degli errori

CALCOLO NUMERICO- A.d'ALESSIO 58

g

23/01/2009

30

(P’,fl(d)) “prossimo” a (P,d)

x prossimo x’P

bencondizionatobencondizionato

x lontano da x’P

malcondizionatomalcondizionato

Dipende solo dal problemaP l di i bi d ll

CALCOLO NUMERICO-A.d'ALESSIO 59

P malcondizionato cambiare modello

problema P

x + 3y= 4x + 3.001 y = 4.001 soluzione x1=(1,1) in R

problema P’

x + 3y= 4x + 3.002 y = 4.001 soluzione x2=(2.5,0.5) in R


perturbazione Completamentediversa!

23/01/2009

31

x1

Errore dati

x2


Errore sol. = dist(x1,x2) >> errore dati

Sistema malcondizionato

Stima quantitativa del condizionamento

Errore relativosoluzione

= K x Errore relativo nei dati

Errore relativo soluzione 10q x 10-t ≅ 10q-t

Errore relativosoluzione ≅ 10q x 10

-t

La soluzione ha al più t q cifre corrette

Indice di condizionamento


t ≅ q soluzione inaffidabile

La soluzione ha al più t-q cifre corrette

23/01/2009

32

AlgoritmoAlgoritmo StabileStabilel’errore di round off nelle operazioni dell’algoritmo

n nn n sisi mplificmplific n l isult t

AlgoritmoAlgoritmo InstabileInstabilel’errore di round off nelle operazioni dell’algoritmo

si mplific “ n m m nt ”mplific “ n m m nt ” n l isult t

nonnon sisi amplificaamplifica nel risultatoil risultato calcolato è la soluzione “esatta”“esatta” didi P’P’


si amplifica enormemente”amplifica enormemente” nel risultatoche differisce sostanzialmente da quello di P’

x2-12.4x +.494=0 sol: x=12.36 y=.03996 in un sistema F. con t=4

x=(12.4+Δ1/2)/2 = 12.36y=(12.4-Δ1/2)/2 = 0.1000 errato!

C ll i t t fiCancellazione catastrofica

usando le formule alternativex=(12.4+Δ1/2)/2 = 12.36y=c/ax =.494/12.36=0.0401 esatto !

algoritmi equivalenti danno soluzioni diverse


algoritmi equivalenti danno soluzioni diversea causa dell’amplificazione

dell’errore di round off

23/01/2009

33

P, d P, d x x

P’, P’, flfl(d )(d )yy

Algoritmo stabile stabile

ProblemaProblemaverovero

Problema nelProblema nelcomputercomputer

z soluzione calcolata di P’ |z-y| “piccolo”“piccolo”

PPbencondizionatobencondizionato

PPmalcondizionatomalcondizionato


x,y “vicini”|x-z| “piccolo”

x,y “lontani”|x-z| “grande”

Un problema può essere malcondizionatoper certi dati e non per altri

Un problema bencondizionato può dare risultati non corretti se risolto con un algoritmo instabile

Se un problema è bencondizionato esisteun algoritmo stabile per risolverlo

Un semplice calcolo può darerisultati errati seeffettuato con un algoritmo instabileA h i di d ff i li i t ti i


Anche errori di round off piccoli, ripetuti inalgoritmi lunghi e complessi possono avereeffetti catastrofici

23/01/2009

34

Fallimento di un Missile PATRIOT

25 02 1991 guerra del Golfo Arabia Saudita Un missile Patriot25.02.1991,guerra del Golfo,Arabia Saudita. Un missile Patriotfallisce l'intercettazione di un missile iracheno SCUD

conseguenze

morte di 28 soldati e oltre 100 feriti

Cause


auseRapporto del General Accaunting Office:

Patriot Missile Defense:Softwareproblem led to System Failure at Dharan

calcolo inaccurato del tempo di caricamentoa causa dell’aritmetica del computer

il tempo, calcolato in decimi di secondo nelcomputer del PATRIOT è

f i ditrasformato in secondi

t in sec.= tx1/10 0.1 in base 2 ha infinite cifre

sistema di calcolo a virgola fissa a 24 bit

troncamento


errore .95 x 10-7

23/01/2009

35

la batteria è stata in funzione 100 oreprima del lancio

al momento dell’intercettazioneerrore = 0.34sec.

Velocità SCUD 1676 m/secVelocità SCUD 1676 m/sec

in 0.34sec percorre 569.84 metri

SCUD fuori traccia del patriot di


f pcirca mezzo chilometro

non è intercettato

Un programma si deve basare suun buon algoritmo

risolvere un problema con“ minima spesa” = efficienza“massima accuratezza”= stabilitàm m

numero di operazioninumero di variabili utilizzate

tempo di calcolospazio di memoria

70CALCOLO NUMERICO-A.d'ALESSIO

spazio di memoria

23/01/2009

36

Costo computazionale di un algoritmo

numero di operazioni richiestecorrelato ad un

i di di l i à di indicatore di velocità di un computerflop = 1 operazione f.p.

velocità = numero di flops al secondo

(mega) Mflops=106 flops(giga) Gflops=109 flops

71CALCOLO NUMERICO-A.d'ALESSIO

(g g ) f p f p(tera) Tflops=1012 flops(peta) Pflops=1015 flops

i computer più veloci (supercomputer)

RoadRunner IBM-Los Alamos1 Pflops (1mil. di miliardi al sec.)

BlueGene/L IBM200 fl200 Tflops

Earth Simulator NEC (Giappone)35 Tflops

BlueFire National Center Atmospheric Research75 Tflops

mutamenti climatici


35 Tflopsprevisioni metereologiche

Lagrange HP Italia (CILEA)22 Tflops

23/01/2009

37

n dimensione del problema

efficienzaefficienza= si valuta al crescere di n mediante

T(n)T(n) = complessità di temponumero di flops più significativep p g

S(n)S(n) = complessità di spazionumero di variabili utilizzate


sistema aritmeticoorganizzazione della memoria compilatoresistema operativosistema operativo

prestazioni dialgoritmo software

il software non è la banale traduzione di un algoritmo


f gin un linguaggio di programmazione

23/01/2009

38

performanceperformance

tempo di esecuzione =tempo di CPU

elapsed timetempo di CPU

impiegato dall’unità centrale per

eseguire il programma

non conta i tempi di attesa

tempo tra il momento in cui il Software

è mandato in esecuzione

ed il suo completamento


pdi input/output

comp tam nto

AffidabilitàAffidabilitàEfficienzaEfficienzaRobustezzaRobustezza

AffidabilitàAffidabilità

Accuratezza Accuratezza Correttezza codiceCorrettezza codiceMisura dell’effettivo Congruenza tra specifiche

di d ff d l i lt ti


errore di round off del programma e risultatiintrodotto nel risultato

23/01/2009

39

Tecniche per analisi di correttezza

Statiche (senza esecuzione)corrispondenza non corretta tra

parametri formali ed attuali utilizzo di variabili non inizializzate, etc

Dinamiche (con esecuzione)rilevabili dal software di base (compilatore..)

radice quadrata negativaindice di array fuori del range fissato,etc

Testing strutturalescelta di dati in modo da soddisfare un


criterio di “copertura”ogni istruzione eseguita almeno una volta

EfficienzaEfficienza

tempo di elaborazione edccup zi n

RobustezzaRobustezza

controllooccupazionedi memoria su

un dato calcolatore

fortemente dipendentedall’ambiente di calcolo

controllodi overflow, underflow,…sui dati di input

gestione dierrori e diagnostica


dall ambiente di calcolo

23/01/2009

40

interfaccia utente-software

Istruzioni d’usoInformazioni sull’organizzazione

internaT f i di iTrasferimento di esperienze

DocumentazioneDocumentazione

InternaInterna Esterna


Linee di commenti manuale d’usoesplicativi di sequenze help on line

Scopobreve descrizione dei problemi risolubili, algoritmo usato, raccomandazioni sull’uso Complessità e

Specificheintestazione della procedura

Lista parametriinput/output

tipo,dimensione e descrizione

accuratezzaeventuali informazioni su

complessità di tempo e spazioaccuratezza del risultato

Esempio d’usoesempio di un semplice

h


Diagnosticadescrizione delle

situazioni di errore previste

programma chiamantecon elenco dei dati di input

e dei risultati

Corso di Calcolo Numerico II profAl df.Alessandra d’Al id’Alessio …unina.stidue.net/Calcolo...

Documents

Transcript of Corso di Calcolo Numerico II profAl df.Alessandra d’Al id’Alessio …unina.stidue.net/Calcolo...