Calcolo numerico/Metodi multipasso

Aumento dell'ordine di convergenza

Per aumentare l'ordine di convergenza del metodo ci sono due strade possibili:

1. metodi multipasso: usano informazioni precedenti ma rimangono lineari, ne sono esempi i metodi di Adam-Bashword (espliciti), di Adam-Multonn (impliciti), di Mistrom, i generalized Simpson (espliciti e impliciti), BDF usate nel caso dei problemi stiff;
2. metodi a un passo non lineari.

Metodi multipasso lineare

I metodi di questo tipo si definiscono come

$\sum _{j=0}^{k}\alpha _{j}y_{n+j}=h\sum _{j=0}^{k}b_{j}f_{n+j}$

dove

$f_{n+j}=f(t_{n+j},y_{n+j})$

$\alpha _{j},\beta _{j}$ sono tali che $\alpha _{k}=1$ e

$\sum |\alpha _{0}|+|\beta _{0}|\neq 0$

Osserviamo che se $\beta _{k}=0$ il metodo è esplicito perché a secondo membro non compare la soluzione da calcolare al passo corrente, mentre diversamente è implicito, e si può pensare di applicare il metodo del punto fisso

$y_{n+k}=h\beta _{k}f(t_{n+k},y_{n+k})+\phi$

Dev'essere $h<1/(|\beta _{k}|*L)$.

Teorema 7.2

Consideriamo i seguenti polinomi

$\rho (z)=\sum _{j=0}^{k}\alpha _{j}z^{j}$

con $z$ variabile complessa (primo polinomio caratteristico).

$\sigma (z)=\sum _{j=0}^{k}\beta _{j}z^{j}$

e consideriamo il metodo

$\sum _{j=0}^{k}\alpha _{j}y_{n+j}=h\phi _{f}(y_{n+k},y_{n+k-1},y_{n},t_{n},h)=h\sum _{j=0}^{k}\beta _{j}f_{n+j}$

