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 dei
fenomeni su base fisico-matematica, ma che mutua alcune caratteristiche
dalle Scienze naturali, quali l’accurata raccolta e analisi di “reperti”
osservativi (tassonomia). Ciò implica la necessità di utilizzare metodi
statistici di indagine. Inoltre l’astronomia-astrofisica moderna è caratterizzata 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.
Metodi matematici dell’astronomia
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
Metodi matematici dell’astronomia
• 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.
Metodi matematici dell’astronomia
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
Metodi matematici dell’astronomia
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”
Metodi matematici dell’astronomia
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)
point spread
function
k(x,y)
risposta strumentale
g(x)
convoluzione
g(x)
F(y)
deconvoluzione
Un’approssimazione F(y) di f(y) si ottiene per deconvoluzione della
risposta strumentale, ma quanto vale F( y )  f ( y ) ?
Metodi matematici dell’astronomia
Metodi matematici dell’astronomia
g(’)
f()
Metodi matematici dell’astronomia
Il problema della deconvoluzione (o inversione dei dati).
x
x
g( x )   k( x , y ) f ( y )dy
caso k(x,y)=1
0
ha soluzione formale
g ( x)   f ( y)dy
0
f ( x)  g ' ( x)
Semplice, ma dà problemi; infatti supponiamo che
g( x )  1  e x  f ( x )  g' ( x )   ex
Se al dato g(x) è sovrapposto rumore di frequenza  si ha:
g o ( x )  g( x )  senx  f o ( x )  e x   cos x
Ne consegue che le 2 soluzioni f(x) e fo(x) differiscono
f ( x )  f o ( x )  max f  f o  max  cos  x  
x
x
Metodi matematici dell’astronomia
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
  
df o 1  f o
f o
d
d 
 
d 
d  
(cos  x   xsen x )   cos  x  
fo
f o  

f


o 
le variazioni relative d/ e d/ sono amplificate dal fattore /fo.
Metodi matematici dell’astronomiainstabilità
instabilità
output
input
=0.8, =0.04,=20
Metodi matematici dell’astronomia
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)+acosy.
Il relativo osservabile sarà
x
g o 2 ( x )  g o1( x )  a  cos ydy g o1( x ) 
0
per cui
f  a cos y  g o 
a

a

senx
senx  g o 
a

che corrisponde al problema visto prima. Una perturbazione in f di grande
ampiezza a può essere smorzata dalla sua alta frequenza ω.
Metodi matematici dell’astronomia
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 osservativi
per esempio flussi in bande diverse dello spettro e.m..
Metodi matematici dell’astronomia
Problemi tipici dell’astronomia
E  esenE 
2
(t   )
P
Equazione di Keplero
N
mi m j
d 2 ri
r  rj , (i  1,2,..., N )
mi 2  G 
3 i
dt
j 1 | ri  r j |
j i
ri (0)  ri 0
.
.
ri (0)  ri 0
(i  1,2,..., N )
Equazioni del moto di N
corpi autogravitanti
sistema di N eq. diff. vett. ordinarie del 20 ordine
 ordine del sistema = 6N
servono 6N condizioni iniziali su posizioni e velocità
Metodi matematici dell’astronomia
Equazioni della struttura stellare
1 dp
M( r )
 G 2
 dr
r
dM
 4r 2
dr
p  p g  pr 

m H
Sistema di 4 eq. diff. + 1 EOS con incognite:
(r). P(r), M(r), T(r), L(r)
kT 
dT
3 L( r ) 1

dr
4 4r 2 T 3
dL
 4r 2 
dr
M ( 0)  0
p ( 0)  0
p( R)  pR
T ( R )  TR
1 4
T
3c
mentre:
= (,T;Xi), = (,T;Xi), = (,T;Xi)
sono funzioni note.
4 condizioni al contorno
Metodi matematici dell’astronomia
2
E  esenE 
( t  )
P
Si vuole ottenere E(t); è necessario quindi
risolvere un’equazione trascendente:
f ( E ; t )  E  esenE 
2
( t  )  0
P
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
f ( 0; t )  
2
( t  )  0
P
f ' ( E ; t )  1  e cos E  0
f ( 2 ; t )  2 

