Equazioni differenziali ordinarie

Alvise Sommariva

Università degli Studi di Padova

Dipartimento di Matematica

9 maggio 2017

Problema di Cauchy

Problema. (Problema di Cauchy ) Si determini la funzione y tale che

\{ \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}) ∈ Ω. Nota. Di seguito supporremo che tale problema abbia una sola soluzione.

Alvise Sommariva  Introduzione  2/ 74

Esempi

Esempio (Equazione differenziale ordinaria)

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

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

Esempio (Sistema di equazioni differenziali ordinarie)

\{ \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)), è 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)).

Alvise Sommariva  Introduzione  3/ 74

Teoremi di Cauchy

Teorema (Teorema di Picard-Lindelöf)

Sia D = [t_{0} − ε, t_{0} + ε] × ℝ. Dato il problema ai valori iniziali

\{ \begin{matrix}y′(t) = f (t, y(t)), t ∈ [t_{0} − ε, t_{0} + ε] \\ y(t_{0}) = y_{0} \end{matrix}

se f : D → ℝ è una funzione lipschitziana in y cioè

|f (t, y_{1}) − f (t, y_{2})| ≤ K∥y_{1} − y_{2}∥, (t, y_{1}), (t, y_{2}) ∈ D

ed è continua in t allora per qualche ε > 0, allora esiste un’unica soluzione y al problema ai valori iniziali sull’intervallo [t_{0} − ε, t_{0} + ε].

Teorema (Cauchy in piccolo, [2, p.7]) Si supponga che

D sia un aperto connesso di ℝ × ℝ^{n}; f : D → ℝ sia una funzione continua in D; (t_{0}, y_{0}) un punto interno a D; la funzione f verifichi la condizione di Lipschitz

∥f (t, y_{1}) − f (t, y_{2})∥ ≤ K∥y_{1} − y_{2}∥, (t, y_{1}), (t, y_{2}) ∈ D. per qualche K ≥ 0.

Allora esiste un unica funzione y definita su un intervallo [t_{0} − α, t_{0} + α], per qualche α > 0 tale che

\{ \begin{matrix}y′(t) = f (t, y(t)), t ∈ [t_{0} − α, t_{0} + α] \\ y(t_{0}) = y_{0} \end{matrix}  (4)

Teorema (Cauchy in grande, [4, p.331]) Sia f : [a, b] × ℝ^{n} → ℝ tale che

f (t, y), t ∈ [a, b], y ∈ ℝ^{n} continua nella prima variabile, soddisfi rispetto alla seconda variabile la condizione di Lipschitz

∥f (t, y_{1}) − f (t, y_{2})∥ ≤ K∥y_{1} − y_{2}∥, t ∈ [a, b], y_{1}, y_{2} ∈ ℝ^{n}.

per qualche norma vettoriale ∥ · ∥ e K ≥ 0. Allora il problema di Cauchy

\{ \begin{matrix}y′(t) = f (t, y(t)), t ∈ [a, b] \\ y(a) = y_{0} \end{matrix}  (5)

ha una e una sola soluzione y in [a, b], per un arbitrario y_{0} ∈ ℝ^{n}.

Nota.

Nel teorema di Cauchy in piccolo si fornisce un teorema di esistenza e unicità locale.

Nel teorema di Cauchy in grande si fornisce un teorema di esistenza e unicità globale.

Nota. (Soluzioni multiple)

Esistono casi in cui il problema di Cauchy ha soluzioni multiple. Ad esempio, il problema

\{ \begin{matrix}y′(x) = 2\sqrt{y(x)}, x > 0 \\ y(x_{0}) = 0 \end{matrix}  (6)

ha per soluzioni, per α ∈ ℝ^{+}

y(x) = \{ \begin{matrix}y(x) = (x − α)^{2}, x ≥ α \\ y(x) = 0, altrimenti. \end{matrix}  (7)

Metodo di Eulero esplicito

Indichiamo con I(x, \overline{x}) il più piccolo intervallo 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}) (8)

cioè

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

Metodo di Eulero esplicito

Metodo (Eulero esplicito, 1768-1770) Tale metodo consiste nell’approssimare y(x_{k+1}) con u_{k+1} dove

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

ove h = x_{k+1} − x_{k}.

Teorema (Errore Eulero esplicito) Se l’unica soluzione y del problema di Cauchy in [a, b] è sufficientemente regolare in [x_{k}, x_{k+1}] ⊆ [a, b] e u_{k} = y(x_{k}) allora

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

per qualche ξ ∈ (x_{k}, x_{k+1}).

Alvise Sommariva  Introduzione  9/ 74

Metodo di Eulero esplicito

Dimostrazione. Essendo

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

u_{k+1} = u_{k} + h · f (x_{k}, u_{k}) = u_{k} + h · f (x_{k}, y(x_{k}))
= y(x_{k}) + y′(x_{k})(x_{k+1} − x_{k}),

(12) e quindi si ha per qualche ξ in (x_{k}, x_{k+1})

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

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

− (y(x_{k}) + y′(x_{k})(x_{k+1} − x_{k}))

Alvise Sommariva  Introduzione  10/ 74

Metodo di Eulero implicito

Similmente al caso di Eulero esplicito, da

y(x) ≈ y(\overline{x}) + f (\overline{x}, y(\overline{x}))(x − \overline{x})

se poniamo 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}), (13)

e quindi

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

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

u_{k+1} = u_{k} + h · f (x_{k+1}, u_{k+1}), u_{0} = y(x_{0})  (15)

Alvise Sommariva  Introduzione  11/ 74

Eulero implicito

Nota. (Risoluzione equazioni nonlineari) Evidentemente ad ogni iterazione si richiede di risolvere un’equazione nonlineare nella variabile z del tipo

z = u_{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.

Osservazione

Nota. (Metodi di Eulero e sistemi di equazioni differenziali) Il problema di Cauchy è definito da

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

dove f è a valori in ℝ^{n} e definita in un sottoinsieme Ω di ℝ × ℝ^{n}, con (x_{0}, y_{0}) ∈ Ω. Nei casi multivariati con n > 1, indipendentemente dalla derivazione dovuta alla espansione di Taylor, i metodi

Eulero esplicito: u_{k+1} = u_{k} + h · f (x_{k}, u_{k}), u_{0} = y(x_{0})

Eulero implicito: u_{k+1} = u_{k} + h · f (x_{k+1}, u_{k+1}), u_{0} = y(x_{0}), possono essere utilizzati, in quanto la sequenza è ben definita.

Alvise Sommariva  Introduzione  13/ 74

Linear Multistep Methods

Definizione (Linear Multistep Methods )

Sono i metodi che approssimano y(x_{k}) con u_{k}, dove x_{k} = x_{0} + kh e

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

Nota. Si osservi che

Sono metodi per cui se voglio calcolare u_{n+1} necessitano u_{n−p}, …, u_{n} (si chiamano metodi a p + 1 passi). Siccome inizialmente si dispone solo di u_{0}, per calcolare u_{1}, …, u_{p}, necessari ad innescare il metodo, si utilizzano di solito altre strategie (come ad esempio applicare più volte metodi a un passo, cioè in cui p = 0, come ad esempio Eulero implicito).

Se b_{−1} = 0 allora il metodo è esplicito, altrimenti è implicito.

Alvise Sommariva  Introduzione  14/ 74

Analisi convergenza

Nota. Osserviamo che il metodo di Eulero esplicito

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

e Eulero implicito

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

hanno la forma

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

ponendo rispettivamente Eulero esplicito: p = 0, a_{0} = 1, b_{−1} = 0, b_{0} = 1, Eulero implicito: p = 0, a_{0} = 1, b_{−1} = 1, b_{0} = 0.

Alvise Sommariva  Introduzione  15/ 74

Analisi convergenza

Definizione (Convergenza )

Supponiamo di analizzare il problema di Cauchy nell’intervallo compatto I = [x_{0}, xfin]. Siano x_{s} = x_{0} + sh, con Nh_{N} = x_{fin} − x_{0},

y^{(h_{N})} = {y(x_{k})} la soluzione esatta di un fissato problema di Cauchy,

u^{(h_{N})} = {u_{k}} dove u_{k} è l’approssimazione di y(x_{k}) fornita nei punti x_{k} da un metodo numerico a p + 1 passi.

Un metodo Linear Multistep a p + 1 passi si dice convergente se qualora i passi iniziali u_{0}^{(h_{N})}, …, u_{p}^{(h_{N})} sono tali che

η(h_{N}) = max_{0≤k≤p}|y^{(h_{N})}(x_{0} + kh) − u_{k}^{(h_{N})}| → 0 per h_{N} → 0

si ha che

∥y^{(h_{N})} − u^{(h_{N})}∥_{∞} → 0 per h_{N} → 0.

Se qualora η(h_{N}) ≤ C_{1}h^{q} implica ∥y^{(h_{N})} − u^{(h_{N})}∥_{∞} ≤ C_{2}h^{q}, con C_{1}, C_{2} indipendenti da h, allora il metodo si dice convergente con ordine q.

Definizione (Errori di troncamento, [2], p.112)

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

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

ove u_{n} ≈ y(t_{n}), y soluzione di (1) in t_{n} = x_{0} + nh, e sia

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

La quantità τ_{n}(h) chiama errore locale di troncamento del metodo. La quantità τ(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  Introduzione  17/ 74

Nota. La consistenza rende conto di come la soluzione y del problema di Cauchy verifichi lo schema discreto del metodo linear multistep.

Nota. In queste note ci interesseremo esclusivamente della convergenza e consistenza di metodi linear multistep, ma esistono definizioni che permettono l’analisi di altri metodi non rientranti in questa famiglia.

Nota. Un teorema dovuto a Lax/Dahlquist, mostra che un metodo è convergente se e solo se consistente e stabile. Un metodo Linear Multistep a p + 1 passi si dice (zero-)stabile per un certo problema di Cauchy nel compatto [x_{0}, x_{fin}] se per ogni ϵ > 0 esiste h tale che se h_{N} ≤ h e

max_{k=0,...,p−1}|u_{k}^{(h_{N})} − y(x_{0} + kh_{N})| ≤ ϵ

allora

∥u_{k}^{(h_{N})} − y(x_{0} + kh_{N})∥_{∞} ≤ Kϵ

con K indipendente da h.

Alvise Sommariva  Introduzione  18/ 74

Di seguito intendiamo mostrare, sotto opportune ipotesi, la convergenza del metodo di Eulero esplicito.

Consideriamo il metodo di Eulero esplicito, con x_{n} = x_{0} + nh per un prefissato passo h = (x_{fin} − x_{0})/N. Sia y soluzione del problema di Cauchy, e

\overline{u}_{n} = y(x_{n−1}) + h · f (x_{n−1}, y(x_{n−1})),
u_{n} = u_{n−1} + h · f (x_{n−1}, u_{n−1}).

Osserviamo che la prima sequenza \overline{u}^{(h)} = \{\overline{u}_{n}\}, ha il generico \overline{u}_{n} ottenuto applicando il metodo di Eulero esplicito partendo dal valore assunto dalla soluzione y in ogni x_{n−1}, la seconda sequenza u^{(h)} = {u_{n}}, ha il generico u_{n} ottenuto applicando il metodo di Eulero esplicito partendo dal valore assunto u_{n−1}, che a priori non coincide con y(x_{n−1}).

Da

e_{n} = y(x_{n}) − u_{n} = (y(x_{n}) − \overline{u}_{n}) + (\overline{u}_{n} − u_{n})  (17)

per la disuguaglianza triangolare

|e_{n}| ≤ |y(x_{n}) − \overline{u}_{n}| + |\overline{u}_{n} − u_{n}|.
Alvise Sommariva  Introduzione  19/ 74

Analisi convergenza

Analizziamo il primo termine |y(x_{n}) − \overline{u}_{n}|. Se la derivata seconda di y esiste ed è continua allora abbiamo per il teorema di Weierstrass

max_{x∈[x_{0},x_{fin}]}|y′′(x)| ≤ M

e quindi per ξ_{n} ∈ (x_{n−1}, x_{n}),

|y(x_{n}) − \overline{u}_{n}| = (h^{2}/2) · |y′′(ξ_{n})| ≤ Mh^{2}/2.

Di conseguenza, h|τ_{n}(h)| = max_{n=0,...,N} |y(x_{n}) − \overline{u}_{n}| ≤ Mh^{2}/2.

Alvise Sommariva  Introduzione  20/ 74

Analisi convergenza Eulero esplicito, caso Lipschitziano

Sketch della dimostrazione (supposto y suff. regolare e f continua e L-lipschitziana):

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

2 1 + (1 + hL) + … + (1 + hL)^{n−1} ≤ \frac{exp(L(x_{fin}−x_{0}))}{hL} e (1 + hL)^{n} ≤ exp(L(x_{fin} − x_{0}));

3 conclusione perchè secondo membro infinitesimo per h → 0 e |e_{0}| → 0 per h → 0.

Analisi convergenza Eulero esplicito, caso Lipschitziano

Se f è L-Lipschitziana (rispetto al secondo argomento) allora

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

e otteniamo, ricordando che e_{n} := y(x_{n}) − u_{n},

|\overline{u}_{n} − u_{n}| ≤ (1 + hL)|e_{n−1}|.

Infatti da

\overline{u}_{n} = u_{n−1} + h · f (x_{n−1}, y(x_{n−1})), u_{n} = u_{n−1} + h · f (x_{n−1}, u_{n−1}) ricaviamo dalla L-lipschitzianità

|\overline{u}_{n} − u_{n}| = |(y(x_{n−1}) + h · f (x_{n−1}, y(x_{n−1}))) − (u_{n−1} + h · f (x_{n−1}, u_{n−1}))|
= |(y(x_{n−1}) − u_{n−1}) + h · (f (x_{n−1}, y(x_{n−1})) − f (x_{n−1}, u_{n−1}))|
≤ |y(x_{n−1}) − u_{n−1}| + h|f (x_{n−1}, y(x_{n−1})) − f (x_{n−1}, u_{n−1}))|
≤ |y(x_{n−1}) − u_{n−1}| + hL|y(x_{n−1}) − u_{n−1}| = (1 + hL)|y(x_{n−1}) − u_{n−1}|
= (1 + hL)e_{n−1}.  (18)
Alvise Sommariva  Introduzione  22/ 74

Analisi convergenza Eulero esplicito, caso Lipschitziano

Da e_{n} := y(x_{n}) − u_{n} = (y(x_{n}) − \overline{u}_{n}) + (\overline{u}_{n} − u_{n}),

|y(x_{n}) − \overline{u}_{n}| ≤ h|τ(h)|,
|\overline{u}_{n} − u_{n}| ≤ (1 + hL)|e_{n−1}|

si deduce

|e_{n}| = |(y(x_{n}) − \overline{u}_{n}) + (\overline{u}_{n} − u_{n})|
≤ |y(x_{n}) − \overline{u}_{n}| + |\overline{u}_{n} − u_{n}|
≤ h|τ(h)| + (1 + hL)|e_{n−1}|  (19)

e quindi, ricordando che e_{n} := y(x_{n}) − u_{n}, |e_{n}| ≤ h|τ(h)| + (1 + hL)|e_{n−1}|.

Alvise Sommariva  Introduzione  23/ 74

Analisi convergenza Eulero esplicito, caso Lipschitziano

Quindi essendo |e_{0}| → 0 per h → 0, |e_{n}| ≤ h|τ(h)| + (1 + hL)|e_{n−1}|,

ricaviamo

|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}|
≤ …
≤ (1 + (1 + hL) + … + (1 + hL)^{n−1})h|τ(h)| + (1 + hL)^{n}|e_{0}|

e cioè

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

Alvise Sommariva  Introduzione  24/ 74

Analisi convergenza Eulero esplicito, caso Lipschitziano

Ricordando che

1 + s + … + s^{k} = \frac{(1 − s^{k+1})}{(1 − s)}

posto s = hL deduciamo che 1 + (1 + hL) + … + (1 + hL)^{n−1} = \frac{1 − (1 + hL)^{n}}{1 − (1 + hL)} = \frac{(1 + hL)^{n} − 1}{hL}. Notiamo che per γ > 0

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

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})) ≤ exp (L(x_{fin} − x_{0})).

