Equazioni differenziali ordinarie.

Alvise Sommariva

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

19 maggio 2014

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/ 61

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

Sia \overline{B}_{ϵ}(y_{0}) la palla chiusa di ℝ^{n} di centro y_{0} e raggio ϵ e sia

S_{δ,ϵ,x_{0},y_{0}} ≡ [x_{0} − δ, x_{0} + δ] × \overline{B}_{ϵ}(y_{0}) ⊆ ℝ × ℝ^{n}, δ, ϵ > 0.

Ricordiamo che la funzione f è L-lipschitziana rispetto la seconda variabile in S_{δ,ϵ,x_{0},y_{0}} se L > 0 e

|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} contenente (x_{0}, y_{0}), f : Ω → ℝ^{n} continua e L-lipschitziana (rispetto alla seconda variabile) in S_{δ,ϵ,x_{0},y_{0}} (per qualche δ, ϵ > 0).

Allora esiste un intervallo aperto contenente x_{0} nel quale è definita \underline{una e una sola} soluzione del problema di Cauchy (1).

Alvise Sommariva  Equazioni differenziali ordinarie.  4/ 61

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 : Ω → ℝ^{n} continua e L-lipschitziana (rispetto la seconda variabile) in Ω.

Allora il problema di Cauchy (1) ha una e una sola soluzione in [x_{0}, x_{0} + a].

Si noti che

▶ Nel teorema di Cauchy in piccolo Ω è aperto, mentre in quello

in grande è chiuso. ▶ Il teorema di Cauchy in piccolo fornisce l’esistenza e unicità

locale della soluzione, mentre quello in grande è globale.

Metodo di Eulero esplicito

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

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}) (4)

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})).  (5)

Metodo di Eulero esplicito

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}), y_{0} = y(x_{0}).  (6)

Essendo

y(x_{k+1}) = y(x_{k}) + y′(x_{k})(x − x_{k}) + \frac{y′′(ξ)(x_{k+1} − x_{k})^{2}}{2}

si osserva immediatamente che applicando il metodo di Eulero esplicito ponendo y_{k} = y(x_{k}), si ha per qualche ξ in (x_{k}, x_{k+1}), in virtù del fatto che y′(x_{k}) = f (x_{k}, y(x_{k}))

y(x_{k+1}) − y_{k+1} = y(x_{k}) + y′(x_{k})(x_{k+1} − x_{k})  (7)
+ \frac{y′′(ξ)(x_{k+1} − x_{k})^{2}}{2} − (y(x_{k}) + y′(x_{k})(x_{k+1} − x_{k}))
= \frac{y′′(ξ)(x_{k+1} − x_{k})^{2}}{2} .

Alvise Sommariva  Equazioni differenziali ordinarie.  7/ 61

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}), (8)

e quindi

y(x_{k+1}) = y(x_{k}) + h · f (x_{k+1}, y(x_{k+1})).  (9)

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}), y_{0} = y(x_{0}).  (10)

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.  8/ 61

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}  (11)

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})  (12)

e Eulero implicito

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

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

Alvise Sommariva  Equazioni differenziali ordinarie.  9/ 61

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 (1) 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.

Se un metodo per la soluzione di problemi di Cauchy ha la forma

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

allora, posto y_{n} = y(t_{n}) il valore della sol. del problema di Cauchy in t_{n} (cf. [1, p.112]), si definisce

τ_{n}(h) = \frac{1}{h} · (y_{n+1} − (∑_{j=0}^{p}a_{j}y_{n−j} + h ∑_{j=−1}^{p}b_{j}f (x_{n−j}, y_{n−j}))).

si chiama errore locale di troncamento del metodo, mentre

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

si chiama errore globale di troncamento del metodo. Un metodo per cui τ(h) tende a 0, quando h → 0 si dice consistente.

Alvise Sommariva  Equazioni differenziali ordinarie.  11/ 61

Analisi convergenza

Nota.

Osserviamo che il metodo di Eulero esplicito

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

e Eulero implicito

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

hanno la forma

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

ponendo rispettivamente

p = 0, a_{0} = 1, b_{−1} = 0, b_{0} = 1

e p = 0, a_{0} = 1, b_{−1} = 1, b_{0} = 0.

Alvise Sommariva  Equazioni differenziali ordinarie.  12/ 61

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, x_{n} = x_{0} + nh e

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}).

e osserviamo che

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

Il metodo converge se entrambi i termini a secondo membro di (14) convergono a 0 per h → 0. Osserviamo che

▶ la prima sequenza u_{n}^{∗} parte dal valore assunto dalla soluzione

in x_{n−1}, ▶ la seconda seq. u_{n} parte da valori che a priori non coincidono

con y_{n−1} = y(x_{n−1}).

Alvise Sommariva  Equazioni differenziali ordinarie.  13/ 61

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

|y_{n} − u_{n}^{∗}| = |y_{n} − (y_{n−1} + h · f (x_{n−1}, y_{n−1}))| = (h^{2}/2) · |y′′(ξ_{n})|

per qualche ξ_{n} ∈ (t_{n−1}, t_{n}). Osserviamo che in questo caso, poichè

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

