Calcolo numerico/Algoritmi per la risoluzione

Metodi a un passo

Dato il problema

$\{\begin{array}{c}{y}^{\prime }(t)=f(t,y)\\ y(t_0)=y_0\end{array}$

supponiamo di volerlo risolvere in $I=(t_0,T)\subset I_0$.

Creiamo una successione di valori $\{y_n\}$ dove $y_n$ approssima $y(t_n)$ e $t_n=t_0+nh,n=(T-t_0)/h$, e $h$ viene chiamato passo di discretizzazione.

Per calcolare $y(t+h)$ conoscendo $y(t)$ consideriamo il rapporto incrementale

$\frac{y(t+h)-y(t)}{h}$

che deve essere approssimato. Per l'approssimazione si utilizza la derivata prima, che in questo caso è un dato del problema. Quindi si pone

$\frac{y(t+h)-y(t)}{h}=f(t,y(t))=y^{\prime }(t)$

ed esplicitando $y(t+h)$:

$y(t+h)=y(t)+h\ast f(t,y(t)).$

Metodo di Eulero esplicito

Otteniamo quindi il seguente metodo detto metodo di Eulero esplicito:

$y_{n+1}=y_n+hf(t_n,y_n)$

con $t_n=t_0+nh$, e $y_n$ approssimazione della soluzione esatta all'istante $t_n$.

$h=t_{n+1}-t_n$

Questo metodo è detto esplicito perché la quantità da calcolare, $y_{n+1}$, non compare nel membro di destra.

Sappiamo che

$y(t_{n+1})=y(t_n)+{\int }_{t_n}^{t_{n+1}}f(\tau ,y(\tau ))d\tau$

e sostituendo alla funzione il polinomio interpolante di grado 0 in $(t_n,f(t_n,y_n))$, cioè la retta $p=f(t_n,y(t_n))$ otteniamo il metodo di Eulero esplicito:

$y(t_{n+1})=y(t_n)+{\int }_{t_n}^{t_{n+1}}f(t_n,y(t_n))d\tau$

$y(t_{n+1})=y(t_n)+f(t_n,y(t_n))(t_{n+1}-t_n)$

$y(t_{n+1})=y(t_n)+hf(t_n,y(t_n))$

Generalizzazione del metodo

Consideriamo la generalizzazione del metodo

$y_{n+1}=y_n+h\mathrm{\Phi }(t_n,y_n,f(t_n,y_n);h)$

e lo confronto con la corrispondente espressione nel continuo (soluzione esatta):

$y(t_{n+1})=y(t_n)+h\ast \mathrm{\Phi }(t_n,y(t_n),f(t_n,y(t_n)));h)$

All'istante $t_{n+1}$ definisco il residuo ${\epsilon }_{n+1}$. Ci si chiede di quanto la soluzione esatta non soddisfa lo schema numerico considerato.

${\epsilon }_{n+1}=h{\tau }_{n+1}$ è funzione del passo di discretizzazione, e si può scrivere

${\tau }_{n+1}(h)={\epsilon }_{n+1}/h$

dove ${\tau }_{n+1}$ è detto errore di troncamento locale. ${\epsilon }_{n+1}$ è la differenza tra la soluzione $y(t_{n+1})$ e la soluzione che si ottiene considerando come $y_n$ il valore esatto $y(t_n)$.

Più precisamente

Per dire che l'errore commesso tende a 0, non basta richiedere che ${\epsilon }_{n+1}\to 0$, altrimenti si elimina la dipendenza del metodo da $\mathrm{\Phi }$.

Consistenza

La condizione di consistenza è la minima condizione che dev'essere soddisfatta affinché il metodo funzioni.

Considerando il metodo di Eulero esplicito:

$\mathrm{\Phi }(t,y,f,h)=f(t,y(t))$

$\tau (t,y(t),h)=1/h\ast [y(t+h)-y(t)-h\mathrm{\Phi }(t,y,h)]$

$=1/h\ast [y(t+h)-y(t)-hf(t,y(t))]$

