Metodo di Calcolo Numerico per Equazioni differenziali Ordinarie
Approssimiamo la soluzione in:
[t 0 , T ]
Suddividiamo in N intervalli mediante i punti:
t 0 , t1 ,..., t N  T
Una volta conosciuta la soluzione in questi punti potremo approssimare la soluzione in
tutto l’intervallo (magari con una interpolazione polinomiale).
Per adesso considereremo una suddivisione uniforme:
tn  t0  n  h
T  t0
h
N
Considereremo un problema ai valori iniziali (di Cauchy):
y'  f (t , y)
y(t 0 )  y0
La f definisce una campo di pendenze.
Metodi a un passo
In generale, studiaremo strategie che rientrano nella categoria dei cosiddetti
metodi a un passo; il nome deriva dal fatto che per calcolare la soluzione
numerica al tempo tn+1 é sufficiente conoscere la soluzione numerica al tempo
n:
y n1  y n  h  Φt n ,y n ;h,f 
Funzione Incrementale
Caratterizza un metodo specifico. Rappresenta
una approssimazione numerica della media
della funzione f tra gli istanti n, n+1.
y n1  y n
 Φt n ,y n ;h,f   d (t n , h)
h
Errore locale di discretizzazione
(di troncamento)
Metodi a un passo
Consistenza: Un metodo a un passo si dice consistente nell’intervallo di
integrazione se d(t,h) é infinitesimo per h tendente a zero. Piú precisamente
esiste una funzione d(h) tale che:
d (t , h)  d (h)  lim d (h)  0
h0
Inoltre un metodo a un passo si dirá di ordine di consistenza p se:
d ( h)  O ( h p )
Convergenza: se la Φ é lipschitziana rispetto a y si ha convergenza.
Richiede anche la stabilitá….vedi propagazione dell’errore….
Metodi a un passo
dy (t ) y (t  h)  y (t )
dy (t )

 y (t  h)  y (t )  h 
dt
h
dt
1. Metodo di Eulero Esplicito
y n 1  y n  h  f (t n , y n )
n  0 ,....,N-1
1
y(t n )  y n  y''(ξ )h 2
2
Sviluppo di Taylor sino al primo ordine.
Metodi a un passo
dy (t ) y (t )  y (t  h)
dy (t )

 y (t )  y (t  h)  h 
dt
h
dt
2. Metodo di Eulero Implicito
y n1  y n  h  f (t n1 , y n1 )
n  0,....,N-1
Ha lo stesso grado di accuratezza del metodo esplicito peró richiede la
risoluzione di una equazione non lineare (comunque in alcuni casi i
metodi impliciti possono presentare dei vantaggi).
Metodi a un passo
Se arrestiamo lo sviluppo di Taylor al secondo ordine, e calcoliamo la
derivata seconda di y, si ottiene: (metodo del secondo ordine)
1 2
y(t 0 )  y 0  h  y'(y 0 )  h  y''(y 0 )  O(h 3 )
2
y '  f (t , y )
y' '  f t  f y y'  f t  f y f
3. Metodo basato sullo Sviluppo di Taylor arresto al 2 ordine:
y n 1

1 2  f(t n ,y n ) f(t n ,y n )
 y n  hf(t n , y n )  h 

f(t n , y n )
2  t
y

n  0,...,N-1
Metodi a un passo
Per costruire un metodo del secondo ordine (come il precedente) senza
dover calcolare e valutare le derivate di f, si puó ragionare in questo
modo: t n  t n 1  h
y p  yn
y p  yn
h
m
 f (t n , y n ) 
 y p  y n  f (t n , y n )