Alvise Sommariva  Introduzione  25/ 74

Analisi convergenza Eulero esplicito, caso Lipschitziano

Da

1 + (1 + hL) + … + (1 + hL)^{n−1} = \frac{1−(1+hL)^{n}}{1−(1+hL)} = \frac{(1+hL)^{n}−1}{hL}
(1 + hL)^{n} ≤ exp (L(x_{fin} − x_{0})),

concludiamo che

1 + (1 + hL) + … + (1 + hL)^{n−1} = \frac{(1 + hL)^{n} − 1}{hL} ≤ \frac{exp (L(x_{fin} − x_{0})) − 1}{hL}
≤ \frac{exp (L(x_{fin} − x_{0}))}{hL}  (21)

e quindi

(1 + (1 + hL) + … + (1 + hL)^{n−1})h|τ(h)| ≤ \frac{exp (L(x_{fin} − x_{0}))}{hL}h|τ(h)|
= exp (L(x_{fin} − x_{0}))|τ(h)| → 0.
Alvise Sommariva  Introduzione  26/ 74

Analisi convergenza Eulero esplicito, caso Lipschitziano

Di conseguenza da

|e_{n}| ≤ (1 + (1 + hL) + … + (1 + hL)^{n−1})h|τ(h)| + (1 + hL)^{n}|e_{0}|
≤ exp (L(x_{fin} − x_{0}))|τ(h)| + exp (L(x_{fin} − x_{0}))|e_{0}|
≤ exp (L(x_{fin} − x_{0}))\frac{Mh^{2}}{2} + exp (L(x_{fin} − x_{0}))|e_{0}|  (22)