2
( t  )  0
P
per il teor. dell’esist. degli zeri c’è 1 sola
radice nell’intervallo ]0,P[
Metodi matematici dell’astronomia
Come ottenere e stimare l’accuratezza nei calcoli numerici
• Dati 2 numeri a e b è logico confrontarli:
Se a e b rappresentano le altezze di 2 persone e b-a=50cm
a  b ha senso dire che a<<b, ma se parliamo di due montagne no!
Lo stesso si può dire per ab e ab.
Il significato del confronto, quindi, dipende dal contesto.
• f ( x )  O( g( x ))
per
x a
significa che
lim
xa
f ()
g( x )
 
Metodi matematici dell’astronomia
Errori assoluti e relativi
Indichiamo d’ora in poi con ã un’approssimazione di a.
Data un’approssimazione ã di a:
• l’errore assoluto è a= ã - a
• l’errore relativo è r= (ã – a)/a (se a0)
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.
Metodi matematici dell’astronomia
Sorgenti d’errore
I risultati numerici sono influenzati da varie sorgenti d’errore: alcune
possono essere ridotte o eliminate, altre no.
Possibili errori:
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 quellin che nascono dal “taglio” di
un’espressione, tipo: sn   ai
invece di s   a

i 1
i 1
i
oppure dall’aprossimare un funzione non lineare con una lineare
oppure dall’approssimare una derivata con un rapporto finito (errore
di discretizzazione);
Metodi matematici dell’astronomia
Arrotondamento (round-off)
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
a=0.08 ha 2 decimali e 1 cifra
b=15.6 ha 1 decimale e 3 cifre
Indichiamo d’ora in poi con ã un’approssimazione di a.
• Siano a=0.235 e ã=0.231; ã ha 2 cifre esatte, e ciò corrisponde a un
errore  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
Metodi matematici dell’astronomia
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.2660.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 migliore
precisione (2 dec. corretti).
Metodi matematici dell’astronomia
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.
Metodi matematici dell’astronomia
Nelle operazione di moltiplicazione e divisione, invece, si sommano i
limiti d’errore relativo (approssimativamente).
Si verifica infatti (per esercizio) che:
(ã1ã2 - a1a2)/ a1a2 = r1 +r2 + r1r2  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.
Metodi matematici dell’astronomia
Cancellazione dei termini
Una causa comune di scarsa accuratezza nel calcolo è la sottrazione di
2 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
y  x1  x2 
y
y

x1  x2
x1  x2
L’errore relativo in y può quindi
essere 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%).
Metodi matematici dell’astronomia
Bisogna quindi cercare di evitare differenze di numeri simili, riscrivendo,
se possibile, le formule. Per esempio, sia da calcolare x    x ,  x
Essa si può riscrivere come:
x   x 

x   x

x   x

x   x
x   x

x   x

x   x
evitando così la cancellazione dei termini.
Dovendo calcolare f(x+)-f(x) un’altra possibilità è quella di sviluppare
f(x) in serie di Taylor:
f(x  δ)  f(x)  f’(x)δ 
1
f ' ' ( x )δ 2  ...
2
Metodi matematici dell’astronomia
La formula generale per la propagazione dell’errore
Sia data una funzione y(x) e si voglia valutare
di x. Ci si chiede qual è l’errore y  y( ~x )  y( x )
y( ~
x)
dove ~x è un’approx.
Una via naturale è quella di calcolare il differenziale (approssimato):
y  dy  y' ( ~
x )x
In generale, se y=y(x1,x2,...,xn) si può stimare
n
y  
i 1
y ~
( x ) xi
xi
n
y  
i 1
y ~
( x )xi
xi
da cui:
Che è la formula generale di propagazione dell’errore.
Metodi matematici dell’astronomia
Derivata numerica
Conoscendo una funzione f(x) in 3 punti (xi-1,xi,xi+1), si possono costruire
3 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) se
si 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
Metodi matematici dell’astronomia
Derivata numerica
di punto di partenza x, esplicitando le derivate prime e semisommandole:
f ( x ) f ( x  h )
f 's ( x ) 
 O( h );
