Equazioni differenziali ordinarie

Alvise Sommariva

Università degli Studi di Padova

Dipartimento di Matematica

3 maggio 2016

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.

Esempi

Esempio

\{ \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

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

Teoremi di Cauchy

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(t_{0}) = y_{0} \end{matrix}  (5)

ha una e una sola soluzione y in [a, b], per un arbitrario y_{0} ∈ ℝ^{n}. Inoltre y dipende con continuità da t_{0} e y_{0}.

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 (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  8/ 60

Dimostrazione. 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 u_{k} = y(x_{k}), u_{k+1} = u_{k} + h · f (x_{k}, u_{k}) = u_{k} + (x_{k+1} − x_{k})y′(x_{k}), in virtù del fatto che y′(x_{k}) = f (x_{k}, y(x_{k})),

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  9/ 60

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

e quindi

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

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

Alvise Sommariva  Introduzione  10/ 60

Eulero implicito

Nota. (Risoluzione equazioni nonlineari) Evidentemente ad ogni iterazione si richiede di risolvere un’eqz. 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}  (15)

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  12/ 60

Analisi convergenza

Definizione (Convergenza )

Supponiamo di analizzare il problema di Cauchy nell’intervallo compatto I = [x_{0}, x_{fin}]. 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 da un metodo numerico nei punti x_{k}.

Un tale metodo si dice convergente se

∥y^{(h_{N})} − u^{(h_{N})}∥_{∞} ≤ C(h_{N})

con C(h_{N}) infinitesima rispetto ad h_{N} quando h_{N} tende a 0 (o equivalentemente N → +∞). Se C(h_{N}) = O(h_{N}^{p}) per qualche p > 0 allora si dice che il metodo converge con ordine p.

Alvise Sommariva  Introduzione  13/ 60

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}. 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/ 60

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

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  Introduzione  15/ 60

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  16/ 60

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

Consideriamo il metodo di Eulero esplicito, con ascisse equispaziate per un prefissato passo h. Sia y soluzione del problema di Cauchy, x_{n} = x_{0} + nh 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}).

Alvise Sommariva  Introduzione  17/ 60

Analisi convergenza

Da

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

il metodo converge se entrambi i termini a secondo membro di (16) convergono a 0 per h → 0, cioè

∥y^{(h)} − \overline{u}^{(h)}∥_{∞} ≤ C_{1}(h) ,

e

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

con C_{1}(h), C_{2}(h) infinitesime rispetto ad h, quando h tende a 0.

Alvise Sommariva  Introduzione  18/ 60

Step 1: mostriamo che

∥y^{(h)} − \overline{u}^{(h)}∥_{∞} ≤ C_{1}(h)

Se la derivata seconda di y esiste ed è continua allora abbiamo per il teorema di Weierstrass

max_{x∈I}|y′′(x)| ≤ M

e quindi

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

in cui ξ_{n} ∈ (x_{n−1}, x_{n}). Visto che M non dipende da n deduciamo che

∥y^{(h)} − \overline{u}^{(h)}∥_{∞} = max_{n}|y(x_{n}) − \overline{u}_{n}| ≤ Mh^{2}/2,

e quindi posto C_{1}(h) = Mh^{2}/2 abbiamo mostrato l’asserto.

Alvise Sommariva  Introduzione  19/ 60

Analisi convergenza

Nota. Si osservi che da

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

deduciamo che, sotto queste ipotesi,

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

e quindi Eulero esplicito risulta consistente di ordine 1.

Alvise Sommariva  Introduzione  20/ 60

Analisi convergenza Eulero esplicito, caso Lipschitziano

Step 2: per ottenere la convergenza dobbiamo provare

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

con C_{2}(h) → 0 per h → 0. Se f è L-Lipschitziana (rispetto al secondo argomento) allora

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

e otteniamo che

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

Alvise Sommariva  Introduzione  21/ 60

Analisi convergenza Eulero esplicito, caso Lipschitziano

Infatti da

\overline{u}_{n} = u_{n−1} + h · f (x_{n−1}, u_{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}.  (17)
Alvise Sommariva  Introduzione  22/ 60

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}| |e_{n}| ≤ hτ(h) + |\overline{u}_{n} − u_{n}| ≤ hτ(h) + (1 + hL)|e_{n−1}|

ricaviamo, essendo e_{0} = 0 poichè il metodo di Eulero esplicito parte da y_{0} = y(x_{0}), che