Il metodo è consistente (cioè l'errore di troncamento locale tende a 0 per $h\to 0$), se valgono le seguenti condizioni:

$\sum _{j=0}^{k}\alpha _{j}=0,\quad \longrightarrow \quad \rho (1)=0$

$\sum _{j=0}^{k}j\alpha _{j}f(t_{n},y(t_{n}))=\phi _{f}(t_{n},t_{n},t_{n}$

($h$ tende a 0 quindi tutti gli istanti tendono $t_{n}$)

con $\phi _{f}=\sum _{j}\beta _{j}f_{n+j}$.

Dimostrazione

Affermiamo che

$y(t_{n+j})-y(t_{n})=jhy'(\xi _{j})=jhf(\xi _{j})$

con $\xi _{j}\in (t_{n},t_{n+j})$ con $j=0,\dots ,k$ per il teorema del valor medio.

Consideriamo l'errore di troncamento:

$\tau =1/h[\sum _{j=0}^{k}\alpha _{j}y(t_{n+j})-h\phi _{f}]$

e in base a quanto appena ricavato per $y_{n+j}$:

$=1/h[\sum _{j=0}^{k}\alpha _{j}*(y(t_{n})+jhy'(\xi _{j}))]-\phi _{f}$

e spezzando le somme

$=1/h*\sum _{j=0}^{k}\alpha _{j}*y(t_{n})+\sum _{j=0}^{k}j\alpha _{j}y'(\xi _{j})-\phi _{f}$

Per $h\to 0$ $y'(\xi _{j})\to y'(t_{n})$ e $\phi _{f}(y'(t),y''(t),\dots )=\phi _{f}(t_{n},t_{n},t_{n})$

Segue quindi che si ha consistenza se

$\sum _{j=1}^{k}\alpha _{j}=0$

e se

$\sum _{j=0}^{k}j\alpha _{j}y'(t_{n})-\phi _{f}=0$

Nella pratica tenendo conto che $y'(t_{n})=f(t_{n},y_{n})$ si impongono le condizioni di consistenza:

$\sum _{j}\alpha _{j}=0$

$\sum _{j}[j\alpha _{j}-\beta _{j}]=0$

$\sum _{j}[j^{2}/2\alpha _{j}-j\beta _{j}]=0$

$\sum _{j}[j^{3}/6\alpha _{j}-j^{2}/2\beta _{j}]=0\ldots$

Relazione tra convergenza e consistenza

Se ho convergenza, allora ho consistenza, ma non vale viceversa. Per avere convergenza, oltre alla consistenza dobbiamo richiedere la 0-stabilità.

Definizione

Consideriamo un problema di Cauchy perturbato

${\begin{cases}z'=f(t,z+\delta (t))\\z(t_{0})=z_{0}+\delta \end{cases}}$

Si dice che il problema è ben posto o totalmente stabile se per qualunque coppia di perturbazioni $(\delta (t),\delta )$ e $(\delta ^{*}(t),\delta ^{*})$ tali che

$\vert \delta (t)-\delta ^{*}(t)\vert <\vert \delta -\delta ^{*}\vert$

segue che, posto $z(t),z^{*}(t)$ le soluzioni dei corrispondenti problemi di Cauchy, esiste una costante $s$ positiva tale che per ogni $t\in I_{0}$,

$\vert z(t)-z^{*}(t)\vert \leq s\varepsilon$

Definizione corrispondente nel discreto: consideriamo un sistema di differenze perturbato.

${\begin{cases}z_{n+1}=z_{n}+h\phi _{f}(t_{n},y_{n},h)+\delta _{n+1})\\z_{0}=y_{0}+\delta _{0}\end{cases}}$

Il metodo è 0-stabile se, date le successioni delle perturbazioni $\{\delta _{n}\}_{n=0}^{N}$, e $\{\delta _{n}^{*}\}_{n=0}^{N}$ che corrispondono a soluzioni $z_{n}$ per $n=0,\dots ,N$ e $z_{n}^{*}$ per $n=0,\dots ,N$, se esiste una costante $s$ positiva ed esiste $h_{0}>0$ tale che scelto $0<H<H_{0}$, segue che la norma della differenza tra le due soluzioni discrete è minore di $s\varepsilon$, dove supponiamo che $\vert \delta _{n}-\delta _{n}^{*}\vert <\varepsilon$ (la differenza rispetto alla definizione precedente è che deve esistere $H_{0}$ tale per cui si abbia la conservazione dell'entità della perturbazione).

Stabilità come condizione necessaria

Questo esempio mostra che è necessario richiedere la stabilità.

Esempio 7.1

Consideriamo il problema di Cauchy

${\begin{cases}y'=0\\y(t_{0})=0\end{cases}}$

Applicando il metodo multipasso:

$\sum _{j=0}^{k}\alpha _{j}y_{n+j}=0$

perché $f$ (e quindi il secondo membro) è identicamente nullo.

$(y_{1},\dots ,y_{k-1})$ sono molto piccoli ma non esattamente nulli.

Questa è un'equazione alle differenze omogenea. Il polinomio

$\rho (z)=\sum _{j=0}^{k}\alpha _{j}z^{j}$

ha radici $z_{1},\dots ,z_{s}$ distinte rispettivamente con molteplicità $\mu _{1},\dots ,\mu _{s}$, allora la soluzione dell'equazione alle differenze è

$y_{n}=\sum _{i=0}^{s}\sum _{j=0}^{\mu _{i-1}}c_{ij}z_{i}^{n}*n^{j}$

(simile alla scrittura delle soluzioni nelle equazioni differenziali nel continuo).

Se $|z_{j}|>1$, il metodo diverge, e lo stesso avviene nel caso in cui il modulo di $z$ è uguale a 1 e ha molteplicità maggiore di 1.

Root condition

Definizione

Consideriamo $\rho (z)$ e chiamiamo $z_{i}$ le sue radici, per $i=1,\dots ,k$. $z_{1}=1$ è detta radice principale, e le altre radici $z_{i}$ per $i=2,\dots ,k$ sono le radici spurie.

Definizione

Il metodo numerico

$\sum _{j}\alpha _{j}y_{n+j}=h\sum _{j}\beta _{j}f_{n+j}$

soddisfa la root condition se tutte le radici $z_{j}$ del polinomio caratteristico sono in modulo $\leq 1$ e quelle per cui vale che $|z_{p}|=1$ hanno molteplicità algebrica 1.

Teorema 7.3

Il metodo numerico

$\sum _{j=0}^{s}\alpha _{j}y_{n+j}=h\sum _{j}\beta _{j}f_{n+j}$

è 0-stabile se e solo se vale la root condition.

Teorema 7.4

Il metodo

$\sum _{j=0}^{s}\alpha _{j}y_{n+j}=h\sum _{j}\beta _{j}f_{n+j}$

è convergente se e solo se è consistente e 0-stabile.

Teorema 7.5 (prima barriera di Dahlquist)

Nessun metodo lineare multipasso a $k$ passi che sia 0-stabile può avere ordine di convergenza maggiore di $k+1$ se $k$ è dispari, e di $k+2$ se $k$ è pari.

Espressione dei polinomi caratteristici nei metodi multipasso

Nei metodi di Adam-Bashword e Adam-Multonn il polinomio caratteristico è della forma

$\rho (z)=z^{k-1}*(z-1)$

e $z=1$ è la radice principale necessaria per la consistenza, e c'è poi la radice 0 con molteplicità $k-1$.

Nel Simpson generalizzato

$\rho (z)=z^{k}-z^{k-2}=z^{k-2}*(z^{2}-1)$

dove le radici sono $z_{1}=1,z_{2}=-1$ con molteplicità 1, e $z=0$ con molteplicità $k-2$.

Le BDF hanno il secondo polinomio caratteristico definito come

$\sigma (z)=\beta ^{k}*z^{k}$

Metodi di Adam-Bashword

I metodi di Adam-Bashword derivano da formule di quadratura. Consideriamo l'equazione integrale

$y(t_{n+1})-y(t_{n})=\int _{t_{n}}^{t_{n+1}}f(t,y(t))\,dt$

Consideriamo il polinomio interpolante $f$ con nodi equispaziati: in particolare, per i metodi espliciti consideriamo il polinomio interpolante di grado $p$ nei nodi $(t_{n},\dots ,t_{n-p})$, per gli Adam-Multonn consideriamo il polinomio di grado $p+1$ nei nodi $t_{n+1},\dots t_{n-p}$ e otteniamo la formula

$y_{n+1}=y_{n}+h\sum _{j=-1}^{p}\beta _{j}f_{n-j}$

Se $\beta _{-1}=0$ il metodo è esplicito, altrimenti è implicito.

${\hbox{Metodi di Adam-Bashword:}}\,y_{n+1}=y_{n}+hf_{n}\quad {\hbox{1 nodo, Eulero esplicito}}$

$y_{n+1}=y_{n}+h/2(3f_{n}-f_{n-1})\quad {\hbox{2 nodi, trapezi esplicito}}$

$y_{n+1}=y_{n}+h/12(23f_{n}-16f_{n-1}+5f_{n-2})\quad {\hbox{3 nodi}}$

Adam-Multon

$y_{n+1}=y_{n}+hf_{n+1}\quad {\hbox{1 nodo, Eulero implicito}}$

$y_{n+1}=y_{n}+h/2(f_{n+1}+f_{n})\quad {\hbox{2 nodi, trapezi implicito}}$

$y_{n+1}=y_{n}+h/12(5f_{n+1}+8f_{n}-f_{n-1})\quad {\hbox{3 nodi}}$

In conclusione vale il seguente teorema.

Teorema 7.6

I metodi espliciti di Adam-Bashword hanno ordine $q$ mentre quelli di Adam-Multonn hanno ordine $p=q+1<math>aeccezionedell'Euleropercui<math>q=1$ e $p=1$./

Tecnica del predictor-corrector

Consideriamo un metodo multipasso implicito, da risolvere con un metodo di punto fisso. Poniamo

$z=y_{n+1}$

e consideriamo la relazione

$z=y_{n}+h\beta _{-1}f(t_{n+1},z)+h\sum _{j=0}^{p}\beta _{j}f_{n-j}$

dove il secondo addendo è noto e ne conosciamo tutte le quantità. $z_{k+1}=g(z_{k})$ con $g=y_{n}+h\beta _{-1}f(t_{n+1},z)+h\sum _{j=0}^{p}\beta _{j}f_{n-j}$, e si dà come innesco $y_{n}$ soluzione all'istante precedente. Si impone $g'<1$ quindi

$|g'(z)|=h\beta _{-1}*|{\frac {\partial f}{\partial z}}|\leq h\beta _{-1}L$

e supponendo di avere le condizioni del teorema di esistenza e unicità

$h\beta _{-1}L<1$

cioè sul passo $h$ otteniamo la condizione $h<|\beta _{-1}|*L$ e la costante di Lipschitz condiziona la convergenza del metodo di punto fisso.

Il metodo è costoso, perché a ogni iterazione del metodo di punto fisso si richiede una valutazione della funzione $f$. Più l'innesco fornito è buono, più velocemente viene soddisfatto il criterio d'arresto, e meno oneroso è il metodo dal punto di vista computazionale.

Metodo predictor-corrector

Il metodo predictor-corrector consiste dei seguenti passi:

1. predictor (${\mathcal {P}}$): $y_{n+1}(x_{0})$ è ottenuto applicando uno schema multipasso esplicito, cioè

$y_{n+1}(0)=y_{n}+h\sum _{j=0}^{\bar {p}}\beta _{j}f_{n-j}$

2. evaluation (${\mathcal {E}}$):

$f_{n+1}=f(t_{n+1},y_{n1+}(0))$

(ora un'espressione approssimativa di $f_{n+1}$ è nota)

3. corrector (${\mathcal {C}}$):

$y_{n+1}=y_{n}+h\beta _{-1}f_{n+1}+h\sum _{j=0}^{p}\beta _{j}f_{n-j}$

I passi ${\mathcal {E}},{\mathcal {C}}$ possono essere iterati nel seguente modo:

${\mathcal {E}}\,f_{n+1}(k)=f(t_{n+1},y_{n+1}(k))$

${\mathcal {C}}\,y_{n+1}(k)=y_{n}+h\beta _{-1}f_{n+1}(k)$

dove $f_{n+1}(k)$ è stato calcolato prima.

Si può quindi considerare anche il metodo:

${\mathcal {P}}({\mathcal {E}}{\mathcal {C}})m$

in cui valutazione e correzione vengono alternate $m$ volte.

Il passo finale è

$f_{n+1}^{m}=f(t_{n+1},y_{n+1}^{m})$

cioè il metodo si prepara per un'eventuale iterazione. Se ${\mathcal {P}}$ ha ordine $p-1$ i dati di innesco aggiuntivi sono calcolabili con un errore di ordine di $o(h^{p})$ e se il corrector è di ordine $p$, allora $({\mathcal {P}}{\mathcal {C}})$ ha ordine $p$.

Esempio 7.2 (esempio di metodo predictor corrector)

Usiamo come predictor l'Eulero esplicito

$y_{n+1}(0)=y_{n}+hf_{n}$

e come corrector il metodo dei trapezi:

$y_{n+1}=y_{n}+h/2(f_{n}+f(t_{n+1},y_{n+1}(0)))$

Questo metodo può anche essere considerato come metodo a un passo.