h
f ( x  h ) f ( x )
f 'd ( x ) 
 O( h );
h
f ( x  h ) f ( x  h )
f 'c ( x ) 
 O( h 2 ).
2h
Metodi matematici dell’astronomia
Integrazione numerica
Supponiamo di dover calcolare
b
I   f ( x )dx
a
dove f(x) è una funzione
continua nell’intervallo a ,b oppure è una funzione campionata per
punti: {f(xi}, i=1,2,...,n}. Se non si riesce a ottenerne la primitiva
F'(x)=f(x), per cui F(b)-F(a)=I, l’unica alternativa è di approssimare
~
I 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 basano
sul fatto che l’integrale definito è il limite della somma finita sn, cioè:
n
I  lim sn  lim  f ( xi ) i x
n 
n 
i 1
Metodi matematici dell’astronomia
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);
Metodi matematici dell’astronomia
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è:
n 1
( f  f i 1 )hi
~
~
I t   f ( x )dx   i
2
i 0
a
b
Tale approssimazione equivale ad aver sostituito (dentro l’integrale)
ad f la sua approssimazione costituita dalla spezzata passante per i
punti (xi,fi), cioè:
n 1


f f
~
f ( x )    f i  i 1 i ( x  xi ) ( x  xi ) ( xi 1  x )
xi 1  xi
i 0 

dove (x) è la funzione di Heaviside:  ( x )  10,,xx00
b
~
Verificare per esercizio che  f ( x )dx  ~I .
t
a
Metodi matematici dell’astronomia
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 la
complessità 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 è:
n 1
~
I r ,c   f i 1 / 2 hi
i 0
(form. rett. centrata).
Se il passo è fisso, hi=h=(b-a)/n, la formula
diviene
n 1
~
I r ,c  h f i 1 / 2
i 0
Metodi matematici dell’astronomia
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
f. rett. sin.
n 1
~
I r ,s   f i hi
i 0
f. rett. des.
Il metodo rettangolare
n 1
~
I r ,d   f i 1hi
i 0
Il metodo trapezoidale
Metodi matematici dell’astronomia
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
( x  x
n
pn ( x )   f
i 0
k
)
k 0
k i
i n
( x  x
i
k 0
k i
k
)
(formula d’interpolazione di
Lagrange; verificare per es. che
pn(xi)= fi, i=0,1,...,n.)
n
che si scrive compattamente come
pn ( x )   f i Li ( x ),
i 0
n
( x  x
avendo posto
Li ( x ) 
k
)
k 0
k i
n
( x  x
i
k 0
k i
(funzione moltiplicatirce di Lagrange).
k
)
Metodi matematici dell’astronomia
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 è:
f ( n1 ) (  )
Rn ( x )  f ( x )  pn ( x ) 
( x  x0 )( x  x1 )    ( x  xn ),
( n  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 nel
generico intervallo [xi,xi+1] semplicemente integrando R1(x) su tale int..
Metodi matematici dell’astronomia
xi1
xi1
xi1
xi1
xi
xi
xi
xi
 i   R1( x )   f ( x ) 
 p1( x )  
f '' ( )
( x  xi )( x  xi 1 )dx
2
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:
f ' ' (  i ) i1
i 
x ( x  xi )( x  xi1 )dx , i  [ xi , xi1 ]
2
i
x
Ponendo x=xi+hit, dove hi=xi+1-xi, si ha
i 
f ' ' ( i ) 3 1
f ' ' ( i ) 3  1 1 
f ' ' ( i ) 3
hi  t( t  1 )dt 
hi     
hi
2
2
3
2
12


0
Ricordando che l’err. della somma è la somma degli errori, l’err. di tronc.
globale è dato da che, se i punti x sono equispaziati x -x =h=(b-a)/n,
i
i+1
i
diventa (applicando il teor. del val.medio)
 t ,g
n 1
1 n 1
3
   i    f ' ' (  i )hi ,
12 i 0
i 0
n 1
 t ,g    i  
i 0
h3 n1
h3 n1 f ' ' (  i )
h3
( b  a )h 2
f
'
'
(

)


n


nf
'
'
(

)


f '' ( )

 n
i
12 i 0
12 i 0
12
12
Metodi matematici dell’astronomia
Quindi la formula trapezoidale è localmente del 30 ordine e globalmente
del 20.
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]
f ( x )  f ( xi 1/2 )  f ' ( xi 1/2 )( x  xi 1/2 ) 
f '' ( )
( x  xi 1/2 )2
2
per cui il resto dell’interpol. di f(x) in [xi,xi+1] rispetto al suo valore in
xi+1/2 è
f '' ( )
R( x )  f ( x )  f ( xi 1/2 )  f ' ( xi 1/2 )( x  xi 1/2 ) 
2
( x  xi 1/2 )2 .
Perciò l’err. di troncamento della form. rettangolare nell’interv. [xi,xi+1] è
xi 1
 t   [ f ( x )  f ( xi 1/2 )] dx  
