Metodi di Eulero esplicito, implicito e Crank-Nicolson

Un'equazione differenziale ordinaria descrive come una grandezza evolve nel tempo: il raffreddamento di un pezzo, la velocità di un veicolo con attrito, la corrente in un circuito. Le equazioni risolubili in forma chiusa (Equazioni differenziali a variabili separabiliy' = a(t) b(y) con a continua su I e b continua su J. Gli zeri ȳ di b danno le soluzioni costanti y = ȳ (integrali singolari). Le altre si trovano dividendo per b(y) e integrando: ∫dy/b(y) = ∫a(t)dt + c, poi si ricava y. Se b è C¹ la soluzione del problema di Cauchy è unica, ma può esistere solo su un intervallo più piccolo di I; se b è solo continua può non essere unica (y' = ∛y).Equazioni differenziali a variabili separabili →, Equazioni differenziali lineari del primo ordiney' + a(t)y = b(t) con a, b continue su un intervallo I. Le soluzioni sono definite su tutto I e il problema di Cauchy ha una sola soluzione. Integrale generale: y = e^(−A(t))[c + ∫e^(A(t)) b(t) dt], con A primitiva di a; per l'omogenea (b = 0) y = c e^(−A(t)). Soluzione del problema con y(t0) = x0: stessa formula con A(t) = ∫ da t0 a t di a(s) ds e integrale da t0 a t. Si dimostra moltiplicando per e^(A(t)).Equazioni differenziali lineari del primo ordine →) sono l'eccezione: già $y'=\sin(ty)$ non ha formula. I metodi numerici costruiscono la soluzione passo dopo passo su una griglia di istanti, partendo dal dato iniziale. Qui si trattano i tre metodi più semplici, che mostrano tutte le idee (consistenza, ordine, metodi impliciti) e sono la base dei solutori usati in pratica.

Prerequisiti: Equazioni differenziali - introduzione e problema di CauchyUn'equazione differenziale ordinaria di ordine n ha come incognita una funzione y(t) e lega y⁽ⁿ⁾ a t, y, y', …, y⁽ⁿ⁻¹⁾. L'integrale generale è l'insieme delle soluzioni, definite su un intervallo; di solito sono infinite e dipendono da n costanti. Il problema di Cauchy aggiunge n condizioni iniziali in un istante t0 (y(t0), …, y⁽ⁿ⁻¹⁾(t0)) e sotto opportune ipotesi ha una sola soluzione. y⁽ⁿ⁾ = f(t) si risolve integrando n volte.Equazioni differenziali - introduzione e problema di Cauchy →, Formula di Taylor con resto di LagrangeSe f ha n+1 derivate continue su [a,b], l'errore tra f e il suo polinomio di Taylor si scrive esattamente come f^(n+1)(c)(x-x0)^(n+1)/(n+1)! per un c opportuno: serve a stimare l'errore su tutto un intervallo e a scegliere il grado del polinomio.Formula di Taylor con resto di Lagrange →, Formule di quadratura - trapezi, Cavalieri-Simpson ed estrapolazione di RichardsonUna formula di quadratura approssima $\int_a^b f,dx$ con $\sum w_if(x_i)$ ottenuta integrando il polinomio che interpola $f$ nei nodi. Trapezio $\frac{b-a}2[f(a)+f(b)]$ (grado di esattezza 1, errore $-\frac{(b-a)^3}{12}f''(\xi)$), punto medio $(b-a)f(m)$ (grado 1, errore $+\frac{(b-a)^3}{24}f''(\xi)$), Cavalieri-Simpson $\frac{b-a}6[f(a)+4f(m)+f(b)]$ (grado 3, errore $-\frac{(b-a)^5}{2880}f^{(4)}(\xi)$). Composite con passo $h=\frac{b-a}n$: trapezi errore $-\frac{(b-a)h^2}{12}f''(\xi)$ (ordine 2, rapporto tra errori 4 se $h$ si dimezza), Simpson errore $-\frac{(b-a)h^4}{180}f^{(4)}(\xi)$ (ordine 4, rapporto 16, $n$ pari). Estrapolazione di Richardson: se $Q(h)=I+Ch^p+\dots$ allora $I\approx\frac{2^pQ(h/2)-Q(h)}{2^p-1}$; per i trapezi $p=2$ e $\frac{4T_{h/2}-T_h}3=S_{h/2}$; ripetuta dà la tabella di Romberg. Stima dell'errore a posteriori: $I-Q(h/2)\approx\frac{Q(h/2)-Q(h)}{2^p-1}$.Formule di quadratura - trapezi, Cavalieri-Simpson ed estrapolazione di Richardson → (i metodi sono quadrature); il caso dei sistemi è in Stabilità assoluta, problemi stiff e sistemi di equazioni differenzialiSul problema test $y'=\lambda y$ ($\mathrm{Re},\lambda<0$, soluzione che decade) i metodi danno $u_{n+1}=R(z),u_n$ con $z=h\lambda$: EE $R=1+z$, EI $R=\frac1{1-z}$, CN $R=\frac{1+z/2}{1-z/2}$. La regione di stabilità assoluta è ${z:|R(z)|<1}$: per EE il disco di centro $-1$ e raggio $1$ (per $\lambda<0$ reale: $h<\frac2{|\lambda|}$), per EI l'esterno del disco di centro $1$ e raggio $1$, per CN il semipiano $\mathrm{Re},z<0$. EI e CN sono A-stabili (nessun limite su $h$ per $\mathrm{Re},\lambda<0$); per $z\to-\infty$ è $R_{EI}\to0$ (smorza) e $R_{CN}\to-1$ (oscilla con smorzamento lento). Un problema è stiff quando ha componenti che decadono a velocità molto diverse (rapporto $\frac{\max|\mathrm{Re}\lambda_i|}{\min|\mathrm{Re}\lambda_i|}$ grande): l'accuratezza permetterebbe $h$ grande, ma la stabilità di un metodo esplicito impone $h<\frac2{\max|\lambda_i|}$. Sistemi $\mathbf y'=A\mathbf y$: EE $\mathbf u_{n+1}=(I+hA)\mathbf u_n$, EI $(I-hA)\mathbf u_{n+1}=\mathbf u_n$, CN $(I-\frac h2A)\mathbf u_{n+1}=(I+\frac h2A)\mathbf u_n$; stabilità se $h\lambda_i$ sta nella regione per ogni autovalore. Oscillatore ($\lambda=\pm i\omega$): $|R_{EE}|=\sqrt{1+\omega^2h^2}>1$ (energia cresce), $|R_{EI}|<1$ (smorza), $|R_{CN}|=1$ (conserva). Equazioni di ordine $m$ si riducono a sistemi del primo ordine.Stabilità assoluta, problemi stiff e sistemi di equazioni differenziali →.

1. Il problema di Cauchy e la sua ben posizione

$$y'(t)=f(t,y(t)),\quad t\in[t_0,T],\qquad y(t_0)=y_0 .$$ Se $f$ è continua e lipschitziana in $y$, $|f(t,y_1)-f(t,y_2)|\le L\,|y_1-y_2|$ (basta $\left|\frac{\partial f}{\partial y}\right|\le L$), la soluzione esiste ed è unica e dipende con continuità dai dati: due soluzioni con dati iniziali $y_0$ e $y_0+\delta$ differiscono al più di $\delta\,e^{L(t-t_0)}$. Il problema è ben posto: piccoli errori sul dato (o accumulati dal metodo) non esplodono in tempo finito, anche se possono crescere esponenzialmente. $L$ ha le dimensioni di $\frac1{\text{tempo}}$ ed è la scala dei tempi caratteristici del problema.

Esempio. $y'=t-y$, $y(0)=1$: $f(t,y)=t-y$ ha $\frac{\partial f}{\partial y}=-1$, $L=1$. L'equazione lineare ha soluzione $y(t)=t-1+2e^{-t}$ (somma della soluzione dell'omogenea $Ce^{-t}$ e di una particolare $t-1$). Servirà come banco di prova perché si conosce la soluzione esatta.

2. Tre modi di discretizzare

Si fissa il passo $h$ e i nodi $t_n=t_0+nh$. Si indica con $u_n$ l'approssimazione di $y(t_n)$. Integrando l'equazione su $[t_n,t_{n+1}]$: $$y(t_{n+1})=y(t_n)+\int_{t_n}^{t_{n+1}}f(t,y(t))\,dt .$$ Il metodo dipende da come si approssima l'integrale.

Eulero esplicito (EE). Rettangolo con l'altezza a sinistra, $\int\approx h\,f(t_n,u_n)$: $$u_{n+1}=u_n+h\,f(t_n,u_n).$$ Eulero implicito (EI). Rettangolo con l'altezza a destra, $\int\approx h\,f(t_{n+1},u_{n+1})$: $$u_{n+1}=u_n+h\,f(t_{n+1},u_{n+1}).$$ Crank-Nicolson (CN). Regola dei trapezi, $\int\approx\frac h2\big[f(t_n,u_n)+f(t_{n+1},u_{n+1})\big]$: $$u_{n+1}=u_n+\frac h2\big[f(t_n,u_n)+f(t_{n+1},u_{n+1})\big].$$

Interpretazione geometrica di EE: da $(t_n,u_n)$ si segue la tangente (pendenza $f(t_n,u_n)$) per un tratto $h$. È lo sviluppo di Taylor $y(t_n+h)\approx y(t_n)+hy'(t_n)$ troncato al primo ordine. Equivalentemente, si sostituisce $y'(t_n)$ con il rapporto incrementale in avanti.

Esplicito e implicito. In EE $u_{n+1}$ compare solo a sinistra: basta valutare $f$. In EI e CN compare anche dentro $f$: a ogni passo si deve risolvere un'equazione per $u_{n+1}$.

• Se $f$ è lineare in $y$, $f(t,y)=-\lambda y+g(t)$, si risolve a mano. EI: $u_{n+1}(1+h\lambda)=u_n+h\,g(t_{n+1})$. CN: $u_{n+1}\big(1+\tfrac{h\lambda}2\big)=u_n\big(1-\tfrac{h\lambda}2\big)+\tfrac h2\big[g(t_n)+g(t_{n+1})\big]$.
• Se $f$ è non lineare, si usa l'iterazione di punto fisso $u^{(k+1)}=u_n+h\,f(t_{n+1},u^{(k)})$, che converge se $hL<1$ (contrazione, si veda Metodo di bisezione e iterazioni di punto fissoBisezione: se f è continua e f(a)f(b) < 0 si dimezza l'intervallo tenendo la metà con il cambio di segno; dopo n passi |c_n − α| ≤ (b−a)/2^(n+1), cioè un bit di precisione per iterazione (convergenza lineare, sempre garantita). Punto fisso x = g(x): se g manda [a,b] in sé con |g'| ≤ L < 1 (contrazione) il punto fisso è unico, l'iterazione x_{k+1} = g(x_k) converge da ogni x_0 e |x_k − α| ≤ L^k/(1−L)·|x_1−x_0|; l'errore si moltiplica per circa g'(α) a ogni passo. Una stessa equazione ha molte riscritture x = g(x), alcune convergenti e altre no.Metodo di bisezione e iterazioni di punto fisso →), oppure il metodo di Newton per $G(u)=u-u_n-h\,f(t_{n+1},u)=0$, partendo da $u_n$ (converge quadraticamente e non richiede $hL<1$, perché $G'(u)=1-h\,\frac{\partial f}{\partial y}$ resta lontano da zero anche con $hL$ grande quando $\frac{\partial f}{\partial y}<0$: per questo si usa per i problemi stiff).

Esempio (a mano). $y'=t-y$, $y(0)=1$, $h=0{,}2$ (i dati esatti sono $y(0{,}2)=0{,}837462$, $y(0{,}4)=0{,}740640$, $y(0{,}6)=0{,}697623$).

EE: $u_{n+1}=u_n+0{,}2\,(t_n-u_n)$.

• $u_1=1+0{,}2(0-1)=0{,}8$; $u_2=0{,}8+0{,}2(0{,}2-0{,}8)=0{,}68$; $u_3=0{,}68+0{,}2(0{,}4-0{,}68)=0{,}624$.

EI: $u_{n+1}=u_n+0{,}2\,(t_{n+1}-u_{n+1})\Rightarrow u_{n+1}=\dfrac{u_n+0{,}2\,t_{n+1}}{1{,}2}$.

• $u_1=\frac{1+0{,}2\cdot0{,}2}{1{,}2}=\frac{1{,}04}{1{,}2}=0{,}866667$; $u_2=\frac{0{,}866667+0{,}08}{1{,}2}=0{,}788889$; $u_3=\frac{0{,}788889+0{,}12}{1{,}2}=0{,}757407$.

CN: $u_{n+1}=u_n+0{,}1\,\big[(t_n-u_n)+(t_{n+1}-u_{n+1})\big]\Rightarrow u_{n+1}=\dfrac{0{,}9\,u_n+0{,}1\,(t_n+t_{n+1})}{1{,}1}$.

• $u_1=\frac{0{,}9+0{,}1\cdot0{,}2}{1{,}1}=\frac{0{,}92}{1{,}1}=0{,}836364$; $u_2=\frac{0{,}9\cdot0{,}836364+0{,}1\cdot0{,}6}{1{,}1}=\frac{0{,}812727}{1{,}1}=0{,}738843$; $u_3=\frac{0{,}9\cdot0{,}738843+0{,}1\cdot1}{1{,}1}=0{,}695417$.

Errori a $t=0{,}6$: EE $0{,}0736$ (per difetto), EI $0{,}0598$ (per eccesso), CN $0{,}0022$. EE e EI sbagliano in direzioni opposte (sono i rettangoli a sinistra e a destra di una funzione monotona), e CN, la loro media in un certo senso, è molto più precisa.

3. Errore di troncamento e consistenza

Definizione. L'errore di troncamento locale $\tau_n(h)$ è il residuo che si ottiene inserendo la soluzione esatta nello schema, diviso per $h$. Il metodo è consistente se $\tau_n(h)\to0$ per $h\to0$; ha ordine $p$ se $\tau(h)=O(h^p)$.

Con Taylor (soluzione regolare $y$, $f(t,y(t))=y'(t)$):

• EE: $\tau_n=\dfrac{y(t_{n+1})-y(t_n)}h-y'(t_n)=\dfrac{h}{2}y''(\xi_n)$ perché $y(t_n+h)=y+hy'+\frac{h^2}2y''(\xi)$. Ordine $1$.
• EI: $\tau_n=\dfrac{y(t_{n+1})-y(t_n)}h-y'(t_{n+1})=-\dfrac h2y''(\xi_n)$ (sviluppando attorno a $t_{n+1}$). Ordine $1$.
• CN: $\tau_n=\dfrac{y(t_{n+1})-y(t_n)}h-\dfrac{y'(t_n)+y'(t_{n+1})}2=-\dfrac{h^2}{12}y'''(\xi_n)$: è l'errore della regola dei trapezi (Formule di quadratura - trapezi, Cavalieri-Simpson ed estrapolazione di RichardsonUna formula di quadratura approssima $\int_a^b f,dx$ con $\sum w_if(x_i)$ ottenuta integrando il polinomio che interpola $f$ nei nodi. Trapezio $\frac{b-a}2[f(a)+f(b)]$ (grado di esattezza 1, errore $-\frac{(b-a)^3}{12}f''(\xi)$), punto medio $(b-a)f(m)$ (grado 1, errore $+\frac{(b-a)^3}{24}f''(\xi)$), Cavalieri-Simpson $\frac{b-a}6[f(a)+4f(m)+f(b)]$ (grado 3, errore $-\frac{(b-a)^5}{2880}f^{(4)}(\xi)$). Composite con passo $h=\frac{b-a}n$: trapezi errore $-\frac{(b-a)h^2}{12}f''(\xi)$ (ordine 2, rapporto tra errori 4 se $h$ si dimezza), Simpson errore $-\frac{(b-a)h^4}{180}f^{(4)}(\xi)$ (ordine 4, rapporto 16, $n$ pari). Estrapolazione di Richardson: se $Q(h)=I+Ch^p+\dots$ allora $I\approx\frac{2^pQ(h/2)-Q(h)}{2^p-1}$; per i trapezi $p=2$ e $\frac{4T_{h/2}-T_h}3=S_{h/2}$; ripetuta dà la tabella di Romberg. Stima dell'errore a posteriori: $I-Q(h/2)\approx\frac{Q(h/2)-Q(h)}{2^p-1}$.Formule di quadratura - trapezi, Cavalieri-Simpson ed estrapolazione di Richardson →) applicato a $y'$, diviso per $h$. Ordine $2$.

4. Convergenza e ordine

L'errore globale è $e_n=y(t_n)-u_n$. Gli errori locali si accumulano a ogni passo, e il problema li amplifica (di un fattore al più $e^{L(t-t_0)}$). Per controllare l'accumulo serve la stabilità, nel senso generale visto in Condizionamento di un problema e stabilità di un algoritmoL'errore finale ha tre origini (dati, troncamento, arrotondamento). Un problema è ben posto se la soluzione esiste, è unica e dipende con continuità dai dati; il numero di condizionamento K = |x f'(x)/f(x)| misura di quanto un errore relativo sui dati si amplifica nel risultato (si perdono circa log10 K cifre). La stabilità riguarda invece l'algoritmo: un algoritmo instabile amplifica gli errori anche su un problema ben condizionato (ricorrenza I_n = 1 − n I_{n−1}); un problema mal condizionato non si salva con nessun algoritmo. Uno schema numerico converge se è consistente e stabile (derivata numerica: troncamento contro arrotondamento).Condizionamento di un problema e stabilità di un algoritmo →.

Definizione (stabilità, o zero-stabilità). Un metodo è stabile se, perturbando lo schema con errori $\delta_n$ a ogni passo (e $\delta_0$ sul dato iniziale), la soluzione numerica perturbata $z_n$ dista da quella non perturbata $u_n$ al più $|z_n-u_n|\le C\big(|\delta_0|+\sum_k h|\delta_k|\big)$, con $C$ indipendente da $h$ su un intervallo fissato $[t_0,T]$.

Per EE, EI e CN con $f$ lipschitziana questa proprietà vale sempre (per $h$ abbastanza piccolo): a ogni passo una differenza tra due soluzioni numeriche cresce al più di un fattore $1+hL$ (o $\frac1{1-hL}$ per EI), e dopo $n\le\frac{T-t_0}h$ passi il fattore totale è $\le e^{L(T-t_0)}$, che non dipende da $h$. È una proprietà del metodo per $h\to0$; la stabilità assoluta (a $h$ fissato) è un'altra cosa e si tratta nella nota sui problemi stiff.

Teorema (convergenza). Se $f$ è lipschitziana in $y$ (costante $L$) e il metodo è consistente di ordine $p$ e stabile, allora $\max_n|e_n|\le C\,h^p$ per $h\to0$, con $C$ dipendente da $L$, dall'intervallo e dalle derivate di $y$.

L'idea è quella del teorema di equivalenza: la soluzione esatta soddisfa lo schema a meno del residuo $\tau_n$, cioè è una soluzione numerica perturbata con $\delta_n=\tau_n$; la stabilità dà $|e_n|\le C\sum_kh|\tau_k|\le C\,(T-t_0)\max|\tau|=O(h^p)$.

Dimostrazione per EE. Per la soluzione esatta $y(t_{n+1})=y(t_n)+h\,f(t_n,y(t_n))+h\,\tau_n$, mentre $u_{n+1}=u_n+h\,f(t_n,u_n)$. Sottraendo, $$e_{n+1}=e_n+h\big[f(t_n,y(t_n))-f(t_n,u_n)\big]+h\tau_n\ \Longrightarrow\ |e_{n+1}|\le(1+hL)|e_n|+h|\tau_n| .$$ Con $|\tau_n|\le\frac{Mh}2$, $M=\max|y''|$, e $e_0=0$, iterando: $$|e_n|\le\frac{Mh}{2}\,h\sum_{k=0}^{n-1}(1+hL)^k=\frac{Mh}{2L}\big[(1+hL)^n-1\big]\le\frac{Mh}{2L}\Big(e^{L(t_n-t_0)}-1\Big),$$ usando $1+hL\le e^{hL}$ e $nh=t_n-t_0$. $\blacksquare$

Nella dimostrazione i due ingredienti sono separati: $|\tau_n|\le\frac{Mh}2$ è la consistenza, il fattore $(1+hL)^n\le e^{L(t_n-t_0)}$, limitato indipendentemente da $h$, è la stabilità.

Esempio. Per $y'=t-y$: $L=1$, $y''=2e^{-t}$, $M=2$. A $t=1$ con $h=0{,}2$: $|e|\le\frac{2\cdot0{,}2}{2}(e-1)=0{,}344$. L'errore vero è $0{,}0804$: la stima è corretta, ma pessimistica di un fattore $4$.

Tabella di convergenza per $y'=t-y$ in $t=1$ ($y(1)=2e^{-1}=0{,}735759$). Il rapporto è tra l'errore con il passo precedente (doppio) e quello attuale:

I rapporti $2$ e $4$ confermano gli ordini $1$ e $2$. Per ottenere un errore di $10^{-4}$ EE richiede circa $h\approx2{,}7\cdot10^{-4}$ (3700 passi), CN solo $h\approx0{,}04$ (25 passi): CN costa di più per passo (un'equazione implicita, qui lineare e quindi gratis) ma molto meno complessivamente.

Stima dell'ordine da tre soluzioni con passi $h,\frac h2,\frac h4$ quando la soluzione esatta non è nota: $2^p\approx\frac{u_h-u_{h/2}}{u_{h/2}-u_{h/4}}$, come per la quadratura.

5. Implementazione

import numpy as np

f = lambda t, y: t - y                      # y' = t - y, y(0) = 1
esatta = lambda t: t - 1 + 2 * np.exp(-t)   # soluzione esatta

def eulero_esplicito(f, t0, y0, h, n):
    t, y = t0, y0
    for _ in range(n):
        y += h * f(t, y); t += h
    return y

def eulero_implicito(t0, y0, h, n):         # f lineare: u_{n+1}(1+h) = u_n + h t_{n+1}
    t, y = t0, y0
    for _ in range(n):
        t += h; y = (y + h * t) / (1 + h)
    return y

def crank_nicolson(t0, y0, h, n):           # u_{n+1}(1+h/2) = u_n(1-h/2) + h/2 (t_n + t_{n+1})
    t, y = t0, y0
    for _ in range(n):
        y = (y * (1 - h / 2) + h / 2 * (t + t + h)) / (1 + h / 2); t += h
    return y

print("tre passi con h = 0.2:", [round(eulero_esplicito(f, 0, 1, 0.2, k), 6) for k in (1, 2, 3)],
      [round(eulero_implicito(0, 1, 0.2, k), 6) for k in (1, 2, 3)],
      [round(crank_nicolson(0, 1, 0.2, k), 6) for k in (1, 2, 3)])
print("esatta                :", [round(float(esatta(0.2 * k)), 6) for k in (1, 2, 3)])
print("\n    h     EE: errore  rapp.   EI: errore  rapp.   CN: errore  rapp.")
old = [None] * 3
for n in (5, 10, 20, 40, 80):
    h = 1 / n
    errs = [abs(esatta(1) - v) for v in (eulero_esplicito(f, 0, 1, h, n), eulero_implicito(0, 1, h, n), crank_nicolson(0, 1, h, n))]
    rapp = [f"{o / e:5.2f}" if o else "    -" for o, e in zip(old, errs)]
    print(f"{h:6.4f}  {errs[0]:.3e} {rapp[0]}   {errs[1]:.3e} {rapp[1]}   {errs[2]:.3e} {rapp[2]}")
    old = errs

Output:

tre passi con h = 0.2: [0.8, 0.68, 0.624] [0.866667, 0.788889, 0.757407] [0.836364, 0.738843, 0.695417]
esatta                : [0.837462, 0.74064, 0.697623]

    h     EE: errore  rapp.   EI: errore  rapp.   CN: errore  rapp.
0.2000  8.040e-02     -   6.800e-02     -   2.463e-03     -
0.1000  3.840e-02  2.09   3.533e-02  1.92   6.138e-04  4.01
0.0500  1.879e-02  2.04   1.802e-02  1.96   1.533e-04  4.00
0.0250  9.294e-03  2.02   9.102e-03  1.98   3.832e-05  4.00
0.0125  4.623e-03  2.01   4.575e-03  1.99   9.580e-06  4.00

In MATLAB ode45 (Runge-Kutta adattivo) e ode23, ode15s (per problemi stiff) risolvono [t,y] = ode45(@(t,y) t-y, [0 1], 1); i metodi di questa nota si scrivono con un ciclo for. In Python scipy.integrate.solve_ivp(f, (0,1), [1]) è l'equivalente.

Errori da evitare

• Valutare $f$ al tempo sbagliato: EE usa $(t_n,u_n)$, EI usa $(t_{n+1},u_{n+1})$, CN entrambi.
• Trattare EI come esplicito: $u_{n+1}$ compare dentro $f$; con $f$ lineare si isola $u_{n+1}$, con $f$ non lineare servono punto fisso o Newton.
• Confondere errore locale e globale: il primo è $O(h^{p+1})$ per passo (o $\tau=O(h^p)$ per unità di tempo), il secondo $O(h^p)$ sull'intervallo; si perde una potenza di $h$ perché si fanno $\sim\frac1h$ passi.
• Credere che l'ordine dica tutto: un metodo di ordine alto può essere instabile per $h$ non abbastanza piccolo (si veda la nota sulla stabilità).
• Applicare la stima $\frac{Mh}{2L}(e^{LT}-1)$ come previsione realistica: è un limite superiore spesso pessimistico, specialmente quando $T$ è grande.
• Dimenticare l'arrotondamento nel rapporto degli errori: i rapporti $2$ e $4$ valgono solo per $h$ abbastanza piccolo e prima che gli errori di macchina dominino.

Esercizi svolti: Esercizio - Eulero esplicito a mano e confronto con la soluzione esatta del raffreddamento, Esercizio - Eulero implicito e Crank-Nicolson su un'equazione lineare con forzante, Esercizio - Eulero implicito e Crank-Nicolson per un'equazione non lineare con Newton, Esercizio - ordine di convergenza dalla tabella degli errori e stima del passo.

Domande d'esame

1. Ricavare i metodi di Eulero esplicito, implicito e Crank-Nicolson dall'integrale di $f$ su $[t_n,t_{n+1}]$. Traccia. Rettangolo a sinistra, a destra, trapezi.
2. Che cosa significa che un metodo ha ordine $p$? Calcolare l'errore di troncamento di Eulero esplicito. Traccia. $\tau=O(h^p)$; Taylor: $\tau_n=\frac h2y''(\xi)$.
3. Dimostrare che l'errore globale di Eulero esplicito è $O(h)$. Traccia. $|e_{n+1}|\le(1+hL)|e_n|+h|\tau|$, $|\tau|\le\frac{Mh}2$, somma geometrica: $|e_n|\le\frac{Mh}{2L}(e^{L(t_n-t_0)}-1)$.
4. Qual è il vantaggio e quale il costo dei metodi impliciti? Traccia. Stabilità migliore (nota seguente); a ogni passo un'equazione non lineare da risolvere (punto fisso con $hL<1$ o Newton).