abbiamo che

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

Se la derivata seconda di y esiste ed è continua, posto M = max_{x∈I} |y′′(x)|, per il metodo di Eulero esplicito,

|τ_{n}(h)| = |y_{n}−u_{n}^{∗}|/h = (1/h)·(h^{2}/2)·|y′′(ξ_{n})| = (h/2)·|y′′(ξ_{n})| ≤ Mh/2.

e quindi

τ(h) = max_{n=0,...}|τ_{n}(h)| ≤ Mh/2.

Quindi sotto queste ipotesi il metodo di Eulero esplicito risulta consistente.

Osserviamo ora che essendo

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

abbiamo

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

Analisi convergenza Eulero esplicito, caso Lipschitziano

Consideriamo ora il termine u_{n}^{∗} − u_{n}. Se f è L-Lipschitziana (rispetto al secondo argomento) allora

|f (x, y_{1}) − f (x, y_{2})| ≤ L|y_{1} − y_{2}|

e si mostra che

|u_{n}^{∗} − u_{n}| ≤ (1 + hL)|e_{n−1}|.

Infatti da

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})

ricaviamo dalla L-lipschitzianità

|u_{n}^{∗} − u_{n}| = |(y_{n−1} + h · f (x_{n−1}, y_{n−1})) − (u_{n−1} + h · f (x_{n−1}, u_{n−1}))|
= |(y_{n−1} − u_{n−1}) + h · (f (x_{n−1}, y_{n−1}) − f (x_{n−1}, u_{n−1}))|
≤ |y_{n−1} − u_{n−1}| + h|f (x_{n−1}, y_{n−1}) − f (x_{n−1}, u_{n−1}))|
≤ |y_{n−1} − u_{n−1}| + hL|y_{n−1} − u_{n−1}| = (1 + hL)|y_{n−1} − u_{n−1}|
Alvise Sommariva  Equazioni differenziali ordinarie.  16/ 61

Analisi convergenza Eulero esplicito, caso Lipschitziano

Da e_{n} := y_{n} − u_{n}, |y_{n} − u_{n}^{∗}| ≤ h|τ(h)|, |u_{n}^{∗} − u_{n}| ≤ (1 + hL)|e_{n−1}|

|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}|

Di conseguenza, ragionando induttivamente abbiamo da 1 + s + … + s^{k} = (1 − s^{k+1})/(1 − s) ed 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)|

Alvise Sommariva  Equazioni differenziali ordinarie.  17/ 61

Analisi convergenza Eulero esplicito, caso Lipschitziano

Notiamo che per γ > 0

exp (γ) = ∑_{k=0}^{∞} \frac{γ^{k}}{k!} ≥ ∑_{k=0}^{1} \frac{γ^{k}}{k!} = 1 + γ  (17)

e quindi (1 + γ)^{n} ≤ (exp(γ))^{n} = exp(nγ) implica che per γ = hL si ha essendo nh = x_{n} − x_{0}

(1 + hL)^{n} ≤ exp (nhL) = exp (L(x_{n} − x_{0})).

Così da

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

si ottiene che essendo x_{0} < x_{n} ≤ x_{fin}

|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 (14), possiamo dire che il metodo di Eulero esplicito converge con ordine 1.

Alvise Sommariva  Equazioni differenziali ordinarie.  18/ 61

Analisi convergenza Eulero esplicito, caso dissipativo

Se invece di essere L-lipschitziana, f è L-dissipativa, cioè

−L ≤ \frac{∂f}{∂y} ≤ 0,

abbiamo da

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})

f (x_{n−1}, u_{n−1}) = f (x_{n−1}, y_{n−1}) + \frac{∂f}{∂y} (ξ)(u_{n−1} − y_{n−1}) che per qualche ξ ∈ I(u_{n−1}, y_{n−1})

u_{n} − u_{n}^{∗} = (u_{n−1} + h · f (x_{n−1}, u_{n−1})) − (y_{n−1} + h · f (x_{n−1}, y_{n−1}))
= (u_{n−1} − y_{n−1}) + h · (f (x_{n−1}, u_{n−1}) − f (x_{n−1}, y_{n−1}))
= (u_{n−1} − y_{n−1}) + h\frac{∂f}{∂y} (ξ)(u_{n−1} − y_{n−1})
= (1 + h\frac{∂f}{∂y} (ξ)) · (u_{n−1} − y_{n−1}).  (18)
Alvise Sommariva  Equazioni differenziali ordinarie.  19/ 61

Analisi convergenza Eulero esplicito, caso dissipativo

Da

u_{n} − u_{n}^{∗} = (1 + h\frac{∂f}{∂y} (ξ)) · (u_{n−1} − y_{n−1}).  (19)

abbiamo che se 0 < h ≤ 2/L (non restrittivo per mostrare la conv. visto che si studia il comportamento per h → 0) allora

−1 = 1 − 2L/L ≤ 1 − hL≤ 1 + h\frac{∂f}{∂y} (ξ) ≤ 1

e quindi

|1 + h\frac{∂f}{∂y} (ξ)| ≤ 1

