ODE

ODE

Metodi Numerici e Calcolo Scientifico

Ultimo aggiornamento: September 30, 2024

1
ODE

Argomenti trattati

1 Introduzione

2 Metodi numerici
Consistenza, 0-stabilità, convergenza
Metodi one-step

3 Assoluta stabilità

2
ODE
Introduzione

Molti fenomeni reali sono modellizati da (sistemi di) equazioni differenziali

ordinarie (ODE), i.e. equazioni in cui sono correlate una funzione incognita e le
sue derivate di ordine fino a m

y (m) (t) = f (t, y (t), y 0 (t), ..., y (m−1) (t))

(ODE di ordine m in forma normale)

Sistema di ODE di ordine 1:

 0

 y1 (t) = f1 (t, y1 (t), y2 (t), . . . , , ym (t))
y20 (t) = f2 (t, y1 (t), y2 (t), . . . , , ym (t))


..



 .
 0
ym (t) = fm (t, y1 (t), y2 (t), . . . , , ym (t))

3
ODE
Introduzione

Esempio: Modello di Malthus

Modello per la dinamica di una popolazione isolata.
p(t): numero di individui di una popolazione isolata al tempo t

p(t + ∆t) = p(t) + Np(t)∆t − Mp(t)∆t

N: tasso di natalità (numero di nati per unità di tempo e unità di popolazione)
M: tasso di mortalità

p(t + ∆t) − p(t)

= (N − M)p(t) = ap(t)
∆t
0
∆t → 0 =⇒ p (t) = ap(t)

p 0 (t) = ap(t) + I
I : immigrazioni

4
ODE
Introduzione

Esempio: Modelli Lotka Volterra

Modelli per la dinamica di due popolazioni che si influenzano a vicenda
 0
p1 (t) = a11 p1 (t) + a12 p1 (t)p2 (t)
p20 (t) = a21 p1 (t)p2 (t) + a22 p2 (t)

a11 , a22 tassi natalità mortalità delle due popolazioni presi singolarmente
(termini malthusiani)
a21 , a12 coefficienti di interazione
Sono possibili 3 classi di modelli di interazione

5
ODE
Introduzione

Esempio: (Segue)
1 Cooperazione: a11 , a22 < 0 (nessuna delle due popolazioni sopravvive se
presa da sola) e a21 , a12 > 0 (ogni popolazione trae beneficio dall’altra)
Es. maschi/femmine di una stessa specie
2 Competizione: a11 , a22 > 0 (ogni popolazione prospera se presa da sola) e
a21 , a12 < 0 (competizione fra popolazione danneggia entrambe)
Es. due popolazioni che si nutrono della stessa cosa
3 Preda-Predatore (es. p1 preda, p2 predatore): a11 > 0 (p1 si accrescerebbe
spontaneamente), a12 < 0 (la presenza di p2 danneggia p1), a22 < 0 (p2 si
estinguerebbe, se presa sola) a21 > 0 (grazie a p1, p2 si accresce)

6
ODE
Introduzione

Equazione di ordine m e sistemi di m equazioni di ordine 1 sono strettamente

collegati.

Osservazione
Con una sostituzione standard, è sempre possibile ricondurre una equazione di
ordine m ad un sistema di m equazioni del primo ordine: posto
z1 (t) = y (t), z2 (t) = y 0 (t), . . . , zm (t) = y (m−1) (t)

y (m) (t) = f (t, y (t), y 0 (t), ..., y (m−1) (t))

m

0

 z1 (t) = z2 (t)

 0
 z (t) = z3 (t)
 2


..
 .
0

zm−1 (t) = zm (t)




 0

zm (t) = f (t, z1 (t), z2 (t), ..., zm (t))

7
ODE
Introduzione

Un problema non è ben posto se non si assegna un dato ”iniziale”

Esempio:

y 0 (t) = y (t)
ha infinite soluzioni della forma y (t) = cet ∀c ∈ R

Definizione (Problema di Cauchy)