visto che il secondo membro non dipende da n ricaviamo

max_{n=0,...,N} |e_{n}| ≤ exp (L(x_{fin} − x_{0}))\frac{Mh^{2}}{2} + exp (L(x_{fin} − x_{0}))|e_{0}| e quindi il metodo di Eulero esplicito risulta convergente, visto che il secondo membro è infinitesimo.

Alvise Sommariva  Introduzione  27/ 74

Metodo di Crank-Nicolson

Supponiamo di dover risolvere il problema di Cauchy

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

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

y(x_{n+1}) = y(x_{n}) + ∫_{x_{n}}^{x_{n+1}} y′(x)dx
≈ y(x_{n}) + (h/2)(y′(x_{n+1}) + y′(x_{n}))
= y(x_{n}) + (h/2)(f (x_{n}, y(x_{n})) + f (x_{n+1}, y(x_{n+1}))).

Nota. Utilizzando rispettivamente le formule dei rettangoli,

∫_{a}^{b} g(x)dx ≈ (b − a)g(a), ∫_{a}^{b} g(x)dx ≈ (b − a)g(b)

si possono ottenere similmente i metodi di Eulero esplicito e implicito.

Alvise Sommariva  Introduzione  28/ 74

Metodo di Crank-Nicolson

Essendo

