Subsections


4.2 METODO DI EULERO

Sommario Il più semplice procedimento per approssimare un sistema dinamico continuo con uno discreto è quello di Eulero. Le soluzioni del sistema discreto così ottenuto approssimano le soluzioni del sistema continuo, in un senso che può essere reso rigoroso con la nozione di convergenza uniforme. La differenza tra le due, cioè l'errore di discretizzazione, è una funzione rapidamente crescente con il passare del tempo, per cui la potenza di calcolo richiesta per ottenere soluzioni accurate è notevole.

Metodo di Eulero nel caso lineare

Una definizione alternativa della funzione esponenziale di variabile reale è:

\begin{displaymath}
\lim_{m\to +\infty} \left(1+\frac{at}{m}\right)^m =e^{at}\ .
\end{displaymath}

Vale l'analogo matriciale: se $A$ è una matrice quadrata $n\times n$:

\begin{displaymath}
\exp(A\,t)=\lim_{m\to +\infty} \left(I +\frac {A\,t}m\right)^m
\end{displaymath}

per ogni $t$ in ${\bf R}$. Basta scrivere la differenza

\begin{displaymath}
\exp(A\,t) - \left(I +\frac {A\,t}m\right)^m = \sum_{k=0}^\...
...1{k!} - \left({m\atop k}\right)\frac 1{m^k}\right]\, A^k\, t^k
\end{displaymath}

(i coefficienti binomiali sono nulli quando $k>m$); la serie a secondo membro ha coefficienti non negativi, e se si passa alle norme si trova la maggiorazione

\begin{displaymath}
\left \Vert \exp(A\,t) - \left (I +\frac {A\,t}m\right )^m\r...
...{a\,\vert t\vert} - \left(1+\frac{a\,\vert t\vert}{m}\right)^m
\end{displaymath}

dove $a=\vert\vert A\vert\vert$ è la norma uniforme della matrice $A$. Nel caso scalare la successione converge uniformemente sull'intervallo $[-t,t]$, con $t$ arbitrario (ma fissato). Per il teorema della convergenza in norma, anche la definizione alternativa di esponenziale nel caso matriciale converge uniformemente su $[-t,t]$. Da questo segue (per una delle proprietà della norma uniforme) che

\begin{displaymath}
\exp(A\,t)\,X_0 = \lim_{m\to +\infty} \left(I +\frac {A\,t}m\right)^m\,X_0
\end{displaymath}

uniformemente per $t$ in ogni intervallo limitato, e per ogni condizione iniziale $X_0$. L'aspetto interessante di questa formula è che al primo membro abbiamo il flusso integrale di un sistema dinamico continuo lineare; a secondo membro sotto segno di limite abbiamo l'analoga rappresentazione di tutte le soluzioni di un sistema dinamico discreto lineare, e precisamente

\begin{displaymath}
\frac{d{X}}{d{t}} = A\, X \hspace{5mm},\hspace{5mm}X_{k+1}=\left(I +\frac {A\,t}m\right)\,X_k\;.
\end{displaymath}

La relazione tra il sistema dinamico continuo e quello discreto ottenuta qua sopra è un caso particolare (lineare) del metodo di Eulero. La motivazione geometrica dell'approssimazione di Eulero può essere apprezzata dalla Figura 4.3. Se si deve calcolare la soluzione con condizione iniziale $X(0)=X_0$ al tempo $t=\delta$, si suddivide l'intervallo $[0,\delta]$ in cui varia $t$ in $m$ sottointervalli, ciascuno di lunghezza pari al passo d'integrazione $h=\delta/m$, con $m$ abbastanza grande. Allora si approssima $X(\delta/m)$ con $X_1=X_0+hA\,X_0$, $X(2\delta/m)$ con $X_2=X_1+hA\,X_1$, eccetera, ad ogni passo rimpiazzando la vera soluzione con la sua approssimazione lineare (cioè con il suo differenziale).

Benché nel sistema dinamico discreto appaiano solo i punti che approssimano la soluzione ai tempi $t=k\delta/m$ con $k$ intero, si può immaginare di descrivere un'approssimazione per tutti i tempi $t\in [0,\delta]$ usando l'interpolazione lineare (cioè ancora lo stesso differenziale) per i punti intermedi. Il poligono di Eulero ottenuto congiungendo i punti $(k\delta/m, X_k)$ è il grafico di una funzione di $t$ che è continua ma non differenziabile, e che approssima la soluzione nel senso che la successione di poligoni di Eulero ottenuti per $m$ crescente tende alla soluzione del sistema dinamico continuo, uniformemente su $[0,\delta]$.