(
y 0 (t) = f (t, y (t)), t ∈ [a, b]
y (a) = y0
f : R × Rm 7→ Rm , y0 ∈ Rm

Sotto opportune ipotesi di regolarità su f , si dimostra esistenza e unicità della

soluzione del Problema di Cauchy

8
ODE
Introduzione

Teorema
Sia f : S = [α, β] × Rm 7→ Rm continua in S. Se ∃L > 0 t.c.

kf (t, y1 ) − f (t, y2 )k ≤ Lky1 − y2 k

∀t ∈ [α, β] e ∀y1 , y2 ∈ Rm allora per ogni a ∈ [α, β] e ∀y0 esiste una e una sola
soluzione del problema di Cauchy su tutto [α, β]

Osservazione

kf (t, y1 ) − f (t, y2 )k ≤ Lky1 − y2 k

si chiama condizione di Lipschitz (rispetto alla seconda variabile).

∂fi
La condizione è certamente soddisfatta se le derivate parziali ∂yj
esistono e
sono continue e limitate (in S).

9
ODE
Introduzione

Esempio: Non unicità della soluzione

(
y 0 (t) = |y (t)|, t ∈ [0, ∞)
p

y (0) = 0
p
f (t, y (t)) = |y (t)| non è Lipschitziana in u = 0, infatti
(
1
∂f (t, y ) 1 √y y >0
= 1
∂y 2 − √−y y <0

che è illimitata in y = 0.
Il Problema di Cauchy ammette infinite soluzioni: y (t) ≡ 0 ma anche
(
0 0≤t<c
y (t) = 1 2
4
(t − c) t ≥c

∀c ≥ 0

10
ODE
Introduzione

11
ODE
Metodi numerici

Metodi Numerici

Introduciamo una griglia di punti sull’intervallo [a, b] (nodi)

a ≡ t0 < t1 < t2 < . . . < tn ≡ b

Per semplicità, consideriamo nodi equidistanti:
b−a
h= = tk+1 − tk ∀k
n
Metodo numerico: formula che calcola l’approssimazione yk ' y (tk ) per
k = 1, ..., n.

12
ODE
Metodi numerici

Esempio: Metodo di Eulero esplicito

Consideriamo l’ODE valutata in un generico nodo tk della griglia:

y 0 (tk ) = f (tk , y (tk ))

Ricordiamo la formula alle differenze finite in avanti per approssimare la

derivata in tk , con h > 0:
y (tk + h) − y (tk ) df (t, y (t))
y 0 (tk ) = − |t=ξ · h, ξ ∈ [tk , tk + h]
h dt

y (tk + h) − y (tk ) df (ξ, y (ξ))

− h = f (tk , y (tk ))
h dt
df (ξ, y (ξ)) 2
h
y (tk+1 ) = y (tk ) + hf (tk , y (tk )) +
dt
Trascurando l’errore, si ottiene lo schema numerico, con yk ' y (tk ):

yk+1 = yk + hf (tk , yk )

Daremo nel seguito una definizione formale degli errori colorati

13
ODE
Metodi numerici

Esercizio proposto
Dato il problema di Cauchy
(
y 0 = −y + y 2 , t ∈ [1, 10]
y (1) = 3

effettuare due passi di integrazione con il metodo di Eulero esplicito con

h = 1/2.

Esercizio proposto
Ripetere con 
00 2
y = −y + y ,
 t ∈ [1, 10]
y (1) = 3
 0

y (1) = 2

14
ODE
Metodi numerici
Consistenza, 0-stabilità, convergenza

Errori associati ai metodi numerici

Ci concentriamo sul metodo di Eulero esplicito, la generalizzazione è

immediata. Siano
yk soluzione numerica nel nodo tk ∀k
y (tk ) soluzione esatta nel nodo tk ∀k
ỹk+1 la soluzione numerica nel nodo tk+1 supponendo di partire dal dato
esatto y (tk ) anzichè dalla soluzione numerica yk :

ỹk+1 = y (tk ) + hf (tk , y (tk ))

15
ODE
Metodi numerici
Consistenza, 0-stabilità, convergenza

Definizione (Errore locale di troncamento τ (h))

È l’errore commesso in un passo, supponendo di partire dalla soluzione esatta:

τ (h) = y (tk+1 ) − ỹk+1 = y (tk+1 ) − y (tk ) − hf (tk , y (tk ))

Nell’esempio su EE, è la parte evidenziata in rosso. Per EE si ha τ (h) = O(h2 ).

Definizione (Errore locale di discretizzazione d(h))

È l’errore introdotto in un passo nella discretizzazione della derivata, ed è
τ (h)
d(h) =
h
Nell’esempio su EE, è la parte evidenziata in blu. Per EE si ha d(h) = O(h).

Definizione (Errore globale ek+1 )

È l’errore complessivo commesso in k passi di integrazione,
ek+1 = y (tk+1 ) − yk+1 e può essere scomposto come

ek+1 = y (tk+1 ) − yk+1 = y (tk+1 ) − ỹk+1 + ỹk+1 − yk+1

| {z } | {z }
err .tronc. propagazione
16
ODE
Metodi numerici
Consistenza, 0-stabilità, convergenza

17
ODE
Metodi numerici
Consistenza, 0-stabilità, convergenza

Definizione (Consistenza)
Un metodo numerico si dice consistente se limh→0 d(h) = 0

Definizione (Ordine (di consistenza))

Un metodo numerico si dice di ordine p se d(h) = O(hp )

Definizione (Convergenza)
Un metodo numerico è convergente in un punto t ∈ [a, b] se data una
suddivisione di [a, t] in N intervalli di ampiezza h = t−a
N
si ha

lim yN = y (t)
N→∞

ovvero se l’errore globale eN tende a 0. Il metodo si dirà convergente in [a, b]

se è convergente ∀t ∈ [a, b].

18
ODE
Metodi numerici
Consistenza, 0-stabilità, convergenza

Osservazione
La sola consistenza non è sufficiente per la convergenza, a causa del termine di
propagazione degli errori. Affinchè un metodo numerico sia convergente
occorre che sia consistente e che garantisca la non propagazione degli errori.

Definizione (0-stabilità)
Un metodo numerico è 0-stabile se ∃K > 0, h̄ tali che, dati due valori iniziali
y0 , ŷ0 , le corrispondenti soluzioni yk , ŷk soddisfano, per h ≤ h̄,
b−a
|yk − ŷk | ≤ K |y0 − ŷ0 | ∀k ≤
h

Teorema (Lax-Richtmeyer)
consistenza+0-stabilità = convergenza

19
ODE
Metodi numerici
Metodi one-step

Metodi one-step

yk+1 = yk + hΦ(tk , yk ; h) espliciti

yk+1 = yk + hΦ(tk , yk , yk+1 ; h) impliciti

Metodi multi-step

yk+1 = yk + hΦ(tk , yk−p , ..., yk−1 , yk ; h) espliciti

yk+1 = yk + hΦ(tk , yk−p , ..., yk−1 , yk , yk+1 ; h) impliciti

20
ODE
Metodi numerici
Metodi one-step

Alcuni metodi one-step

1 Eulero esplicito yk+1 = yk + hf (tk , yk )

2 Eulero implicito yk+1 = yk + hf (tk+1 , yk+1 )
3 Trapezi yk+1 = yk + h2 (f (tk , yk ) + f (tk+1 , yk+1 ))
4 Heun yk+1 = yk + h2 (f (tk , yk ) + f (tk+1 , yk + hf (tk , yk )))
5 Eulero modificato yk+1 = yk + hf (tk + h2 , yk + h2 f (tk , yk ))

Osservazione
I metodi one-step sono 0-stabili, quindi se sono consistenti sono convergenti.

21
ODE
Metodi numerici
Metodi one-step

Approfondimento: una prospettiva diversa

Z tk+1 Z tk+1
y 0 (t) = f (t, y (t)) y 0 (t) dt = f (t, y (t)) dt
tk tk
Z tk+1
y (tk+1 ) − y (tk ) = f (t, y (t)) dt
tk

Applico una formula di quadratura all’integrale a secondo membro: Rettangolo:

y (tk+1 ) − y (tk ) ' hf (tk , y (tk ))

y (tk+1 ) − y (tk ) ' hf (tk+1 , y (tk+1 ))

Trapezi:
h
y (tk+1 ) − y (tk ) ' (f (tk , y (tk )) + f (tk+1 , y (tk+1 )))
2

22
ODE
Metodi numerici
Metodi one-step

Metodi di Runge-Kutta

Esiste una famiglia di metodi detta Metodi di Runge Kutta a s stadi.

Metodi ad uno stadio espliciti:

(
yk+1 = yk + ha1 κ1
κ1 = f (tk + b1 h, yk )

Esempio: Eulero esplicito

yk+1 = yk + hf (tk , yk )

a1 = 1, b1 = 0

23
ODE
Metodi numerici
Metodi one-step

Metodi di Runge-Kutta

Metodi a due stadi espliciti (s=2):


yk+1 = yk + h(a1 κ1 + a2 κ2 )

κ1 = f (tk + b1 h, yk )

κ = f (t + b h, y + hc κ )
2 k 2 k 21 1

Esempio: Heun
 
1 1
yk+1 = yk + h f (tk , yk ) + f (tk + h, yk + hf (tk , yk ))
2 2
1
a1 = a2 = , b1 = 0, b2 = 1, c21 = 1
2
Esempio: Eulero modificato
1 1
yk+1 = yk + hf (tk + h, yk + hf (tk , yk ))
2 2
1 1
a1 = 0, a2 = 1, b1 = 0, b2 = , c21 =
2 2

24
ODE
Metodi numerici
Metodi one-step

Metodi di Runge-Kutta

Più in generale, metodi a s stadi espliciti:

s
X
yk+1 = yk + h a i κi
i=1

i−1
X
κi = f (tk + bi h, yk + h cij κj ) i = 1, ..., s
j=1

a ∈ Rs , b ∈ Rs , c ∈ Rs×s
# stadi = # valutazioni di funzione

25
ODE
Metodi numerici
Metodi one-step

Tableau di Butcher

b1
b2 c21
b ∈ Rs c ∈ Rs×s .. ..
a T ∈ Rs . .
bs cs1 ... cs,s−1
a1 ... as
s
X s
X
ai = 1, bi = cij ∀i = 1, .., s
i=1 j=1

26
ODE
Metodi numerici
Metodi one-step

Ordine di alcuni metodi:

1 Eulero I/E: p = 1
2 Trapezi: p = 2
3 Heun: p = 2
4 Runge-Kutta: variabile, ma p ≤ s:
p max = s, s = 1, 2, 3, 4
p max = s − 1, s = 5, 6, 7
p max = s − 2, s = 8, 9
p max ≤ s − 3, s ≥ 10

27
ODE
Metodi numerici
Metodi one-step

Esempio: Confronto sugli ordini

(
y0 = y t ∈ (0, T ]
y (0) = 1
Soluzione esatta: y (t) = exp(t)
h |yNEE − y (tN )| |yNH − y (tN )| |yNRK4 − y (tN )|
1
2
4.6828e-01 7.7657e-02 9.3564e-04
1
4
2.7688e-01 (1.7) 2.3426e-02 (3.3) 7.1889e-05 (13.0)
1
8
1.5250e-01 (1.8) 6.4406e-03 (3.6) 4.9840e-06 (14.4)
1
16
8.0353e-02 (1.9) 1.6883e-03 (3.8) 3.2812e-07 (15.2)
1
32
4.1292e-02 (1.9) 4.3215e-04 (3.9) 2.1048e-08 (15.6)
1
64
2.0937e-02 (2.0) 1.0932e-04 (3.9) 1.3327e-09 (15.8)

28
ODE
Assoluta stabilità

Assoluta stabilità

Convergenza garantisce errore piccolo se h → 0.

Ma quando applico nella pratica un metodo, h è ben diverso da 0!
Tutti i valori di h sono accettabili?
Esempio:
(
y 0 = −10y t ∈ (0, T ]
y (0) = 1
Soluzione esatta: y (t) = exp(−10t)
Eulero esplicito:
yk+1 = yk + h(−10yk ) = (1 − 10h)yk

29
ODE
Assoluta stabilità

Soluzione esatta e−10t e soluzione numerica, Eulero esplicito, h = 0.15

30
ODE
Assoluta stabilità

Soluzione esatta e−10t e soluzione numerica, Eulero esplicito, h = 0.3

31
ODE
Assoluta stabilità

Soluzione esatta e−10t e soluzione numerica, Eulero esplicito, h = 0.5

32
ODE
Assoluta stabilità

Soluzione esatta e−10t e soluzione numerica, Eulero implicito, h = 0.15

33
ODE
Assoluta stabilità

Soluzione esatta e−10t e soluzione numerica, Eulero implicito, h = 0.3

34
ODE
Assoluta stabilità

Soluzione esatta e−10t e soluzione numerica, Eulero implicito, h = 0.5

35
ODE
Assoluta stabilità

Spiegazione qualitativa

Con il metodo di EE i risultati numerici sono accettabili solo se h è

inferiore a una soglia ĥ. Il metodo è allora stabile perché gli errori prodotti
a ogni passo di integrazione non si propagano.
Se h è troppo grande, i risultati non sono accettabili. Il metodo è instabile
perché gli errori prodotti si amplificano passo dopo passo.
La situazione non è limitata solo al problema/metodo specifico:
in generale, per i metodi espliciti è presente una soglia ĥ, che dipende sia
dal metodo usato che dal problema, e solo per h sufficientemente piccolo il
metodo è stabile.
Metodi impliciti come EI e Trapezi non hanno queste restrizioni su h ed è
quindi possibile prendere passi h più grandi. Lo svantaggio di tali metodi è
che è necessario risolvere equazioni (tipicamente nonlineari) ad ogni passo.

36
ODE
Assoluta stabilità

Regione di assoluta stabilità

Formalizziamo, partendo dal caso di una sola equazione.

Problema test o problema modello:
(
y 0 = λy , <λ < 0
y (t0 ) = y0

Soluzione esatta: y (t) = y0 exp(λ(t − t0 ))

lim y (t) = 0
t→∞

(problema detto asintoticamente stabile)

Per h fissato, il metodo numerico riproduce lo stesso comportamento
quando il numero di passi tende a ∞, i.e.

lim yk = 0 ?
k→∞

(metodo detto (assolutamente) stabile)

37
ODE
Assoluta stabilità

Cerchiamo di stabilire se un metodo numerico fissato è assolutamente stabile

per un dato problema.
Applicando un metodo numerico one-step al problema test y 0 = λy , i.e. con

f (t, y ) = λy

ci riconduciamo ad una espressione del tipo

yk+1 = F(hλ)yk

(si ricordi l’esempio precedente: λ = −10 e yk+1 = (1 − 10h)yk ).

38
ODE
Assoluta stabilità

Esempio: Eulero esplicito

yk+1 = yk + hf (tk , yk )

yk+1 = yk + hλyk = (1 + hλ)yk

F(hλ) = 1 + hλ

39
ODE
Assoluta stabilità

Esempio: Eulero implicito

yk+1 = yk + hf (tk+1 , yk+1 )

yk+1 = yk + hλyk+1

1
yk+1 = yk
1 − hλ
1
F(hλ) =
1 − hλ

40
ODE
Assoluta stabilità

Esempio: Heun
h
yk+1 = yk + (f (tk , yk ) + f (tk+1 , yk + hf (tk , yk )))
2
h
yk+1 = yk + (λyk + λ(yk + hλyk ))
2
1
yk+1 = yk + hλyk + (hλ)2 yk
2
1
F(hλ) = 1 + hλ + (hλ)2
2

Esempio: RK esplicito a s stadi di ordine s

1
F(hλ) = 1 + hλ + . . . + (hλ)s
s!

41
ODE
Assoluta stabilità

Allora,
yk+1 = F(hλ)yk = F(hλ)2 yk−1 = . . . = F(hλ)k+1 y0

lim yk = 0 ⇔ |F (hλ)| < 1

k→∞

Regione di assoluta stabilità del metodo numerico:

Ra = {hλ ∈ C : |F(hλ)| < 1}

42
ODE
Assoluta stabilità

Esempi

1 Eulero esplicito |1 + hλ| < 1

1
2 Eulero implicito | 1−hλ |<1
3 Heun |1 + hλ + 21 (hλ)2 | < 1

43
ODE
Assoluta stabilità

Stabilità di sistemi

Nel caso di un sistema di ODE lineari a coefficienti costanti

y 0 (t) = Ay (t)

gli autovalori di A sono i responsabili della stabilità.

Osservazione
Se A è diagonalizzabile e λi sono i suoi autovalori con autovettori vi , per
i = 1, ..., m si ha

y (t) = c1 eλ1 (t−t0 ) v1 + . . . + cm eλm (t−t0 ) vm

Se <λi < 0 per i = 1, ..., m il problema è quindi asintoticamente stabile.

44
ODE
Assoluta stabilità

Affinché il metodo sia assolutamente stabile occorre che ∀λi si abbia hλi ∈ Ra
(in pratica l’autovalore più negativo fornisce la restrizione sul passo)

Eventuali termini sorgenti g (t)

y 0 (t) = Ay (t) + g (t)

non hanno alcun ruolo nella stabilità del problema.

45
ODE
Assoluta stabilità

Problemi stiff

Definizione
Un sistema

y 0 (t) = Ay (t) + g (t)

si dice stiff in I = [t0 , t0 + L] se
1 eventuali autovalori con parte reale positiva hanno <(λi )L non grande
2 esiste almeno un autovalore con parte reale negativa e <(λi )L << −1

maxi |<(λi )|L: grado di stiffness

46
ODE
Assoluta stabilità

A causa dell’autovalore del punto 2, se Ra piccola =⇒ h molto piccolo (rispetto

a L)

Esempio:
 0

 y1 (t) = y2 (t)
y 0 (t) = −100y (t) − 101y (t)

2 1 2
√
y1 (0) = 7



y2 (0) = π
in t ∈ [0, 100]  
0 1
A=
−100 −101
λ1 (A) = −100, λ2 (A) = −1
2
Se usiamo per es. EE h ≤ 100 = 0.02 almeno 5000 passi di integrazione per
coprire l’intervallo [0, 100].

47
ODE
Assoluta stabilità

Esempio:
 0
y1 (t) = y2 (t)

y 0 (t) = −10y (t) − 11y (t)

2 1 2
√


y1 (0) = 7

y2 (0) = π
in t ∈ [0, 10000]  
0 1
A=
−10 −11
λ1 (A) = −10, λ2 (A) = −1
2
Con EE h ≤ 10 = 0.2 almeno 50000 passi di integrazione per coprire
l’intervallo [0, 10000].

48
ODE
Assoluta stabilità

Ma quanto contribuisce veramente l’autovalore più piccolo sulla soluzione

complessiva?

Esempio:
 0
y1 (t) = y2 (t)

y 0 (t) = −1000y (t) − 1001y (t)

2 1 2


y1 (0) = 2
y2 (0) = −1

 
0 1
A=
−1000 −1001
λ1 (A) = −1, λ2 (A) = −1000
2 1
Con EE h < 1000
= 500

49
ODE
Assoluta stabilità

Osservazione

y (t) = c1 e−(t−t0 ) v1 + c2 e−1000(t−t0 ) v2

| {z }
transiente

Il contributo del termine corrispondente all’autovalore −1000 si “spenge”

subito, eppure continua fino alla fine a vincolare il passo!

50
ODE
Assoluta stabilità

51
ODE
Assoluta stabilità

Osservazione (Generalizzazione)
Se il problema non è lineare, o non è assolutamente stabile, è necessario
introdurre diversi concetti di stabilità, volti comunque a garantire che la
soluzione abbia lo stesso comportamento qualitativo della soluzione esatta.

52