y(x_{n+1}) ≈ y(x_{n}) + (h/2)(f (x_{n}, y(x_{n})) + f (x_{n+1}, y(x_{n+1}))) (24)

si introduce il metodo di Crank-Nicolson (1947): Metodo (Crank-Nicolson o dei trapezi ) Tale metodo consiste nell’approssimare y(x_{n+1}) con u_{n+1} dove

u_{n+1} = u_{n} + (h/2)(f (t_{n}, u_{n}) + f (t_{n+1}, u_{n+1}))  (25)

con u_{0} = y_{0}. Si osservi che ad ogni iterazione, essendo il metodo implicito, bisogna risolvere una equazione nonlineare.

Alvise Sommariva  Introduzione  29/ 74

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

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à è composta dagli hλ, per cui il metodo numerico con passo h è tale da avere lo stesso comportamento asintotico del problema di Cauchy (25) relativo al parametro λ (con ℜ(λ) < 0).

Alvise Sommariva  Introduzione  30/ 74

Assoluta stabilità. Una nota.

Nota. Stabilire la stabilità per il problema di Cauchy

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

è in generale complicato. Se invece consideriamo

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

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

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

con

c = y_{0}

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

Alvise Sommariva  Introduzione  31/ 74

Assoluta stabilità. Una nota.