h/2
h/2
2
Coordinate di P:
h
x p  t n 1 / 2  t n 
2
h
y p  y n  f (t n , y n )
2
Soluzione
reale
yp
P
y(t n )
t n 1 / 2
tn
h/2
h
t n 1
Metodi a un passo
4. Metodo di Eulero Modificato:
h
h
y n 1  y n  h  f(t n  , y n  f (t n , y n ))
2
2
n  0,...,N-1
In pratica, abbiamo approssimato la soluzione reale nell’ intervallo t n , t n  h / 2 con la
retta tangente in tn per stimare il valore della funzione al tempo tn+h/2 (punto P); infine
abbiamo approssimato l’incremento della soluzione in t n , t n  h attraverso la pendenza
in P (si cerca di migliorare la stima della retta secante tra due instanti di integrazione).
Questo metodo ha una accuratezza del secondo ordine.
Metodi a un passo
Come abbiamo appena visto per aumentare l’accuratezza il metodo di Eulero
esplicito basta migliorare l’approssimazione della secante tra due instanti di
integrazione: un’altra idea potrebbe essere l’utilizzo di una media tra le
pendenze ai tempi tn e tn+h.
5. Metodo di Heun (differente del Metodo dei Trapezi):
1
y n 1  y n  h   f(t n , y n )  f(t n  h, y n  hf (t n , y n )) 
2
n  0,...,N-1
Per valutare la funzione al passo tn+h si é utilizzato una passo dell’Eulero
esplicito si é considarato un’incremento lineare con pendenza definita al tempo
tn. Differisce dal metodo dei Trapezi poiché quest’ultimo é implicito (vedi
dopo).
Metodi a un passo: Metodi Runge-Kutta
6. Metodo Runge-kutta a M stadi (livelli):
M
y n 1  y n  h    i k i n  0 ,...,N-1
i 1
P
k i  f (t n   i h, y n  h    ij k j )
j 1
Generalizzazione delle osservazioni
utilizzate per arrivare all’ Eulero
modificato e alla formula di Heun.
P  i 1
Metodo Esplicito
PM
Metodo Imsplicito
Pi
Metodo Semi-imsplicito
Metodi a un passo: Metodi Runge-Kutta
La denominazione di esplicito, implicito, o semi-implicito dipende dalla minore
o maggiore facilitá nel derivare i vari ki : in un caso si ricaveranno in cascata,
mentre nell’altro si dovrá risolver un sistema di equazioni. Possiamo giungere
ai metodi Runge-kutta in altro modo; infatti in generale possiamo utlizzare
varie formule di quadratura per ottenere metodi giá conosciuti:
h  f (t n , y n )
y n 1  y n  
t n 1
tn
f (t , y (t )) dt
Eulero esplicito
h  f (t n1 , y n1 )
Eulero implicito
h
 f (t n , y n )  f (t n1 , y n1 )
2
7. Metodo dei Trapezi
Dunque sfruttando delle formule di
quadratura che utilizzano i punti:
tn  i h
i  1,..., M
(ricavato dalla omonima formula
di quadratura)
Metodi a un passo: Metodi Runge-Kutta
Dunque sfruttando altre formule di quadratura
che utilizzano i nodi :
tn  i h
i  1,..., M
Possiamo ottenere formule del tipo:
M
y n 1  y n  h    i f (t n   i h,y (t n   i h))
i 1
che dipendono dai corrispondenti valori incogniti:
y(t n   i h)
Il problema viene risolto approssimando a loro volta
l’integrale seguente con una fomula di quadratura:
( y n  y (t n ))
Metodi a un passo: Metodi Runge-Kutta
La formula di quadratura in questo caso potrá utlizzare
tutti o solo alcuni nodi (diaciamo P):
y n 1 (t n   i h)  y (t n )  
t n  i h
tn
f (t , y (t )) dt
P
 h    ij f (t n   i h, y n 1 (t n   i h))
j 1
Chiaramente se P=i-1 sará un metodo implicito; inoltre P al massimo
sará uguale a M. Se imponiamo che la formula di quadratura sia
esatta almeno per funzioni costanti per ogni i troviamo la relazione:
P
 i    ij
i  1,..., M
j 1
Unendo i varii procendimenti di quadratura troviamo le formule di
Runge-Kutta:
Metodi a un passo: Metodi Runge-Kutta
M
y n 1  y n  h    i f (t n   i h, Yi ) n  0 ,...,N-1
i 1
P
Yi  y n  h    ij f (t n   i h, Y j )
j 1
Queste formule possono essere rappresentate in forma sintentica
con una tabella di coefficienti (consideriamo il caso piú generale
P=M):
………
………
Vettore delle
ascisse
 1 11…………… 1M
Matrice dei
coefficienti B
 M  M 1 …………… MM
1 …………… M
Vettore W dei Pesi
Metodi a un passo: Metodi Runge-Kutta
Se la matrice B é triangolare inferiore il metodo sará esplicito e gli Yi
si calcolano facilmente in cascata. In questo caso la condizione sugli
alfa impone: 1  0 . Se include anche la diagonale sará semiesplicito
e la soluzione sará ricorsiva, mentre con B piena il metodo si dice
implicito e richiede la risoluzione di un sistema non lineare.
0
0
Eulero Esplicito
A un livello esplicito
Eulero Implicito
A un livello implicito
1
1
1
1
0
1
0
0
1/2 1/2
1/2 1/2
Trapezi
A due livelli semi-implicito
Metodi a un passo: Metodi Runge-Kutta
0
0
0
1
1
0
Heun
1/2 1/2
0
0
0
1/2 1/2 0
1
0
Eulero Modificato
Metodi a un passo: Metodi Runge-Kutta
Il massimo ordine di accuratezza p(M) raggiungibile con un metodo a
M livelli varia in questo modo:
Metodi espliciti
M
p(M)
1
1
2
2
3
3
4
4
5
4
6
5
7
5
8
6
Metodi impliciti
p(M)=2M
Metodi a un passo: Metodi Runge-Kutta
Troviamo i coeff. per i metodi a due stadi espliciti imponendo un certo
grado di accuratezza:
y n 1  y n  h  1 k1   2 k 2 
k1  f (t n , y n )
k 2  f (t n  h, y n  h  k1 )
Sviluppando k2 secondo Taylor arrestando al primo ordine:
 f
 f
