4.4 METODI DI INTEGRAZIONE NUMERICA

Sommario Quando un sistema dinamico non è integrabile mediante un'espressione analitica esplicita, il che si verifica nella stragrande maggioranza dei casi, può essere necessario usare un procedimento che fornisca un'approssimazione della soluzione. Un metodo di integrazione numerica è quindi un algoritmo che definisce un sistema dinamico discreto le cui soluzioni $X_k$ approssimano la soluzione del sistema dinamico continuo su una successione di tempi, per esempio $X_k\simeq X(kh)$. Esistono moltissimi metodi di integrazione numerica, e qui presentiamo soltanto i più semplici, commentando i rispettivi vantaggi e svantaggi.

Metodi simplettici a scorrimento

Come si è visto nella dimostrazione del teorema della mappa standard, gli scorrimenti sono mappe conservative del piano in se stesso. Perciò, se si vuole in qualche modo approssimare un sistema continuo conservativo con un sistema discreto pure conservativo, si possono cercare dei metodi di approssimazione che consistano in una sequenza di scorrimenti. Questi metodi vanno sotto il nome di metodi simplettici a scorrimento. Per esempio la mappa standard è un metodo simplettico a due scorrimenti.

Per integrare il sistema dinamico derivato da un sistema newtoniano di dimensione 1