Sia y la soluzione di

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

e y_{ϵ} quella di

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

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

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

la cui soluzione è z_{ϵ} = ϵ · exp (λx). 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à.

Alvise Sommariva  Introduzione  32/ 74

Assoluta stabilità: Eulero esplicito

Definizione (Problema stiff) 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).

Alvise Sommariva  Introduzione  33/ 74

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_{n} è l’approssimazione della soluzione in x_{n} = x_{0} + nh fornita da un metodo numerico a passo h, si desidera sia u_{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/|λ|.

Nel caso di problemi stiff, cioè con ℜ(λ) ≪ 0, si deve scegliere un passo h molto piccolo affinchè la soluzione di Eulero esplicito tenda a 0.

Regione di assoluta stabilità di Eulero esplicito

Figura : In magenta: regione di assoluta stabilità di Eulero esplicito: disco unitario centrato in (−1, 0).

Alvise Sommariva  Introduzione  37/ 74

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, ragionando per induzione, che

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

e quindi che

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

Alvise Sommariva  Introduzione  38/ 74

Assoluta stabilità: Eulero implicito

Essendo ℜ(λ) < 0, si vede facilmente che |\frac{1}{1−hλ}| < 1 per qualsiasi h. Infatti, λ = a + ib con a < 0, ciò è vero se e solo se

\frac{1}{|1 − hλ|} < 1 ⇔ 1 < |1 − hλ| = \sqrt{(1 − a)^{2} + b^{2}}

ovvero, elevando ai quadrati entrambi i membri, se e solo se

1 < (1 − a)^{2} + b^{2}

ovviamente verificata in quanto a < 0, a, b ∈ ℝ. La regione di ass. stabilità è tutto il semipiano negativo ℜ(hλ) < 0.

Nota. Questa proprietà suggerisce di applicare Eulero implicito invece di Eulero esplicito per risolvere numericamente un problema di Cauchy di tipo stiff.

Regione di assoluta stabilità di Eulero implicito

Figura : In magenta: regione di assoluta stabilità di Eulero implicito, semipiano negativo ℜ(hλ) < 0.

Alvise Sommariva  Introduzione  40/ 74

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  Introduzione  41/ 74

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  Introduzione  42/ 74

Regione di assoluta stabilità di Crank-Nicolson

Figura : In magenta: regione di assoluta stabilità di Crank-Nicolson, semipiano negativo ℜ(hλ) < 0.

Alvise Sommariva  Introduzione  43/ 74

Metodi linear multistep (LM)

Posto h > 0, x_{n} = x_{0} + nh, un metodo linear multistep è 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, … (34)

noti i valori u_{k} ≈ y(x_{k}) per k < p e supposto a_{p}, b_{p} ≠ 0.

Nota. Non è difficile vedere che

Eulero esplicito (a_{0} = 1, b_{0} = 1), Eulero implicito (a_{0} = 1, b_{−1} = 1), Crank-Nicolson (a_{0} = 1, b_{−1} = 1/2, b_{0} = 1/2)

sono metodi che hanno questa struttura.

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}} y′(x) dx = 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, (1883)); Se m = 0, γ = 1 si hanno metodi impliciti (Adams-Moulton, (1926)).