da cui

|u_{n} − u_{n}^{∗}| = |1 + h\frac{∂f}{∂y} (ξ)| · |u_{n−1} − y_{n−1}|
≤ |u_{n−1} − y_{n−1}|.  (20)
Alvise Sommariva  Equazioni differenziali ordinarie.  20/ 61

Analisi convergenza Eulero esplicito, caso dissipativo

Da e_{n} := y_{n} − u_{n}, |y_{n} − u_{n}^{∗}| ≤ h|τ(h)|, |u_{n}^{∗} − u_{n}| ≤ |e_{n−1}|

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

Di conseguenza, ragionando induttivamente abbiamo che essendo |e_{0}| = 0, nh ≤ Nh = xfin − x_{0}, |τ(h)| ≤ Mh/2

|e_{n}| ≤ nh|τ(h)| + |e_{0}|
≤ (xfin − x_{0})|τ(h)| ≤ (xfin − x_{0})Mh/2

e di conseguenza possiamo dire che il metodo di Eulero esplicito converge con ordine 1.

Facoltativo. Analisi convergenza Eulero implicito, caso Lipschitziano

Consideriamo il metodo di Eulero implicito

u_{n} = u_{n−1} + hf (t_{n}, u_{n}), u_{0} = y_{0}.

Per prima cosa calcoliamo l’errore globale di troncamento. Se g ∈ C^{2}(a, b), per qualche ξ ∈ (a, b) si ha

∫_{a}^{b} g(x)dx = (b − a)g(b) − h^{2}g^{(1)}(ξ)/2

e quindi essendo h = x_{n} − x_{n−1}, se la soluzione è suff. regolare

y_{n} = y_{n−1} + ∫_{x_{n−1}}^{x_{n}} y′(x)dx = y_{n−1} + hy′(x_{n}) − h^{2}y′′(ξ)/2
= y_{n−1} + hf (x_{n}, y_{n}) − h^{2}y′′(ξ)/2  (21)

e quindi

τ_{n}(h) = (1/h) · (h^{2}y′′(ξ)/2) = hy′′(ξ)/2.
Alvise Sommariva  Equazioni differenziali ordinarie.  22/ 61

Facoltativo. Analisi convergenza Eulero implicito, caso Lipschitziano

Ricordiamo che y_{n} = y(t_{n}) è la soluzione del problema di Cauchy mentre

u_{n} = u_{n−1} + hf (t_{n}, u_{n}), u_{0} = y_{0}.

Definiamo

u_{n}^{∗} := y_{n−1} + hf (t_{n}, y_{n}), u_{0}^{∗} = y_{0}.

Come in precedenza

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

e quindi per la disuguaglianza triangolare

|y_{n} − u_{n}| ≤ |y_{n} − u_{n}^{∗}| + |u_{n}^{∗} − u_{n}|.

Studiamo separatamente i termini al secondo membro. Osserviamo che

|y_{n} − u_{n}^{∗}| = |y_{n} − (y_{n−1} + hf (t_{n}, y_{n}))| = h|τ_{n}(h)|

Alvise Sommariva  Equazioni differenziali ordinarie.  23/ 61

Facoltativo. Analisi convergenza Eulero implicito, caso Lipschitziano

Studiamo |u_{n}^{∗} − u_{n}| ricordando che

u_{n}^{∗} := y_{n−1} + hf (t_{n}, y_{n}), u_{0}^{∗} = y_{0},

u_{n} := u_{n−1} + hf (t_{n}, u_{n}), u_{0} = y_{0}.

Di conseguenza, dalla lipschitzianità (rispetto la seconda variabile)

|u_{n}^{∗} − u_{n}| = |y_{n−1} + hf (t_{n}, y_{n}) − u_{n−1} − hf (t_{n}, u_{n})|
= |y_{n−1} − u_{n−1} + hf (t_{n}, y_{n}) − hf (t_{n}, u_{n})|
≤ |y_{n−1} − u_{n−1}| + h|f (t_{n}, y_{n}) − f (t_{n}, u_{n})|
≤ |y_{n−1} − u_{n−1}| + hL|y_{n} − u_{n}|  (22)

In definitiva, essendo

|y_{n} − u_{n}| ≤ |y_{n} − u_{n}^{∗}| + |u_{n}^{∗} − u_{n}|,
|y_{n} − u_{n}^{∗}| = |y_{n} − (y_{n−1} + hf (t_{n}, y_{n}))| = h|τ_{n}(h)|,
|u_{n}^{∗} − u_{n}| ≤ |y_{n−1} − u_{n−1}| + hL|y_{n} − u_{n}|,

deduciamo che

|y_{n} − u_{n}| ≤ h|τ_{n}(h)| + |y_{n−1} − u_{n−1}| + hL|y_{n} − u_{n}|.

Visto che siamo interessati ad analizzare l’errore per h → 0, non è restrittivo supporre hL < 1 e portando hL|y_{n} − u_{n}| a primo membro e dividendo ambo i membri per 1 − hL > 0 abbiamo, posto e_{n} := y_{n} − u_{n},