xi
xi 1
 f'( x
i 1/2
xi
)( x  xi 1/2 )dx 
xi1

xi
f '' ( )
( x  xi 1/2 )2 dx
2
Metodi matematici dell’astronomia
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
f ' ' ( i ) 3 1
f ' ' ( i ) 3  1 1 1 
f ' ' ( i ) 3 1
f ' ' ( i ) 3
1
t 
h  ( t  )2 dt 
h    
h

h ,
2
2
2
3
4
2
2
12
24


0
per cui l’err. di troncamento globale è
 t ,g
( b  a )h 2

f ' ' (  ).
24
Metodi matematici dell’astronomia
Metodo di Simpson
Risponde all’esigenza di avere una formula di quadratura di ordine
elevato 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.
xi 1
Scritta la formula  f ( x )dx  haf
i 1
xi 1
 bf i  cf i 1 
si tratta di trovare i coefficienti
a,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
f ( x )  f ( 0 )  f ' ( 0 )x 
che integrata in [x-1,x1] dà
1
1
f ' ' ( 0 )x 2  f ' ' ' ( 0 )x 3  O( h 4 ),
2
3!
Metodi matematici dell’astronomia
h

h
x

2
 
f ( x )dx  2hf ( 0 )  f ' ( 0 ) 
2
h

h
3
x
1
f '' ( 0 )
2
3



h

