Equazioni differenziali ordinarie

Alvise Sommariva

Università degli Studi di Padova Dipartimento di Matematica Pura e Applicata

16 gennaio 2012

Problema di Cauchy

Si consideri il problema di Cauchy

\{ \begin{matrix}y′(x) = f (x, y(x)), x ≥ x_{0} \\ y(x_{0}) = y_{0} \end{matrix}  (1)

dove f è a valori in ℝ^{n} e definita in un sottoinsieme Ω di ℝ × ℝ^{n}, con (x_{0}, y_{0}) ∈ Ω. Di seguito supporremo che tale problema abbia una sola soluzione.

Alvise Sommariva  Equazioni differenziali ordinarie.  2/ 28

Esempi

Esempio 1.

\{ \begin{matrix}y′(x) = y(x), x ≥ 0 \\ y(0) = 1 \end{matrix}  (2)

la cui soluzione è exp (x). E’ un problema di Cauchy con f (x, y(x)) = y(x).

Esempio 2.

\{ \begin{matrix}y′_{1}(x) = −y_{2}(x), x ≥ 0 \\ y′_{2}(x) = y_{1}(x), x ≥ 0 \\ y_{1}(0) = 1, y_{2}(0) = 0 \end{matrix}  (3)

la cui soluzione è y(x) = (y_{1}(x), y_{2}(x)) = (cos (x), sin (x)). E’ un problema di Cauchy con y = (y_{1}, y_{2}), x_{0} = 0, y(x_{0}) = (1, 0) e

f (x, y(x)) = f (x, y_{1}(x), y_{2}(x)) = (−y_{2}(x), y_{1}(x)).

Teorema di Cauchy in piccolo

A tal proposito supponiamo che per qualche δ, ϵ > 0 si abbia

S_{δ,ϵ,x_{0},y_{0}} ≡ [x_{0} − δ, x_{0} + δ] × \overline{B}_{ϵ}(y_{0}) ⊆ Ω,

e che la funzione sia L-lipschitziana in S_{δ,ϵ,x_{0},y_{0}} cioè esista L > 0 tale che

|f (x, z_{1})−f (x, z_{2})| ≤ L|z_{1}−z_{2}|, ∀x ∈ [x_{0}−δ, x_{0}+δ], ∀z_{1}, z_{2} ∈ \overline{B}_{ϵ}(y_{0}).

Vale il seguente teorema di esistenza ed unicità (in piccolo)

Teorema. Siano Ω un aperto di ℝ × ℝ^{n}, f una funzione di Ω in ℝ^{n} e (x_{0}, y_{0}) un punto di Ω. Se f è continua e valgono le condizioni di L-lipschitzianità in S_{δ,ϵ,x_{0},y_{0}} (per qualche δ, ϵ > 0) allora esiste un intervallo aperto contenente x_{0} nel quale è definita una e una sola soluzione del problema di Cauchy (9).

Teorema di Cauchy in grande

Accanto a tale teorema esiste quello di esistenza ed unicità (in grande)

Teorema. Siano Ω = [x_{0}, x_{0} + a] × ℝ^{n}, f una funzione di Ω in ℝ^{n}. Se f è continua e valgono le condizioni di L-lipschitzianità in Ω allora il problema di Cauchy (9) ha una e una sola soluzione in [x_{0}, x_{0} + a].

Metodo di Eulero esplicito

Indichiamo con I(x, \overline{x}) il più piccolo intv. aperto contenente x e \overline{x}. Assumendo che la soluzione sia sufficientemente regolare, abbiamo dalla formula di Taylor per x ≈ \overline{x}, ξ ∈ I(x, \overline{x}) e (9)

y(x) = y(\overline{x}) + y′(\overline{x})(x − \overline{x}) + \frac{y′′(ξ)(x − \overline{x})^{2}}{2}
≈ y(\overline{x}) + y′(\overline{x})(x − \overline{x}) = y(\overline{x}) + f (\overline{x}, y(\overline{x}))(x − \overline{x})

Di conseguenza se si desidera calcolare la soluzione nei punti x_{k+1} = x_{0} + kh = x_{k} + h con k > 0, ponendo x = x_{k+1}, \overline{x} = x_{k}

y(x_{k+1}) ≈ y(x_{k}) + h · f (x_{k}, y(x_{k})).  (4)

Il metodo di Eulero esplicito consiste nell’approssimare y(x_{k+1}) con y_{k+1} definito da

y_{k+1} = y_{k} + h · f (x_{k}, y_{k})  (5)

ove y_{0} = y(x_{0}).

Alvise Sommariva  Equazioni differenziali ordinarie.  6/ 28

Metodo di Eulero implicito

Similmente al caso di Eulero esplicito, se poniamo invece x = x_{k}, \overline{x} = x_{k+1} abbiamo

y(x_{k}) ≈ y(x_{k+1}) + f (x_{k+1}, y(x_{k+1}))(x_{k} − x_{k+1}), (6)

e quindi

y(x_{k+1}) ≈ y(x_{k}) + h · f (x_{k+1}, y(x_{k+1})).  (7)

Il metodo di Eulero implicito consiste nell’approssimare y(x_{k+1}) con y_{k+1} definito da

y_{k+1} = y_{k} + h · f (x_{k+1}, y_{k+1}),  (8)

ove y_{0} = y(x_{0}).

Evidentemente ad ogni iterazione si richiede di risolvere un’eqz. nonlineare nella variabile z del tipo z = y_{k} + h · f (x_{k+1}, z), la cui soluzione può essere calcolata utilizzando ad esempio col metodo di Newton o del punto fisso o delle secanti.

Alvise Sommariva  Equazioni differenziali ordinarie.  7/ 28

Osservazione

Il problema di Cauchy è definito da

\{ \begin{matrix}y′(x) = f (x, y(x)), x ≥ x_{0} \\ y(x_{0}) = y_{0} \end{matrix}  (9)

dove f è a valori in ℝ^{n} e definita in un sottoinsieme Ω di ℝ × ℝ^{n}, con (x_{0}, y_{0}) ∈ Ω. Quindi a priori i metodi di Eulero esplicito

y_{k+1} = y_{k} + h · f (x_{k}, y_{k}), y_{0} = y(x_{0})  (10)

e Eulero implicito

y_{k+1} = y_{k} + h · f (x_{k+1}, y_{k+1}), y_{0} = y(x_{0})  (11)

possono essere utilizzati per sistemi di equazioni differenziali, in quanto la sequenza è ben definita.

Alvise Sommariva  Equazioni differenziali ordinarie.  8/ 28

Analisi convergenza

Supponiamo di analizzare il problema di Cauchy nell’intervallo compatto I = [x_{0}, xfin]. Sia y la soluzione esatta di un fissato problema di Cauchy (9) e u^{(h)} l’approssimazione fornita da un metodo numerico, campionando la soluzione nei punti x_{s} = x_{0} + sh, con Nh = x_{fin} − x_{0}.

Un tale metodo si dice convergente se

∀n, ∥y − u^{(h)}∥_{∞} ≤ C(h)

con C(h) infinitesimo rispetto ad h quando h tende a 0. Se

C(h) = O(h^{p})

per qualche p > 0 allora si dice che il metodo converge con ordine p.

Analisi convergenza, I

Consideriamo il metodo di Eulero esplicito, con ascisse equispaziate per un prefissato passo h. Sia y_{n} = y(x_{n}), con y sol. del problema di Cauchy, e

e_{n} = y_{n} − u_{n} = (y_{n} − u_{n}^{∗}) + (u_{n}^{∗} − u_{n})  (12)

dove

u_{n}^{∗} = y_{n−1} + h · f (x_{n−1}, y_{n−1}),
u_{n} = u_{n−1} + h · f (x_{n−1}, u_{n−1}).

Il metodo converge se entrambi i termini a secondo membro di (12) convergono a 0 per h → 0. Osserviamo che la prima sequenza parte dal valore assunto dalla soluzione in x_{n}, mentre la seconda da valori che a priori non coincidono con y(x_{n}).

Analisi convergenza, II

Se la derivata seconda di y esiste ed è continua allora si dimostra che per quanto concerne la sequenza di Eulero esplicito,

y_{n} − u_{n}^{∗} = (h^{2}/2) · y′′(ξ_{n})

per qualche ξ_{n} ∈ (t_{n−1}, t_{n}). La quantità

τ_{n}(h) = (y_{n} − u_{n}^{∗})/h

si chiama errore locale di troncamento del metodo, mentre

τ(h) = max_{n=0,...}|τ_{n}(h)|

si chiama errore globale di troncamento del metodo.

Alvise Sommariva  Equazioni differenziali ordinarie.  11/ 28

Analisi convergenza, III

Sotto queste condizioni, per il metodo di Eulero espl., posto M = max_{x∈I} |y′′(x)| abbiamo τ(h) → 0 per h → 0 in quanto 0 ≤ τ(h) ≤ Mh/2.

In generale, un metodo per cui τ(h) → 0 per h → 0 si dice consistente. Da (12) abbiamo

|e_{n}| = |y_{n} − u_{n}| ≤ |y_{n} − u_{n}^{∗}| + |u_{n}^{∗} − u_{n}| ≤ h|τ(h)| + |u_{n}^{∗} − u_{n}| (13)

Consideriamo ora il termine u_{n}^{∗} − u_{n}. Se f è L-Lips. (rispetto al secondo argomento) si prova che |u_{n}^{∗} − u_{n}| ≤ (1 + hL)|e_{n−1}| e

|e_{n}| ≤ h|τ(h)| + (1 + hL)|e_{n−1}|
≤ h|τ(h)| + (1 + hL)(h|τ(h)| + (1 + hL)|e_{n−2}|)
= (1 + (1 + hL))h|τ(h)| + (1 + hL)^{2}|e_{n−2}|

Analisi convergenza, IV

Di conseguenza, ragionando induttivamente abbiamo da 1 + s + … + s^{k} = (1 − s^{k+1})/(1 − s) che essendo |e_{0}| = 0,

|e_{n}| ≤ (1 + (1 + hL) + … + (1 + hL)^{n−1})h|τ(h)| + (1 + hL)^{n}|e_{0}|
≤ \frac{1 − (1 + hL)^{n}}{1 − (1 + hL)}h|τ(h)| + (1 + hL)^{n}|e_{0}| ≤ \frac{(1 + hL)^{n} − 1}{L}|τ(h)|

Osservato che (1 + hL)^{n} ≤ exp (nhL) e che nh = x_{n} − x_{0} , si ottiene

|e_{n}| ≤ \frac{exp (L(x_{n} − x_{0})) − 1}{L}τ(h) ≤ \frac{exp (L(x_{fin} − x_{0})) − 1}{L} · \frac{Mh}{2}

e di conseguenza, da (12), possiamo dire che il metodo di Eulero esplicito converge con ordine 1.

Metodi linear multistep

Posto h > 0, x_{n} = x_{0} + nh ed n ≥ 0, un metodo si dice linear multistep (abbr. LM) a p + 1 passi, se è del tipo

u_{n+1} = ∑_{j=0}^{p}a_{j}u_{n−j} + h ∑_{j=−1}^{p}b_{j}f (x_{n−j}, u_{n−j}), n = p, p + 1, … (14)

noti i valori u_{k} ≈ y(x_{k}) per k < p e supposto a_{p}, b_{p} ≠ 0. Si dimostra che un metodo LM è consistente se e solo se

∑_{j=0}^{p}a_{j} = 1, − ∑_{j=0}^{p}ja_{j} + ∑_{j=−1}^{p}b_{j} = 1

e inoltre che se y ∈ C^{q+1} per qualche q ≥ 1, il metodo ha ordine q se e solo se è consistente e

∑_{j=0}^{p}(−j)^{i}a_{j} + i ∑_{j=−1}^{p}(−j)^{i−1}b_{j} = 1, i = 2, …, q.

Metodi linear multistep: Adams

Da y′(x) = f (x, y(x)) e y(x_{n+1}) = y(x_{n−m}) + ∫_{x_{n−m}}^{x_{n+1}} y′(x) dx ricaviamo che

y(x_{n+1}) = y(x_{n−m}) + ∫_{x_{n−m}}^{x_{n+1}} f (x, y(x)) dx.

Se è nota una approssimazione della soluzione nei punti x_{n+γ}, …, x_{n+γ−p} e se P_{p}(x) è il polinomio che interpola le coppie (x_{n+γ−k},f (x_{n+γ−k}, u_{n+γ−k})), si ha

y(x_{n+1}) = y(x_{n−m}) + ∫_{x_{n−m}}^{x_{n+1}} P_{p}(x) dx

e quindi facilmente un metodo per quadratura numerica dell’integrale. Se m = 0, γ = 0 si hanno metodi espliciti (Adams-Bashforth), se m = 0, γ = 1 si hanno metodi impliciti (Adams-Moulton).

Alvise Sommariva  Equazioni differenziali ordinarie.  15/ 28

Esempi LM

▶ Adams-Bashforth (p=0): u_{n+1} = u_{n} + hf_{n},

τ_{n+1} = (h/2)y^{(2)}(ξ_{n}). ▶ Adams-Bashforth (p=1): u_{n+1} = u_{n} + \frac{h}{2}(3f_{n} − f_{n−1}),

τ_{n+1} = (5h^{2}/12)y^{(3)}(ξ_{n}). ▶ Adams-Bashforth (p=2):

u_{n+1} = u_{n} + \frac{h}{12}(23f_{n} − 16f_{n−1} + 5f_{n−1}), τ_{n+1} = (3h^{3}/8)y^{(4)}(ξ_{n}).

▶ Adams-Moulton (p=0): u_{n+1} = u_{n} + hf_{n+1},

τ_{n+1} = (−h/2)y^{(2)}(ξ_{n}). ▶ Adams-Moulton (p=1): u_{n+1} = u_{n} + \frac{h}{2}(f_{n+1} − f_{n}),

τ_{n+1} = (−h^{2}/12)y^{(3)}(ξ_{n}). ▶ Adams-Moulton (p=2): u_{n+1} = u_{n} + \frac{h}{12}(5f_{n+1} + 8f_{n} − f_{n−1}),

τ_{n+1} = (−h^{3}/24)y^{(4)}(ξ_{n}).

Convergenza linear LM

Si supponga che i dati iniziali siano calcolati in modo che

η(h) = max_{i=0,...,p}|y(x_{i}) − u_{i}| → 0, per h → 0

e che il metodo sia consistente. Supponiamo inoltre che

a_{j} ≥ 0, j = 0, …, p

mentre il passo di discretizzazione h sia tale h ≤ 1/(2c), c = L ∑_{j=−1}^{p}|b_{j}|,

dove L è la costante di Lipschitz di f . Allora il metodo LM converge ed inoltre esistono C_{1}, C_{2} positive tali che per ogni n

|y(x_{n}) − u_{n}| ≤ C_{1}η(h) + C_{2}τ(h).
Alvise Sommariva  Equazioni differenziali ordinarie.  17/ 28

Metodi di Eulero in Matlab

Codice di Eulero esplicito, eseguito da Atkinson, Han, Steward:

f u n c t i o n [ t , y ] = eulero_esplicito ( t0 , y0 , t_end , h , fcn )

n = f i x ( ( t_end−t0 ) /h ) +1;
t = l i n s p a c e ( t0 , t0+(n−1)∗h , n ) ’ ;
y = z e r o s ( n , 1 ) ;
y ( 1 ) = y0 ;
f o r i = 2 : n
y ( i ) = y ( i−1) + h∗ f e v a l ( fcn , t ( i−1) , y ( i−1) ) ;

end

Metodi di Eulero in Matlab

Codice di Eulero implicito, eseguito da Atkinson, Han, Steward:

f u n c t i o n [ t , y]= eulero_implicito ( t0 , y0 , t_end , h , fcn , tol )

n = f i x ( ( t_end−t0 ) /h ) + 1 ;
t = l i n s p a c e ( t0 , t0+(n−1)∗h , n ) ’ ;
y = z e r o s ( n , 1 ) ; y ( 1 ) = y0 ;
f o r i=2:n
yt1 = y ( i−1) + h∗ f e v a l ( fcn , t ( i−1) , y ( i−1) ) ;
count = 0 ; d i f f = 1 ;
w h i l e d i f f > tol && count < 10
yt2 = y ( i−1) + h∗ f e v a l ( fcn , t ( i ) , yt1 ) ;
d i f f = abs ( yt2−yt1 ) ; yt1 = yt2 ;

  count = count +1;

end

i f count >= 10
d i s p ( ’ Not c o n v e r g i n g a f t e r 10 s t e p s at t = ’ )
f p r i n t f ( ’ %5.2 f \n ’ , t ( i ) )

end

y ( i ) = yt2 ;

end

Alvise Sommariva  Equazioni differenziali ordinarie.  19/ 28

Esercizio

Esercizio: ▶ Implementare i metodi LM della tabella precedente. ▶ Per ognuno di questi, approssimare la soluzione dell’equazione

test dell’esempio 1, nell’intervallo [0, 10], in N punti equispaziati x_{n}, con N = 2, 4, 8, 16, 32, descrivendone l’errore ∥y_{n} − u_{n}∥_{∞}, dove y_{n} = y(x_{n}) e u_{n} l’approssimazione fornita dal metodo numerico.

Runge Kutta, I

Consideriamo un metodo numerico del tipo

u_{n+1} = u_{n} + F (x_{n}, u_{n}; h)

con

F (x, y; h) = γ_{1}f (x, y)γ_{2}f (x + αh, y + βhf (x, y))

e determiniamo i parametri γ_{1}, γ_{2}, α, β così da ottenere un metodo del second’ordine. Osserviamo che questo implica che τ(h) sia O(h^{2}). Per ottenere questo risultato usiamo la formula di Taylor bivariata. Denotate con f_{x}, f_{y} le derivate parziali rispetto al primo e secondo argomento di f , abbiamo

f (x+αh, y+βhf (x, y)) = f (x, y)+αhf_{x}(x, y)+βhf_{y}(x, y)f (x, y)+O(h^{2})

Runge Kutta, II

Da f (x+αh, y+βhf (x, y)) = f (x, y)+αhf_{x}(x, y)+βhf_{y}(x, y)f (x, y)+O(h^{2}) abbiamo che

F (x, y; h) = γ_{1}f (x, y) + γ_{2}f (x + αh, y + βhf (x, y))
= γ_{1}f (x, y) + γ_{2}(f (x, y) + αhf_{x}(x, y) + βhf_{y}(x, y)f (x, y))
+O(h^{2})

Facilmente, dalla formula di Taylor

y(x + h) − y(x) = hy^{(1)}(x) + (h^{2}/2)y^{(2)}(x) + O(h^{3})
= hf (x, y) + (h^{2}/2)(f_{x}(x, y) + f_{y}(x, y)f (x, y)) + O(h^{3}) (15)