Metodi matematici dellastronomia Che cosa sono i Metodi matematici dellastronomia? In che cosa sono...

Metodi matematici dell’astronomia

Che cosa sono i Metodi matematici dell’astronomia?

In che cosa sono differenti da quelli della fisica?

L’astronomia è una scienza che si basa sull’analisi quantitativa deifenomeni su base fisico-matematica, ma che mutua alcune caratteristichedalle Scienze naturali, quali l’accurata raccolta e analisi di “reperti” osservativi (tassonomia). Ciò implica la necessità di utilizzare metodistatistici di indagine. Inoltre l’astronomia-astrofisica moderna è carat-terizzata dal dover interpretare quantitativamente fenomeni complessi (sovrapporsi di processi fisici interconnessi e su un ampio intervallo di scale spazio-temporali) il che implica l’utilizzo di metodi matematici e, principalmente, numerici specifici.


L’astronomia è l’esempio massimo di scienza osservativa, che è diverso da sperimentale. Ciò significa che ci si limita a raccogliere dati (output) provenienti da sorgenti lontane.

Se si osserva un fenomeno particolare (per es. l’esplosione di una SN):

La SN 1987a

prima dopo


• se ne possono raccogliere i dati e ipotizzare un’interpretazione che non può, ovviamente, essere confermata in laboratorio:

• si può solo sperare che ci siano un numero sufficiente di altre esplosioni osservabili che forniscano i dati utili alle ipotesi interpretative teoriche. La possibilità di azione del ricercatore è solo quella di potenziare i mezzi di raccolta e analisi dei dati osservativi.


Ammesso di avere statistica sufficiente che cosa si può imparare come fisica?In astronomia i processi fisici sconosciuti alla base del flusso radiativo osservato sono convoluti attraverso la struttura della sorgente, il mezzo interstellare, l’atmosfera terrestre e infine il ricettore e il rivelatore.

Lo schema è quindi

input Struttura, ISM, ecc. output

Quindi si tratta di risolvere un problema di inversione.

A parte i problemi tecnici (matematici) implicati, la difficoltà intrinseca dell’astronomia di ottenere buoni dati osservativi rende il processo di falsificazione delle ipotesi (al fine dell’ottenimento dell’unica interpretazione esatta) difficile. Infatti è possibile che modelli fisici anche significativamente diversi portino a “osservabili” indistinguibili entro


l’errore osservativo. In altro linguaggio: il test “pratico” cui si sottopongono le varie teorie è troppo debole per selezionarne una sola (o poche).

Situazione “auspicabile”

Situazione “reale”


In un esperimento ideale la risposta dello strumento allo stimolo (sconosciuto) è univoca, nella realtà tale risposta strumentale è “sbrodolata” su un insieme (tramite la point-spread function)

Stimolo f(y)

risposta strumentale g(x)

point spread function k(x,y)convoluzione

F(y) g(x)

deconvoluzione

Un’approssimazione F(y) di f(y) si ottiene per deconvoluzione della risposta strumentale, ma quanto vale ?)y(f)y(F


g(’) f()


Il problema della deconvoluzione (o inversione dei dati).

caso k(x,y)=1 x

dyyfxg0

)()(

ha soluzione formale )(')( xgxf

Semplice, ma dà problemi; infatti supponiamo che xx e)x('g)x(fe)x(g 1

Se al dato g(x) è sovrapposto rumore di frequenza si ha:

xcose)x(fxsen)x(g)x(g xoo

Ne consegue che le 2 soluzioni f(x) e fo(x) differiscono

xcosmaxffmax)x(f)x(fx

ox

o

x

dy)y(f)y,x(k)x(g0


che significa che l’errore sulla soluzione ha ampiezza

Il problema è quindi instabile a variazioni di alta frequenza in g(x), infatti la soluzione fo è tale che nel suo errore relativo

d

xcosd

)xxsenx(cosf

df

df

ff

df

o

oo

oo

o 1

le variazioni relative d/ e d/ sono amplificate dal fattore /fo.

Metodi matematici dell’astronomiainstabilità

=0.8, =0.04,=20

output input

instabilità


E quindi? Invece di deconvolvere (procedura “all’indietro”)si può effettuare il cosiddetto “model fitting” (procedura “in avanti”)cioè fare un’ipotesi su f(y), convolvere e valutare il risultato g(x).

Supponiamo di avere un’ipotesi f1(y) e di perturbarla di modo che f2(y)=f1(y)+acosy.

Il relativo osservabile sarà

xsena

)x(gydycosa)x(g)x(g o

x

oo

10

12

per cui

a

gxsena

gycosaf oo

che corrisponde al problema visto prima. Una perturbazione in f di grande ampiezza a può essere smorzata dalla sua alta frequenza ω.


Conclusione: 1) la deconvoluzione è instabile; 2) il model fitting (convoluzione) dà risultati illusori.

La via giusta è affinare la tecnica di deconvoluzione (regolarizzando),utilizzando, se possibile, vari kernels, cioè un insieme di dati osservativiper esempio flussi in bande diverse dello spettro e.m..


Problemi tipici dell’astronomia

)(2

tP

esenEE Equazione di Keplero

),...,2,1(,||132

2

Nimm

Gdt

dm

N

ijj

jiji

jiii

rrrr

r Equazioni del moto di N corpi autogravitanti

sistema di N eq. diff. vett. ordinarie del 20 ordine ordine del sistema = 6N

),...,2,1(

)0(

)0(

0

..

0

Ni

ii

ii

rr

rr

servono 6N condizioni iniziali su posizioni e velocità


2

32

4

2

2

4

1

44

3

3

1

4

1

rdr

dLTr

)r(L

dr

dT

Tc

kTm

ppp

rdr

dM

r

)r(MG

dr

dp

Hrg

R

R

TRT

pRp

p

M

)(

)(

0)0(

0)0(

Sistema di 4 eq. diff. + 1 EOS con incognite: (r). P(r), M(r), T(r), L(r)

mentre: = (,T;Xi), = (,T;Xi), = (,T;Xi) sono funzioni note.

4 condizioni al contorno

Equazioni della struttura stellare


)t(P

esenEE

2 Si vuole ottenere E(t); è necessario quindi risolvere un’equazione trascendente:

02

)t(P

esenEE)t;E(f

dove: E (incognita) è l’anomalia eccentrica, e<1 è l’eccentricità della traiettoria (ellittica), t (il tempo) è un parametro tale che 0 t- P dove è l’istante di passaggio al perielio e P il periodo dell’orbita.L’equazione non ha soluzione esplicita, per cui bisogna cercarne un’approssimazione numerica. Nell’intervallo 0<t-<P (cioè 0<E<2) la soluzione esiste unica, infatti la funzione è continua e

02

0 )t(P

)t;(f 0

222 )t(

P)t;(f

01 Ecose)t;E('f per il teor. dell’esist. degli zeri c’è 1 solaradice nell’intervallo ]0,P[


• Dati 2 numeri a e b è logico confrontarli:

Come ottenere e stimare l’accuratezza nei calcoli numerici

ba

Il significato del confronto, quindi, dipende dal contesto.

Se a e b rappresentano le altezze di 2 persone e b-a=50cmha senso dire che a<<b, ma se parliamo di due montagne no!Lo stesso si può dire per ab e ab.

))x(g(O)x(f per ax significa che )x(g

()flim

ax•


Errori assoluti e relativi

Data un’approssimazione ã di a: • l’errore assoluto è a= ã - a • l’errore relativo è r= (ã – a)/a (se a0)

La conoscenza esatta dell’errore è, ovviamente, di solito impossibile, per cui si cerca di averne una stima, o meglio una limitazione superiore 0 tale che |a| .

Con la notazione a= ã si intende |a |=|ã - a| . In molti casi con tale notazione ha il senso di deviazione standard o altra misura statistica di errore.

Indichiamo d’ora in poi con ã un’approssimazione di a.


Sorgenti d’errore

I risultati numerici sono influenzati da varie sorgenti d’errore: alcune possono essere ridotte o eliminate, altre no.

1) nei dati di input;2) dovuti a semplificazioni nel modello matematico;3) di arrotondamento durante il calcolo (se la macchina gestisce fino a s cifre, un prodotto, che avrebbe 2s o 2s-1 cifre, viene troncato a s);4) di troncamento: sono quelli che nascono dal “taglio” di un’espressione, tipo: invece di oppure dall’aprossimare un funzione non lineare con una lineare oppure dall’approssimare una derivata con un rapporto finito (errore di discretizzazione);5) “umani”: imprevedibili.