Alvise Sommariva  Introduzione  45/ 74

Esempi LM: Adams-Bashforth

Adams-Bashforth (p=0): u_{n+1} = u_{n} + hf_{n}, τ_{n+1} = (h^{2}/2)y^{(2)}(ξ_{n}). Adams-Bashforth (p=1): u_{n+1} = u_{n} + \frac{h}{2} (3f_{n} − f_{n−1}), τ_{n+1} = (5h^{3}/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^{4}/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^{5}/720)y^{(5)}(ξ_{n}). Adams-Bashforth (p=4): u_{n+1} = u_{n} + \frac{h}{720} (1901f_{n} − 2774f_{n−1} + 2616f_{n−2} − 1274f_{n−3} + 251f_{n−4}), τ_{n+1} = (95h^{6}/2888)y^{(5)}(ξ_{n}).

dove τ_{n+1} è l’errore locale di troncamento

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}/2)y^{(2)}(ξ_{n}). Adams-Moulton (p=1): u_{n+1} = u_{n} + \frac{h}{2} (f_{n+1} + f_{n}), τ_{n+1} = (−h^{3}/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^{4}/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^{5}/720)y^{(5)}(ξ_{n}). Adams-Moulton (p=4): u_{n+1} = u_{n} + \frac{h}{720} (251f_{n+1} + 646f_{n} − 264f_{n−1} + 106f_{n−2} − 19f_{n−3}), τ_{n+1} = (−3h^{6}/160)y^{(6)}(ξ_{n}).

Nota. Si noti che sono effettivamente metodi impliciti a p passi.

Consistenza LM

Teorema Un metodo linear multistep

è consistente se e solo se

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

se y ∈ C^{q+1} per q ≥ 1, ha ordine di consistenza 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.

Alvise Sommariva  Introduzione  48/ 74

Convergenza LM

Teorema (Convergenza LM (cf.[1, p.116])) 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;

il metodo sia consistente; a_{j} ≥ 0, j = 0, …, p; il passo di discretizzazione h sia tale che h ≤ 1/(2c), dove c = L ∑_{j=−1}^{p} |b_{j}| e 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).

Se la soluzione è m + 1 volte derivabile con continuità, il metodo ha ordine di consistenza m e gli errori iniziali soddisfano η(h) = O(h^{m}) allora il metodo è di ordine m.

Alvise Sommariva  Introduzione  49/ 74

Nota. Nelle ipotesi del teorema precedente, per ricavare la convergenza chiedavamo:

i dati iniziali siano calcolati in modo che

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

il metodo sia consistente;

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

il passo di discretizzazione h sia tale che h ≤ 1/(2c), dove c = L ∑_{j=−1}^{p} |b_{j}| e L è la costante di Lipschitz di f .

Se prendiamo un metodo di Adams-Bashforth o Adams-Moulton di quelli esposti, sono tutti consistenti con ordine p + 1. Inoltre a_{0} = 1, a_{1} = 0, …, a_{p} = 0. Di conseguenza tutti i metodi esposti sono convergenti. Se la soluzione è p + 2 volte derivabile con continuità, il metodo ha ordine di consistenza p + 1 e gli errori iniziali soddisfano η(h) = O(h^{m}) allora il metodo è di ordine p + 1.