|e_{n}| = |y_{n} − u_{n}| ≤ (1/(1 − hL))(h|τ_{n}(h)| + |y_{n−1} − u_{n−1}|)
= \frac{h|τ_{n}(h)|}{1 − hL} + \frac{|e_{n−1}|}{1 − hL}  (23)
Alvise Sommariva  Equazioni differenziali ordinarie.  25/ 61

Con ragionamenti simili a quelli utilizzati per mostrare la convergenza di Eulero esplicito, abbiamo per s = 1/(1 − hL), posto τ(h) = max_{n} |τ_{n}(h)|

|e_{n}| ≤ sh|τ_{n}(h)| + s|e_{n−1}|
≤ sh|τ_{n}(h)| + s(sh|τ_{n−1}(h)| + s|e_{n−2}|)
≤ …
≤ h ∑_{k=1}^{n}s^{k}|τ_{n−k+1}(h)| + s^{n}|e_{0}|

≤ h(∑_{k=1}^{n}s^{k})|τ(h)| + s^{n}|e_{0}|  (24)

Facoltativo. Analisi convergenza Eulero implicito, caso Lipschitziano

Notiamo che se hL < 1 allora s = 1/(1 − hL) > 1 e che

∑_{k=1}^{n}s^{k} = \frac{s^{n+1} − 1}{s − 1} ≤ \frac{s^{n+1}}{s − 1}

Ora

s − 1 = \frac{1}{1 − hL} − 1 = \frac{hL}{1 − hL} > 0  (25)

e da (1 + x)^{k} ≤ exp(kx) per x > 0, abbiamo

s^{n+1} = (1 + (s − 1))^{n+1} ≤ exp((n + 1)(s − 1))
= exp(\frac{hL(n + 1)}{1 − hL}) = exp(\frac{L(x_{n+1} − x_{0})}{1 − hL})  (26)

e quindi

∑_{k=1}^{n}s^{k} ≤ \frac{(1 − hL) exp(\frac{L(x_{n+1}−x_{0})}{1−hL})}{hL}

Alvise Sommariva  Equazioni differenziali ordinarie.  27/ 61

Facoltativo. Analisi convergenza Eulero implicito, caso Lipschitziano

Così ricapitolando

|e_{n}| ≤ h(∑_{k=1}^{n}s^{k})|τ(h)| + s^{n}|e_{0}|,

∑_{k=1}^{n}s^{k} ≤ \frac{(1 − hL) exp(\frac{L(x_{n+1}−x_{0})}{1−hL})}{hL}

s^{n} ≤ exp(\frac{L(x_{n} − x_{0})}{1 − hL}),

implica finalmente

|e_{n}| ≤ h\frac{(1 − hL) exp(\frac{L(x_{n+1}−x_{0})}{1−hL})}{hL}|τ(h)| + exp(\frac{L(x_{n} − x_{0})}{1 − hL})|e_{0}|.

Alvise Sommariva  Equazioni differenziali ordinarie.  28/ 61

Notiamo che τ_{n}(h) = h|y′′(ξ)|/2 e che se y′′ è continua per il teorema di Weierstrass |y′′| ha massimo nel compatto [x_{0}, xfin], diciamo M. Di conseguenza

τ(h) = max(τ_{n}(h)) ≤ Mh/2.

Ricordando che |e_{0}| = 0, che xfin = x_{0} + Nh abbiamo

|e_{n}| ≤ h\frac{(1 − hL) exp(\frac{L(x_{n+1}−x_{0})}{1−hL})}{hL}|τ(h)| + exp(\frac{L(x_{n} − x_{0})}{1 − hL})|e_{0}|
≤ h\frac{(1 − hL) exp(\frac{L(x_{n+1}−x_{0})}{1−hL})}{hL} \frac{Mh}{2} + exp(\frac{L(xfin − x_{0})}{1 − hL})|e_{0}|
≤ \frac{(1 − hL) exp(\frac{L(x_{N}−x_{0})+hL}{1−hL})}{L} \frac{Mh}{2}  (27)

e quindi |e_{n}| → 0 per h → 0, per cui Eulero implicito converge.

Alvise Sommariva  Equazioni differenziali ordinarie.  29/ 61

Metodo di Crank-Nicolson

Supponiamo di dover risolvere il problema di Cauchy