n

iin as

1

1i

ias

Possibili errori:


Arrotondamento (round-off)

a=0.08 ha 2 decimali e 1 cifrab=15.6 ha 1 decimale e 3 cifre

• Siano a=0.235 e ã=0.231; ã ha 2 cifre esatte, e ciò corrisponde a unerrore a=ã - a=-0.004=-0.4×10-2. In generale, se |a|0.5 ×10-t, l’approssimazione ha t decimali corretti (qui 2 cifre significative).

• Se a=0.001 e ã=0.002, =0.001=0.1×10-2 2 decimali corretti (e nessuna cifra significativa).

• Il numero a=0.0654±0.0003 ha 3 decimali corretti e 2 cifre significative

Il numero di cifre di un numero è quello che si ottiene escludendo gli zeri all’inizio: a=0.0078 ha 2 cifre a=7.8×10-3

Indichiamo d’ora in poi con ã un’approssimazione di a.


Ci sono 2 modi di arrotondare un numero a un dato numero (t) di decimali:

1) tagliare al t-esimo decimale: a=58.9978345 ã=58.997 (t=3)

2) arrotondare (realmente) il numero in modo da lasciare il t-esimo decimale inalterato se la parte del numero che resta alla sua destra è < 0.5 ×10-t, aumentandolo di 1 altrimenti.

Naturalmente l’arrotondamento causa un errore:per esempio a=0.2660.003 ha 2 decimali corretti (e 2 cifre significative)arrotondandolo a 2 decimali si ha ã=0.27 il cui secondo decimale non ècorretto. In questo caso, il taglio al secondo decimale porta a una miglioreprecisione (2 dec. corretti).


Propagazione d’errore

D. Se il limite d’errore (assoluto) di a1>0 è 1 e quello di a2>0 è 2 , qual è il limite d’errore di a1- a2 ?

R. Poichè a1 = ã1 ± 1 e a2 = ã2 ± 2, si ha ã1 - 1 - (ã2 + 2) a1- a2 ã1 + 1 – (ã2 - 2), cioè ã1 - ã2 - (1 + 2) a1- a2 ã1 – ã2 + (1 +2), quindi a1 -a2 = ã1 - ã2 ± (1 +2).

Analogamente si vede che a1+ a2 = ã1 + ã2 ± (1 +2). Quindi: nelle operazioni di addizione e sottrazione i limiti di errore assoluto si sommano.


Nelle operazione di moltiplicazione e divisione, invece, si sommano ilimiti d’errore relativo (approssimativamente).

Si verifica infatti (per esercizio) che:

(ã1ã2 - a1a2)/ a1a2 = r1 +r2 + r1r2 r1 +r2, se |r1|<<1 e |r2|<<1

e

(ã1/ã2 – a1/a2)/ (a1/a2) = (1+ r1)/(1+ r2 ) – 1 = (r1- r2 )/(1+ r2 )r1 - r2 , se |r1|<<1 e |r2|<<1

Perciò nella moltiplicazione e divisione i limiti d’errore relativo si sommano e sottraggono, rispettivamente.


Cancellazione dei termini

Una causa comune di scarsa accuratezza nel calcolo è la sottrazione di2 numeri molto simili (per cui la differenza è << dei numeri stessi).Questo problema è detto della cancellazione dei termini (term cancellation)

Infatti, dati 2 numeri x1 e x2 affetti da errori x1 e x2, e ponendo y= x1- x2

21

2121 xx

xx

y

yxxy

L’errore relativo in y può quindiessere molto grande

Esempio: x1=6.3456±½·10-4 e x2 =6.3455±½·10-4 y=0.0001 ±10-4,che equivale a |y/y| 10-4/ 10-4 (stima d’errore relativo del 100%).


Bisogna quindi cercare di evitare differenze di numeri simili, riscrivendo,se possibile, le formule. Per esempio, sia da calcolare Essa si può riscrivere come:

xxxx

xx

xx

xxxxxx

x,xx

evitando così la cancellazione dei termini.

Dovendo calcolare f(x+)-f(x) un’altra possibilità è quella di sviluppare f(x) in serie di Taylor:

...δ)x(''ff ’(x)δf(x)δ)f(x 2

2

1


La formula generale per la propagazione dell’errore

Sia data una funzione y(x) e si voglia valutare )x~(y dove x~ è un’approx.di x. Ci si chiede qual è l’errore )x(y)x~(yy

Una via naturale è quella di calcolare il differenziale (approssimato):

x)x~('ydyy

In generale, se y=y(x1,x2,...,xn) si può stimare

n

ii

i

x)x~(x

yy

1da cui:

n

ii

i

x)x~(x

yy

1Che è la formula generale di propagazione dell’errore.


Derivata numerica

Conoscendo una funzione f(x) in 3 punti (xi-1,xi,xi+1), si possono costruire3 rette: a) passante per (xi-1,fi-1) e (xi,fi); b) per (xi,fi) e (xi+1,fi+1); c) per(xi-1,fi-1) e (xi+1,fi+1).

Una ragionevole approssimazione della derivata, f‘ (x), di f(x) all’interno dell’intervallo [xi-1,xi[ può quindi essere data dal coefficiente angolare della retta a) passante per tali punti, così come da quello della retta c) sesi vuole approssimare f‘ (x) in ]xi,xi+1].

E per f‘ (xi)? E’ possibile usare l’approssimazione a) (sinistra), c)(destra) ma anche b) (centrale). Si verifica che, se xi-1,xi,xi+1 sono spaziati con passo h, allora l’approssimazione centrale è del 20 ordine in h mentre le altre due sono del 10 ordine. Si ha infatti, esprimendo f(x+h) e f(x-h) con la f. di Taylor