$$\left\{\begin{array}{lcl} {\displaystyle \dot x} & {\displ... ...} & {\displaystyle=} &{\displaystyle f(x)} \end{array}\right.$$

usiamo una sequenza di tre scorrimenti, alternando l'aggiornamento della $x$ e della $y$:

$$\left\{\begin{array}{lcl} {\displaystyle x_M} & {\displays... ...yle=} &{\displaystyle x_M+c_2\,h\, y'} \end{array}\right. \;.$$

Cerchiamo quindi di determinare i tre coefficienti $c_1,d_1,c_2$ in modo che il metodo abbia ordine di troncamento 2, cioè con errore di troncamento locale $O(h^3)$. Per ottenere questo, basta imporre che risolva esattamente il caso $f(x)=2$, con soluzione per esempio $x=t^2,\;y=2t$:

$$\left\{\begin{array}{lcl} {\displaystyle x_M} & {\displays... ...} & {\displaystyle=} &{\displaystyle 1} \end{array}\right.\;.$$

Esiste una sola soluzione $c_1=c_2=1/2\;, d_1=1$ che ha la notevole proprietà aggiuntiva di essere reversibile, cioè di dare $(x,y)$ a partire da $(x',y')$ se il metodo viene applicato con passo $-h$; per avere questa proprietà occorre che $c_1=c_2$.

$$\left\{\begin{array}{lcl} {\displaystyle x_M} & {\displays... ...style=} &{\displaystyle x_M+(h/2)\, y'} \end{array}\right.\;.$$

Per una trattazione più ampia dei metodi di integrazione simplettici a scorrimento si veda [Yoshida 93].

Esempio:

Figura 4.7: Il metodo simplettico a tre scorrimenti e di ordine due applicato al pendolo nonlineare conservativo, con passo h=1. Rispetto alla mappa standard usata con lo stesso passo, le soluzioni si muovono vicino a curve più simili alle linee di livello dell'integrale dell'energia, con l'eccezione della regione attorno alla separatrice.


Applichiamo il metodo simplettico e reversibile a tre scorrimenti (di ordine 2) al pendolo nonlineare conservativo $\ddot x= -\sin(x)$;

$$\left\{\begin{array}{lcl} {\displaystyle x_M} & {\displays... ...laystyle=} &{\displaystyle x_M+(h/2)\, y'} \end{array}\right.$$


Per confrontare questa approssimazione con quella della mappa standard del pendolo, si veda la Figura 4.7, in cui abbiamo usato lo stesso valore $h=1$ della Figura 4.6.

Esercizio Nell'esempio precedente, mostrare che i punti fissi sono ancora $(x,y)=(0,0)$ e $(x,y)=(\pi,0)$. (Soluzione)

Metodi di Runge-Kutta

Una classe importante di metodi di integrazione numerica va sotto il nome di Runge e Kutta. Dato il sistema dinamico

$$\dot X= F(X)\;,$$

un metodo di Runge-Kutta ad $s$ passi intermedi è un sistema dinamico discreto la cui mappa $X\mapsto X'$ è della forma

$$X'=X+ h\sum_{i=1}^s b_i F(Z_i) = R_h(X)$$

dove $h$ è il passo e $Z_i, i=1,\ldots ,s$ sono approssimazioni di valori intermedi ottenuti con le formule

$$Z_i= X+ h\sum_{j=1}^s a_{ij}F(Z_j)\;.$$

Se la matrice $A=[a_{ij}]$ è tale che $a_{ij}=0$ per $j\geq i$, il metodo è esplicito ed i valori intermedi $Z_j$ si possono calcolare in sequenza, altrimenti il metodo è implicito e le equazioni per gli $Z_j$ vanno risolte ad ogni passo. La soluzione di questo sistema di equazioni esiste ed è unica per $h$ abbastanza piccolo; infatti si tratta di un'equazione del tipo del punto unito, e per $h$ piccolo il secondo membro è una contrazione (in un intorno di $X$).

Lo stesso metodo si potrebbe applicare anche ad un sistema di equazioni differenziali non autonomo

$$\dot X=F(t,X)$$

con il procedimento di omogeneizzazione, cioè aggiungendo un'equazione per la variabile $x_0=t$:

$$\left\{\begin{array}{lcl} {\displaystyle \dot x_0} & {\dis... ...isplaystyle=} &{\displaystyle F(x_0,X)} \end{array}\right.\;.$$

Questo significa usare la formula modificata per il calcolo dei valori intermedi

$$Z_i= X+ h\sum_{j=1}^s a_{ij}F(t_j,Z_j)\;,$$

dove i tempi intermedi sono $t_i=t+hc_i$ dati dai coefficienti

$$c_i=\sum_{j=1}^s a_{ij}\;.$$

Lo scopo di un metodo di Runge-Kutta, come di ogni metodo di integrazione numerica, è quello di approssimare $X(t+h)$ con $X'$, supponendo che $X=X(t)$. Il principale requisito per un metodo di integrazione numerica accurato è quello di avere un ordine di troncamento più grande possibile. L'ordine è $q$ se la differenza tra l'approssimazione discreta $R^h$ e la soluzione esatta $X(t+h)$ è infinitesimo di ordine superiore a $q$ rispetto al passo $h$ per $h\to 0$:

$$\vert X(t+h)-R^h[X(t)]\vert={\cal o}(h^q)\;.$$

Per imporre che il metodo sia di ordine $q$ si può richiedere che risolva esattamente le equazioni differenziali che hanno come soluzioni dei polinomi di grado fino a $q$: in altri termini le equazioni

$$\frac{d{x}}{d{t}}=kt^{k-1}$$

devono essere risolte esattamente. Applicando un metodo di Runge-Kutta a queste equazioni si trovano delle condizioni sui coefficienti $c_i,b_i, a_{ij}$. Se $x(t)=t^k$ è la soluzione esatta,

$$x(h)=x(0)+h\sum_{i=1}^s b_i\,k\,t_i^{k-1}= h\sum_{i=1}^s b_i\,k\,(h\,c_i)^{k-1}$$

$$h^k=h^k \sum_{i=1}^s b_i\,k\,c_i^{k-1}\Longrightarrow \sum_{i=1}^s c_i^{k-1}b_i=\frac 1k$$

Si noti che se i $c_i$ sono già stati scelti, questo è un sistema lineare con incognite $b_i$ e come matrice dei coefficienti la matrice di Vandermonde $V=[v_{ij}]=[c_i^{k-1}]$; il sistema è quadrato se si richiede che l'ordine sia $s$, e in tal caso ha soluzione unica se i $c_i$ sono diversi tra loro (perché il determinante della matrice di Vandermonde $s\times s$ è diverso da zero).

Vediamo allora alcune soluzioni di queste equazioni: per esempio, se cerchiamo un metodo di ordine 2 a 2 passi intermedi, abbiamo un sistema di due equazioni in quattro incognite

$$\left\{\begin{array}{lcl} {\displaystyle b_1+b_2} & {\disp... ...isplaystyle=} &{\displaystyle \frac{1}{2}} \end{array}\right.$$

che ha soluzione per ogni coppia di valori diversi di $c_1,c_2$ arbitrariamente scelti; per esempio per $c_1=0,\;c_2=1$ abbiamo $b_1=b_2=1/2$. Bisogna poi soddisfare alle condizioni per cui le $c_i$ sono somma delle righe di $A$; l'unica soluzione che dia un metodo esplicito (con $a_{11}=a_{22}=a_{12}=0$) è $a_{21}=1$. In conclusione un metodo di Runge-Kutta esplicito di ordine 2 è dato dalle formule

$$\left\{\begin{array}{lcl} {\displaystyle Z_1} & {\displays... ... {\displaystyle=} &{\displaystyle X+hF(X)} \end{array}\right.$$

cioè

$$X'=X+\frac h2 \left[ F(X)+ F(X+h\,F(X))\right]\;.$$

Se invece cerchiamo un metodo ad un solo passo intermedio e di ordine 2, è chiaro che non potrà essere esplicito (avendo una matrice $A$ di tipo $1\times 1$, questa non può essere nulla sulla diagonale, altrimenti è nulla e si ottiene un metodo di ordine 1). Abbiamo allora due equazioni, nelle due sole incognite $b_1,\,c_1$:

$$\left\{\begin{array}{lcl} {\displaystyle b_1\,c_1^0} & {\d... ...isplaystyle=} &{\displaystyle \frac 12} \end{array}\right.\;.$$

La soluzione in questo caso è unica, con $a_{11}=1/2$; il metodo è implicito, nel senso che il valore intermedio appare nell'equazione del punto unito:

$$Z_1=X+\frac h2\,F\left(t+\frac h2, Z_1\right)$$

e solo dopo aver risolto questa (per esempio per approssimazioni successive) si può calcolare il nuovo punto:

$$X'=X+h\,F\left(t+\frac h2,Z_1\right)\ .$$

Il metodo di Runge-Kutta implicito ad un passo intermedio ha notevoli proprietà qualitative, delle quali la più importante è quella del teorema seguente.

Il metodo di Runge-Kutta implicito ad 1 passo intermedio e di ordine 2 è conservativo, nel senso che se applicato ad un sistema newtoniano ad un grado di libertà fornisce un sistema dinamico discreto conservativo.

Dimostrazione omessa.

Figura 4.8: L'equazione del pendolo nonlineare conservativo integrata con il metodo di Runge-Kutta implicito di ordine 2 ad un solo passo intermedio, e con un passo molto lungo (h=1). Benché il risultato non sia quantitativamente molto accurato, l'andamento generale è molto simile a quello del sistema continuo, con una sorprendentemente piccola regione caotica vicino alla separatrice.

Esempio:

Applichiamo il metodo di Runge-Kutta implicito di ordine 2 e ad un solo passo intermedio al pendolo nonlineare conservativo $\ddot x= -\sin(x)$. Ponendo come al solito $y=\dot x$, l'equazione del punto unito per il valore intermedio $(x_M,y_M)$ è


$$\left[\begin{array}{c}{x_M}\\ {y_M}\end{array}\right]= \le... ...ft[\begin{array}{c}{y_M}\\ {-\sin(x_M)}\end{array}\right]\;;$$


una volta trovato il punto intermedio (con il metodo delle approssimazioni successive), un passo di lunghezza $h$ è approssimato dal sistema dinamico discreto:

$$\left[\begin{array}{c}{x'}\\ {y'}\end{array}\right]=\left[... ...ft[\begin{array}{c}{y_M}\\ {-\sin(x_M)}\end{array}\right]\;.$$



Poiché il sistema dinamico così ottenuto è conservativo, le proprietà di stabilità delle approssimazioni lineari sono conservate; anche alcune proprietà qualitative, come la presenza di curve invarianti, sono conservate, e la regione caotica che circonda la separatrice è di area molto ridotta, persino per valori di $h$ che con altri metodi di discretizzazione danno una larga regione caotica, per esempio $h=1$ (Figura 4.8).

La ragione per cui questo metodo è particolarmente efficace nel generare una piccola regione caotica attorno alla separatrice ha a che fare con un'altra proprietà aggiuntiva rispetto a quella di essere conservativo: i metodi di questa classe (di ordine $2s$ ad $s$ passi intermedi, come il metodo di Runge-Kutta-Gauss) hanno la proprietà di preservare esattamente le forme quadratiche che siano integrali primi. Quindi nell'intorno del punto di sella nonlineare il comportamento delle soluzioni è descritto in modo molto accurato, perché nell'equazione linearizzata ci sarebbe un integrale primo quadratico.

Costruire metodi con ordine di troncamento superiore al secondo non è tanto facile, ma si conoscono molte soluzioni. Prendiamo per esempio l'ordine 4: il metodo di Runge-Kutta classico (il primo ad essere stato scoperto, ed ancora largamente in uso) è di ordine 4 con 4 passi intermedi, e con la seguente scelta dei tempi intermedi: $c_1=0, c_2=c_3= 1/2,c_4=1$.

Figura 4.9: Il pendolo nonlineare conservativo integrato con il metodo di Runge-Kutta classico, di ordine 4 e a 4 passi intermedi, con un passo abbastanza corto (h=0.2). Il risultato appare molto simile alla soluzione esatta, ma se l'integrazione continuasse per un tempo molto lungo l'effetto dissipativo introdotto dalla discretizzazione diventerebbe evidente, per esempio osservando l'andamento nel tempo dell'integrale primo dell'energia.

Da qui il sistema di 4 equazioni lineari nelle 4 incognite $b_i$

$$\left\{\begin{array}{lcl} {\displaystyle b_1+(b_2+b_3)+b_4... ...4} & {\displaystyle=} &{\displaystyle 1/4} \end{array}\right.$$

con soluzione unica: $b_1=b_4=1/6, b_2=b_3=1/3$. Si definisce così il metodo di Runge-Kutta classico, che può essere espresso in termini dei valori intermedi del secondo membro

$$\left\{\begin{array}{lcl} {\displaystyle K_1} & {\displayst... ...array}\right. \ \ , \ \ X'=X+\frac h6\,(K_1+K_4+2K_2+2K_3) \ .$$

Un metodo del quarto ordine, anche se non è conservativo, può essere usato con un passo abbastanza piccolo per ottenere soluzioni che siano quantitativamente molto simili a quelle esatte (Figura 4.9). Bisogna però fare attenzione a non usare il risultato dell'integrazione numerica per dedurne delle proprietà qualitative; per esempio, in un sistema newtoniano si introduce una ``dissipazione numerica''.

A titolo di esempio di un metodo ancora più accurato, che però richiede un numero maggiore di valutazioni del secondo membro per ogni passo, diamo i coefficienti dell'unico metodo implicito di ordine 4 a due passi intermedi, detto metodo di Runge-Kutta-Gauss:

$$\left\{\begin{array}{lcl} {\displaystyle c_1} & {\displays... ...sqrt{3})/12}\\ {(3+2\sqrt{3})/12}&{1/4}\end{array}\right]\;.$$

Anche il metodo di Gauss è una discretizzazione conservativa, ed ha proprietà di stabilità molto forti: per esempio, fornisce una mappa stabile se usato per un sistema lineare stabile. Per una trattazione completa delle proprietà note dei metodi di integrazione numerica di Runge-Kutta si veda [Butcher 87].

Esercizio Per esercitarsi sui metodi di integrazione numerica occorre usare un calcolatore. Dato un sistema newtoniano $\ddot x= f(x)$, e scelto un linguaggio di programmazione con capacità grafiche, si scrivano e si collaudino (osservando le orbite rappresentate graficamente) programmi che eseguano le seguenti discretizzazioni:

a)

mappa standard

b)

metodo simplettico a tre scorrimenti

c)

Runge-Kutta-Gauss implicito a un passo intermedio

d)

Runge-Kutta classico.

(Soluzione)