f 
f 
2
  O(h 2 )
k 2  f  h   
 k1   O(h )  f  h   
 f
y 
y 
 t
 t
y n 1
 f
f 
  O(h 3 )
 y n  h(1   2 ) f  h  2  
 f
y 
 t
2
Qui risulta al
secondo…
Metodi a un passo: Metodi Runge-Kutta
Confrontando la formula precedente con il reale sviluppo di Taylor
della formala arrestato al secondo ordine, troviamo le condizioni sui
coefficienti:
 f
f 

  O(h 3 )
y n 1  y n  h  f  h   f
y 
 t
 f
f 
2
  O(h 3 )
y n 1  y n  h(1   2 ) f  h  2 
 f
y 
 t
2
I coeff. dovranno compiere:
1   2  1
1
2
1
2   
2
2   
Queste condizioni sono rispettate da Eulero Modificato
e da Heun. Inoltre si vede α=β come avevamo giá
detto.
Analisi dell’errore: Metodi Runge-Kutta
L’errore di troncamento locale (errore locale di discretizzazione)
indica l’errore all’integrare un passo tra due instanti di tempo…ora
questo errore si propagherá al passo successivo sommandosi al
seguente errore di integrazione.
y n1  y n
 Φt n ,y n ;h,f   d (t n , h)
h
Errore locale di discretizzazione
E’ utile definire ora l’errore di troncamento locale al passo nesimo:
Valori esatti !!
 (t n , h)  y(t n1)  y(t n )  h  Φt n ,y(t n )
In pratica l’errore che commetto integrando dal passo n a n+1,
supponendo perfettamente conosciuta la soluzione al tempo n.
Analisi dell’errore: Metodi Runge-Kutta
E’ chiaro che all’errore di trocamento locale al passo n, dovremo
sommare l’errore commesso dalle precedenti integrazioni:
Soluzione esatta
En  y n  y(t n )  E PR  E LOC  E PR   (t n , h)
Errore accumulato totale
Errore di troncamento locale
Valore ottenuto
dalle varie
integrazioni
Errore di propagazione
Consistenza
Stabilitá
E PR  y n1   y(t n )  h  (t n , y(t n ))
Valore esatto
Analisi dell’errore: Metodi Runge-Kutta
In generale per definire un metodo convergente si richiede che
l’errore di troncamento locale tenda a zero al decrescere il passo di
integrazione h (consistenza) e che l’errore di propagazione non si
amplifichi passo dopo passo (stabilitá) .
Convergenza = Consistenza + Stabilitá
Se h é infinitesimo lo é anche
l’errore di troncamento
Gli errori non si amplificano al
propagarsi
Analisi dell’errore: Metodi Runge-Kutta
 Diciamo che la consistenza é una condizione statica, che suppone la
convergenza al diminuire il passo di integrazione h (ovvero che la nostra
approssimazione migliori con un passo h piú piccolo).
La stabilitá controlla la dinamica del nostro modello, in modo tale che
errori successivi non portino a approssimazioni assolutamente erronee.
Stabilitá: Metodi Runge-Kutta
In generale e’ difficile analizzare la stabilitá di un metodo, per questo
ci limitiamo a una classe particolare di equazioni differenziali test:
y '  f (t , y )    y
y ( 0)  y 0
    i
La cui soluzione generale é :
y (t )  y 0  et  (cos( t )  i sin( t ))
Noi conisidereremo ovviamente le soluzioni stabili con alfa<0.
Stabilitá: Metodi Runge-Kutta
Applichiamo ora per esempio il metodo di Eulero esplicito:
y n1  y n  h  f (t n , y n )  y n  h    y n  (1  h )  y n
Questa equazione alle differenze é stabile se:
1  h  1
Regione di assoluta stabilitá per il
metodo di Eulero esplicito: Il passo h
deve essere sufficientemente piccolo
(dato un λ).
hα
    i
-2
-1
hβ
Stabilitá: Metodi Runge-Kutta
Altre regioni di assolutá stabilitá:
1
1
1  h
1  h / 2
1
1  h / 2
Eulero Implicito
Trapezi
Stabilitá: Metodi Runge-Kutta
 Se le regioni di assoluta stabilitá contengono il semipiano α<0 allora il
metodo si dice incondizionatamente stabile o assolutamente stabile
poiché risulta stabile per tutti i λ della equazione test stabili, e per ogni
passo h.
 I Metodi impliciti risultano migliori se si analizza la stabilitá.
Pur essendo l’equazione test un caso particolare puó servire per
studiare almeno localmente equazioni piú generali. Infatti intorno a un
punto (tn,yn) possiamo linearizzare rispetto a y:
f (t , y )
y '  f (t , y ) 
y y
y
f (t , y )

y