Metodi matematici dell’astronomia

1. 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 troncamento globale 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: che corrisponde al fatto che, globalmente, un metodo localmente di ordine p diventa di ordine p-1. Naturalmente, non conoscendo la soluzione esatta y(x) la maggiorazione d’errore va fatta utilizzando il fatto che y=f(x,y), per cui y= f (x,y)

2. Metodi matematici dell’astronomia che permette una stima d’errore effet- tuando la maggiorazione di f’(x,y) come funzione 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. di y in xn+1= xn+h come 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

3. Metodi matematici dell’astronomia 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 per cui: 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:

4. Metodi matematici dell’astronomia Poichè, come si può verificare, si ha y(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 L’idea dei metodi di R-K consiste nell’ approssimare l’integrale usando i dati di- sponibili. Ad es., il metodo di Heun si

5. Metodi matematici dell’astronomia ricava immediatamente se f dipende solo da x. In tal caso, infatti, si può approssimare l’integrale con la formula del trapezio: che è appunto la formula di Heun per f che non dipende da y. L’errore della form. di Heun sarebbe lo stesso della formula trapezoidale se si conoscesse f(xn+1,y(xn+1)) da mettere nell’appross. dell’integrale. Poichè invece si usa f(xn+1,y(1)(xn+1)), che ha un errore locale: (dove y(xn)=f(xn,yn)) ecco che la f. di Heun ha errore che contiene tutte le potenze 2 di h, mentre quella trapezoidale contiene solo le potenze pari maggiori o uguali a 2.

6. Metodi matematici dell’astronomia Similmente, il metodo di R-K del 40 ordine si ricava immediatamente se f=f(x) approssimando l’integrale con la formula di Simpson, considerando anche il punto di mezzo tra xn e xn+1, xn+1/2xn+h/2: 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).

7. Metodi matematici dell’astronomia Metodi impliciti (predictor-corrector) Il più semplice metodo implicito è quello trapezoidale (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 Eulero modificato”). Il metodo è chiaramente implicito, in quanto yn+1 appare come argomento di 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

8. Metodi matematici dell’astronomia 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. con l’usuale significato dei simboli di norma matriciale, derivate di vettori, ecc.. La convergenza è tanto più rapida quanto più piccola è la norma 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.

9. Metodi matematici dell’astronomia La fine del procedimento iterativo può avvenire quando è soddisfatta una condizione d’errore 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 necessarie per 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

10. Metodi matematici dell’astronomia (Adams-Bashforth) che nel corrector (Adams-Moulton). Le espressioni sono Predictor (A-B) Corrector (A-M) 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+1ottenuto col predictor va messo in fn+1 nel corrector). Il metodo risulta quindi globalmente del 40 ordine.

11. Metodi matematici dell’astronomia 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 spesso in Fisica e astronomia (le equazioni del moto sono di quel tipo, dove a primo membroc’è 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) e della derivata prima (sempre al 20 ordine e sempre con la combinazione

12. Metodi matematici dell’astronomia lineare di 2 sviluppi di Taylor in avanti e indietro): Ne risulta il metodo (chiamato metodo centrale esplicito alle differenze) Che non può essere utilizzato finchè non si elimina y-1 dall’espress. Alle differenze 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

13. Metodi matematici dell’astronomia Problemi stiff Alcuni problemi differenziali sono tali da essere intrinsecamente difficili da risolvere numericamente in maniera affidabile. Questi problemi sono detti stiff (rigidi, difficili). Vediamo con alcuni esempli. • L’eq. y=100y ha soluzioneesatta y(x)=c1e10x+c2e-10x. L’esponenziale • crescente è assente quando le c.i. sono y(0)=1 e y(0)=-10. In tal caso la • c1=0, c2=1 e la sol. è y(x)=e-10x. Applicando al problema detto i metodi • Numerici precedentemente visti si verifica però che la soluzione dopo un • po’ 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 corrisponde • alla generale c. lineare dei 2 esponenziali.

14. Metodi matematici dell’astronomia Per capire meglio la cosa si suggerisce per es. di studiare il problema perturbato 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 temporali Multiple” nella soluzione, come si può vedere con l’esempio del sistema: La sol. del sistema si ottiene ponendo u=2y-z e v=-y+z, sostituendo e sommando m. a m. e moltiplicando la 2a eq. per 2 e poi sommando m. a m.

15. Metodi matematici dell’astronomia Il sistema si disaccoppia in cioè: Nelle sol. ci sono due componenti esponenziali che decadono entrambe ma con “tempi di decadimento” molto diversi, una delle 2 avendo un “tempo di decadimento” 1000 volte più grande. Questo vuol dire che un metodo numerico esplicito per essere accettabile deve usare un passo che Sia in grado di seguire la soluzione più rapidamente variabile. Poichè una ragionevole scelta del passo si ottiene richiedendo che tra n e n+1 l’incremento relativo di y sia inferiore a una costante  prefissata: ecco che nell’esempio sopra dato il passo risulta hn+1=min(1,1/1000)

16. Metodi matematici dell’astronomia 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 metodi Impliciti. 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ò si usa h grande la soluzione può essere molto poco accurata anche se tende correttamente a zero per x grande). Anche il metodo del trapezio (implici-

17. Metodi matematici dell’astronomia to del 20 ordine) è stabile, se applicato all’eq. sopra scritta. Infatti dà 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 a sistemi y=f(x,y), dopo linearizzazione di f(x,y). Controllo del passo d’integrazione Un’integrazione accurata e sufficientemente rapida richiede controllo e modifica del passo h. Una possibilità è il controllo e modifica precedente-

18. Metodi matematici dell’astronomia 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 variazione relativa di y tra yne 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 Si può fermare l’operazione di dimezzamento del passo quando si ha ln<h.