In particolare, nelle ipotesi richieste,

il metodo di Eulero esplicito, avendo ordine di consistenza 1, è convergente con ordine di convergenza 1;

il metodo di Eulero implicito, avendo ordine di consistenza 1, è convergente con ordine di convergenza 1;

il metodo di Crank-Nicolson, avendo ordine di consistenza 2, è convergente con ordine di convergenza 2.

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.

Parte facoltativa

Gli argomenti che seguono sono facoltativi. Servono solo a dare una immagine più completa sul tema dell’analisi numerica per la risoluzione di problemi di Cauchy.

Analisi convergenza Eulero esplicito, caso dissipativo

Se f non è L-Lipschitziana ma è L-dissipativa, cioè

−L ≤ \frac{∂f}{∂y} (ξ) ≤ 0, ξ ∈ Ω := (x_{0}, x_{fin}) × ℝ

abbiamo da \overline{u}_{n} = y(x_{n−1}) + h · f (x_{n−1}, y(x_{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(x_{n−1})) + \frac{∂f}{∂y}(ξ)(u_{n−1} − y(x_{n−1})),

che per qualche ξ ∈ (x_{n−1}, x_{n}) × ℝ ⊆ Ω

u_{n} − \overline{u}_{n} = (u_{n−1} + h · f (x_{n−1}, u_{n−1})) − (y(x_{n−1}) + h · f (x_{n−1}, y(x_{n−1})))
= (u_{n−1} − y(x_{n−1})) + h · (f (x_{n−1}, u_{n−1}) − f (x_{n−1}, y(x_{n−1})))
= (u_{n−1} − y(x_{n−1})) + h \frac{∂f}{∂y} (ξ)(u_{n−1} − y(x_{n−1}))
= (1 + h \frac{∂f}{∂y} (ξ)) · (u_{n−1} − y(x_{n−1})).
Alvise Sommariva  Introduzione  54/ 74

Analisi convergenza Eulero esplicito, caso dissipativo

Da

u_{n} − \overline{u}_{n} = (1 + h\frac{∂f}{∂y} (ξ)) · (u_{n−1} − y(x_{n−1})).  (35)

abbiamo che se 0 < h ≤ 2/L (non restrittivo per mostrare la conv. visto che si studia il comportamento per h → 0) allora da −L ≤ \frac{∂f}{∂y} ≤ 0,

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

e quindi |1 + h\frac{∂f}{∂y}(ξ)| ≤ 1 da cui

|u_{n} − \overline{u}_{n}| = |1 + h\frac{∂f}{∂y} (ξ)| · |u_{n−1} − y(x_{n−1})|
≤ |u_{n−1} − y(x_{n−1})| = |e_{n−1}|.  (36)

Analisi convergenza Eulero esplicito, caso dissipativo

Da e_{n} := y(x_{n}) − u_{n}, |y(x_{n}) − \overline{u}_{n}| ≤ h|τ(h)|, |\overline{u}_{n} − u_{n}| ≤ |e_{n−1}|, |e_{n}| ≤ hτ(h) + |\overline{u}_{n} − u_{n}| ≤ h|τ(h)| + |e_{n−1}|

ricaviamo, essendo h|τ(h)| > 0, e_{0} = 0, n ≤ N, Nh = xfin − x_{0}

|\overline{u}_{n} − u_{n}| ≤ |e_{n−1}| ≤ h|τ(h)| + |e_{n−1}|
≤ h|τ(h)| + (h|τ(h)| + |e_{n−2}|)
= 2h|τ(h)| + |e_{n−2}|
≤ 2h|τ(h)| + (h|τ(h)| + |e_{n−3}|)
= 3h|τ(h)| + |e_{n−3}| ≤ …
≤ nh|τ(h)| + |e_{0}| = nh|τ(h)|
≤ Nh|τ(h)| = (xfin − x_{0})|τ(h)| ≤ \frac{(xfin − x_{0})Mh}{2}.

Alvise Sommariva  Introduzione  56/ 74

Analisi convergenza Eulero esplicito, caso dissipativo

Di conseguenza, visto che (xfin − x_{0}) non dipende da n, essendo per ogni n

|\overline{u}_{n} − u_{n}| ≤ \frac{(xfin − x_{0})Mh}{2}

abbiamo

∥u^{(h)} − \overline{u}^{(h)}∥_{∞} = max_{n}|u_{n} − \overline{u}_{n}| ≤ \frac{(xfin − x_{0})Mh}{2}.

Posto C_{2}^{∗}(h) = (xfin − x_{0})Mh/2 abbiamo che

∥u^{(h)} − \overline{u}^{(h)}∥_{∞} ≤ C_{2}^{∗}(h)

con C_{2}^{∗}(h) → 0 qualora h → 0 e quindi il metodo di Eulero esplicito, è convergente anche nel caso in cui f sia L-dissipativa.

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} = u_{n−1} + ∫_{x_{n−1}}^{x_{n}} y′(x)dx = u_{n−1} + hy′(x_{n}) − h^{2}y′′(ξ)/2
= u_{n−1} + hf (x_{n}, u_{n}) − h^{2}y′′(ξ)/2  (37)

e quindi

τ_{n}(h) = (1/h) · (h^{2}y′′(ξ)/2) = hy′′(ξ)/2.
Alvise Sommariva  Introduzione  58/ 74

Analisi convergenza Eulero implicito, caso Lipschitziano

Ricordiamo che u_{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

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

Come in precedenza

u_{n} − u_{n} = (u_{n} − \overline{u}_{n}) + (\overline{u}_{n} − u_{n})

e quindi per la disuguaglianza triangolare

|u_{n} − u_{n}| ≤ |u_{n} − \overline{u}_{n}| + |\overline{u}_{n} − u_{n}|.

Studiamo separatamente i termini al secondo membro. Osserviamo che

|u_{n} − \overline{u}_{n}| = |u_{n} − (u_{n−1} + hf (t_{n}, y_{n}))| = h|τ_{n}(h)|

Alvise Sommariva  Introduzione  59/ 74

Analisi convergenza Eulero implicito, caso Lipschitziano

Studiamo |\overline{u}_{n} − u_{n}| ricordando che \overline{u}_{n} := u_{n−1} + hf (t_{n}, y_{n}), \overline{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)

|\overline{u}_{n} − u_{n}| = |u_{n−1} + hf (t_{n}, y_{n}) − u_{n−1} − hf (t_{n}, u_{n})|
= |u_{n−1} − u_{n−1} + hf (t_{n}, y_{n}) − hf (t_{n}, u_{n})|
≤ |u_{n−1} − u_{n−1}| + h|f (t_{n}, y_{n}) − f (t_{n}, u_{n})|
≤ |u_{n−1} − u_{n−1}| + hL|y_{n} − u_{n}|  (38)
Alvise Sommariva  Introduzione  60/ 74

Analisi convergenza Eulero implicito, caso Lipschitziano

In definitiva, essendo

|y(x_{n}) − u_{n}| ≤ |y(x_{n}) − \overline{u}_{n}| + |\overline{u}_{n} − u_{n}|,

|y(x_{n}) − \overline{u}_{n}| = |y(x_{n}) − (y(x_{n−1}) + hf (t_{n}, y_{n}))| = h|τ_{n}(h)|,

|\overline{u}_{n} − u_{n}| ≤ |y(x_{n−1}) − u_{n−1}| + hL|y(x_{n}) − u_{n}|,

deduciamo che

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

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

|e_{n}| = |u_{n} − u_{n}| ≤ (1/(1 − hL))(h|τ_{n}(h)| + |y(x_{n−1}) − u_{n−1}|)
= \frac{h|τ_{n}(h)|}{1 − hL} + \frac{|e_{n−1}|}{1 − hL}  (39)
Alvise Sommariva  Introduzione  61/ 74

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

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

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

e quindi

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

Alvise Sommariva  Introduzione  63/ 74

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  Introduzione  64/ 74

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

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

Alvise Sommariva  Introduzione  65/ 74

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  Introduzione  67/ 74

Runge-Kutta

Consideriamo un metodo numerico, che diremo di Runge-Kutta di ordine 2, 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}).  (44)

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

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})  (46)
Alvise Sommariva  Introduzione  70/ 74

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

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 (43), (44), (45), necessita

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

Alvise Sommariva  Introduzione  71/ 74

Runge-Kutta

Vediamo alcuni metodi in cui

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

Metodo di Heun (α = 1), (scoperto nel 1900):

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  Introduzione  72/ 74

La scelta particolare di un maggior numero di vincoli permette, con qualche fatica, di calcolare un metodo di ordine 4 (scoperto da Kutta nel 1901). 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})  (48)

si ottiene un metodo del quart’ordine.

Alvise Sommariva  Introduzione  73/ 74

Bibliografia

K. Atkinson e W. Han, Elementary Numerical Analysis, Wiley, (2004). 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. W. Gautschi, Numerical Analysis, Birkhaüser, second edition, 2012.