e sviluppando con Taylor $y(t+h)$ e tenendo conto che $f(t,y(t)=y^{\prime }(t)$:

$=1/h\ast [y(t)+h{y}^{\prime }(t)+h^2/2y^{(2)}(\xi )-y(t)-h{y}^{\prime }(t)]$

e in seguito a cancellazioni:

$=1/h\ast h^2/2y^{(2)}(\xi )=o(h)$

e il metodo di Eulero esplicito è consistente di ordine $p=1$.

Si ricava anche

$\tau (h,y,f)=h/2y^{(2)}(\xi )\text{formula 1}$

Convergenza nell'intervallo $[t_{0}T]$

Per $h\to 0$, l'insieme di valori $t_n=t_0+nh$ tende a ricoprire l'intervallo in questione.

Chiediamo che $y_n$ converga alla soluzione esatta $y(t)$ per ogni $t\in I$.

Il problema è convergente di ordine $p$ se

$C=o(h^p)$

Disegniamo un grafico con il tempo in ascissa e la condizione $y_0$ in ordinata. Tenendo conto che $t=t_0+nh$, segue che al diminuire del passo $h$, sono necessari più passi per raggiungere un fissato istante di tempo $t$. Non si può scegliere come condizione di convergenza

$\underset{h_i\to 0}{lim}{y}_n(h_i)=y_0$

perché questo tiene conto solo della condizione iniziale e non dà importanza al fatto che dopo un fissato numero di passi, a differenti scelte di $h$ corrispondono istanti $t$ diversi. Bisogna quindi fissare $t$, e fare contemporaneamente

$\underset{h_i\to 0,n\to \mathrm{\infty }}{lim}{y}_n(h_i)=y_0$

Abbiamo un metodo convergente se per ogni problema ai valori iniziali

$\underset{h\to 0,t=t_0+nh}{lim}y(t)=y_0$

per ogni successione $\{y_n\}$ generata a partire da $y_0$.

Analisi di convergenza

Dato il metodo

$y_{n+1}=y_n+hf(t_n,y_n)$

Consideriamo l'errore all'istante $t_{n+1}$, e definiamo

$e_{n+1}=y(t_{n+1})-y_{n+1}=y(t_{n+1})-y_{n+1}^{\ast }+y_{n+1}^{\ast }-y_{n+1}$

dove $y_{n+1}^{\ast }$ è il valore che otteniamo se all'istante $t_{n+1}$ prendiamo come dato $y_n$ il valore esatto $y(t_n)$, cioè

$y_{n+1}^{\ast }=y(t_n)+hf(t_n,y(t_n))$

e rappresenta l'errore locale, di consistenza.

Possiamo sostituire i primi due addendi con l'errore di troncamento moltiplicato per $h$, quindi

$E_{n+1}=h{\tau }_{n+1}+y(t_n)+hf(t_n,y(t_n))-y_{n+1}$

e sostituendo l'espressione di $y_{n+1}=y_n+hf(t_n,y_n)$:

$E_{n+1}=h{\tau }_{n+1}+y(t_n)+hf(t_n,y(t_n))-y_n-hf(t_n,y_n)$

e siccome $E_n=y(t_n)-y_n$,

$E_{n+1}\le h{\tau }_{n+1}+E_n+h\ast f(t_n,y(t_n))-h\ast f(t_n,y_n)$

Passando ai moduli

$|E_{n+1}|\le h\ast |{\tau }_{n+1}(h)|+|E_n|+h\ast (|f(t_n,y(t_n))|-|f(t_n,y_n)|)$

Poniamo ${\tau }_h=\underset{n}{max}\{{\tau }_{n+1}(h)\}$ e teniamo conto della lipschitzianità di $f$, allora otteniamo

$\le h|{\tau }_h|+|E_n|+Lh(y(t_n)-y_n)\le h{\tau }_h+|E_n|+Lh\ast E_n\le h{\tau }_h+(1+L_h)|E_n|$

e avendo trovato la relazione ricorsiva si prosegue fino al passo 0:

$\le h{\tau }_h+(1+Lh)\ast (h{\tau }_h+(1+Lh)|E_{n-1}|)$

e siccome $E_0\le h{\tau }_h$ in conclusione si ha

$|E_{n+1}\le h{\tau }_h\ast (1+(1+Lh)+\dots +(1+Lh{)}^n)$

$\le h{\tau }_h\sum _{k=0}^n(1+LH{)}^k$

$=\frac{(1+Lh{)}^{n+1}-1}{Lh}h{\tau }_h$

$=\frac{(1+Lh{)}^{n+1}-1}{L}{\tau }_h$

Siccome in generale $1+x\le e^x$, si ha

$\le \frac{(e^{hL}{)}^{n+1}-1}{L}\ast {\tau }_h$

Siccome $t_{n+1}-t_0=h$ si ha

$=\frac{e^{L(t_{n+1}-t_0)}-1}{L}\ast {\tau }_h$

Se $Y\in C_2$, nel caso di Eulero esplicito si ha per la formula 1

$\tau =h/2f^{(2)}(\xi )$

quindi ${\tau }_h\le Mh/2$ dove $M=max|f^{(2)}(\xi )|$.

Allora segue che

$E_{n+1}\le \frac{e^{L(t_{n+1}-t_0)}-1}{L}\ast M/2h$

per ogni $n\ge 0$.

La quantità tende a 0 per $h\to 0$, con lo stesso ordine 1 dell'errore di troncamento locale.

Metodo di Eulero implicito

Il metodo di Eulero implicito è dato da

$y_{n+1}=y_n+hf(t_{n+1},y_{n+1})$

Consideriamo l'errore di troncamento

${\tau }_h(y,t,n)=1/h[y(t+h)-y(t)-h\ast f(t+h,y(t+h))]$

e tenendo conto della definizione del problema di Cauchy:

$=1/h[y(t+h)-y(t)-h{y}^{\prime }(t+h)]$

$=1/h[y(t)+h{y}^{\prime }(t)+h^2/2y^{(2)}(\xi )-y(t)-h(y^{\prime }(t)+h\ast y^{(2)}(\eta ))]$

$=1/h\ast h^2(c_1{y}^{″}(\xi )+c_2{y}^{″}(\eta )$

e abbiamo un ordine di consistenza $h=1$.

Metodo dei trapezi

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

A destra ci sono due valutazioni di $f$ e l'unica incognita è $y_{n+1}$.

Esercizio 7.1

Verificare, usando Taylor, che l'ordine di consistenza di questo metodo è 2.

Dim. L'errore di troncamento si definisce come

$\tau =1/h\ast [y(t+h)-y(t)-h/2(f(t_n,y_n)+f(t_{n+1},y_{n+1})]$

$\tau =1/h\ast [y(t+h)-y(t)-h/2y^{\prime }(t)-h/2y^{\prime }(t+h)]$

$\tau =1/h\ast [y(t)+h{y}^{\prime }(t)+h^2/2y^{″}(t)+h^3/3y^{(3)}(\xi )-y(t)-h/2y^{\prime }(t)-h/2y^{\prime }(t)-h/2y^{″}(t)h-h/2y^{(3)}(\eta )h^2/2]$

$\tau =1/h\ast [h^3(a{y}^{(3)}(\xi )-b{y}^{(3)}(\eta ))]=o(h^2)$

Considerando il metodo di Eulero implicito:

$y_{n+1}=y_n+hf(t_{n+1},y_{n+1})$

si ha un'equazione non lineare e può essere trattata in due modi.

1. Definiamo il metodo di punto fisso:

$y_{k+1}=y_k+hf(t_0+nh,y_k),$

Affinché il metodo converga si richiede che $|g^{\prime }(z)|<1$. In questo caso, con $z=y_{k+1}$ si ha:

$g^{\prime }(z)=h\frac{\partial f}{\partial z}<HL$

quindi la costante di Lipschitz influenza la scelta del passo di discretizzazione. Se la costante di Lipschitz è grande, bisogna scegliere un passo di discretizzazione piccolo per poter avere convergenza.

2. Applichiamo il metodo di Newton. Poniamo $z=y_{n+1}$ e cerchiamo lo zero di

$\varphi (z)=z-y_n+hf(t_{n+1},z)$

con il metodo

$z_{k+1}=z_k-\frac{\varphi (z_k)}{{\varphi }^{\prime }(z_k)}$

con vettore d'innesco $z_0=y_n$.

${\varphi }^{\prime }(z)=1-h\frac{\partial f}{\partial z}(t_{n+1})$

Come test d'arresto si ha il test dell'incremento.