Calcolo numerico/Algoritmi per la risoluzione

Metodi a un passo

Dato il problema

${\begin{cases}y'(t)=f(t,y)\\y(t_{0})=y_{0}\end{cases}}$

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,\quad 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'(t)$

ed esplicitando $y(t+h)$:

$y(t+h)=y(t)+h*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\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*\Phi (t_{n},y(t_{n}),f(t_{n},y(t_{n})));h)$

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

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

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

dove $\tau _{n+1}$ è detto errore di troncamento locale. $\varepsilon _{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

Definizione

Si definisce errore di troncamento locale (LTE) la quantità

$\tau (t,y,h)={\frac {y(t+h)-y(t)}{h}}-\Phi (t,y,f,h)$

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

Consistenza

Definizione

Il metodo è consistente con il problema ai valori iniziali se vale che

$\lim _{h\to 0}\tau (t,y,h)=0$

con $t=t_{0}+nh$. Se per ogni istante $t\in I_{0}$ la soluzione $y(t)$ è tale che $\tau (t,y,h)=o(h^{p})$, allora il metodo è consistente di ordine $p$.

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

Considerando il metodo di Eulero esplicito:

$\Phi (t,y,f,h)=f(t,y(t))$

$\tau (t,y(t),h)=1/h*[y(t+h)-y(t)-h\Phi (t,y,h)]$

$=1/h*[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'(t)$:

$=1/h*[y(t)+hy'(t)+h^{2}/2y^{(2)}(\xi )-y(t)-hy'(t)]$

e in seguito a cancellazioni:

$=1/h*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 )\quad {\hbox{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$.

Definizione

Un metodo si dice convergente se per ogni problema ai valori iniziali tale che valga il teorema di esistenza e unicità, è verificata la seguente proprietà: per ogni $n$,

$|y_{n}-y(t_{n})|\leq C(h)$

dove

$\lim _{h\to 0}C(h)=0$

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

$\lim _{h_{i}\to 0}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

$\lim _{h_{i}\to 0,n\to \infty }y_{n}(h_{i})=y_{0}$

Abbiamo un metodo convergente se per ogni problema ai valori iniziali

$\lim _{h\to 0,t=t_{0}+nh}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}^{*}+y_{n+1}^{*}-y_{n+1}$

dove $y_{n+1}^{*}$ è 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}^{*}=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}\leq h\tau _{n+1}+E_{n}+h*f(t_{n},y(t_{n}))-h*f(t_{n},y_{n})$

Passando ai moduli

$|E_{n+1}|\leq h*|\tau _{n+1}(h)|+|E_{n}|+h*(|f(t_{n},y(t_{n}))|-|f(t_{n},y_{n})|)$

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

$\leq h|\tau _{h}|+|E_{n}|+Lh(y(t_{n})-y_{n})\leq h\tau _{h}+|E_{n}|+Lh*E_{n}\leq h\tau _{h}+(1+L_{h})|E_{n}|$

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

$\leq h\tau _{h}+(1+Lh)*(h\tau _{h}+(1+Lh)|E_{n-1}|)$

e siccome $E_{0}\leq h\tau _{h}$ in conclusione si ha

$|E_{n+1}\leq h\tau _{h}*(1+(1+Lh)+\dots +(1+Lh)^{n})$

$\leq 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\leq e^{x}$, si ha

$\leq {\frac {(e^{hL})^{n+1}-1}{L}}*\tau _{h}$

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

$={\frac {e^{L(t_{n+1}-t_{0})}-1}{L}}*\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}\leq Mh/2$ dove $M=\max |f^{(2)}(\xi )|$.

Allora segue che

$E_{n+1}\leq {\frac {e^{L(t_{n+1}-t_{0})}-1}{L}}*M/2h$

per ogni $n\geq 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*f(t+h,y(t+h))]$

e tenendo conto della definizione del problema di Cauchy:

$=1/h[y(t+h)-y(t)-hy'(t+h)]$

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

$=1/h*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*[y(t+h)-y(t)-h/2(f(t_{n},y_{n})+f(t_{n+1},y_{n+1})]$

$\tau =1/h*[y(t+h)-y(t)-h/2y'(t)-h/2y'(t+h)]$

$\tau =1/h*[y(t)+hy'(t)+h^{2}/2y''(t)+h^{3}/3y^{(3)}(\xi )-y(t)-h/2y'(t)-h/2y'(t)-h/2y''(t)h-h/2y^{(3)}(\eta )h^{2}/2]$

$\tau =1/h*[h^{3}(ay^{(3)}(\xi )-by^{(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'(z)|<1$. In questo caso, con $z=y_{k+1}$ si ha:

$g'(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

$\phi (z)=z-y_{n}+hf(t_{n+1},z)$

con il metodo

$z_{k+1}=z_{k}-{\frac {\phi (z_{k})}{\phi '(z_{k})}}$

con vettore d'innesco $z_{0}=y_{n}$.

$\phi '(z)=1-h{\frac {\partial f}{\partial z}}(t_{n+1})$

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