Derivata numerica

).h(Oh

)hx(f)hx(f)x('f

);h(Oh

)x(f)hx(f)x('f

);h(Oh

)hx(f)x(f)x('f

c

d

s

2

2

di punto di partenza x, esplicitando le derivate prime e semisommandole:


Integrazione numerica

Supponiamo di dover calcolare b

a

dx)x(fI dove f(x) è una funzione

continua nell’intervallo b,a

Se non si riesce a ottenerne la primitiva

F'(x)=f(x), per cui F(b)-F(a)=I, l’unica alternativa è di approssimareI tramite una valutazione numerica. I

~

Le formule di approssimazione si chiamano “di quadratura” perchè si tratta di valutare l’area sottesa dalla funzione, cioè di quadrarla.

Tutti i metodi per ottenere un’approssimazione di I si basanosul fatto che l’integrale definito è il limite della somma finita sn, cioè:

n

iii

nn

nx)x(flimslimI

1

oppure è una funzione campionata per

punti: {f(xi}, i=1,2,...,n}.


dove n è il numero di intervalli ix in cui si è suddiviso [a,b], entro ognuno dei quali si è scelto xi.

Ci limiteremo a studiare 4 metodi:1) il metodo trapezoidale,2) il metodo rettangolare,3) il metodo di Simpson,4) il metodo Montecarlo.

Il metodo trapezoidale

E’ forse il più intuitivo; si basa su:i) considerare una suddivisione di [a,b] in n intervalli, il generico dei quali ha estremi xi e xi+1 e ampiezza hi , definendo quindi una griglia di n+1 punti (mesh-points): x0=a , xi+1=xi+ hi , i=0,1,...,n-1 (naturalmente xn=xn-1+hi=b);


ii) valutare la funzione sui punti griglia, ottenendo l’insieme di n+1 valutazioni {f(xi)fi, i=0,...,n}, di modo da avere il campionamento della funzione dato dalle coppie (xi, fi);iii) stimare I come somma delle aree dei trapezi di base fi e fi+1 e altezza hi, cioè:

1

0

1

2

n

i

iiib

a

t

