# 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)
}Metodi numerici per ODE
Appunti di calcolo numerico
Soluzioni Numeriche per Problemi a Valore Iniziale (IVP)
In molti casi, trovare una soluzione analitica per i problemi a valore iniziale (IVP) è impraticabile a causa della complessità degli integrali o delle funzioni coinvolte. Pertanto, si ricorre a metodi numerici per approssimare le soluzioni. Questi metodi sono applicabili non solo a singoli IVP, ma anche a sistemi di equazioni differenziali ordinarie (ODE). Dato che le equazioni differenziali di ordine superiore possono essere riformulate come sistemi di equazioni del primo ordine, i metodi numerici per sistemi del primo ordine sono centrali nella risoluzione delle ODE.
Passo Temporale
Su un intervallo temporale I = [x_0, x_0 + T], suddividiamo l’intervallo in n sottointervalli, o passi temporali. Ogni sottointervallo I_k = [x_{k-1}, x_k] ha una dimensione del passo h_k = x_k - x_{k-1}, e si garantisce che x_n = x_0 + T. Se i passi temporali sono uniformi, si ha h_k = h per tutti i k. La scelta del passo h è cruciale, poiché influenza la precisione e la stabilità del metodo numerico.
Metodi Numerici per IVP: Categorie Principali
Metodi a Un Passo
Nei metodi a un passo, la soluzione al tempo x_{k+1} dipende solo dalla soluzione calcolata al tempo x_k. In altre parole, l’approssimazione y_{k+1} si basa unicamente sull’informazione relativa a y_k. Il più semplice tra questi è il metodo di Eulero.
Metodi a Più Passi
I metodi a più passi utilizzano informazioni provenienti da diversi passi precedenti per calcolare la nuova approssimazione y_{k+1}. Questi metodi sono più complessi e solitamente offrono una maggiore precisione rispetto ai metodi a un passo, poiché tengono conto di una parte maggiore della storia della soluzione. Un esempio noto è il metodo di Adams-Bashforth.
Metodi Espliciti
Un metodo esplicito (o integrazione in avanti) calcola la soluzione al passo successivo usando solo i valori già noti. A partire dalla condizione iniziale y_0 al tempo x_0, la sequenza delle soluzioni viene determinata tramite una relazione del tipo:
y_{k+1} = y_k + h_k F(x_k, y_k),
dove F è una funzione che rappresenta l’incremento della soluzione al passo k. I metodi espliciti sono computazionalmente meno costosi e semplici da implementare, ma possono richiedere passi temporali piccoli per mantenere la stabilità della soluzione.
Metodi Impliciti
Al contrario, nei metodi impliciti (o integrazione all’indietro), la funzione F dipende dal valore della soluzione al tempo x_{k+1}, il che porta alla necessità di risolvere un’equazione per y_{k+1} a ogni passo:
y_{k+1} = y_k + h_k F(x_{k+1}, y_{k+1}).
Questi metodi sono computazionalmente più intensivi, poiché richiedono la risoluzione iterativa di sistemi di equazioni, ma offrono una maggiore stabilità, soprattutto per problemi rigidi. Un esempio è il metodo di Eulero implicito.
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 Taylor. Consideriamo un problema di valore iniziale:
\begin{cases} y'(x) = f(x,y), \\ y(x_0) = y_0. \end{cases}
L’espansione di Taylor di y(x) attorno al punto iniziale x_0 è:
y(x_0 + h) = y_0 + h y'(x_0) + O(h^2),
dove
y'(x_0) = f(x_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(x_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 x_{k+1} anziché in x_k. L’equazione risultante è:
y_{k+1} = y_k + h f(x_{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.
Convergenza del Metodo di Eulero Esplicito
Consideriamo la soluzione esatta y(x_k) al punto x_k e introduciamo l’errore
\epsilon_k = y(x_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.
Sviluppo di Taylor e Relazione con il Metodo di Eulero
L’espansione di Taylor della soluzione esatta y(x) attorno a x_k è:
y(x_{k+1}) = y(x_k + h) = y(x_k) + y'(x_k) h + \frac{y''(\zeta_k)}{2} h^2,
dove \zeta_k \in [x_k, x_{k+1}] è un punto intermedio.
Il metodo di Eulero esplicito approssima la soluzione al passo successivo come:
y_{k+1} = y_k + h f(x_k, y_k).
Sottraendo l’espressione di y_{k+1} dall’espansione di Taylor di y(x_{k+1}), otteniamo:
\begin{aligned} y(x_{k+1}) - y_{k+1} &= y(x_k) - y_k \\[1em] &+h \big(f(x_k, y(x_k)) - f(x_k, y_k)\big)\\[1em] &+ \frac{y''(\zeta_k)}{2} h^2. \end{aligned}
Questo può essere riscritto come:
\epsilon_{k+1} = \epsilon_k + h \big(f(x_k, y(x_k)) - f(x_k, y_k)\big) + \frac{y''(\zeta_k)}{2} h^2.
Stima dell’Errore
Assumendo che f(x,y) sia lipschitziana rispetto a y, ovvero:
|f(x_k, y(x_k)) - f(x_k, y_k)| \leq L |y(x_k) - y_k| = L |\epsilon_k|,
e che esista una costante M_2 tale che
|y''(x)| \leq M_2 \quad \text{per ogni } x \in [a,b],
possiamo stimare l’errore al passo successivo come:
|\epsilon_{k+1}| \leq (1 + hL) |\epsilon_k| + \frac{M_2}{2} h^2. \tag{1}
usando questa diseguaglianza possiamo dimostrale la convegenza del metodo di Eulero.
Lemma 1 Sia
E_h = \max_{k=0}^{n} |\epsilon_k|
l’errore massimo e
h = \frac{b - a}{n}
il passo temporale. Per il metodo di Eulero esplicito si ha la seguente stima:
E_h \leq C h, \qquad C = \frac{2L e^{(b-a)}}{L M_2}
dove L è la costante di Lipschitz di f(x,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 \frac{h^2}{2} nella stima dell’errore data in Equazione 1, otteniamo che l’errore si propaga con la seguente limitazione:
|\varepsilon_{k+1}| \leq A |\varepsilon_k| + B,
dove \varepsilon_0 = y(x_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 = \frac{A^{k+1} - 1}{A - 1}.
Di conseguenza, abbiamo:
E_h = \max |\varepsilon_k| \leq \frac{A^n - 1}{A - 1} B \leq \frac{A^n}{A - 1} B.
Sostituendo A = 1 + hL e B = M_2 \frac{h^2}{2}, otteniamo:
E_h \leq \frac{A^n}{A - 1} B = \frac{(1+hL)^n}{hL} M_2 \frac{h^2}{2}.
Utilizzando la disuguaglianza valida per z \geq 0:
1+z \leq e^z, \tag{2}
risultante dalla serie di Taylor, otteniamo:
E_h \leq \frac{(e^{hL})^n}{2L} M_2 h = \frac{e^{nhL}}{2L} M_2 h = \left(\frac{e^{(b-a)L}}{2L} M_2\right) h.
Convergenza del Metodo di Eulero Implicito
Analogamente al metodo esplicito, consideriamo la soluzione esatta y(x_k) al punto x_k e l’errore \epsilon_k = y(x_k) - y_k.
Sviluppo di Taylor e Relazione con il Metodo di Eulero Implicito
L’espansione di Taylor della soluzione esatta y(x) attorno a x_{k+1} è:
y(x_{k+1}) = y(x_k) + y'(x_{k+1}) h - \frac{y''(\zeta_k)}{2} h^2,
dove \zeta_k \in [x_k, x_{k+1}].
Nel metodo di Eulero implicito, l’approssimazione al passo successivo si basa su:
y_{k+1} = y_k + h f(x_{k+1}, y_{k+1}),
dove la funzione f(x, y) è valutata al passo successivo (x_{k+1}, y_{k+1}).
Sottraendo l’espressione di y_{k+1} dalla serie di Taylor, otteniamo:
\begin{aligned} y(x_{k+1}) - y_{k+1} & = y(x_k) - y_k \\[1em] & + h \big( f(x_{k+1}, y(x_{k+1})) - f(x_{k+1}, y_{k+1}) \big) \\[1em] & - \frac{y''(\zeta_k)}{2} h^2, \end{aligned}
che può essere riscritta come:
\epsilon_{k+1} = \epsilon_k + h \big( f(x_{k+1}, y(x_{k+1})) - f(x_{k+1}, y_{k+1}) \big) - \frac{y''(\zeta_k)}{2} h^2.
Stima dell’Errore
Assumiamo che f(x, y) sia lipschitziana rispetto a y, cioè:
|f(x_{k+1}, y(x_{k+1})) - f(x_{k+1}, y_{k+1})| \leq L |y(x_{k+1}) - y_{k+1}| = L |\epsilon_{k+1}|,
e che esista una costante M_2 tale che |y''(x)| \leq M_2 per ogni x \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}| + \frac{M_2}{2} h^2.
Riorganizzando i termini, otteniamo:
|\epsilon_{k+1}| \leq \frac{|\epsilon_k| + \frac{M_2}{2} h^2}{1 - hL}.
Da questa disuguaglianza possiamo dimostrare la convergenza del metodo di Eulero implicito.
Lemma 2 (Convergenza del Metodo di Eulero Implicito) Sia
E_h = \max_{k=0}^{n} |\epsilon_k|
l’errore massimo e sia
h = \frac{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(x, 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 \frac{|\epsilon_k| + \frac{M_2}{2} h^2}{1 - hL}.
Poniamo
A = \frac{1}{1 - hL} \quad\textrm{e}\quad B = A\frac{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 \frac{A^{n+1} - 1}{A - 1} B.
Assumendo h sufficientemente piccolo affinché Lh\leq 1/2, possiamo scrivere:
\begin{aligned} A &= \frac{1}{1 - hL} = \sum_{k=0}^\infty (hL)^k \\ &= 1 + hL \sum_{k=1}^\infty (hL)^k \leq 1 + 2hL, \\ A-1 &= \frac{1}{1 - hL} - 1 \geq hL, \end{aligned}
e quindi, usando anche Equazione 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\frac{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'(x) = f(x,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( x_k, y_k; h )
dove \Phi( x, y; h ) è la funzione che rappresenta il passo di avanzamento. Ad esempio, nel caso di Eulero esplicito:
\Phi( x, y; h ) = y + h f(x,y)
Funzione di Avanzamento per Altri Schemi Numerici
Vediamo ora come scrivere la funzione di avanzamento \Phi( x, 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 + \frac{h}{2}f(x_k,y_k) \\ y_{k+1} & = y_k + h f\left(x_k+\frac{h}{2},\; y_{k+1/2}\right) \end{aligned}
Questo schema può essere riarrangiato in una forma compatibile con \Phi( x, y; h ):
y_{k+1} = y_k + h f\left(x_k+\frac{h}{2},\; y_k + \frac{h}{2}f(x_k,y_k)\right)
per cui la funzione di avanzamento per Collatz diventa
\Phi( x, y; h ) = y + h f\left(x+\frac{h}{2},\; y + \frac{h}{2}f(x,y)\right)
Metodo di Heun
Il metodo di Heun è definito da:
\begin{aligned} \widetilde{y}_{k+1} & = y_k + hf(x_k,y_k) \\ y_{k+1} & = y_k + \frac{h}{2}\left( f(x_k,y_k) + f(x_{k+1},\widetilde{y}_{k+1}) \right) \end{aligned}
Riarrangiando anche questo schema, otteniamo:
y_{k+1} = y_k + \frac{h}{2}\Big( f(x_k,y_k) + f\big(x_k+h,y_k + hf(x_k,y_k)\big) \Big)
per cui la funzione di avanzamento per Heun diventa
\Phi( x, y; h ) = y + \frac{h}{2}\Big(f(x,y) + f\big(x+h,y+ hf(x,y)\big)\Big)
Propagazione dell’errore
Sostituendo la soluzione esatta nella funzione di avanzamento, otteniamo:
y(x+h) = \Phi( x, y(x) ) + \color{blue} \tau(x,h)
dove \color{blue} \tau(x,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(x,h) &= y(x+h)-\Phi( x, y(x), h ) \\ & = {\color{red}y(x+h)}-{\color{blue}y(x)-h f(x,y(x))} \\ & = {\color{red}\cancel{y(x)}+\cancel{y'(x)h}+y''(\zeta)h^2/2} -{\color{blue}\cancel{y(x)}-\cancel{h y'(x)}} \\ &= \dfrac{h^2}{2}y''(\zeta) \end{aligned}
dove \zeta \in [x, x+h].
Errore di Troncamento per Collatz
Per il metodo di Collatz:
\begin{aligned} \tau(x,h) &= y(x+h)-\Phi( x, y(x), h ) \\ & = {\color{red}y(x+h)}-{\color{blue}y(x) - h f\left(x+\frac{h}{2},\; y(x) + \frac{h}{2}f(x,y(x))\right)} \\ & = {\color{red}\cancel{y(x)}+y'(x)h+y''(x)h^2/2+y'''(\zeta)h^3/6} \\ &-{\color{blue}\cancel{y(x)}-h\left(f(x,y(x))+ \frac{h}{2}\frac{\partial f(x,y(x))}{\partial x}+ \frac{h}{2}\frac{\partial f(x,y(x))}{\partial x}f(x,y(x))+ \mathcal{O}(h^2) \right)} \\ & = {\color{red}\cancel{f(x,y(x))h}+\frac{h^2}{2}\dfrac{\mathrm{d}}{\mathrm{d}x}f(x,y(x))+y'''(\zeta)h^3/6} \\ &{\color{blue}-\cancel{hf(x,y(x))}- \frac{h^2}{2}\frac{\partial f(x,y(x))}{\partial x}- \frac{h^2}{2}\frac{\partial f(x,y(x))}{\partial x}f(x,y(x))+ \mathcal{O}(h^3)} \\ & = {\color{red}\frac{h^2}{2} \cancel{\left( \frac{\partial f(x,y(x))}{\partial x}+ \frac{\partial f(x,y(x))}{\partial x}f(x,y(x)) \right)} +y'''(\zeta)h^3/6} \\ &{\color{blue} -\frac{h^2}{2} \cancel{\left(\frac{\partial f(x,y(x))}{\partial x} +\frac{\partial f(x,y(x))}{\partial x}f(x,y(x)) \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(x, y) è sufficientemente regolare, l’errore locale di troncamento soddisferà:
\color{blue} |\tau(x,h)| \;\leq\; C\, h^{p+1} \tag{3}
ove p \geq 1 dipende dal metodo utilizzato.
Convergenza di un Metodo a un Passo
Consideriamo un metodo generico a un passo: y_{k+1} = \Phi(x_k,y_k;h)
Sfruttando l’errore locale di troncamento:
y(x_{k+1}) = \Phi(x_k,y(x_k);h)+\tau(x_k,h)
possiamo scrivere l’errore come:
y(x_{k+1})-y_{k+1}= \Phi(x_k,y(x_k);h) -\Phi(x_k,y_k;h) +\tau(x_k,h) \tag{4}
Se la funzione \Phi(x, y; h) è Lipschitziana rispetto a y, con costante vicino a 1 del tipo 1+hL come segue: |\Phi(x,y;h)-\Phi(x,z;h)| \leq (1+hL) |y-z|
Utilizzando questa proprietà nell’equazione dell’errore Equazione 4 con l’ipotesi sull’errore locale di troncamento Equazione 3, otteniamo:
|\epsilon_{k+1}|\leq (1+hL)|\epsilon_k| + C h^{p+1}
Questa disuguaglianza è analoga a quella trovata per il metodo di Eulero esplicito Equazione 1, e replicando il ragionamento nel Lemma 1 sulla convergenza di Eulero esplicito, otteniamo che l’errore massimo per un metodo a un passo soddisfa:
E_h \leq \widetilde{C}\,h^p
dove \widetilde{C} è una costante che dipende dalle proprietà della funzione f(x, y) e dalla regolarità della soluzione.
Generalizzazione (metodi impliciti)
Un metodo implicito a un passo può essere scritto genericamente come:
y_0 = a, \qquad y_{k+1} = \Phi( x_k, y_k, y_{k+1}; h )
dove \Phi(x, 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( x, y, z; h ) = z + h f(x,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(x+h) e la previsione del metodo numerico. L’errore di troncamento per un metodo implicito è dato da:
\tau(x,h) = y(x+h)-\Phi( x, y(x), y(x+h); h )
Errore di Troncamento per Eulero Implicito
Per il metodo di Eulero implicito, l’errore locale di troncamento si calcola come segue: \begin{aligned} \tau(x,h) &= y(x+h)-\Phi( x, y(x), y(x+h), h ) \\ & = {\color{red}y(x+h)}-{\color{blue}y(x)-h f(x+h,y(x+h))} \\ & = {\color{red}y(x+h)}-{\color{blue}y(x)-hy'(x+h)} \\ & = {\color{red}\cancel{y(x+h)}}-{\color{blue}(\cancel{y(x+h)}-h \cancel{y'(x+h)}+h^2/2 y''(\zeta))-h\cancel{y'(x+h)}} \\ &=-\dfrac{h^2}{2}y''(\zeta) \end{aligned}
dove \zeta \in [x, x+h]. Notiamo che l’errore locale di troncamento per Eulero implicito è:
\tau(x,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-Nicolson, basato sulla regola dei trapezi, è un metodo implicito che può essere scritto come: y_{k+1} = y_k + \dfrac{h}{2}\Big( f(x_k,y_k)+ f(x_{k+1},y_{k+1}) \Big)
per cui la funzione di avanzamento prende la forma
\Phi( x, y, z; h ) = y + \dfrac{h}{2}\Big(f(x,y)+f(x+h,z)\Big)
L’errore locale di troncamento per questo metodo si calcola come: \begin{aligned} \tau(x,h) &= y(x+h)-\Phi( x, y(x), y(x+h), h ) \\ &= y(x+h) - y(x) - \dfrac{h}{2}\Big(f(x,y(x))+f(x+h,y(x+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: Schema Simile al Metodo Esplicito
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.
Confronto tra la soluzione esatta e quella numerica
Consideriamo l’equazione che lega la soluzione esatta y(x_{k+1}) alla soluzione numerica y_{k+1}, includendo l’errore locale di troncamento:
y(x_{k + 1}) = \Phi(x_k, y(x_k), y(x_{k + 1}); h) + \tau_k(h).
Sottraendo l’equazione che descrive il metodo numerico:
y_{k + 1} = \Phi(x_k, y_k, y_{k + 1}; h),
otteniamo:
y(x_{k + 1}) - y_{k + 1} = \Phi(x_k, y(x_k), y(x_{k + 1}); h) - \Phi(x_k, y_k, y_{k + 1}; h) + \tau_k(h).
Uso della condizione di Lipschitz
Supponiamo che la funzione \Phi(x, y, z; h) sia Lipschitziana rispetto a y e z con costanti L_y e L_z. Applicando la condizione di Lipschitz:
|\Phi(x,y,z;h) - \Phi(x,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(x,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.
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(x_k)|
al variare del passo h segue la relazione: E_h \approx C h^p. \tag{5}
dove:
- y_k è la soluzione numerica al passo x_k
- y(x_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{6}
Verifica Grafica
Riscrivendo E_h della Equazione 5 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 x) 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 formula Equazione 6
# 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)
)