|\overline{u}_{n} − u_{n}| ≤ (1 + hL)|e_{n−1}| ≤ 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}|
= (1 + (1 + hL) + … + (1 + hL)^{n−1})h|τ(h)|.  (18)
Alvise Sommariva  Introduzione  23/ 60

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

|\overline{u}_{n} − u_{n}| ≤ (1 + (1 + hL) + … + (1 + hL)^{n−1})h|τ(h)| concludiamo

|\overline{u}_{n} − u_{n}| ≤ (1 + (1 + hL) + … + (1 + hL)^{n−1})h|τ(h)|
≤ \frac{(1 + hL)^{n} − 1}{hL}h|τ(h)|
= \frac{(1 + hL)^{n} − 1}{L}|τ(h)|
Alvise Sommariva  Introduzione  24/ 60

Analisi convergenza Eulero esplicito, caso Lipschitziano

Notiamo che per γ > 0

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

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ì, essendo x_{0} < x_{n} ≤ x_{fin}, τ(n) ≤ Mh/2,

|u_{n} − \overline{u}_{n}| ≤ \frac{(1 + hL)^{n} − 1}{L}|τ(h)| ≤ \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, visto che \frac{exp (L(x_{fin}−x_{0}))−1}{L} · \frac{Mh}{2} non dipende da n, ∥u^{(h)} − \overline{u}^{(h)}∥_{∞} = max_{n}|u_{n} − \overline{u}_{n}| ≤ \frac{exp (L(x_{fin} − x_{0})) − 1}{L} · \frac{Mh}{2}.

Alvise Sommariva  Introduzione  25/ 60

Analisi convergenza Eulero esplicito, caso Lipschitziano

Posto

C_{2}(h) = \frac{exp (L(x_{fin} − x_{0})) − 1}{L} · \frac{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, in virtù di quanto visto nello step 1 e nello step 2, è convergente.

Analisi convergenza Eulero esplicito, caso dissipativo

Step 2 (variante): 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  27/ 60

Analisi convergenza Eulero esplicito, caso dissipativo

Da

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

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

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  29/ 60

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, in virtù di quanto visto nello step 1 e nello step 2 (variante), è convergente anche nel caso in cui f sia L-dissipativa.

Alvise Sommariva  Introduzione  30/ 60

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

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(t_{n}) + ∫_{x_{n}}^{x_{n+1}} y′(x)dx = y(x_{n}) + ∫_{x_{n}}^{x_{n+1}} f (x, y(x))dx
≈ y(x_{n}) + (h/2)(f (x_{n}, y(x_{n})) + f (x_{n+1}, y(x_{n+1})))

in quanto, posto G_{y}(x) = f (x, y(x)), dalla formula dei trapezi,

∫_{x_{n}}^{x_{n+1}} f (x, y(x))dx = ∫_{x_{n}}^{x_{n+1}} G_{y}(x)dx ≈ (h/2)(G_{y}(x_{n}) + G_{y}(x_{n+1}))
= (h/2)(f (x_{n}, y(x_{n})) + f (x_{n+1}, y(x_{n+1}))).
Alvise Sommariva  Introduzione  31/ 60

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

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

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

Alvise Sommariva  Introduzione  32/ 60

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

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  33/ 60

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/|λ|.

Definizione (Problema stiff)

Un problema di Cauchy

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

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.

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

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.

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  Introduzione  38/ 60

Regione di assoluta stabilità di Eulero implicito

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

Alvise Sommariva  Introduzione  39/ 60

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

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  40/ 60

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  41/ 60

Regione di assoluta stabilità di Crank-Nicolson

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

Alvise Sommariva  Introduzione  42/ 60

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

è in generale complicato. Se invece consideriamo

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

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

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

con

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

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

Alvise Sommariva  Introduzione  43/ 60

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

e y_{ϵ} quella di

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

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

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

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  44/ 60

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, … (33)

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

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 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  45/ 60

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.

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

Alvise Sommariva  Introduzione  47/ 60

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

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

Facoltativo. 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  Introduzione  50/ 60

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

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

Facoltativo. 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  53/ 60

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

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

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})  (36)
Alvise Sommariva  Introduzione  56/ 60

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

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 (34), (35), (36), necessita

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

Alvise Sommariva  Introduzione  57/ 60

Facoltativo. 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  58/ 60

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

si ottiene un metodo del quart’ordine.

Alvise Sommariva  Introduzione  59/ 60

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.