\{ \begin{matrix}y′(t) = f (t, y(t)), t ≥ t_{0} \\ y(t_{0}) = y_{0} \end{matrix}  (28)

Osserviamo che se t_{n+1} = t_{n} + h allora ricordando la formula del trapezio per il calcolo di integrali

y(t_{n+1}) = y(t_{n}) + ∫_{t_{n}}^{t_{n+1}} y′(t)dt
= y(t_{n}) + ∫_{t_{n}}^{t_{n+1}} f (t, y(t))dt
≈ y(t_{n}) + (h/2)(f (t_{n}, y(t_{n})) + f (t_{n+1}, y(t_{n+1})))
Alvise Sommariva  Equazioni differenziali ordinarie.  30/ 61

Metodo di Crank-Nicolson

Da

y(t_{n+1}) ≈ y(t_{n}) + (h/2)(f (t_{n}, y(t_{n})) + f (t_{n+1}, y(t_{n+1}))) (29)

si introduce il metodo implicito (detto di Crank-Nicolson o dei dei trapezi)

u_{n+1} ≈ u_{n} + (h/2)(f (t_{n}, u_{n}) + f (t_{n+1}, u_{n+1}))  (30)

per risolvere numericamente il problema di Cauchy (28).

Si osservi che ad ogni iterazione, essendo il metodo implicito, bisogna risolvere una equazione nonlineare.

Assoluta stabilitá

Nell’ambito delle equazioni differenziali ordinarie, esistono vari criteri di stabilitá. Un classico problema è quello di vedere se un metodo è assolutamente stabile. Definito il problema di Cauchy

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

per un certo λ ∈ ℂ con ℜ(λ) < 0, visto che l’unica soluzione è y(x) = exp (λx) si cerca di definire il passo h cosicchè il metodo numerico per la soluzione del problema di Cauchy (1) abbia lo stesso comportamento asintotico di exp (λx).

La regione di assoluta stabilità e’ composta dagli hλ, per cui il metodo numerico con passo h è tale da avere lo stesso comportamento asintotico del problema di Cauchy (31) relativo al parametro λ (con ℜ(λ) < 0).

Assoluta stabilitá

Ricordiamo a tal proposito che dalla formula di Eulero, se z = a + ib = ℜ(z) + ℑ(z) allora

exp (z) = exp (a) · (cos (b) + i sin (b)).

Quindi se ℜ(λ) < 0, visto che | cos (ℑ(λx)) + i sin (ℑ(λx))| = 1 e x > 0 abbiamo

| exp (λx)| = |exp(ℜ(λx))|| cos (ℑ(λx)) + i sin ((ℑ(λx)))|
= |exp (ℜ(λx))| → 0, per x → ∞.

Visto il comportamento asintotico di exp (λx), se u^{(h)}(x_{n}) è l’approssimazione della soluzione in x_{n} fornita da un metodo numerico a passo h, si desidera sia u^{(h)}(x_{n}) → 0 per n → +∞.

Assoluta stabilitá: Eulero esplicito

Nel caso del metodo di Eulero esplicito,

u_{n} = u_{n−1} + hf (x_{n−1}, u_{n−1}) = (1 + hλ)u_{n−1}, u_{0} = 1.

Si verifica facilmente che questa equazione alle differenze (di ordine 1) ha quale unica soluzione

u_{n} = (1 + hλ)^{n}

e che

u_{n} → 0 se e solo se |1 + hλ| < 1

il che significa che hλ deve stare nel disco del piano complesso di centro −1 e raggio 1.

Assoluta stabilitá: Eulero esplicito

Di conseguenza, fissato λ, il metodo risulta stabile se e solo se si sceglie un passo sufficientemente piccolo, cioè minore di 1/|λ|.

Un problema di Cauchy

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

si dice stiff, quando ℜ(λ) ≪ 0 (cioè ℜ(λ) è molto grande in valore assoluto).

In questo caso, si deve scegliere un passo h molto piccolo affinchè la soluzione di Eulero esplicito tenda a 0.

Assoluta stabilitá: Eulero implicito

Nel caso del metodo di Eulero implicito,

u_{n} = u_{n−1} + hf (x_{n}, u_{n}) = u_{n−1} + hλu_{n}, u_{0} = 1.

Portando a primo membro hλu_{n}, dividendo i membri per 1 − hλ

u_{n} = \frac{u_{n−1}}{1 − hλ}.

Si verifica che

u_{n} = \frac{1}{(1 − hλ)^{n}}

e che

u_{n} → 0 se e solo se \frac{1}{|1 − hλ|} < 1.

Essendo ℜ(λ) < 0, si vede facilmente che |\frac{1}{1−hλ}| < 1 per qualsiasi h. Questa proprietà suggerisce di applicare Eulero impl. invece di Eulero espl. per risolvere num. un problema di Cauchy.

La regione di ass. stabilità è tutto il semipiano negativo ℜ(hλ) < 0.

Alvise Sommariva  Equazioni differenziali ordinarie.  36/ 61

Assoluta stabilitá: metodo di Crank-Nicolson

Nel caso del metodo di Crank-Nicolson, da f (x, y) = λy

u_{n} = u_{n−1} + (h/2)(f (x_{n}, u_{n}) + f (x_{n−1}, u_{n−1}))
= u_{n−1} + (h/2)λ(u_{n} + u_{n−1}), u_{0} = 1  (33)

Quindi,

(1 − \frac{hλ}{2})u_{n} = u_{n−1} + \frac{hλ}{2}u_{n}

cioè

\frac{2 − hλ}{2} · u_{n} = \frac{2 + hλ}{2} · u_{n−1}

e di conseguenza

u_{n} = \frac{2 + hλ}{2 − hλ}u_{n−1}

Alvise Sommariva  Equazioni differenziali ordinarie.  37/ 61

Assoluta stabilitá: metodo di Crank-Nicolson

Da

u_{n} = \frac{2 + hλ}{2 − hλ}u_{n−1}

si verifica facilmente che questa equazione alle differenze ha quale unica soluzione

u_{n} = \frac{(2 + hλ)^{n}}{(2 − hλ)^{n}}

e che

u_{n} → 0 se e solo se \frac{|2 + hλ|}{|2 − hλ|} < 1.

Da ℜ(λ) < 0, |\frac{2+hλ}{2−hλ}| < 1 per ogni h. Infatti, se λ = a + ib, a < 0

\frac{|2 + hλ|}{|2 − hλ|} = \frac{|2 + h(a + ib)|}{|2 − h(a + ib)|} = \frac{|(2 + ha) + ib|}{|(2 − ha) + ib|} = \begin{matrix}√ & \overline{(2 + ha)^{2} + b^{2}} \\ √ & \overline{(2 − ha)^{2} + b^{2}} \end{matrix} < 1

in quanto (2 + ha)^{2} < (2 − ha)^{2}, qualora a < 0. La regione di ass. stabilità è tutto il semipiano negativo ℜ(hλ) < 0.

Alvise Sommariva  Equazioni differenziali ordinarie.  38/ 61

Facoltativo: Assoluta stabilitá. Una nota

Nota.

Stabilire la stabilità per il problema di Cauchy

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

è in generale complicato. Se invece consideriamo

\{ \begin{matrix}y′(x) = λy(x) + g(x), x ≥ x_{0} \\ y(x_{0}) = y_{0} \end{matrix}  (35)

la questione è più semplice. Ricordiamo che la soluzione di (??) è

y(x) = c · exp (λx) + ∫_{x_{0}}^{x} exp (λ(x − t))g(t)dt

con

c = exp (−λx_{0})y(x_{0})

(cf. [?, p.369]).

Alvise Sommariva  Equazioni differenziali ordinarie.  39/ 61

Facoltativo: Assoluta stabilitá. Una nota

Sia y la soluzione di

\{ \begin{matrix}y′(x) = λy(x) + g(x), x ≥ x_{0} \\ y(x_{0}) = y_{0} \end{matrix}  (36)

e y_{ϵ} quella di

\{ \begin{matrix}y′(x) = λy(x) + g(x), x ≥ x_{0} \\ y(x_{0}) = y_{0} + ϵ \end{matrix}  (37)

Posto z_{ϵ} = y_{ϵ} − y, da (??), (??),

\{ \begin{matrix}z′_{ϵ}(x) = λz_{ϵ}(x), x ≥ x_{0} \\ y(x_{0}) = ϵ \end{matrix}  (38)

la cui soluzione è z_{ϵ} = ϵ · exp (λx).

Alvise Sommariva  Equazioni differenziali ordinarie.  40/ 61

Facoltativo: Assoluta stabilitá. Una nota

Di solito nelle applicazioni ci si interessa al caso λ < 0 o complesso con parte reale negativa. In questi casi z_{ϵ} → 0 per n → ∞ e quindi l’effetto delle perturbazioni si annulla per valori grandi di x.

Si desidera che il metodo numerico goda delle stesse proprietà e la risposta la si ha nuovamente dallo studio della regione di stabilità .

Metodi linear multistep (LM)

Posto h > 0, x_{n} = x_{0} + nh ed n ≥ 0, un metodo si dice linear multistep (abbr. LM), 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, … (39)

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 (LM)

Nota.

Si osservi che per metodi a più passi non essendo noti valori approssimati della soluzione nei punti x_{0}, …, x_{p}, usualmente li si fornisce utilizzando un metodo a meno passi, come ad esempio Eulero implicito.

Nota.

I metodi già proposti sono di questo tipo:

▶ Eulero esplicito: a_{0} = 1, b_{−1} = 0, b_{0} = 1;

▶ Eulero implicito: a_{0} = 1, b_{−1} = 1, b_{0} = 0;

▶ Eulero esplicito: a_{0} = 1, b_{−1} = 1/2, b_{0} = 1/2; Applicando i teoremi precedenti si vede, come noto che

▶ Eulero esplicito ha ordine di consistenza 1,

▶ Eulero implicito ha ordine di consistenza 1,

▶ il metodo dei trapezi ha ordine di consistenza uguale a 2.

Alvise Sommariva  Equazioni differenziali ordinarie.  43/ 61

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).

Esempi LM: Adams-Bashforth

▶ 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−2}), τ_{n+1} = (3h^{3}/8)y^{(4)}(ξ_{n}).

▶ Adams-Bashforth (p=3):

u_{n+1} = u_{n} + \frac{h}{24}(55f_{n} − 59f_{n−1} + 37f_{n−2} − 9f_{n−3}), τ_{n+1} = (251h^{4}/720)y^{(5)}(ξ_{n}).

Nota.

Si noti che sono effettivamente metodi espliciti a p + 1 passi.

Esempi LM: Adams-Moulton

▶ 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}). ▶ Adams-Moulton (p=3):

