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:
N
p
p
h N ( p)
h
t   y ( n )  N
p! n1
p!
y
( p)
n 1
N
( n )
hp (b  a ) ( p)

y ( )
p! h
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)
Metodi matematici dell’astronomia
f ' ( x, y ) 
df f f
f f


y' 

f
dx x y
x y
che permette una stima d’errore effettuando la maggiorazione di f’(x,y) com
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


1
yn 1  yn  h f ( xn , yn )  f ( xn 1 , y ( 1 ) n 1 ) ,
2
dove al posto di yn+1 (che renderebbe il
metodo implicito) c’è una sua approx. y(1)n+
data da una passo di Eulero in avanti
Metodi matematici dell’astronomia
yn( 11)  yn  hf ( xn , yn ).
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
(2)
n 1
(1)
(2)
y

y
n

1
n 1
 yn  hf ( xn 1 , y n 1 ) per cui: yn 1 
.
2
(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:
Metodi matematici dell’astronomia
k1  hf ( xn , yn ),
1
1
k 2  hf ( xn  h , yn  k1 ),
2
2
1
1
k3  hf ( xn  h , yn  k 2 ),
2
2
k 4  hf ( xn  h , yn  k3 ),
yn 1  yn 
Poichè, come si può verificare, si ha
y(x;h)-y(x)=c4(x)h4+c5(x)h5+···+
Il metodo è del 40 ordine.
1
k1  2k2  2k3  k4 .
6
Interpretazione euristica delle formule di R-K
Nell’intervallo [xn,xn+1], dove xn+1 xn+h, la sol. esatta dell’eq. diff. darebbe
L’idea dei metodi di R-K consiste nell’
y( x  h )  y( x )   f ( x , y( x ))dx;
approssimare l’integrale usando i dati disponibili. Ad es., il metodo di Heun si
xn  h
n
n
xn
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:
xn  h

xn
f ( x )dx 
f ( x n )  f ( xn  h )
h,
2
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:

y( k ) ( xn ) k
y( xn1 )  yn1  y( xn1 )  yn  y' ( xn )h  
h ,
k!
k 2
(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.
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:
y( xn 1 )  y( xn ) 
xn1

xn
1
h

f ( x )dx   f ( xn )  4 f ( xn  )  f ( xn  h ) 
6
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).
Metodi matematici dell’astronomia
Metodi impliciti (predictor-corrector)
Il più semplice metodo implicito è quello trapezoidale
1
yn 1  yn  h f ( xn , yn )  f ( xn 1 , yn 1 ),
2
(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 considerazioni generali sui m. iterativi risulta spontaneo l’utilizzo di un
metodo iterativo tipo
Metodi matematici dell’astronomia
1
1
1
yn( k11 )  F ( xn , xn 1 , yn ; yn( k1) )  yn  hf ( xn , yn )  hf ( xn 1 , yn( k1) )  un  hf ( xn 1 , yn( k1) ),
2
2
2
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.
1 f
h
1
2 yn
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.
Metodi matematici dell’astronomia
La fine del procedimento iterativo può avvenire quando è soddisfatta una
condizione d’errore
yn( k11 )  yn( k1)   ,
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
Metodi matematici dell’astronomia
(Adams-Bashforth) che nel
corrector (Adams-Moulton).
Le espressioni sono
Predictor (A-B)
yn 1  yn 
1
h55 f n  59 f n 1  37 f n  2  9 f n 3   O( h 5 ).
24
Corrector (A-M)
yn 1  yn 
1
h9 f n 1  19 f n  5 f n 1  f n  2   O( h 5 ).
24
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.
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 pass
attraverso la consueta riscrittura come sistema di eq. diff. del 10 ordine.
Si possono usare, però, anche approssimazioni dirette (alle differenze) dell
derivata seconda come quella (ottenibile dalla somma m. a m. di uno
sviluppo di Taylor per yn+1 e per yn-1)
y' ' ( xn ) 
yn 1  yn 1  2 yn
,
h2
e della derivata prima (sempre al 20 ordine e sempre con la combinazione
Metodi matematici dell’astronomia
lineare di 2 sviluppi di Taylor in avanti e indietro):
y' ( xn ) 
yn 1  yn 1
.
2h
Ne risulta il metodo (chiamato metodo centrale esplicito alle differenze)
yn 1  2 yn  yn1  h 2 f n ,
y0   , y1  y1  2h .
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 esprimendo y-1 tramite la prima relazione del metodo scritta per n=0, ottenend
yn1  2 yn  yn1  h 2 f n ,
1
y0   , y1  y0  h  h 2 f 0 .
2
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.
a) L’eq. y=100y ha soluzione esatta 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.
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:
u'  998u  1998v
v'  999u  1999v


u( 0 )  1
v( 0 )  0.
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.
Metodi matematici dell’astronomia
 y'   y ,

 z'  1000 z ,
Il sistema si disaccoppia in
cioè:
u( x )  2e  x  e 1000x ,

v( x )  e  x  e 1000x .
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:
hn1  
yn
,
f ( xn , yn )
ecco che nell’esempio sopra dato il passo risulta hn+1=min(1,1/1000)
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-
Metodi matematici dell’astronomia
to del 20 ordine) è stabile, se applicato all’eq. sopra scritta. Infatti dà
1
1  hc
2 y
yn 1 
n,
1
1  hc
2
e quindi yn+1 converge a zero.
Tutte queste considerazioni permangono valide per sistemi di eq. diff. lineari 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-
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 variazion
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 ch
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
ln  y( x  h )  y
(2)
n 1

yn( 21)  yn( 11)
2 1
p
h
 cn  
2
p 1
,
Si può fermare l’operazione di dimezzamento del passo quando si ha
ln<h.