h)ff(dx)x(f

~I~

Tale approssimazione equivale ad aver sostituito (dentro l’integrale)ad f la sua approssimazione costituita dalla spezzata passante per ipunti (xi,fi), cioè:

)xx()xx()xx(xx

fff)x(f

~ii

n

ii

ii

iii

1

1

0 1

1

dove (x) è la funzione di Heaviside:

01

00

x,

x,)x(

Verificare per esercizio che .I~

dx)x(f~b

a

t


Quando è possibile, è comodo usare un passo fisso hi=h=(b-a)/n, di modo che i punti xi son definiti come x0=a, xi= x0+ih (i=1,2,...,n),il che semplifica l’espressione delle formule di quadratura (e ne riduce lacomplessità computazionale).

Il metodo rettangolare

Si basa sull’approssimazione di f(x) con una funzione a gradini tale che fi=f(xi+1/2) nel generico [xi,xi+1], dove xi+1/2 è il punto di mezzo dell’intervallo. L’approssimazione di I che ne risulta è:

1

021

n

ii/ic,r hfI

~

(form. rett. centrata).Se il passo è fisso, hi=h=(b-a)/n, la formula diviene

1

021

n

i/ic,r fhI

~


Il metodo rettangolare Il metodo trapezoidale

E’ intuitivamente chiaro (vedi figura in basso a sin.) che prendere la funzione nel punto di mezzo xi+1/2 invece che nel punto xi (form. rettangolare sinistra) o xi+1 (form. rett. destra) dà un’approx migliore

1

0

n

iiis,r hfI

~f. rett. des.

1

01

n

iiid,r hfI

~f. rett. sin.


L’errore di troncamento nella formula trapezoidale e rettangolare

Data una funzione f(x) che assume valori fi su una griglia di n+1 punti xi (i=0,1,...,n) esiste un solo polinomio, pn(x), di grado n che passa per i punti (xi ,fi), la cui espressione è:

n

in

ikk

ki

n

ikk

k

in

)xx(

)xx(

f)x(p0

0

0 (formula d’interpolazione diLagrange; verificare per es. che pn(xi)= fi, i=0,1,...,n.)

),x(Lf)x(p i

n

iin

0

che si scrive compattamente come

avendo posto

n

ikk

ki

n

ikk

k

i

)xx(

)xx(

)x(L

0

0

(funzione moltiplicatirce di Lagrange).


L’espressione di Lagrange si ottiene dall’espressione polinomiale pn(x)=c0+c1(x-x0)+c2(x-x0)(x-x1)+···+cn (x-x0)(x-x1) ···(x-xn),determinando (ricorsivamente) le costanti ck tramite l’imposizione che pn(x) assuma i valori fk in xk, k=0,1,...,n.

Se f(x) è continua insieme, almeno, alla (n+1)-esima derivata nell’interv.Int(x,x0,x1,...,xn) (che è per def. il minimo interv. contenente x,x0,x1,...,xn),si può verificare che il resto (errore) della formula di Lagrange è:

),xx()xx)(xx()!n(

)(f)x(p)x(f)x(R n

)n(

nn

10

1

1

dove è un punto (incognito) in int(x,x0,x1,...,xn). Si noti che l’espr. dataè simile al resto di uno sviluppo di Taylor (uguale, se x0=x1=···=xn).

Possiamo ora valutare l’err. di troncamento della formula trapezoidale nelgenerico intervallo [xi,xi+1] semplicemente integrando R1(x) su tale int..


11 11

111 2

i

i

i

i

i

i

i

i

x

x

ii

x

x

x

x

x

x

i dx)xx)(xx()(''f

)x(p)x(f)x(R

in cui [xi,xi+1] dipende da x. Poichè (x-xi)(x-xi+1)0 in [xi,xi+1], si può applicare il teor. della media:

33

0

3

122

1

3

1

21

2

1

ii

ii

ii

i h)(''f

h)(''f

dt)t(th)(''f

Ponendo x=xi+hit, dove hi=xi+1-xi, si ha

1

112

i

i

x

x

iiiiii

i ]x,x[,dx)xx)(xx()(''f

Ricordando che l’err. della somma è la somma degli errori, l’err. di tronc. globale è dato da

,h)(''fn

iii

n

iig,t

1

0

31

0 12

1 )(''f

h)ab()(''nf

h

n

)(''fn

h)(''f

h n

i

in

ii

n

iig,t

12121212

231

0

31

0

31

0

che, se i punti xi sono equispaziati xi+1-xi=h=(b-a)/n,diventa (applicando il teor. del val.medio)


Es. Valutare l’errore di troncamento della formula rettangolare centrata.

Si può sviluppare f(x) con la formula di Taylor attorno al punto xi+1/2 nell’intervallo [xi,xi+1]

2

2)xx(

)(''f)xx)(x('f)x(f)x(f 1/2i1/2i1/2i1/2i

per cui il resto dell’interpol. di f(x) in [xi,xi+1] rispetto al suo valore in xi+1/2 è

.)xx()(''f

)xx)(x('f)x(f)x(f)x(R 1/2i1/2i1/2i1/2i2

2

Perciò l’err. di troncamento della form. rettangolare nell’interv. [xi,xi+1] è

1i

i

i

i

i

i

x

x

1/2i

x

x

1/2i1/2i

x

x1/2it dx)xx(

)(''fdx)xx)(x('fdx)]x(f)x(f[ 2

2

11

Quindi la formula trapezoidale è localmente del 30 ordine e globalmentedel 20.


Facendo la sostituzione x=xi+ht=xi+1/2 –h/2+ht si verifica che il 10 dei 2 integrali è nullo, mentre il 20, usando il t. della media integrale, porta a

,h)(''f

h)(''f

h)(''f

dt)t(h)(''f iiii

t333

0

23

2412

1

22

1

4

1

3

1

22

1

2

1

).(''fh)ab(

g,t 24

2per cui l’err. di troncamento globale è


Metodo di Simpson

Risponde all’esigenza di avere una formula di quadratura di ordineelevato e semplice da utilizzare al contempo.In pratica, si tratta di trovare una formula che, dati 3 punti xi-1,xi,xi+1 (xi=xi-1+h) sia esatta per polinomi di più alto grado possibile.

1i

1i

x

x

1ii1i cfbfafhdx)x(fScritta la formula si tratta di trovare i coefficientia,b,c

Per comodità poniamo i=0 e x0=0 e sviluppiamo f(x) in serie di McLaurin, per cui, per |x|h, si ha

),h(Ox)('''f!

x)(''fx)('f)(f)x(f 432 03

10

2

100

che integrata in [x-1,x1] dà


,h)h(O)('''f!

)(''f)('f)(hfdx)x(fxxx

h

h

h

h

h

h

h

h

2

4

03

13

02

12

002 4

432

in cui i termini di potenza pari sono nulli, per simmetria, e l’ultimo termine è O(h5).

Sviluppando in serie di Taylor in “avanti” e “indietro” attorno a x, si ha

f(x+h)=f(x)+f '(x)h+(1/2)f '' (x)h2+(1/3!)f '''(x)h3+O(h4), f(x-h)=f(x)-f ‘(x)h+(1/2)f ''(x)h2-(1/3!)f '''(x)h3+O((-h)4),

dove in O(h4) e O((-h)4) appare la derivata quarta, per cui l’espressione è esatta per polinomi fino al 30 grado. Si può quindi esplicitare f ''(x) sommando m. a m. e dividendo per h2:

).h(O)x(f)hx(f)hx(fh

)x(''f 22

21


Valutando l’espressione precedente in x=0 e inserendola nell’integrale si ha

),h(Offfh)h(Offfh

hfdx)x(fh

h

5110

50110 3

1

3

1

3

42

32

),h(Offfh

dx)x(f iii

x

x

i

i

511 4

3

1

1

Da cui la formula generale, nell’intervallo [xi-1,xi+1]

esatta per polinomi fino al 30 grado, che è la formula di Simpson.Si verifica che l’errore di troncamento locale è

11

5

90 ii)iv( x,x),(f

h

per cui quello globale è b,a),(f)ab(

h)(f

nh)(f

hi

)iv()iv(

pariii

)iv(g,t

18029090

455

(il fattore 2 a dividere viene dal fatto che h=(b-a)/n, per cui ci sono n/2


intervalli [xi-1,xi+1] in [a,b]). Il metodo è quindi di ordine 4.

Metodo Montecarlo

La base del metodo Montecarlo è l’approssimazione di un integrale multidimensionale nel modo seguente

Vn

VfVVfd)(f fV

rr 3

Dove <f> è la deviazione standard della media, <f>, di f (r) valutata prendendo n punti distribuiti a caso nel dominio di integrazione, di volume V. Tale deviazione standard è uguale alla deviazione standard di f , cioè ((f-f)2)1/2, diviso la radice del numero, n, dei punti usati per il campionamento:

n

ff

n

ff

nf

222


In una dimensione, la valutazione Montecarlo dell’int. di f(x) sull’interv.[a,b] si effettua distribuendo a caso n punti (xi,yi) nel dominio rettangolareA di base [a,b] e altezza max f(x), stimando IA come

)x(f)ab(n

)x(k)A(area

n

)x(kdx)x(fI max

n

iib

a

n

ii

11

dove il rapporto che moltiplica l’area di A è quello fra il numero n* di punti del campionamento che stanno entro l’area I e il numero totale dipunti. Definendo la funzione k(x) come

n

ii* n)x(kn

1

0

),x(fy,

),x(fy,)x(k

ii

iii 1

0risulta

I punti (xi,yi) sono indicati;nel caso in fig. xmax=b


Poichè l’ordine del metodo è h1/2 (la stima d’errore converge a zero conn-1/2) questo metodo è conveniente solo in più dimensioni, essendo facile da implementare.

RADICI DI EQUAZIONI NON LINEARI

Per quanto riguarda le equazioni algebriche, ricordiamo il teorema diAbel-Ruffini: non è possibile esprimere sotto forma di radicali le soluzionidi equazioni algebriche superiori al 40 grado. Il problema di trovare approssimazioni numeriche alla soluzione di un’eq. del tipo f(x)=0 riguarda quindi non solo espressioni trascendenti di f(x) ma anche il caso incui f(x) è un polinomio Pn(x) con n>4. Esamineremo alcuni metodi classici per trovare tali approssimazioni.Banalmente, la prima approssimazione possibile è quella che si ottiene dauna tabulazione della funzione, di dato passo h, a partire da un punto x0 .


E’ chiaro infatti che se si opera una tabulazione f(x0+ih), i=0,1,2,... e si haa un certo punto f(x0+ih)·f(x0+(i+1)h)<0, una radice di f(x) è certo in ]x0+ih,x0+(i+1)h[. Si può quindi prendere come approssimazione della radice incognita la quantità, di modo che

h)i(xx~2

10

2

hx~x x

Si può procedere analogamente per cercare altre radici. Si può al contempo ridurre l’errore semplicemente riducendo il passo quando si trova l’intervallo [xi,xi+1] entro cui è la radice, ritabulando lì dentro conh/2 o meno.Con l’accresciuta potenza dei calcolatori tale metodo naif è valido perchè ovviamente semplicissimo da implementare.

Vediamo ora qualche metodo più sofisticato.


Metodo della bisezione

Supponiamo che f(x) sia continua e che a0 e b0 siano tali che f(a0)·f(b0)<0(se f(a0)· f(b0)=0 allora o a0 è radice o lo è b0 o lo sono entrambe) per cui

.b,ax,)x(f:x 000 Ipotizziamo, inoltre, che la radice sia semplice,cioè

.)x('f 0 È possibile costruire una successione di intervalli {In} taleche In+1 In tutti contenenti la radice, per cui:

.xIlim nn

Lo schema di costruzione degli {In}è chiaro dalla figura


Si procede così: sia I0=[a0,b0], con f(a0)<0 e f(b0)>0. Prendiamo il punto di mezzo di I0, m0= a0+(b0-a0)/2=(a0+b0)/2; I1 si costruisce così:

se f(m0)<0 allora I1=[m0,b0], se f(m0)>0 allora I1=[a0,m0],

.)m(fse]m,a[

,)m(fse]b,m[]b,a[I

kkk

kkkkkk 0

0

111

111e così via, per cui il generico Ik è

Chiaramente f(ak)<0, f(bk)>0 (se no si è trovata la radice ak o bk) e.k,Ix k Dopo n passi si ha mis(In)=(bn-an)=(bn-1-an-1)/2=(bn-2-an-2)/2/2=

=···=2-n(b0-a0)= 2-n mis (I0), per cui

,)I(mislim nn

0

da cui, prendendo mn come stima di x , si ha

)ab()I(mis,mx )1n(nnnn 002

2

1

.xIlim nn


La convergenza è quindi certa, ma lenta: poiché 10-12-3.3, ci vogliono, inmedia, 3,3 passi per guadagnare una cifra decimale di precisione.Si noti che la convergenza non dipende da f(x) poiché il metodo fa uso solo di sign(f(x)). Ciò é un vantaggio da una parte, ma costituisce alcontempo il limite del metodo perché non se ne può sveltire la convergenza non utilizzando proprietà della f(x) nè delle sue derivate.

Metodo della tangente (o di Newton-Raphson)

Questo metodo si basa su una stima iniziale, x0, della radice e sull’approssimazione locale della funzione tramite la retta tangente in (x0f(x0)). L’intersezione di tale retta con l’asse x dà una stima successiva x1, e così via (il metodo é iterativo). Il procedimento per determinare x1 è:


i) Si scrive l’espress. della tangente alla curva in (x0, f(x0))

),x('fm),xx(m)x(fy 000

ii) se ne trova l’intersezione x1 con l’asse x (cioè si trova la radice di y(x))

,)x('f

)x(fxx

0

001

iii) si generalizza il procedimento, ottenendo la successione {xn}

),x(Fhx)x('f

)x(fxx nnn

n

nnn 1

iv) si possono fermare le iterazioni quando |hn|=|xn+1-xn|<, dove èla tolleranza prescelta.

Il met. di N-R è un esempio di metodo iterativo (vedi iii)).


Esaminiamo le proprietà di convergenza del metodo diNewton-Raphson, ipotizzando f(x)C2[a,b]. _ _ _Se x è radice semplice, allora f‘(x)0 e f '(x)0 in tutto un intervallo I(x). _Se n=xn-x è l’errore all’n-esimo passo si ha

),x,xint(,)xx)((''f)xx)(x('f)x(f)x(f nnnnnnn 2

2

10

cioè2

2

1)xx(

)x('f

)(''fx

)x('f

)x(fx n

n

n

n

nn

1nx

21 2

1)xx(

)x('f

)(''fxx n

n

nn

per cui

,)x('f

)(''f)xx(

)x('f

)(''fxx n

n

nn

n

nnn

2211 2

1

2

1 _quindi per xnx: .

)x('f

)x(''f

)x('f

)x(''f

n

n

n

n

2

1

2

121

che vuol dire che, vicino alla radice, l’errore (per passo) scala col quadratodell’errore al passo prec. (n+1n

2): il metodo è detto “del 20 ordine”.


Il concetto di ordine di un metodo iterativo si generalizza così: _ _• Sia {xn} una successione che converge a x e sia n=xn-x; se esiste un numero p e una costante C0 tali che

Clim pn

n

n

1

allora si dice che p è l’ordine di convergenza della successione e C è la costante asintotica di errore.


Metodo iterativi: generalità

Un metodo iterativo è quello per cui xn+1 è dato da una funzione dipendente da una m-pla di valori precedenti xn cioè tale che: xn+1=(xn,xn-1,xn-2,...,xn-m+1) dove è chiamata funzione d’iterazione.

i) Nel met. di Newton-Raphson si ha (xn)=f (xn)/f´(xn) (m=1);ii) nel met. della secante, si ha (xn,xn-1)=f (xn)(xn- xn-1)/(f(xn)-f(xn-1) (m=2).

Ovviamente la teoria più semplice si sviluppa per m=1: xn+1= (xn) (iterazione a 1 punto). _Supponiamo che {xn} converga a x. Se è continua allora

),x()x(limxlimx nn

nn

1 per cui x é un punto fisso di , cioè é

radice di F(x)x-(x). Quindi, data la funz. f(x) di cui si cercano le radici, si può provare a scriverla come f(x)=x-(x) per definire il metodo iterativoxn+1= f(xn)=xn-(xn) augurandosi che converga al punto fisso di (x), che


Es. Sia f(x)=x3-x-5. E’ possibile scrivere f(x)=0 in vari modi nella forma x=(x): i) x=x3-5=1(x); ii) x=(x+5)1/3=2(x); iii) x=5/(x2-1)= 3(x). Si possono quindi definire 3 metodi iterativi: xn+1=1(xn); xn+1=2(xn); xn+1=3(xn). Se convergono, allora convergono alla radice di f(x), ma non è detto che convergano.

Si può verificare che la convergenza è assicurata se |´(x)|<1 in un intorno della radice contenente x0 e x1 (vedi figura). Se |´(x)|1laconvergenza c’é solo in casi molto particolari, anche se x0 é già moltovicino alla radice.

che sarebbe radice di f(x).


Vale infatti il teorema: _Ipotesi: (x) é continua, ha un punto fisso x, è derivabile e’(x) m<1 in xx:xJ

allora _ _ Tesi: Per tutti gli x0J, i) xnJ, n=0,1,2,... ii) lim xn=x , iii) x è il solo p. fisso di (x) in J.

Dimostrazione: Dimostriamo per induzione che

i) xnJ, n. Si ha per ip. che x0 J, inoltre se xn-1J allora xnJ, infatti:

xxxxmxx)()x()x(xx nnnnnn 1111


_La tesi ii), cioè che xn converge a x , si ottiene in virtù della catena di diseguaglianze

nnnnn mxxm...xxmxxmxx 02

21

da cui ne scende .mlimxxlim n

nn

n0

_ _Che x sia l’unico punto fisso in J dipende dal fatto che, se y è un altropunto fisso in J, allora si avrebbe

,yxyx)('yx 0

cioè un assurdo, da cui ne segue che il punto fisso è unico.


Ordine dei metodi iterativi

I metodi iterativi a un punto sono di solito del prim’ordine. Se però il puntofisso è tale che la funz. iteratrice (x) ha tutte le prime p-1 derivate nulle in tal punto mentre la p-esima è diversa da zero, allora il metodo è di ordine p, come vediamo.

pn

)p(

n1n )xx(!p

)()x()x(x

Sviluppando (x) in serie di Taylor si ha

da cui si deduce ,!p

)x(lim

)p(

pn

n

n

1 cioè che il metodo è di ordine p.

E’ possibile verificare che il metodo di N-R è almeno del second’ordine perchè la funz. iteratrice (x)=x-f(x)/f(x) ha derivata prima nulla nel punto fisso.


Il second’ordine, nel metodo di N-R, è pagato dal fatto che a ogni passobisogna valutare sia f(x) che la sua derivata. E’ possibile verificare che,in generale, si possono costruire m. iterativi di ordine elevato. Essi sono, però, più dispendiosi, perchè per avere un metodo di ordine p bisogna valutare la funzione e le sue prime p-1 derivate.

Stima dell’errore

Sia ,..,,n,)x~(x~ nnn 2101 la successione dei valori “calcolati” col

metodo iterativo (inficiati quindi da errore di valutazione della funzionee dall’errore di arrotondamento). Sottraendo m. a m. l’eguaglianza

)x(x si ottiene ),x()x~(xx~ nnn 1 da cui, usando il teorema della

media, ),x,x(I,xx)('xx~ nnnnn 1 che porta, dopo sottraz. ad ambo


i membri della quantità ,xx~)x~(' nn 1 a .)x~x~)((')xx~()(' nnnnnn 111

Se si ha |(n)|m<1 e |n|< n, allora ne segue:

.m

x~x~m

m

)('x~x~

)('

)('xx~ n

nnn

nnn

n

nn

1111 111

In questa maggiorazione dell’errore all’ n+1-esimo passo il primotermine dell’ultimo m. stima l’errore di troncamento, il secondo quello dicalcolo delle funzioni implicate. E’ chiaro che per n grande il primotermine è trascurabile rispetto al secondo che quindi è da un certo puntodominante e impedisce una diminuzione d’errore pur al diminuire dell’errore di troncamento.


Accuratezza raggiungibile

Se xn è un’approssimazione di una radice semplice, , di f(x), allora dal teorema della media f(xn)=(x-)f(), I(xn,)=J. Se ne deduce la stima d’erroreindipendente dal metodo

.Jx,m)x('f,m

)x(fx n

n

L’accuratezza entro cui si può determinare è limitata dall’errore nellavalutazione di f(x); se infatti ,)x(),x()x(f)x(f

~nnnn allora, poichè

al meglio si ha ,)x(f~

n 0 il valore esatto di f in xn sarà f(xn) per cui,

se f(x) non varia molto attorno ad , si ha dalle diseguaglianze di sopra:,

mxn

dove = / f (). Quindi il miglior limite d’errore per


qualsiasi metodo è , che si chiama accuratezza massima raggiungibile per la radice . Se f () è molto piccolo allora l’accuratezza è bassa: il problema del calcolo di è mal-condizionato.


Equazioni differenziali e sistemi di equazioni differenziali

Le equazioni dl moto di un oggetto in un campo gravitazionale (pianeti,sonde interplanetarie, satelliti, stelle nelle galassie, ecc.) sono equazioni differenziali ordinarie dove la variabile indipendente è il tempo (t) e quelledipendenti le coordinate spaziali (x,y,z). Il problema differenziale completocomprende anche l’insieme delle condizioni iniziali su posizioni evelocità (le equazioni sono infatti del 20 ordine, poichè riguardano leaccelerazioni). I sistemi di eq. differenziali della meccanica costituisconouna classe di problemi detti, appunto, “ai valori iniziali”, per distinguerlida un’altra classe di problemi differenziali, che sono quelli “ai limiti” o“ai valori al contorno”. Un esempio tipico, in astronomia, di quest’ultimacategoria di problemi differenziali è il sistema d’equazioni che regolanola struttura stellare in equilibrio sferico, viste in una delle prime lezioni.


Le eq. differenziali di cui parleremo sono esclusivamente ordinarie, cioècoinvolgono solo derivate rispetto a una sola variabile indipendente.Un’espressione tipica di sistema di eq. diff. ordinarie con condizioni iniziali è:

00 yy

yfy

)(

),x('

dove y è un vettore a n componenti, come la funzione vettoriale f (x,y).

In generale, la trattazione numerica di eq. differenziali di ordine maggiore di uno richiede la loro riscrittura in termini di sistemi di eq. diff. del prim’ordine utilizzando il fatto che un’eq. diff. di ordine n si puo’ scrivere sotto forma di n eq.diff. del 10 ordine e scrivendo l’equivalente insieme di c. in. sulle nuove funzioni incognite.


Famiglie di soluzioni dell’eq. diff. y=f(x,y) al variare della condizioneiniziale y(0); notare l’evoluzione dell’errore globale, che può portarela soluzione relativa a y(0)=c1 verso quella relativa a y(0)=c2.


Vediamo come esempio l’eq. diff. d’ordine n: F(x,y,y´,y´´,...,y(n))=0 con le c.in. y(0)=y0,y´(0)=y0´,...,y(n-1)(0)= y0

(n-1). La riduzione al 10 ordine s’ottiene introducendo n variabili ausiliarie (u1,u2,...,un) u nel modo seguente: u1=y,u2=y´,...,un=y(n), per cui si ottiene il sistema:

)n(n

'

n

n'n

'

'

y)(u

y)(u

y)(u

)u,...,u,u,x(F

uu

uu

uu

10

02

01

21

1

32

21

0

0

0

0

Per avere una trattazione numerica dell’eq. diff., cioè una soluzione approssimata, bisogna anche esplicitare la n-esima variabile un comeun=f(x,u1,u2,...,un-1), per il motivo che vedremo fra poco.


Ci sono molti metodi per risolvere numericamente eq. diff. e sistemi di eq. diff. ordinarie. Tratteremo qui metodi per la soluzione di problemi aivalori iniziali e non ai limiti.

Metodo di Eulero

Il metodo di Eulero (dell’800) è il più semplice e si basa sull’approssimazione della derivata prima della funzione y(x) tramite unosviluppo di Taylor:

1) “in avanti”: y´(x0)=[y(x0+h)-y(x0)]/h+O(h) che, eguagliando a f(x0,y0),

dà yn+1=yn+hf(xn,yn)+O(h2),

dove l’errore (locale) è O(h2)=(1/2)y´´(n)h2 e si può stimare maggiorando y´´(x)=f´(x,y) nell’intervallo [xn,xn+1] (xn è x0 e xn+1 è x0+h) (questo é il m. di Eulero in avanti).


Il metodo di E. in avanti è esplicito, a un sol punto (one step) e localmen-te del 20 ordine (e quindi globalmente del 10).

2) “all’indietro”: y(x0-h)=y(x0)-y´(x0)h+(1/2)y´´()h2, da cui y´(x0)=[y(x0)-y(x0-h)]/h+O(h) che, eguagliando a f(x0,y0), dà

yn=yn-1+hf(xn,yn)+O(h2).

Anche questo metodo è a un sol punto e localmente del 20 ordine (e

quindi globalmente del 10), però è implicito (l’incognita yn appare anche

come argomento di f).

Entrambi questi metodi sono molto semplici da implementare, hanno losvantaggio di una scarsa precisione e, quello all’indietro (impicito), di richiedere la soluzione di un’equazione a ogni passo, che può essere complicata se f è non-lineare.


Un semplice miglioramento del metodo di Eulero si ottiene utilizzando2 sviluppi di Taylor su 3 punti (xn-1,xn,xn+1) spaziati di h:

yn+1=yn+y(xn)h+(1/2)y(xn)h2+O(h3), yn-1=yn-y(xn)h+(1/2)y(xn)h2+O(h3);

sottraendo la seconda dalla prima rel. si ha: y(xn)=(yn+1-yn-1)/2h+O(h2),per cui, eguagliando a f(xn, yn):

yn+1=yn-1+2hf(xn, yn)+O(h3) (metodo del punto di mezzo o leap-frog).Tale metodo è esplicito, a due punti (multi-step), localmente del 30 ordine (quindi globalmente del 20). Poichè il calcolo di yn+1 richiede laconoscenza di 2 valutazioni precedenti, yn-1 e yn, il metodo ha bisogno diuna valutazione indipendente di y1, oltre alla condizione iniziale y0, per essere utilizzato. Normalmente tale inizializzazione avviene con un passo


di Eulero. Il leap-frog mostra “debole instabilità” come si può vedere con l’esempio dell’equazione

10 )(y

y'y

la cui sol. è y=e-x. L’instabilità ha una possibile cura col leap-frog modificato: si inizializza il leap-frog con un passo di E., poi si usa illeap-frog fino a x+H=x+Nh (N pari; H “macropasso”). Il valore yN

così ottenuto si corregge così:

)y,x(hfyyy NNNNN 12

1

Per poi proseguire per altri N passi col leap-frog inizializzato da Ny

e ).y,x(hfyy NNNN 11 Poichè l’errore è y(x;h)-y(x)=c1(x)h2+c2(x)h4+···

è possibile effettuare l’usuale estrapolazione di Richardson con pk=2k.


Soluzione di y´=-y, y(0)=1 col leap-frog con h=0.1 e inizializzato con 2metodi diversi: y1=0.9 (cerchi neri); y1=0.85 (cerchi bianchi). La sol. esatta in x=0.1 è y(0.1)=0.9048 con 4 dec..Nel secondo caso si nota l’insorgere della debole instabilità.


Soluzione di y´=-y, y(0)=1 con vari metodil leap-frog con h=0.1


Proseguendo nella tratt dei sist. di eq. diff. si sottolinea che ci riferiremo per semplicità a una singola eq. diff. ma tutto ciò che si dirà vale anche per sistemi (basta sostituire il simbolo di vettore a y e f(x,y)).

Per quanto riguarda l’errore, una maggiorazione dell’errore di troncamentoglobale si può fare in analogia con quanto visto per la formula di quadratura numerica del trapezio e del rettangolo, sommando gli errori locali fatti a ogni passo:

)(yh

)ab(

!p

h

N

)(yN

!p

h)(y

!p

h )p(p

N

nn

)p(pN

nn

)p(p

t

1

1

che corrisponde al fatto che, globalmente, un metodo localmente diordine p diventa di ordine p-1.Naturalmente, non conoscendo la soluzione esatta y(x) la maggiorazioned’errore va fatta utilizzando il fatto che y=f(x,y), per cui y= f (x,y)


fy

f

x

f'y

y

f

x

f

dx

df)y,x('f

che permette una stima d’errore effet-tuando la maggiorazione di f (x,y) comefunzione di 2 variabili.

Metodi di Runge-Kutta

Sono metodi che utilizzano valutazioni della funzione f(x,y) in un insieme di punti entro l’intervallo xn, xn+h. Il più semplice di tali metodi è il Metodo di Heun (metodo di R-K del 20 ordine). Si basa su un’approssimazione trapezoidale esplicita, cioè sull’ottenere l’approx. diy in xn+1= xn+h come

,)y,x(f)y,x(fhyy)(

nnnnnn

1

111 2

1 dove al posto di yn+1 (che renderebbe il

metodo implicito) c’è una sua approx. y(1)n+1

data da una passo di Eulero in avanti


).y,x(hfyy nnn)(

n 1

1

Risulta chiaramente che tale metodo equivale alla valutazione di yn+1 come media aritmetica di 2 stime avanzate: y(1)

n+1, appunto, e y(2)n+1

definita come un passo di Eulero semi-implicito:

)y,x(hfyy)(

nnn)(

n

1

1121 per cui: .

yyy

n)(

n)(

n 21

21

1

1

Si verifica che y(x;h)-y(x)=c2(x)h2+c3(x)h3+···+, quindi il metodo è di 20

ordine.

Il più usato tra i metodi di R-K è quello del 40 ordine, definito dalla sequenza di calcoli:


.kkkkyy

),ky,hx(hfk

),ky,hx(hfk

),ky,hx(hfk

),y,x(hfk

nn

nn

nn

nn

nn

43211

34

23

12

1

226

1

2

1

2

12

1

2

1

Poichè, come si può verificare, si hay(x;h)-y(x)=c4(x)h4+c5(x)h5+···+Il metodo è del 40 ordine.

Interpretazione euristica delle formule di R-K

Nell’intervallo [xn,xn+1], dove xn+1 xn+h, la sol. esatta dell’eq. diff. dareb-be

hx

x

nn

n

n

;dx))x(y,x(f)x(y)hx(y L’idea dei metodi di R-K consiste nell’approssimare l’integrale usando i dati di-sponibili. Ad es., il metodo di Heun si


ricava immediatamente se f dipende solo da x. In tal caso, infatti, si puòapprossimare l’integrale con la formula del trapezio:

hx

x

nnn

n

,h)hx(f)x(f

dx)x(f2

che è appunto la formula di Heun per f che non dipende da y. L’errore dellaform. di Heun sarebbe lo stesso della formula trapezoidale se si conoscessef(xn+1,y(xn+1)) da mettere nell’appross. dell’integrale. Poichè invece si usaf(xn+1,y(1)(xn+1)), che ha un errore locale:

2111

k

kn)k(

nnnnn ,h!k

)x(yh)x('yy)x(yy)x(y

(dove y(xn)=f(xn,yn)) ecco che la f. di Heun ha errore che contiene tutte lepotenze 2 di h, mentre quella trapezoidale contiene solo le potenze parimaggiori o uguali a 2.


Similmente, il metodo di R-K del 40 ordine si ricava immediatamente se f=f(x) approssimando l’integrale con la formula di Simpson, considerandoanche il punto di mezzo tra xn e xn+1, xn+1/2xn+h/2:

1

24

6

11

n

n

x

x

nnnnn )hx(f)h

x(f)x(fdx)x(f)x(y)x(y

che è proprio l’espress. di R-K, tenuto conto che k2=k3 poichè f dipende solo da x. L’errore globale è quindi del 40 ordine, come nel metodo di S.,anche nel caso generale f=f(x,y).


Metodi impliciti (predictor-corrector)

Il più semplice metodo implicito è quello trapezoidale

,)y,x(f)y,x(fhyy nnnnnn 111 2

1

(si noti che tale metodo corrisponde alla media aritmetica fra un passo di Eulero in avanti e uno indietro, da cui il nome, anche, di m. “di Euleromodificato”).Il metodo è chiaramente implicito, in quanto yn+1 appare come argomentodi f(x,y); l’espressione è quindi del tipo yn+1=F(xn,xn+1,yn,yn+1). Se, quindi, f è una funzione non-lineare si tratta di risolvere un’eq. (o un sist. d’eq.) non-lineare a ogni passo d’integrazione. Ricordando le consi-derazioni generali sui m. iterativi risulta spontaneo l’utilizzo di un metodo iterativo tipo


),y,x(hfu)y,x(hf)y,x(hfy)y;y,x,x(Fy )k(nnn

)k(nnnnn

)k(nnnn

)k(n 111111

11 2

1

2

1

2

1

in cui l’indice iterativo è in realtà un apice (k).E’ possibile verificare che un criterio suff. per la convergenza del m. iterativo nel caso di un sistema è, in analogia col caso della singola eq.

12

1

n

hy

f

con l’usuale significato dei simboli di norma matriciale, derivate di vettori, ecc.. La convergenza è tanto più rapida quanto più piccola è lanorma della matrice delle derivate della funz. vett. f(x,y).Una scelta iniziale y(0)

n+1 valida e spontanea è quella di un passo di Eulero esplicito: y(0)

n+1= y(0)n+hf(xn,yn). La scelta iniziale si chiama

“predittore” (predictor) e la correzione iterativa “correttore” (corrector),per cui il m. implicito si chiama predictor-corrector.


La fine del procedimento iterativo può avvenire quando è soddisfatta unacondizione d’errore

,yy )k(n

)k(n

1

11

con >0 prescelto, oppure prefissando un numero massimo, kmax , di iterazioni (la cosa migliore è la combinazione dei 2 criteri). Si noti che il metodo di Heun corrisponde alla scelta kmax=1.E’ ovvio che un buon predictor riduce il numero di iterazioni necessarieper arrivare a una buona approssimazione di yn+1.

Metodo di Adams-Bashforth-Moulton

Il più noto, e usato, m. predictor-corrector è quello di Adams-Bashforth-Moulton. E’ del 50 ordine localmente sia nel predictor (Adams-Bashforth)


che nel corrector (Adams-Moulton).

Le espressioni sono

Predictor (A-B)

Corrector (A-M)

).h(Offffhyy nnnnnn5

3211 937595524

1

).h(Offffhyy nnnnnn5

2111 519924

1

Nelle espressioni di sopra fn=f(xn,yn), ecc.. Il predictor serve chiaramente a evitare che il corrector sia una complicata espressione implicita per yn+1 (yn+1 ottenuto col predictor va messo in fn+1 nel corrector).Il metodo risulta quindi globalmente del 40 ordine.


Un metodo alle differenze per un’eq. diff. del 20 ordine

Eq. della forma y=f(x,y), con le c.i. y(a)=, y(a)= s’incontrano spessoin Fisica e Astronomia (le equazioni del moto sono di quel tipo, dove aprimo membro c’è l’accelerazione e l’espressione a secondo membro èla legge di forza ). Naturamente una possibilità di soluzione numerica Passa attraverso la consueta riscrittura come sistema di eq. diff. del 10 ordine.Si possono usare, però, anche approssimazioni dirette (alle differenze) della derivata seconda come quella (ottenibile dalla somma m. a m. di uno sviluppo di Taylor per yn+1 e per yn-1)

,h

yyy)x(''y nnn

n 211 2

e della derivata prima (sempre al 20 ordine e sempre con la combinazione


lineare di 2 sviluppi di Taylor in avanti e indietro):

.h

yy)x('y nn

n 211

Ne risulta il metodo (chiamato metodo centrale esplicito alle differenze)

.hyy,y

,fhyyy nnnn

2

2

110

211

che non può essere utilizzato finchè non si elimina y-1 dall’espress. alledifferenze della c.i. sulla derivata. Tale eliminazione si può fare espri-mendo y-1 tramite la prima relazione del metodo scritta per n=0, ottenendo

.fhhyy,y

,fhyyy nnnn

02

010

211

2

1

2


Problemi stiff

Alcuni problemi differenziali sono tali da essere intrinsecamente difficilida risolvere numericamente in maniera affidabile. Questi problemi sonodetti stiff (rigidi, difficili). Vediamone alcuni esempli.

a) L’eq. y=100y ha soluzione esatta y(x)=c1e10x+c2e-10x. L’esponenzialecrescente è assente quando le c.i. sono y(0)=1 e y(0)=-10. In tal casoc1=0, c2=1 e la sol. è y(x)=e-10x. Applicando al problema detto i metodinumerici precedentemente visti si verifica però che la soluzione dopo unpo’ invece di convergere a zero esplode positivamente o negativamente con andamento e10x, come se c1 fosse diverso da zero. Il motivo è la transizione, per errore di arrotondamento, dalla soluzione“esatta” corrispondente alle c.i. date a una adiacente che corrispondealla generale c. lineare dei 2 esponenziali.


Per capire meglio la cosa si suggerisce per es. di studiare il problemaperturbato y=100y, y(0)=1 e y(0)=-10+, la cui sol. è

y(x)=( /20)e10x+(1- /20) e-10x,

che corrisponde, in pratica, all’effettiva soluzione numerica del problema.

Un altro tipico problema stiff è quello della presenza di “scale temporalimultiple” nella soluzione, come si può vedere con l’esempio del sistema:

.)(v

)(u

vu'v

vu'u

00

10

1999999

1998998

La sol. del sistema si ottiene ponendo u=2y-z e v=-y+z, sostituendo esommando m. a m. e moltiplicando la 2a eq. per 2 e poi sommando m. a m.


Il sistema si disaccoppia in

,z'z

,y'y

1000

cioè:

.ee)x(v

,ee)x(uxx

xx

1000

10002

Nelle sol. ci sono due componenti esponenziali che decadono entrambema con “tempi di decadimento” molto diversi, una delle 2 avendo un “tempo di decadimento” 1000 volte più grande. Questo vuol dire che unmetodo numerico esplicito per essere accettabile deve usare un passo chesia in grado di seguire la soluzione più rapidamente variabile.Poichè una ragionevole scelta del passo si ottiene richiedendo che tra ne n+1 l’incremento relativo di y sia inferiore a una costante prefissata:

,)y,x(f

yh

nn

nn 1

ecco che nell’esempio sopra dato il passo risulta hn+1=min(1,1/1000)


con ovvio sovraccarico di calcoli, che implica lungo tempo d’attesa e aggravio nell’errore accumulato di arrotondamento e troncamento.

Qunado ci sono problemi stiff, instabili, è opportuno ricorrere a metodiimpliciti.Vediamo, infatti, che il metodo esplicito di Eulero soffre di instabilitàper h grande nel caso di un’ eq. tipo y=-cy, c>0.In tal caso il m. di E. esplicito dà: yn+1=(1-ch)yn. Tale metodo diverge se |1-ch|>1, cioè (essendo c e h>0) se ch>2 h>2/c,mentre la sol. esatta converge a zero.

Il metodo di E. implicito applicato all’eq. y=-cy dà invece yn+1=yn/(1+ch)che risulta stabile perchè converge a zero anche se h è grande (se però siusa h grande la soluzione può essere molto poco accurata anche se tendecorrettamente a zero per x grande). Anche il metodo del trapezio (implici-


Un’integrazione accurata e sufficientemente rapida richiede controllo emodifica del passo h. Una possibilità è il controllo e modifica precedente-

to del 20 ordine) è stabile, se applicato all’eq. sopra scritta. Infatti dà

,nn yhc

hcy

21

1

21

1

1

e quindi yn+1 converge a zero. Tutte queste considerazioni permangono valide per sistemi di eq. diff. li-neari del tipo y=-Cy, dove C è una matrice definita positiva, e anche persistemi y=f(x,y), dopo linearizzazione di f(x,y).

Controllo del passo d’integrazione


mente visti, che consistevano nell’ utilizzare (per avanzare la soluzione da xn a xn+1) un passo, hn+1, che fosse abbastanza piccolo da limitare la Variazione relativa di y tra yn e yn+1 ottenibile con un passo di Eulero in avanti.Tale metodo è grossolano ma semplice da implementare.Un controllo forse migliore è quello che viene dalla scelta di un passo h che limiti l’errore per passo a un valore prefissato.Esso si basa su una scelta di h e sulle valutazioni y(1)

n+1 e y(2)n+1 a

xn+1=xn+h ottenute, rispettivamente, con un passo h e con 2 passi h/2.Ricordando l’espressione che dà origine all’estrapolazione di Richardson

,h

cyy

y)hx(ylp

np

)(n

)(n)(

nn

111

212

1 212

si può fermare l’operazione di dimezzamento del passo quando si haln<h.

Metodi matematici dellastronomia Che cosa sono i Metodi matematici dellastronomia? In che cosa sono...

Documents

Transcript of Metodi matematici dellastronomia Che cosa sono i Metodi matematici dellastronomia? In che cosa sono...