u_{n+1} = u_{n} + \frac{h}{24}(9f_{n+1} + 19f_{n} − 5f_{n−1} + f_{n−2}), τ_{n+1} = (−19h^{4}/720)y^{(5)}(ξ_{n}).

Nota.

Si noti che sono effettivamente metodi impliciti a p passi.

Convergenza 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.  47/ 61

Barriere LM

Definizione

Un metodo numerico è detto A-stabile se la sua regione di stabilità assoluta contiene tutto il semipiano negativo.

Valgono le seguenti barriere di stabilità (Dahlquist, 1963)

Teorema

Nessun metodo LM esplicito è A-stabile.

Teorema

Nessun metodo LM implicito di ordine maggiore di 2 è A-stabile.

Predictor-corrector

In generale un metodo LM implicito ha proprietà di stabilità migliori rispetto ad un LM esplicito e questo ne suggerisce l’utilizzo.

Purtroppo ad ogni iterazione bisogna risolvere una equazione nonlineare.

Un approccio comunemente utilizzato è quello del predictor-corrector. Consiste nell’utilizzare il metodo di punto fisso con una certa scelta del punto iniziale.

Nel caso di un metodo di Adams implicito (detto corrector), si utilizza ad esempio un metodo di Adams esplicito (detto predictor).

