Metodi matematici dell’astronomia

Metodi matematici dell’astronomia

Equazioni differenziali e sistemi di equazioni differenziali(continua)

Si ricorda innanzitutto che, pur riferendoci per semplicità a una singola eq.diff., tutto ciò che si dice 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 n)(

nnnnn 11

11 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 n)(

nn)(

n 11

121 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(yL’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 nelcorrector (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 chiaramentea 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 membroc’è l’accelerazione e l’espressione a secondo membro èla legge di forza ). Naturamente una possibilità di soluzione numerica passaattraverso la consueta riscrittura come sistema di eq. diff. del 10 ordine.Si possono usare, però, anche approssimazioni dirette (alle differenze) delladerivata 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). Vediamo con 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 caso lac1=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,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 asistemi y=f(x,y), dopo linearizzazione di f(x,y).

Controllo del passo d’integrazione


mente vista, che consisteva nell’ utilizzare (per avanzare la soluzione da xn

a xn+1) un passo, hn+1, che fosse abbastanza piccolo da limitare la variazionerelativa 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 chelimiti 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 dell’astronomia

Documents

Transcript of Metodi matematici dell’astronomia