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}$$

$$\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})$$

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

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

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

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

Osservo che se $\beta _{k}=0$$\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$$

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

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

Teorema 7.2

Considero i seguenti polinomi

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

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

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

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

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

e considero 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}$$

$$\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$$h\to 0$) se valgono le seguenti condizioni:

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

   $$\sum _{j=0}^{k}\alpha _{j}=0,\quad \longrightarrow \quad \rho (1)=0$$
2. $$\sum _{j=0}^{k}j\alpha _{j}f(t_{n},y(t_{n}))=\phi _{f}(t_{n},t_{n},t_{n}$$

   $$\sum _{j=0}^{k}j\alpha _{j}f(t_{n},y(t_{n}))=\phi _{f}(t_{n},t_{n},t_{n}$$
   ($h$$h$ tende a 0 quindi tutti gli istanti tendono a $t_{n}$$t_{n}$)

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

 

Dimostrazione

Affermo che

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

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

con $\xi _{j}\in (t_{n},t_{n+j})$$\xi _{j}\in (t_{n},t_{n+j})$ con $j=0,\dots ,k$$j=0,\dots ,k$ per il teorema del valor medio.
Considero l'errore di troncamento:

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

$$\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}$$y_{n+j}$:

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

$$=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}$$

$$=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$$h\to 0$$y'(\xi _{j})\to y'(t_{n})$$y'(\xi _{j})\to y'(t_{n})$ e $\phi _{f}(y'(t),y''(t),\dots )=\phi _{f}(t_{n},t_{n},t_{n})$$\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$$

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

e se

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

$$\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})$$y'(t_{n})=f(t_{n},y_{n})$ si impongono le condizioni di consistenza:

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

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

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

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

$$\sum _{j}[j^{2}/2\alpha _{j}-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$$

$$\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 devo richiedere la 0-stabilità.

Definizione 7.4

Considero un problema di Cauchy perturbato

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

$${\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 )$$(\delta (t),\delta )$ e $(\delta ^{*}(t),\delta ^{*})$$(\delta ^{*}(t),\delta ^{*})$ tali che

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

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

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

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

$$\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}}$$

$${\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}$$\{\delta _{n}\}_{n=0}^{N}$, e $\{\delta _{n}^{*}\}_{n=0}^{N}$$\{\delta _{n}^{*}\}_{n=0}^{N}$ che corrispondono a soluzioni $z_{n}$$z_{n}$ per $n=0,\dots ,N$$n=0,\dots ,N$ e $z_{n}^{*}$$z_{n}^{*}$ per $n=0,\dots ,N$$n=0,\dots ,N$, se esiste una costante $s$$s$ positiva ed esiste $h_{0}>0$$h_{0}>0$ tale che scelto $0<H<H_{0}$$0<H<H_{0}$, segue che la norma della differenza tra le due soluzioni discrete è minore di $s\varepsilon$$s\varepsilon$, dove suppongo che $\vert \delta _{n}-\delta _{n}^{*}\vert <\varepsilon$$\vert \delta _{n}-\delta _{n}^{*}\vert <\varepsilon$. (la differenza rispetto alla definizione precedente è che deve esistere $H_{0}$$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

Considero il problema di Cauchy

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

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

Applicando il metodo multipasso:

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

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

perché $f$$f$ (e quindi il secondo membro) è identicamente nullo.
$(y_{1},\dots ,y_{k-1})$$(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}$$

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

ha radici $z_{1},\dots ,z_{s}$$z_1, \dots, z_s$ distinte rispettivamente con molteplicità $\mu _{1},\dots ,\mu _{s}$$\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}$$

$$y_n = \sum_{i=0}^s \sum_{j=0}^{\mu_{i-1}} c_{ij} z_i^n*n^j$$

(simile alla scrittura delle soluzioni nell'equazioni differenziali nel continuo)
Se $|z_{j}|>1$$|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$$1$.

 

Root condition

Definizione 7.5

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

 

Definizione 7.6 (root condition)

Il metodo numerico

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

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

soddisfa la root condition se tute le radici $z_{j}$$z_{j}$ del polinomio caratteristico sono in modulo $\leq 1$$\leq 1$ e quelle per cui vale che $|z_{p}|=1$$|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}$$

$$\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}$$

$$\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$$k$ passi che sia 0-stabile può avere ordine di convergenza maggiore di $k+1$$k+1$ se $k$$k$ è dispari, e di $k+2$$k+2$ se $k$$k$ è pari.

 

Espressione dei polinomi caratteristici nei metodi multipasso

1. Nei metodi di Adam-Bashword e Adam-Multonn il polinomio caratteristico è della forma
   $$\rho (z)=z^{k-1}*(z-1)$$

   $$\rho (z)=z^{k-1}*(z-1)$$
   e $z=1$$z=1$ è la radice principale necessaria per la consistenza, e c'è poi la radice 0 con molteplicità $k-1$$k-1$.
2. Nel Simpson generalizzato
   $$\rho (z)=z^{k}-z^{k-2}=z^{k-2}*(z^{2}-1)$$

   $$\rho (z)=z^{k}-z^{k-2}=z^{k-2}*(z^{2}-1)$$
   dove le radici sono $z_{1}=1,z_{2}=-1$$z_{1}=1,z_{2}=-1$ con molteplicità 1, e $z=0$$z=0$ con molteplicità $k-2$$k-2$.
3. Le BDF hanno il secondo polinomio caratteristico definito come
   $$\sigma (z)=\beta ^{k}*z^{k}$$

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

Metodi di Adam-Bashword

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

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

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

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

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

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

Se $\beta _{-1}=0$$\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}}$$

$${\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/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}}$$

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

$${\hbox{Adam-Multon}}$$

$${\hbox{Adam-Multon}}$$

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

$$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/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}}$$

$$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$$q$ mentre quelli di Adam-Multonn hanno ordine $p=q+1$$p=q+1$ ad eccezione dell'eulero per cui $q=1$$q=1$ e $p=1$$p=1$.

 

Tecnica del predictor-corrector

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

$$z=y_{n+1}$$

$$z=y_{n+1}$$

e considero la relazione

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

$$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 conosco tutte le quantità. $z_{k+1}=g(z_{k})$$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}$$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}$$y_n$ soluzione all'istante precedente. Si impone $g'<1$$g'<1$ quindi

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

$$|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$$

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

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

Il metodo è costoso, perché ad ogni iterazione del metodo di punto fisso si richiede una valutazione della funzione $f$$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}}$${\mathcal {P}}$): $y_{n+1}(x_{0})$$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}$$

   $$y_{n+1}(0)=y_{n}+h\sum _{j=0}^{\bar {p}}\beta _{j}f_{n-j}$$
2. evaluation (${\mathcal {E}}$${\mathcal {E}}$):
   $$f_{n+1}=f(t_{n+1},y_{n1+}(0))$$

   $$f_{n+1}=f(t_{n+1},y_{n1+}(0))$$
   (ora un'espressione approssimativa di $f_{n+1}$$f_{n+1}$ è nota)
3. Corrector (${\mathcal {C}}$$\mathcal C$):
   $$y_{n+1}=y_{n}+h\beta _{-1}f_{n+1}+h\sum _{j=0}^{p}\beta _{j}f_{n-j}$$

   $$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}}$${\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 {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)$$

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

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

Si può quindi considerare anche il metodo:

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

$${\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})$$

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

cioè il metodo si prepara per un'eventuale iterazione.

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

Esempio 7.2 (esempio di metodo predictor corrector)

Uso come predictor l'Eulero esplicito

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

$$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)))$$

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