Esempio Predictor-corrector

Consideriamo i seguenti metodi di Adams (porremo f_{i} = f (x_{i}, u_{i})). Predictor: u_{i+1} = u_{i} + \frac{h}{2}(3f_{i} − f_{i−1}) Corrector: u_{i+1} = u_{i} + \frac{h}{2}(f_{i+1} + f_{i}).

Posto f_{i}^{(η)} = f (x_{i}, u_{i}^{(η)}) abbiamo

P:  u_{i+1}^{(0)} = u_{i}^{(1)} + \frac{h}{2}(3f_{i}^{(1)} − f_{i−1}^{(1)})
E:  f_{i+1}^{(0)} = f (x_{i+1}, u_{i+1}^{(0)})
C:  u_{i+1}^{(1)} = u_{i}^{(1)} + \frac{h}{2}(f_{i+1}^{(0)} + f_{i}^{(1)})
E:  f_{i+1}^{(1)} = f (x_{i+1}, u_{i+1}^{(1)})
Alvise Sommariva  Equazioni differenziali ordinarie.  50/ 61

Runge-Kutta

Consideriamo un metodo numerico che definisca una sequenza u_{n} t.c.

y(x_{n+1}) ≈ u_{n+1} = u_{n} + hF (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. Ricordiamo che questo significa 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})

Per definizione di errore locale di troncamento:

τ(x, h) = (y(x + h) − u(x + h))/h

con u(x + h) calcolata dal metodo qualora u(x) = y(x). 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})  (40)

Runge-Kutta

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}) = ϕ_{1}(x, h) + O(h^{3})  (41)

Ma è pure

u(x + h) − y(x) = hF (x, y; h) = h(γ_{1}f (x, y) + γ_{2}(f (x, y)

+ αhf_{x}(x, y) + βhf_{y}(x, y)f (x, y))) + O(h^{3})
= h(γ_{1} + γ_{2})f (x, y) + h^{2}γ_{2}(αf_{x}(x, y) + βf_{y}(x, y)f (x, y)))
+ O(h^{3}) = ϕ_{2}(x, h, γ_{1}, γ_{2}, α, β) + O(h^{3})  (42)
Alvise Sommariva  Equazioni differenziali ordinarie.  53/ 61

Se

ϕ_{1}(x, h) ≡ ϕ_{2}(x, h, γ_{1}, γ_{2}, α, β)

allora da

y(x + h) − y(x) = ϕ_{1}(x, h) + O(h^{3})
u(x + h) − y(x) = ϕ_{2}(x, h, γ_{1}, γ_{2}, α, β) + O(h^{3})

sottraendo membro a membro

y(x + h) − u(x + h) = (ϕ_{1}(x, h) − ϕ_{2}(x, h, γ_{1}, γ_{2}, α, β)) + O(h^{3})
= O(h^{3})  (43)

e quindi si determina un metodo avente ordine di consistenza 2, poichè τ_{n}(h) = O(h^{2}) per ogni n e quindi τ(h) = O(h^{2}). Affinchè ϕ_{1}(x, h) = ϕ_{2}(x, h, γ_{1}, γ_{2}, α, β), per confronto tra (35), (36), (37), necessita

γ_{1} + γ_{2} = 1, γ_{2}α = 1/2, γ_{2}β = 1/2.

Alvise Sommariva  Equazioni differenziali ordinarie.  54/ 61

Runge-Kutta

Vediamo alcuni metodi in cui

γ_{1} + γ_{2} = 1, γ_{2}α = 1/2, γ_{2}β = 1/2.

Metodo di Heun (α = 1):

u_{n+1} = u_{n} + \frac{h}{2}(f (x_{n}, u_{n}) + f (x_{n} + h, u_{n} + hf (x_{n}, u_{n}))). Metodo di Eulero modificato (α = 1/2):

u_{n+1} = u_{n} + hf (x_{n} + \frac{h}{2}, u_{n} + \frac{h}{2}f (x_{n}, u_{n})).