h
4
x
1
f ''' ( 0 )
3!
4



h
 O( h 4 )2h ,
h
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:
f '' ( x ) 
1
 f ( x  h )  f ( x  h )  2 f ( x )  O( h 2 ).
2
h
Metodi matematici dell’astronomia
Valutando l’espressione precedente in x=0 e inserendola nell’integrale
si ha
h

f ( x )dx  2hf 0 
h
h
 f1  f 1  2 f 0   O( h5 )  h 4 f 0  1 f 1  1 f1   O( h5 ),
3
3
3 
3
Da cui la formula generale, nell’intervallo [xi-1,xi+1]
xi1

f ( x )dx 
xi1
h
 fi 1  4 f i  f i 1   O( h5 ),
3
esatta per polinomi fino al 30 grado, che è la formula di Simpson.
Si verifica che l’errore di troncamento locale è
h5 ( iv )

f (  ),  xi 1 , xi 1 
90
per cui quello globale è
 t ,g
h 5 ( iv )
h 5 n ( iv )
h4

f ( i ) 
f ( ) 
( b  a ) f ( iv ) (  i ),  a ,b
90 2
180
i pari 90
(il fattore 2 a dividere viene dal fatto che h=(b-a)/n, per cui ci sono n/2
Metodi matematici dell’astronomia
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

f ( r )d 3r  f V   f V  f V 

n
V
V
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:

f


n

f 
f
n

2

f2  f
n
2
Metodi matematici dell’astronomia
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 rettangolare
A di base [a,b] e altezza max f(x), stimando IA come
n
b
I   f ( x )dx 
a
 k( x )
i
i 1
n
n
area( A ) 
 k( x )
i
i 1
n
( b  a ) f ( xmax )
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 di
punti. Definendo la funzione k(x) come
0, yi  f ( xi ),
k( xi )  
1, yi  f ( xi ),
n
risulta
0  n*   k ( xi )  n
i 1
I punti (xi,yi) sono indicati;
nel caso in fig. xmax=b
Metodi matematici dell’astronomia
Poichè l’ordine del metodo è h1/2 (la stima d’errore converge a zero con
n-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 di
Abel-Ruffini: non è possibile esprimere sotto forma di radicali le soluzioni
di 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 in
cui 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 da
una tabulazione della funzione, di dato passo h, a partire da un punto x0 .
Metodi matematici dell’astronomia
E’ chiaro infatti che se si opera una tabulazione f(x0+ih), i=0,1,2,... e si ha
a 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
x
1
~
x  x0  ( i  )h
2
h
x~
x
2
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 con
h/2 o meno.
Con l’accresciuta potenza dei calcolatori tale metodo naif è valido
perchè ovviamente semplicissimo da implementare.
Vediamo ora qualche metodo più sofisticato.
Metodi matematici dell’astronomia
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
x : f ( x )  0, x  a0 ,b0 . Ipotizziamo, inoltre, che la radice sia semplice,cioè
f ' ( x )  0.
È possibile costruire una successione di intervalli {In} tale
che In+1 In tutti contenenti la radice, per cui:
lim I n  x .
n
Lo schema di
costruzione degli {In}
è chiaro dalla figura
Metodi matematici dell’astronomia
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],
[ mk 1 ,bk 1 ] se f ( mk 1 )  0,
I k  [ ak ,bk ]  
[ ak 1 , mk 1 ] se f ( mk 1 )  0.
e così via, per cui il generico Ik è
Chiaramente f(ak)<0, f(bk)>0 (se no si è trovata la radice ak o bk) e
x  I k ,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
lim mis( I n )  0 ,
n 

lim I n  x .
n
da cui, prendendo mn come stima di x , si ha
x  mn   n , n 
1
mis( I n )  2 ( n 1 ) ( b0  a0 )
2
Metodi matematici dell’astronomia
La convergenza è quindi certa, ma lenta: poiché 10-12-3.3, ci vogliono, in
media, 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 al
contempo 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 è:
Metodi matematici dell’astronomia
i) Si scrive l’espress. della tangente alla curva in (x0, f(x0))
y  f ( x0 )  m( x  x0 ), m  f ' ( x0 ),
ii) se ne trova l’intersezione x1 con l’asse x (cioè si trova la radice di y(x))
x1  x0 
f ( x0 )
,
f ' ( x0 )
iii) si generalizza il procedimento, ottenendo la successione {xn}
xn 1  xn 
f ( xn )
 xn  hn  F ( xn ),
f ' ( xn )
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)).
Metodi matematici dell’astronomia
Esaminiamo le proprietà di convergenza del metodo di
Newton-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
1
0  f ( x )  f ( xn )  f ' ( xn )( x  xn )  x f ' ' (  n )( x  xn )2 , n  int( x , xn ),
2
n 1
cioè
xn 
 n 1  xn 1  x 
f ( xn )
1 f ' ' ( n )
x
( x  x n )2
f ' ( xn )
2 f ' ( xn )
1 f ' ' (n )
1 f ' ' ( n ) 2
( x  x n )2 
n ,
2 f ' ( xn )
2 f ' ( xn )

xn 1  x 
1 f ' ' ( n )
( x  x n )2
2 f ' ( xn )
_ 
quindi per xnx: 
n 1
2
n

per cui
1 f ' ' ( xn ) 1 f ' ' ( x )

.
2 f ' ( xn ) 2 f ' ( x )
che vuol dire che, vicino alla radice, l’errore (per passo) scala col quadrato
dell’errore al passo prec. (n+1n2): il metodo è detto “del 20 ordine”.
Metodi matematici dell’astronomia
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 C0 tali che
 n1
lim
C
n   p
n
allora si dice che p è l’ordine di convergenza della successione e C è la
costante asintotica di errore.
Metodi matematici dell’astronomia
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  lim xn 1  lim  ( xn )   ( x ),
n 
n 
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 iterativo
xn+1= f(xn)=xn-(xn) augurandosi che converga al punto fisso di (x), che
Metodi matematici dell’astronomia
che sarebbe radice di f(x).
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)|1la
convergenza c’é solo in casi molto particolari, anche se x0 é già molto
vicino alla radice.
Metodi matematici dell’astronomia
Vale infatti il teorema:
_
Ipotesi: (x) é continua, ha un punto fisso x, è derivabile e
’(x) m<1 in J  x : x  x  
allora
_
_
Tesi: Per tutti gli x0J, i) xnJ, n=0,1,2,... ii) lim xn=x , iii) x è il solo p.
fisso di (x) in J.
Dimostrazione: Dimostriamo per induzione che
i) xnJ, n. Si ha per ip. che x0 J, inoltre se xn-1J allora xnJ,
infatti:
xn  x  ( xn1 )  ( x )   (  n ) xn1  x  m xn1  x  xn1  x  
Metodi matematici dell’astronomia
_
La tesi ii), cioè che xn converge a x , si ottiene in virtù della catena di
diseguaglianze
xn  x  m xn1  x  m2 xn2  x  ...  mn x0  x  mn 
da cui ne scende
lim xn  x   lim mn  0.
n
n
_
_
Che x sia l’unico punto fisso in J dipende dal fatto che, se y è un altro
punto fisso in J, allora si avrebbe
0  x  y  ' (  ) x  y  x  y ,
cioè un assurdo, da cui ne segue che il punto fisso è unico.
Metodi matematici dell’astronomia
Ordine dei metodi iterativi
I metodi iterativi a un punto sono di solito del prim’ordine. Se però il punto
fisso è 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.
Sviluppando (x) in serie di Taylor si ha
da cui si deduce
 n 1  ( p ) ( x )
lim

,
n   p
p
!
n
xn1   ( xn )   ( x ) 
 ( p )(  )
p!
( xn  x ) p
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.
Metodi matematici dell’astronomia
Il second’ordine, nel metodo di N-R, è pagato dal fatto che a ogni passo
bisogna 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
~
xn 1   ( ~
xn )   n , n  0,1,2,..
la successione dei valori “calcolati” col
metodo iterativo (inficiati quindi da errore di valutazione della funzione
e dall’errore di arrotondamento). Sottraendo m. a m. l’eguaglianza
x  ( x )
media,
si ottiene
~
xn1  x   ( ~
xn )   n   ( x ),
~
xn1  x  ' (  n )xn  x    n ,  I ( xn , x ),
da cui, usando il teorema della
che porta, dopo sottraz. ad ambo
Metodi matematici dell’astronomia
i membri della quantità ' ( ~xn )~xn1  x , a 1  ' (  n )( ~xn1  x )  ' (  n )( ~xn  ~xn1 )   n .
Se si ha |(n)|m<1 e |n|< n, allora ne segue:
' (  n ) ~ ~
n

m ~ ~
~
xn1  x 
xn  xn1 

xn  xn1  n .
1  ' (  n )
1  ' (  n ) 1  m
1 m
In questa maggiorazione dell’errore all’ n+1-esimo passo il primo
termine dell’ultimo m. stima l’errore di troncamento, il secondo quello di
calcolo delle funzioni implicate. E’ chiaro che per n grande il primo
termine è trascurabile rispetto al secondo che quindi è da un certo punto
dominante e impedisce una diminuzione d’errore pur al diminuire dell’errore di troncamento.
Metodi matematici dell’astronomia
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
indipendente dal metodo
xn   
f ( xn )
m
, f ' ( x )  m, x  J .
L’accuratezza entro cui si può determinare  è limitata dall’errore nella
valutazione di f(x); se infatti ~f ( x )  f ( x )   ( x ), ( x )   , allora, poichè
n
al meglio si ha
~
f ( xn )  0,
n
n
n
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:

x      , dove = / f (). Quindi il miglior limite d’errore per
m
n

Metodi matematici dell’astronomia
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.
Metodi matematici dell’astronomia
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 quell
dipendenti le coordinate spaziali (x,y,z). Il problema differenziale completo
comprende anche l’insieme delle condizioni iniziali su posizioni e
velocità (le equazioni sono infatti del 20 ordine, poichè riguardano le
accelerazioni). I sistemi di eq. differenziali della meccanica costituiscono
una classe di problemi detti, appunto, “ai valori iniziali”, per distinguerli
da un’altra classe di problemi differenziali, che sono quelli “ai limiti” o
“ai valori al contorno”. Un esempio tipico, in astronomia, di quest’ultima
categoria di problemi differenziali è il sistema d’equazioni che regolano
la struttura stellare in equilibrio sferico, viste in una delle prime lezioni.
Metodi matematici dell’astronomia
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 è:
 y'  f ( x , y )

 y( 0 )  y0
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.
Metodi matematici dell’astronomia
Famiglie di soluzioni dell’eq. diff. y=f(x,y) al variare della condizione
iniziale y(0); notare l’evoluzione dell’errore globale, che può portare
la soluzione relativa a y(0)=c1 verso quella relativa a y(0)=c2.
Metodi matematici dell’astronomia
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:
u1'  u2
 '
u2  u3
     

u'n 1  un

 F ( x ,u1 ,u2 ,...,un )  0
u ( 0 )  y
0
 1
u2 ( 0 )  y0'

     
u ( 0 )  y ( n 1 )
0
 n
Per avere una trattazione numerica dell’eq. diff., cioè una soluzione
approssimata, bisogna anche esplicitare la n-esima variabile un come
un=f(x,u1,u2,...,un-1), per il motivo che vedremo fra poco.
Metodi matematici dell’astronomia
Ci sono molti metodi per risolvere numericamente eq. diff. e sistemi di
eq. diff. ordinarie. Tratteremo qui metodi per la soluzione di problemi ai
valori 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 uno
sviluppo 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).
Metodi matematici dell’astronomia
Il metodo di E. in avanti è esplicito, a un sol punto (one step) e localmente 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 lo
svantaggio 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.
Metodi matematici dell’astronomia
Un semplice miglioramento del metodo di Eulero si ottiene utilizzando
2 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 la
conoscenza di 2 valutazioni precedenti, yn-1 e yn, il metodo ha bisogno di
una valutazione indipendente di y1, oltre alla condizione iniziale y0, per
essere utilizzato. Normalmente tale inizializzazione avviene con un passo
Metodi matematici dell’astronomia
di Eulero. Il leap-frog mostra “debole instabilità” come si può vedere
con l’esempio dell’equazione
 y'   y

 y( 0 )  1
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 il
leap-frog fino a x+H=x+Nh (N pari; H “macropasso”). Il valore yN
così ottenuto si corregge così:
ŷ N 
1
 y N  y N 1  hf ( xN , y N )
2
Per poi proseguire per altri N passi col leap-frog inizializzato da ŷ N
e ŷN 1  ŷN 1  hf ( xN , ŷN ). Poichè l’errore è y(x;h)-y(x)=c1(x)h2+c2(x)h4+···
è possibile effettuare l’usuale estrapolazione di Richardson con pk=2k.
Metodi matematici dell’astronomia
Soluzione di y´=-y, y(0)=1 col leap-frog con h=0.1 e inizializzato con 2
metodi 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à.
Metodi matematici dell’astronomia
Soluzione di y´=-y, y(0)=1 con vari metodil leap-frog con h=0.1
Metodi matematici dell’astronomia
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 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
N
p
h
h
 t   y( p )(  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)
1
yn 1  yn  h f ( xn , yn )  f ( xn 1 , y 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
 yn  hf ( xn1 , yn1 )
(1)
y ( 1 ) n 1  y ( 2 ) n 1
.
per cui: yn 1 
2
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
y( x  h )  y( x )   f ( x , y( x ))dx; L’idea dei metodi di R-K consiste nell’
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 (Adams-Bashforth)
Metodi matematici dell’astronomia
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 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)
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, ottenendo
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). Vediamone 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
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 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ò 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 per
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 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
ln  y( x  h )  yn( 21) 
yn( 21)  yn( 11)
2 p 1
h
 cn  
2
p 1
,
si può fermare l’operazione di dimezzamento del passo quando si ha
ln<h.