Figura 4.3: Il poligono di Eulero: ad ogni passo temporale si approssima la soluzione del sistema dinamico continuo con il suo differenziale nel punto precedentemente calcolato; in questo esempio la soluzione dell'equazione y'=y è approssimata con il metodo di Eulero usando 2,4,6 ed 8 passi.
\begin{figure}{\centerline{\epsfig{figure=figures/figeulero.ps,height=8cm}}}
\end{figure}

L'errore di troncamento locale commesso in un passo di lunghezza $h>0$ del metodo di Eulero al posto della soluzione esatta si può maggiorare usando il confronto in norma con l'esponenziale:

\begin{displaymath}
\vert\vert\exp(Ah)-(I+Ah)\vert\vert\leq e^{a\,h}-(1+ah)= \fr...
...2\,h^2}2\,
e^{a\,\theta\,h}\leq \frac{a^2\,h^2}2\, e^{a\,h}\;,
\end{displaymath}

dove $a=\vert\vert A\vert\vert$ e $0<\theta<1$; il penultimo passaggio è il resto di Lagrange della formula di Taylor per la funzione esponenziale di variabile reale.

Esempio:


Esercizio Discretizzare, con il metodo di Eulero, il sistema dinamico

\begin{displaymath}\dot x = \omega\,y \hspace{5mm},\hspace{5mm}\dot y = \omega \,x\end{displaymath}

e discutere l'andamento delle soluzioni per il tempo discreto che tende a $\pm\infty$. (Soluzione)

Convergenza del metodo di Eulero

Supponiamo di dover approssimare la soluzione di un sistema dinamico nonlineare $\dot X = F(X)$ con condizione iniziale $X(0)=X_0$. Il metodo di Eulero consiste nell'approssimare la soluzione al tempo $h$ con l'approssimazione lineare:

\begin{displaymath}
X(h)\simeq X_1 = X_0 + h\,F(X_0)
\end{displaymath}

La soluzione a tempi successivi può essere approssimata ripetendo il procedimento, cioè con la successione definita per ricorrenza:

\begin{displaymath}\left\{\begin{array}{lcl}
{\displaystyle X_0} & {\displays...
...ystyle=} &{\displaystyle X_k + h\, F(X_k)}
\end{array}\right.
\end{displaymath}

Se si desidera avere delle informazioni sui valori assunti da $X(t)$ per valori di $t$ intermedi tra $kh$ e $(k+1)h$ si può ricorrere ad un'interpolazione, per esempio lineare ottenendo ancora il poligono di Eulero. L'errore di troncamento locale è sempre infinitesimo del secondo ordine rispetto ad $h$:

\begin{displaymath}
X(t+h)= X(t) + h\,\dot X(t) + O(h^2) \ .
\end{displaymath}

L'errore di troncamento accumulato è la differenza tra la successione definita per ricorrenza dal metodo di Eulero e la soluzione del sistema dinamico con le stesse condizioni iniziali:

\begin{displaymath}
\Delta^{(k)}=X(k\,h)-X_k \ .
\end{displaymath}

Se il campo vettoriale $F(X)$ è lipschitziano di costante $L$ e limitato in modulo dalla costante $M$, cioè se

\begin{displaymath}
\vert F(X)-F(Y)\vert\leq L\,\vert X-Y\vert \hspace{5mm},\hspace{5mm}\vert F(X)\vert \leq M
\end{displaymath}

per ogni $X,Y$ in ${\bf R}^n$, allora l'errore accumulato $\Delta^{(k)}$ dopo il tempo $t=k\,h$ soddisfa alla diseguaglianza

\begin{displaymath}
\vert\Delta^{(k)}\vert=\delta_k\leq M\,h\,\left [e^{\sqrt{n}\, L\, t}-1 \right]\;.
\end{displaymath}

Dimostrazione:

 C.D.D.


La maggiorazione per l'errore di troncamento accumulato che si ricava in generale, cioè per ogni possibile sistema dinamico, si rivela pessimistica in casi specifici. Che l'accumulazione dell'errore possa essere esponenziale è inevitabile, ma in molti casi l'esponente sarà più piccolo della costante di Lipschitz, e sarà piuttosto legato al massimo esponente di Lyapounov della soluzione che si cerca di approssimare.

Andrea Milani 2009-06-01