Alvise Sommariva  Equazioni differenziali ordinarie.  55/ 61

Runge-Kutta

La scelta particolare di un maggior numero di vincoli permette, con qualche fatica, di calcolare un metodo di ordine 4. Posto

u_{n+1} = u_{n} + h · \frac{f_{1} + 2f_{2} + 2f_{3} + f_{4}}{6}

ove

f_{1} = f (t_{n}, u_{n})
f_{2} = f (t_{n} + \frac{h}{2}, u_{n} + \frac{h · f_{1}}{2} )
f_{3} = f (t_{n} + \frac{h}{2}, u_{n} + \frac{h}{2}f_{2})
f_{4} = f (t_{n} + h, u_{n} + h · f_{3})  (44)

si ottiene un metodo del quart’ordine.

Alvise Sommariva  Equazioni differenziali ordinarie.  56/ 61

Metodi di Eulero esplicito in Matlab

Di seguito descriviamo i codici Matlab di Eulero esplicito e implicito, eseguiti 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

Metodo di Eulero implicito in Matlab

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

Metodo di Eulero implicito in Matlab

Osserviamo che

▶ la routine in questione risolve l’equazione del metodo di

Eulero implicito

z = u_{n} + hf (x_{n+1}, z)

tramite l’applicazione del metodo di punto fisso. ▶ quale valore iniziale del metodo di punto fisso si sceglie il

valore fornito da Eulero esplicito.

Esercizio 1. Eulero esplicito

Si risolva utilizzando il metodo di Eulero esplicito con passi h uguali rispettivamente a 0.2, 0.1, 0.05

\{ \begin{matrix}y′(x) = (cos (y(x)))^{2}, 0 ≤ x ≤ 10 \\ y(0) = 0 \end{matrix}  (45)

la cui soluzione è Y (x) = tan^{−1}(x).

Si calcolino (o plottino in scala semi-logaritmica)

▶ l’errore assoluto rispetto la soluzione esatta; ▶ l’errore relativo rispetto la soluzione esatta.

Per selezionati valori di x, ad esempio x = 10, si osservi il rapporto per cui l’errore decresce quando h è dimezzato.

Esercizio 2. Crank-Nicolson

Modificando opportunamente il metodo di Eulero implicito, si implementi il metodo di Crank-Nicolson.

Alvise Sommariva  Equazioni differenziali ordinarie.  61/ 61

Esercizio 3. Assoluta stabilitá

Approssimare per λ = −100 il valore assunto dalla soluzione \overline{y} del problema di Cauchy

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

nel punto x = 100.

A tal proposito

▶ Si utilizzano i metodi di Eulero esplicito e Eulero implicito con

passi h = 0.1, h = 0.05, h = 0.02, h = 0.01, h = 0.001. ▶ Al variare di h verificare quando |1 + hλ| < 1. ▶ Si osservi che la mancata convergenza a 0 di Eulero implicito

è dovuta al metodo di Eulero esplicito come predittore (si veda il codice di Eulero esplicito). Cosa succede se invece la si risolve con il merodo di Newton?

Esercizio 4. Metodi di Adams

▶ Implementare i metodi di Adams-Bashforth e Moulton per

p = 3. ▶ Si consideri il problema di Cauchy

\{ \begin{matrix}y′(x) = (cos (y(x)))^{2}, 0 ≤ x ≤ 10 \\ y(0) = 0 \end{matrix}  (47)

la cui soluzione è Y (x) = tan^{−1}(x).

Approssimare con tali metodi la soluzione di (41) nell’intervallo [0, 10], in N punti equispaziati x_{n}, con N = 2, 4, 8, 16, 32, descrivendone con un plot l’errore assoluto ∥y_{n} − u_{n}∥_{∞} (dove y_{n} = y(x_{n}) e u_{n} l’approssimazione fornita dal metodo numerico).

Quali valori iniziali si considerino quelli propri della soluzione esatta Y (x) = tan^{−1}(x) nei punti richiesti.

Alvise Sommariva  Equazioni differenziali ordinarie.  63/ 61

Esercizio Facoltativo: Confronto di metodi.

Si risolva utilizzando il metodo di Eulero esplicito, implicito e quello dei trapezi, con passi h uguali rispettivamente a 0.5, 0.1, 0.01

\{ \begin{matrix}y′(x) = λy(x) + (1 − λ) cos(t) − (1 + λ)sin(t), 0 ≤ x ≤ 5 \\ y(0) = 1 \end{matrix}

(48)

Sapendo che la soluzione è Y (x) = sin(x) + cos(x) (indipendentemente da λ), si calcoli l’errore assoluto rispetto la soluzione esatta nei punti x = 1, x = 2, x = 3, x = 4, x = 5.

Bibliografia

K.E. Atkinson, W. Han, D.E. Steward Numerical Solution of Ordinary Differential Equations, Wiley, (2009). V. Comincioli, Analisi Numerica, metodi modelli applicazioni, Mc Graw-Hill, 1990.

Alvise Sommariva  Equazioni differenziali ordinarie.  65/ 61