# Funzione rhs della ODE (y' = f(t, y))
f <- function(t, y) { -2 * y + t^2 }
# Soluzione esatta del problema per confronto
exact_solution <- function(t,y0) {
t*(t-1)/2 + 1/4 + exp(-2*t)*(y0 - 1/4)
}28 Soluzione numerica
In molti casi, trovare una soluzione analitica di un problema a valore iniziale (IVP) è impraticabile, a causa della complessità degli integrali o delle funzioni coinvolte. Si ricorre quindi a metodi numerici per approssimare la soluzione. Questi metodi si applicano sia ai singoli IVP sia ai sistemi di equazioni differenziali ordinarie (ODE). Poiché le equazioni differenziali di ordine superiore possono essere riformulate come sistemi del primo ordine, i metodi numerici per questi sistemi sono centrali nella risoluzione delle ODE.
28.1 Passo Temporale
Sia I = [t_0, t_0 + T] un intervallo temporale, che suddividiamo in n sottointervalli, detti anche passi temporali. Ogni sottointervallo I_k = [t_{k-1}, t_k] ha ampiezza h_k = t_k - t_{k-1}, con t_n = t_0 + T. Se i passi sono uniformi, h_k = h per ogni k. La scelta di h è cruciale, poiché influisce sulla precisione e sulla stabilità del metodo numerico.
28.2 Principali metodi numerici per IVP
Metodi a Un Passo
Nei metodi a un passo, la soluzione al tempo t_{k+1} dipende soltanto da quella calcolata al tempo t_k. In altre parole, l’approssimazione y_{k+1} si basa unicamente sull’informazione relativa a y_k. Il più semplice di questi metodi è quello di Eulero1.
Metodi a Più Passi
I metodi a più passi utilizzano le informazioni di diversi passi precedenti per calcolare la nuova approssimazione y_{k+1}. Sono più complessi e spesso più accurati dei metodi a un passo, poiché tengono conto di una parte più ampia della storia della soluzione. Un esempio noto è il metodo di Adams-Bashforth.
Metodi Espliciti
Un metodo esplicito (o di integrazione in avanti) calcola la soluzione al passo successivo usando soltanto valori già noti. A partire dalla condizione iniziale y_0 al tempo t_0, le approssimazioni successive sono determinate da una relazione del tipo:
y_{k+1} = y_k + h_k F(t_k, y_k),
dove F rappresenta l’incremento della soluzione al passo k. I metodi espliciti sono meno costosi dal punto di vista computazionale e semplici da implementare, ma possono richiedere passi temporali piccoli per garantire la stabilità della soluzione.
Metodi Impliciti
Nei metodi impliciti (o di integrazione all’indietro), invece, la funzione F dipende dal valore della soluzione al tempo t_{k+1}. A ogni passo occorre quindi risolvere un’equazione per determinare y_{k+1}:
y_{k+1} = y_k + h_k F(t_{k+1}, y_{k+1}).
Questi metodi hanno un costo computazionale maggiore, poiché richiedono la risoluzione di sistemi di equazioni, spesso in modo iterativo; possono però offrire una maggiore stabilità, soprattutto per i problemi rigidi. Un esempio è il metodo di Eulero implicito.
28.3 Il Metodo di Eulero
Il metodo di Eulero è uno dei metodi numerici più semplici per la risoluzione degli IVP. Si tratta di un metodo a un passo che approssima la soluzione tramite una forma semplificata dell’espansione in serie di Taylor2. Consideriamo un problema di valore iniziale:
\begin{cases} y'(t) = f(t,y), \\ y(t_0) = y_0. \end{cases}
L’espansione di Taylor di y(t) attorno al punto iniziale t_0 è:
y(t_0 + h) = y_0 + h y'(t_0) + O(h^2),
dove
y'(t_0) = f(t_0, y_0).
Metodo di Eulero Esplicito
Trascurando i termini di ordine superiore, otteniamo la formula del metodo di Eulero esplicito:
y_{k+1} = y_k + h f(t_k, y_k).
Metodo di Eulero Implicito
Nel metodo di Eulero implicito, si usa un’analisi simile, ma si considera l’approssimazione alla derivata in t_{k+1} anziché in t_k. L’equazione risultante è:
y_{k+1} = y_k + h f(t_{k+1}, y_{k+1}),
che richiede la risoluzione implicita di y_{k+1} a ogni passo, rendendo il metodo più stabile ma computazionalmente più costoso.
28.4 Convergenza Metodo di Eulero
Convergenza Metodo di Eulero Esplicito
Consideriamo la soluzione esatta y(t_k) al punto t_k e introduciamo l’errore
\epsilon_k = y(t_k) - y_k,
che rappresenta la differenza tra la soluzione esatta e quella approssimata al passo k. Per analizzare l’errore del metodo di Eulero esplicito, confrontiamo l’espansione di Taylor della soluzione esatta con il passo di avanzamento del metodo.
Taylor e metodo di Eulero
L’espansione di Taylor della soluzione esatta y(t) attorno a t_k è:
y(t_{k+1}) = y(t_k + h) = y(t_k) + y'(t_k) h + \dfrac{y''(\zeta_k)}{2} h^2,
dove \zeta_k \in [t_k, t_{k+1}] è un punto intermedio.
Il metodo di Eulero esplicito approssima la soluzione al passo successivo come:
y_{k+1} = y_k + h f(t_k, y_k).
Sottraendo l’espressione di y_{k+1} dall’espansione di Taylor di y(t_{k+1}), otteniamo:
\begin{aligned} y(t_{k+1}) - y_{k+1} &= y(t_k) - y_k \\[1em] &+h \big(f(t_k, y(t_k)) - f(t_k, y_k)\big)\\[1em] &+ \dfrac{y''(\zeta_k)}{2} h^2. \end{aligned}
Questo può essere riscritto come:
\epsilon_{k+1} = \epsilon_k + h \big(f(t_k, y(t_k)) - f(t_k, y_k)\big) + \dfrac{y''(\zeta_k)}{2} h^2.
Stima dell’Errore
Assumendo che f(t,y) sia lipschitziana rispetto a y, ovvero:
|f(t_k, y(t_k)) - f(t_k, y_k)| \leq L |y(t_k) - y_k| = L |\epsilon_k|,
e che esista una costante M_2 tale che
|y''(t)| \leq M_2 \quad \text{per ogni } t \in [a,b],
possiamo stimare l’errore al passo successivo come:
|\epsilon_{k+1}| \leq (1 + hL) |\epsilon_k| + \dfrac{M_2}{2} h^2. \tag{28.1}
usando questa diseguaglianza possiamo dimostrare la convergenza del metodo di Eulero.
Lemma 28.1 Sia
E_h = \max_{k=0}^{n} |\epsilon_k|
l’errore massimo e
h = \dfrac{b - a}{n}
il passo temporale. Per il metodo di Eulero esplicito si ha la seguente stima:
E_h \leq C h, \qquad C = \dfrac{2L e^{(b-a)}}{L M_2}
dove L è la costante di Lipschitz3 di f(t,y) rispetto a y, e M_2 è un limite superiore per il modulo della derivata seconda della soluzione.
Dimostrazione. Ponendo A = 1 + hL e B = M_2 \dfrac{h^2}{2} nella stima dell’errore data in 28.1, otteniamo che l’errore si propaga con la seguente limitazione:
|\varepsilon_{k+1}| \leq A |\varepsilon_k| + B,
dove \varepsilon_0 = y(t_0) - y_0 = 0. Utilizzando la ricorrenza, possiamo scrivere:
\begin{aligned} |\varepsilon_1| & \leq A |\varepsilon_0| + B = B, \\ |\varepsilon_2| & \leq A |\varepsilon_1| + B = AB + B = (A + 1)B, \\ |\varepsilon_3| & \leq A |\varepsilon_2| + B = A(A + 1)B + B = (A^2 + A + 1)B, \\ |\varepsilon_4| & \leq A |\varepsilon_3| + B = A(A^2 + A + 1)B + B = (A^3 + A^2 + A + 1)B, \\ & \vdots \\ |\varepsilon_{k+1}| & \leq (A^k + A^{k-1} + \cdots + A^2 + A + 1)B. \end{aligned}
La sommatoria A^k + A^{k-1} + \cdots + A^2 + A + 1 è una serie geometrica, che può essere espressa come:
A^k + A^{k-1} + \cdots + A^2 + A + 1 = \dfrac{A^{k+1} - 1}{A - 1}.
Di conseguenza, abbiamo:
E_h = \max |\varepsilon_k| \leq \dfrac{A^n - 1}{A - 1} B \leq \dfrac{A^n}{A - 1} B.
Sostituendo A = 1 + hL e B = M_2 \dfrac{h^2}{2}, otteniamo:
E_h \leq \dfrac{A^n}{A - 1} B = \dfrac{(1+hL)^n}{hL} M_2 \dfrac{h^2}{2}.
Utilizzando la disuguaglianza valida per z \geq 0:
1+z \leq e^z, \tag{28.2}
risultante dalla serie di Taylor, otteniamo:
E_h \leq \dfrac{(e^{hL})^n}{2L} M_2 h = \dfrac{e^{nhL}}{2L} M_2 h = \left(\dfrac{e^{(b-a)L}}{2L} M_2\right) h.
Convergenza Metodo di Eulero Implicito
Analogamente al metodo esplicito, consideriamo la soluzione esatta y(t_k) al punto t_k e l’errore \epsilon_k = y(t_k) - y_k.
Taylor e metodo di Eulero Implicito
L’espansione di Taylor della soluzione esatta y(t) attorno a t_{k+1} è:
y(t_{k+1}) = y(t_k) + y'(t_{k+1}) h - \dfrac{y''(\zeta_k)}{2} h^2,
dove \zeta_k \in [t_k, t_{k+1}].
Nel metodo di Eulero implicito, l’approssimazione al passo successivo si basa su:
y_{k+1} = y_k + h f(t_{k+1}, y_{k+1}),
dove la funzione f(t, y) è valutata al passo successivo (t_{k+1}, y_{k+1}).
Sottraendo l’espressione di y_{k+1} dalla serie di Taylor, otteniamo:
\begin{aligned} y(t_{k+1}) - y_{k+1} & = y(t_k) - y_k \\[1em] & + h \big( f(t_{k+1}, y(t_{k+1})) - f(t_{k+1}, y_{k+1}) \big) \\[1em] & - \dfrac{y''(\zeta_k)}{2} h^2, \end{aligned}
che può essere riscritta come:
\epsilon_{k+1} = \epsilon_k + h \big( f(t_{k+1}, y(t_{k+1})) - f(t_{k+1}, y_{k+1}) \big) - \dfrac{y''(\zeta_k)}{2} h^2.
Stima dell’Errore
Assumiamo che f(t, y) sia lipschitziana rispetto a y, cioè:
|f(t_{k+1}, y(t_{k+1})) - f(t_{k+1}, y_{k+1})| \leq L |y(t_{k+1}) - y_{k+1}| = L |\epsilon_{k+1}|,
e che esista una costante M_2 tale che |y''(t)| \leq M_2 per ogni t \in [a,b].
Otteniamo quindi la seguente stima per l’errore al passo successivo:
|\epsilon_{k+1}| \leq |\epsilon_k| + h L |\epsilon_{k+1}| + \dfrac{M_2}{2} h^2.
Riorganizzando i termini, otteniamo:
|\epsilon_{k+1}| \leq \dfrac{|\epsilon_k| + \dfrac{M_2}{2} h^2}{1 - hL}.
Da questa disuguaglianza possiamo dimostrare la convergenza del metodo di Eulero implicito.
Lemma 28.2 (Convergenza di Eulero implicito) Sia
E_h = \max_{k=0}^{n} |\epsilon_k|
l’errore massimo e sia
h = \dfrac{b - a}{n}
il passo temporale. Per il metodo di Eulero implicito si ha la seguente stima:
E_h \leq C h, \qquad C = \dfrac{\mathrm{e}^{2(b-a)L}}{L} M_2,
dove L è la costante di Lipschitz di f(t, y) rispetto a y e M_2 è una stima per il massimo modulo della derivata seconda della soluzione.
Dimostrazione. Riprendiamo la stima dell’errore per |\epsilon_{k+1}|:
|\epsilon_{k+1}| \leq \dfrac{|\epsilon_k| + \dfrac{M_2}{2} h^2}{1 - hL}.
Poniamo
A = \dfrac{1}{1 - hL} \quad\text{e}\quad B = A\dfrac{M_2}{2} h^2
in modo da riscrivere la relazione come:
|\epsilon_{k+1}| \leq A |\epsilon_k| + B.
Facendo lo stesso schema del Lemma sulla convergenza del metodo di Eulero esplicito, otteniamo:
E_h \leq \dfrac{A^{n+1} - 1}{A - 1} B.
Assumendo h sufficientemente piccolo affinché Lh\leq 1/2, possiamo scrivere:
\begin{aligned} A &= \dfrac{1}{1 - hL} = \sum_{k=0}^\infty (hL)^k \\ &= 1 + hL \sum_{k=1}^\infty (hL)^k \leq 1 + 2hL, \\ A-1 &= \dfrac{1}{1 - hL} - 1 \geq hL, \end{aligned}
e quindi, usando anche 28.2:
\dfrac{A^n - 1}{A - 1} \leq \dfrac{\left(1 + 2hL\right)^n - 1}{hL} \leq\dfrac{\mathrm{e}^{2nhL} - 1}{hL} \leq\dfrac{\mathrm{e}^{2(b-a)L}}{hL}.
Poiché A \leq 2 per hL\leq 1/2, si ha
B = A\dfrac{M_2}{2} h^2 \leq M_2 h^2,
e unendo le disuguaglianze, otteniamo la stima cercata. \square
Generalizzazione
Il metodo per dimostrare la convergenza di schemi numerici, come quelli di Eulero esplicito e implicito, può essere esteso a metodi generici sotto certe ipotesi. In questa sezione, vediamo come formalizzare un metodo numerico generico (sia esplicito che implicito) e dimostrare la sua convergenza.
Consideriamo l’equazione differenziale ordinaria (ODE):
\begin{cases} y'(t) = f(t,y), \\ y(a) = a. \end{cases}
Un metodo esplicito a un passo può essere scritto genericamente come:
y_0 = a, \qquad y_{k+1} = \Phi( t_k, y_k; h )
dove \Phi( t, y; h ) è la funzione che rappresenta il passo di avanzamento. Ad esempio, nel caso di Eulero esplicito:
\Phi( t, y; h ) = y + h f(t,y)
Funzione di avanzamento per altri schemi
Vediamo ora come scrivere la funzione di avanzamento \Phi( t, y; h ) per altri schemi numerici.
Metodo di Collatz
Il metodo di Collatz può essere descritto come segue:
\begin{aligned} y_{k+1/2} & = y_k + \dfrac{h}{2}f(t_k,y_k) \\ y_{k+1} & = y_k + h f\left(t_k+\dfrac{h}{2},\; y_{k+1/2}\right) \end{aligned}
Questo schema può essere riarrangiato in una forma compatibile con \Phi( t, y; h ):
y_{k+1} = y_k + h f\left(t_k+\dfrac{h}{2},\; y_k + \dfrac{h}{2}f(t_k,y_k)\right)
per cui la funzione di avanzamento per Collatz diventa
\Phi( t, y; h ) = y + h f\left(t+\dfrac{h}{2},\; y + \dfrac{h}{2}f(t,y)\right)
Metodo di Heun
Il metodo di Heun è definito da:
\begin{aligned} \widetilde{y}_{k+1} & = y_k + hf(t_k,y_k) \\ y_{k+1} & = y_k + \dfrac{h}{2}\left( f(t_k,y_k) + f(t_{k+1},\widetilde{y}_{k+1}) \right) \end{aligned}
Riarrangiando anche questo schema, otteniamo:
y_{k+1} = y_k + \dfrac{h}{2}\Big( f(t_k,y_k) + f\big(t_k+h,y_k + hf(t_k,y_k)\big) \Big)
per cui la funzione di avanzamento per Heun diventa
\Phi( t, y; h ) = y + \dfrac{h}{2}\Big(f(t,y) + f\big(t+h,y+ hf(t,y)\big)\Big)
Propagazione dell’errore
Sostituendo la soluzione esatta nella funzione di avanzamento, otteniamo:
y(t+h) = \Phi( t, y(t) ) + \color{blue} \tau(t,h)
dove \color{blue} \tau(t,h) è l’errore locale di troncamento. Ora calcoliamo l’errore di troncamento per alcuni metodi.
Errore di Troncamento per Eulero
Per Eulero esplicito, abbiamo: \begin{aligned} \tau(t,h) &= y(t+h)-\Phi( t, y(t), h ) \\ & = {\color{red}y(t+h)}-{\color{blue}y(t)-h f(t,y(t))} \\ & = {\color{red}\cancel{y(t)}+\cancel{y'(t)h}+y''(\zeta)h^2/2} -{\color{blue}\cancel{y(t)}-\cancel{h y'(t)}} \\ &= \dfrac{h^2}{2}y''(\zeta) \end{aligned}
dove \zeta \in [t, t+h].
Errore di Troncamento per Collatz
Per il metodo di Collatz:
\begin{aligned} \tau(t,h) &= y(t+h)-\Phi( t, y(t), h ) \\ & = {\color{red}y(t+h)}-{\color{blue}y(t) - h f\left(t+\dfrac{h}{2},\; y(t) + \dfrac{h}{2}f(t,y(t))\right)} \\ & = {\color{red}\cancel{y(t)}+y'(t)h+y''(t)h^2/2+y'''(\zeta)h^3/6} \\ &-{\color{blue}\cancel{y(t)}-h\left(f(t,y(t))+ \dfrac{h}{2}\dfrac{\partial f(t,y(t))}{\partial t}+ \dfrac{h}{2}\dfrac{\partial f(t,y(t))}{\partial t}f(t,y(t))+ \mathcal{O}(h^2) \right)} \\ & = {\color{red}\cancel{f(t,y(t))h}+\dfrac{h^2}{2}\dfrac{\mathrm{d}}{\mathrm{d}t}f(t,y(t))+y'''(\zeta)h^3/6} \\ &{\color{blue}-\cancel{hf(t,y(t))}- \dfrac{h^2}{2}\dfrac{\partial f(t,y(t))}{\partial t}- \dfrac{h^2}{2}\dfrac{\partial f(t,y(t))}{\partial t}f(t,y(t))+ \mathcal{O}(h^3)} \\ & = {\color{red}\dfrac{h^2}{2} \cancel{\left( \dfrac{\partial f(t,y(t))}{\partial t}+ \dfrac{\partial f(t,y(t))}{\partial t}f(t,y(t)) \right)} +y'''(\zeta)h^3/6} \\ &{\color{blue} -\dfrac{h^2}{2} \cancel{\left(\dfrac{\partial f(t,y(t))}{\partial t} +\dfrac{\partial f(t,y(t))}{\partial t}f(t,y(t)) \right)}+ \mathcal{O}(h^3)} \\ &= \mathcal{O}(h^3) \end{aligned}
Analogamente, si può calcolare l’errore per altri metodi, ad esempio Heun, che risulta essere \mathcal{O}(h^2). In generale, se la soluzione esatta non esplode e la funzione f(t, y) è sufficientemente regolare, l’errore locale di troncamento soddisferà:
\color{blue} |\tau(t,h)| \;\leq\; C\, h^{p+1} \tag{28.3}
ove p \geq 1 dipende dal metodo utilizzato.
Convergenza dei metodi a un passo
I conti precedenti si ripetono quasi identici per ogni metodo a un passo. Conviene allora enunciarli una volta sola, in forma generale, e ricavarne una stima esplicita dell’errore globale.
Scriviamo un generico metodo a un passo nella forma incrementale
y_{k+1} = y_k + h\, \Psi(t_k, y_k; h) \tag{28.4}
dove \Psi è detta funzione di incremento. Ad esempio \Psi = f(t_k,y_k) per Eulero esplicito e \Psi = f\big(t_k+\tfrac{h}{2},\,y_k+\tfrac{h}{2}f(t_k,y_k)\big) per Collatz.
Definizione 28.1 (Difetto locale e ordine) Detta y(t) la soluzione esatta, il difetto locale è
d_{k+1} = y(t_{k+1}) - y(t_k) - h\, \Psi\big(t_k, y(t_k); h\big)
Il metodo ha ordine p se, su ogni intervallo limitato,
|d_{k+1}| \leq C\, h^{p+1}
ed è consistente se |d_{k+1}|/h \to 0 per h\to 0.
Si noti la relazione con l’errore locale di troncamento usato sopra: d_{k+1} = \tau(t_k,h). L’errore globale è invece
\epsilon_k = y(t_k) - y_k
La distinzione è essenziale: un difetto \mathcal{O}(h^{p+1}) viene prodotto ad ognuno dei circa T/h passi, quindi ci si attende un errore globale \mathcal{O}(h^{p}) purché le perturbazioni non vengano amplificate troppo. Il teorema seguente rende preciso questo “purché”.
Teorema 28.1 (Convergenza dei metodi a un passo) Si supponga che, per 0<h\leq h_0, la funzione di incremento sia uniformemente lipschitziana rispetto alla seconda variabile,
\big| \Psi(t,u;h) - \Psi(t,v;h) \big| \leq L_{\Psi}\, |u-v| \tag{28.5}
e che il difetto locale soddisfi |d_{k+1}|\leq C h^{p+1}. Allora, su un intervallo fissato [t_0,t_0+T], l’errore globale soddisfa
\boxed{\; |\epsilon_k| \leq \dfrac{C}{L_{\Psi}} \Big( e^{L_{\Psi}(t_k-t_0)} - 1 \Big) h^{p} \;} \tag{28.6}
quando L_{\Psi}>0, e |\epsilon_k| \leq C\,T\,h^{p} nel caso limite L_{\Psi}=0. In particolare il metodo converge con ordine p.
Dimostrazione. Sottraendo il passo numerico 28.4 dal passo esatto si ottiene
\epsilon_{k+1} = \epsilon_k + h\Big[ \Psi\big(t_k,y(t_k);h\big) - \Psi\big(t_k,y_k;h\big) \Big] + d_{k+1}
e quindi, per la 28.5,
|\epsilon_{k+1}| \leq (1 + h L_{\Psi})\, |\epsilon_k| + C h^{p+1} \tag{28.7}
Poiché \epsilon_0 = 0, iterando la 28.7 si ricava
|\epsilon_k| \leq C h^{p+1} \sum_{j=0}^{k-1} (1+hL_{\Psi})^{j}
Per L_{\Psi}>0 la somma geometrica vale \big[(1+hL_{\Psi})^{k}-1\big]/(hL_{\Psi}), da cui
\begin{aligned} |\epsilon_k| &\leq \dfrac{C}{L_{\Psi}} \Big[ (1+hL_{\Psi})^{k} - 1 \Big] h^{p} \\ &\leq \dfrac{C}{L_{\Psi}} \Big( e^{L_{\Psi} k h} - 1 \Big) h^{p} = \dfrac{C}{L_{\Psi}} \Big( e^{L_{\Psi}(t_k-t_0)} - 1 \Big) h^{p} \end{aligned}
dove si è usato 1+\xi \leq e^{\xi} per \xi\geq 0. Se L_{\Psi}=0 la somma vale k e si ottiene direttamente |\epsilon_k|\leq C k h^{p+1} \leq C T h^{p}. \square
Il teorema è il meccanismo elementare dietro la regola
\text{errore locale } \mathcal{O}(h^{p+1}) \quad\Longrightarrow\quad \text{errore globale } \mathcal{O}(h^{p})
e mostra anche a quale prezzo: la costante cresce come e^{L_{\Psi}T}, quindi la stima è utile su intervalli finiti e si degrada esponenzialmente al crescere di T. Si confronti con la 27.7: è la stessa struttura, ed è la stessa ragione.
Generalizzazione (metodi impliciti)
Un metodo implicito a un passo può essere scritto genericamente come:
y_0 = a, \qquad y_{k+1} = \Phi( t_k, y_k, y_{k+1}; h )
dove \Phi(t, y, z; h) è la funzione di avanzamento che coinvolge sia il punto noto y_k che il punto incognito y_{k+1}. Ad esempio, per il metodo di Eulero implicito, la funzione di avanzamento assume la forma:
\Phi( t, y, z; h ) = z + h f(t,z)
dove z = y_{k+1}, ovvero il nuovo valore che stiamo cercando.
Errore Locale di Troncamento
Per definire la precisione del metodo, consideriamo l’errore locale di troncamento. Questo rappresenta la differenza tra la soluzione esatta al passo successivo y(t+h) e la previsione del metodo numerico. L’errore di troncamento per un metodo implicito è dato da:
\tau(t,h) = y(t+h)-\Phi( t, y(t), y(t+h); h )
Errore di troncamento di Eulero implicito
Per il metodo di Eulero implicito, l’errore locale di troncamento si calcola come segue: \begin{aligned} \tau(t,h) &= y(t+h)-\Phi( t, y(t), y(t+h), h ) \\ & = {\color{red}y(t+h)}-{\color{blue}y(t)-h f(t+h,y(t+h))} \\ & = {\color{red}y(t+h)}-{\color{blue}y(t)-hy'(t+h)} \\ & = {\color{red}\cancel{y(t+h)}}-{\color{blue}(\cancel{y(t+h)}-h \cancel{y'(t+h)}+h^2/2 y''(\zeta))-h\cancel{y'(t+h)}} \\ &=-\dfrac{h^2}{2}y''(\zeta) \end{aligned}
dove \zeta \in [t, t+h]. Notiamo che l’errore locale di troncamento per Eulero implicito è:
\tau(t,h) = -\dfrac{h^2}{2}y''(\zeta)
che è dello stesso ordine dell’errore per il metodo di Eulero esplicito, ovvero \mathcal{O}(h^2).
Metodo di Crank-Nicolson
Il metodo di Crank-Nicolson4, basato sulla regola dei trapezi, è un metodo implicito che può essere scritto come: y_{k+1} = y_k + \dfrac{h}{2}\Big( f(t_k,y_k)+ f(t_{k+1},y_{k+1}) \Big)
per cui la funzione di avanzamento prende la forma
\Phi( t, y, z; h ) = y + \dfrac{h}{2}\Big(f(t,y)+f(t+h,z)\Big)
L’errore locale di troncamento per questo metodo si calcola come: \begin{aligned} \tau(t,h) &= y(t+h)-\Phi( t, y(t), y(t+h), h ) \\ &= y(t+h) - y(t) - \dfrac{h}{2}\Big(f(t,y(t))+f(t+h,y(t+h))\Big) \\ & = \mathcal{O}(h^3) \end{aligned}
che si ottiene dopo una lunga espansione di Taylor in più variabili e che qui è omessa. L’errore di troncamento è di ordine \mathcal{O}(h^3), più preciso rispetto a Eulero implicito.
Dimostrazione della convergenza
Analogamente al caso esplicito, possiamo dimostrare la convergenza dei metodi impliciti utilizzando l’errore locale di troncamento e le proprietà di Lipschitz delle funzioni di avanzamento.
Errore locale di troncamento
Per i metodi impliciti, come nel caso di Eulero esplicito, l’errore locale di troncamento \tau(h) è della forma:
|\tau(h)| \leq C h^{p + 1},
dove p dipende dal metodo utilizzato. Per Eulero implicito abbiamo p = 1, mentre per Crank-Nicolson p = 2.
Soluzione esatta e soluzione numerica
Consideriamo l’equazione che lega la soluzione esatta y(t_{k+1}) alla soluzione numerica y_{k+1}, includendo l’errore locale di troncamento:
y(t_{k + 1}) = \Phi(t_k, y(t_k), y(t_{k + 1}); h) + \tau_k(h).
Sottraendo l’equazione che descrive il metodo numerico:
y_{k + 1} = \Phi(t_k, y_k, y_{k + 1}; h),
otteniamo:
y(t_{k + 1}) - y_{k + 1} = \Phi(t_k, y(t_k), y(t_{k + 1}); h) - \Phi(t_k, y_k, y_{k + 1}; h) + \tau_k(h).
Uso della condizione di Lipschitz
Supponiamo che la funzione \Phi(t, y, z; h) sia Lipschitziana rispetto a y e z con costanti L_y e L_z. Applicando la condizione di Lipschitz:
|\Phi(t,y,z;h) - \Phi(t,y',z';h)| \leq L_y |y - y'| + L_z |z - z'|,
otteniamo una relazione dell’errore ricorsivo simile a quella del metodo esplicito:
|\epsilon_{k + 1}| \leq (1 + h L_z) |\epsilon_k| + C h^{p + 1}.
Conclusione
Ripetendo il ragionamento induttivo utilizzato per Eulero esplicito, possiamo dedurre che, sotto opportune ipotesi di regolarità su f(t,y) e con una condizione di Lipschitz, l’errore globale soddisfa:
E_h \leq C' h^p,
dove p dipende dal metodo scelto. Ad esempio, per Eulero implicito abbiamo p = 1, mentre per Crank-Nicolson p = 2.
Abbiamo mostrato che la dimostrazione di convergenza per i metodi impliciti può seguire un approccio simile a quello utilizzato per i metodi espliciti. La chiave sta nell’analisi dell’errore locale di troncamento e nell’uso della condizione di Lipschitz sulla funzione di avanzamento. I metodi impliciti, come Eulero implicito e Crank-Nicolson, presentano vantaggi in termini di stabilità e precisione, soprattutto per sistemi stiff, sebbene richiedano la risoluzione di equazioni non lineari a ogni passo.
28.5 Verifica Numerica dell’Ordine
Quando si sviluppa o implementa un metodo numerico, è fondamentale verificarne il corretto comportamento rispetto all’ordine teorico di convergenza. Questo controllo garantisce che il metodo funzioni correttamente e rispetti le aspettative teoriche.
Analisi Teorica
Un metodo numerico si dice di ordine p se l’errore massimo
E_h = \max_{k=0}^n |y_k-y(t_k)|
al variare del passo h segue la relazione: E_h \approx C h^p. \tag{28.8}
dove:
- y_k è la soluzione numerica al passo t_k
- y(t_k) è la soluzione esatta nello stesso punto,
- C è una costante che dipende dal problema ma non da h
- h è il passo del di avanzamento metodo.
Stima dell’Ordine Numerico
Utilizzando soluzioni calcolate con passi h_1 e h_2, e indicando i rispettivi errori come E_{h_1} e E_{h_2}, possiamo calcolare il rapporto:
\dfrac{E_{h_1}}{E_{h_2}} \approx \dfrac{C h_1^p}{C h_2^p} = \left(\dfrac{h_1}{h_2}\right)^p
Prendendo il logaritmo si ottiene una relazione lineare tra il rapporto degli errori e quello dei passi:
\log \dfrac{E_{h_1}}{E_{h_2}}\approx p \log\dfrac{h_1}{h_2},
Questa formula permette di stimare l’ordine p:
\boxed{ \quad p \approx \dfrac{\log\dfrac{E_{h_1}}{E_{h_2}}} {\log\dfrac{h_1}{h_2}},\quad } \tag{28.9}
Verifica Grafica
Riscrivendo E_h della 28.8 in forma logaritmica:
\log E_h \approx \log\left(C h^p\right)= \log C + p \log h
si nota che rappresentando graficamente \log h (asse t) e \log E_h (asse y) si ottiene una retta con:
- intercetta pari a \log C
- La pendenza stimata p rappresenta l’ordine numerico del metodo. Se il valore coincide con l’ordine teorico, il metodo è implementato correttamente.
Una buona linearità tra \log h e \log E_h conferma che il comportamento del metodo è coerente con la stima della convergenza.
Esempio numerico
Per verificare numericamente l’ordine di un metodo, è utile calcolare gli errori per una serie di valori di h e rappresentarli graficamente.
ODE con soluzione esatta
metodo di Eulero
# Implementazione del metodo di Eulero
euler_method <- function(f, t0, y0, h, t_end) {
t <- seq(t0, t_end, by = h)
n <- length(t)
y <- numeric(n)
y[1] <- y0
for (i in 1:(n - 1)) {
y[i + 1] <- y[i] + h * f(t[i], y[i])
}
data.frame(t = t, y = y)
}calcolo la soluzione per vari valori di h
# Parametri iniziali
t0 <- 0
y0 <- 0.25 # Soluzione iniziale esatta: y(0) = 1/4
t_end <- 1
h_values <- c(0.1, 0.05, 0.025, 0.01, 0.005, 0.0025 )
# Calcolo dell'errore massimo per ogni passo
errors <- numeric(length(h_values))
for (i in seq_along(h_values)) {
h <- h_values[i]
solution <- euler_method(f, t0, y0, h, t_end)
exact <- exact_solution(solution$t,y0)
errors[i] <- max(abs(solution$y - exact)) # Errore massimo
}stampo lista degli errori in funzione di h
# Tabella degli errori
error_table <- data.frame(h = h_values, Error = errors)
print(error_table) h Error
1 0.1000 0.0223156454
2 0.0500 0.0109802918
3 0.0250 0.0054467990
4 0.0100 0.0021684511
5 0.0050 0.0010825254
6 0.0025 0.0005408387Stimo errore con 28.9
# Calcolo dell'ordine di convergenza
log_h <- log(h_values)
log_E <- log(errors)
fit <- lm(log_E ~ log_h) # Regressione lineare
p_estimate <- coef(fit)[2] # Stima dell'ordine
cat("Stima dell'ordine di convergenza p:", p_estimate, "\n")Stima dell'ordine di convergenza p: 1.007675 errore in log-log plot
# Plot dell'errore
plot(
log_h, log_E,
type = "b", col = "blue", pch = 19,
xlab = "log(h)", ylab = "log(E_h)",
main = "Verifica dell'Ordine di Convergenza"
)
abline(fit, col = "red", lwd = 2)
legend(
"bottomright",
legend = c("Dati", "Fit Lineare"),
col = c("blue", "red"),
pch = c(19, NA), lty = c(NA, 1)
)
In sintesi
- Un metodo a un passo si scrive y_{k+1}=y_k+h\Psi(t_k,y_k;h); tutto il capitolo studia il legame fra il difetto locale \mathcal{O}(h^{p+1}) e l’errore globale \mathcal{O}(h^{p}).
- Il Teorema 28.1 dà la stima esplicita \dfrac{C}{L_\Psi}(e^{L_\Psi(t_k-t_0)}-1)h^{p}: la lipschitzianità di f è l’unica ipotesi sostanziale.
- Eulero esplicito ha ordine 1 ed è il caso elementare; Collatz (punto medio) e Heun salgono a ordine 2 con due valutazioni di f.
- Eulero implicito e Crank-Nicolson richiedono di risolvere un’equazione a ogni passo, ma hanno proprietà di stabilità molto migliori (Capitolo 31).
- L’ordine si verifica numericamente dal rapporto degli errori a passi dimezzati: p\approx\log(E_{h_1}/E_{h_2})/\log(h_1/h_2).