4.3 DISCRETIZZAZIONE CONSERVATIVA

Sommario Equazioni differenziali di ordine superiore possono essere discretizzate utilizzando diverse approssimazioni per le derivate non solo prime ma anche di ordine superiore. Alcune di queste approssimazioni, come quella che è equivalente al metodo di Eulero, hanno però delle proprietà qualitative indesiderabili, per esempio modificano le proprietà di stabilità delle soluzioni. Altre approssimazioni sono conservative, cioè se applicate ad un sistema dinamico continuo conservativo producono un sistema dinamico discreto ugualmente conservativo. Non esiste però un metodo che consenta di preservare sempre gli integrali primi del sistema continuo, quindi le discretizzazioni vanno sempre interpretate come approssimazioni anche nel senso qualitativo, perché riproducono in modo imperfetto proprietà qualitative come l'integrabilità.

Discretizzazione di ordine superiore

Per approssimare le soluzioni di un'equazione differenziale di ordine superiore al primo con un sistema dinamico discreto si possono seguire due approcci: il più semplice è quello di ridursi ad un sistema dinamico e poi usare una discretizzazione per quest'ultimo, per esempio il metodo di Eulero $S\simeq I+ h\,D$ (oppure altri più accurati come sarà descritto nella Sezione 4.4).

Può essere più conveniente tentare di scrivere direttamente una formula di approssimazione per le derivate di ordine superiore. Queste approssimazioni possono essere scritte in termini di un'algebra di operatori.

Definizione:

Data una funzione $f(t)$ definita per ogni $t\in {\bf R}$, gli operatori spostamento

$$(Sf)(t)=f(t+h)\;,$$


differenza in avanti

$$(\Delta_+ f)(t)=f(t+h)-f(t)\;,$$


ed identità:

$$(If)(t)=f(t)$$


mandano lo spazio di tali funzioni in sé. Valgono le relazioni:

$$\Delta_+= S-I\hspace{5mm},\hspace{5mm}I\,S=S\,I=S\hspace{5mm},\hspace{5mm}I\,\Delta_+=\Delta_+\,I=\Delta_+$$


e tutti questi operatori commutano, in particolare

$$\Delta_+\,S= S\, \Delta_+$$

Segue immediatamente dalle proprietà di commutazione che si può applicare la formula del binomio di Newton per calcolare le potenze di $\Delta_+$: la differenza k-esima in avanti è

$$\Delta_+^k = \sum_{j=0}^k \left({k \atop j}\right)\,S^j\, I^{k-j}\, (-1)^{k-j}\;;$$

per esempio

$$\Delta_+^2= S^2 -2\,S + I\;.$$

Queste formule possono essere impiegate direttamente in approssimazioni discrete delle derivate di ordine superiore, del tipo

$$\frac{\Delta_+}h=D+ O(h)\hspace{5mm},\hspace{5mm}\frac{\Delta_+^2}{h^2}=D^2 +O(h)$$

dove $D f(t)=df/dt(t)$ è l'operatore derivata.

Esempio:

Riprendiamo l'oscillatore armonico come equazione differenziale del secondo ordine


$$\frac{d^2{x}}{d{t}^2} = -\omega^2\,x$$


ed usiamo direttamente la discretizzazione della derivata seconda

$$\frac {\Delta_+^2}{h^2}\,x_k=-\omega^2\,x_k\Leftrightarrow x_{k+2}-2x_{k+1} +x_k= -\omega^2\,h^2\,x_k\;.$$



Si ottiene l'equazione alle differenze finite lineare del secondo ordine


$$x_{k+2}=2x_{k+1}+ \left(-1-\omega^2\,h^2\right)\,x_k \Longle... ...A\;\left[\begin{array}{c}{x_k}\\ {x_{k+1}}\end{array}\right]$$




$$A=\left[\begin{array}{cc}{0}&{1}\\ {-1-\omega^2\,h^2}&{2}\end{array}\right]\;,$$


la cui soluzione si calcola usando gli autovalori di $A$:

$$det[A-\lambda\,I]=\lambda^2 -2\,\lambda+\left(1+\omega^2\,h^2\right)=0 \Longrightarrow \lambda = 1\pm J\,\omega\,h\;.$$


La discretizzazione diretta della derivata seconda fatta con le differenze in avanti ed il metodo di Eulero applicato al corrispondente sistema dinamico hanno gli stessi autovalori, quindi sono lo stesso sistema dinamico discreto, visto in un diverso sistema di riferimento. Perciò hanno anche lo stesso inconveniente: la soluzione del sistema dinamico discreto diverge, seguendo un cammino a spirale il cui raggio cresce con $k$ come $(1+\omega^2\,h^2)^{k/2}$.

La matrice $A$ può essere scritta nella forma canonica


$$B=\left[\begin{array}{cc}{1}&{\omega\,h}\\ {-\omega\,h}&{1}\end{array}\right]\;.$$


Il cambiamento di coordinate che riduce alla forma canonica si può ottenere con il solito metodo a base di autovettori, ma anche con un altro ragionamento. Per definire la successione $x_k$ occorrono due condizioni iniziali, $x_0$ ed $x_1$. Invece, se si adotta il punto di vista che risulta dalla discretizzazione con il metodo di Eulero dell'equazione ridotta al primo ordine, la condizione iniziale è $(x_0,y_0)$, con $\dot x(0)=\omega \,y_0$. Adottando la discretizzazione $\Delta_+/h\simeq D$ si trova

$$x_1=x_0 +\Delta_+x_0\simeq x_0 + h\,\omega\,y_0\;,$$


e in generale, per ogni $k$, la relazione:

$$x_{k+1}=x_k +\Delta_+x_k\simeq x_k + h\,\omega\,y_k$$


per cui il cambiamento di coordinate cercato è

$$\left[\begin{array}{c}{x_k}\\ {x_{k+1}}\end{array}\right] ... ...gin{array}{cc}{1}&{0}\\ {1}&{\omega\,h}\end{array}\right]\;.$$


Si può verificare che $C^{-1}\,A\, C=B$, cioè che la matrice $C$ è associata al cambiamento di coordinate che riduce $A$ in forma canonica reale.

Esercizio Se si usasse una discretizzazione basata sulla differenza all'indietro:

$$\Delta_-f(t)=f(t)-f(t-h) =I-S^{-1}\hspace{5mm},\hspace{5mm}\... ... \hspace{5mm},\hspace{5mm} \frac{\Delta_-^2}{h^2}\simeq D^2\;,$$

quali sarebbero le proprietà delle soluzioni della discretizzazione dell'oscillatore armonico? Quale sarebbe il limite della soluzione $x_k$ per $k\to+\infty$? (Soluzione)

Differenze centrali

Per equazioni differenziali che contengono solo le derivate di ordine pari si può trovare un procedimento di approssimazione che ha proprietà molto migliori, in particolare che non introduce moltiplicatori di Lyapounov spurii come il metodo di Eulero.

Definizione:

Data una funzione $f(t)$ definita per ogni $t$, l'operatore differenza centrale è definito da:

$$\Delta_0 f(t)=f(t+\frac h2) -f(t-\frac h2)$$

Usando le regole di commutazione e il binomio di Newton si verifica che le potenze pari di $\Delta_0$ non usano ``mezzi passi'', quindi sono definite come operatori su di una successione; per esempio,

$$\Delta_0^2=S-2\,I +S^{-1}\;.$$

Le differenze centrali sono quindi utili per definire una discretizzazione di equazioni senza derivate di ordine dispari, per esempio usando

$$\frac{\Delta_0^2}{h^2}=D^2 + O(h)\;.$$

Esempio:

La discretizzazione dell'oscillatore armonico usando le differenze centrali è

$$\Delta_0^2 \, x_k=-\omega^2\, h^2\,x_k \Leftrightarrow x_{k+1}=2\, x_k -x_{k-1}-\omega^2\,h^2\,x_k\;,$$


da cui il sistema dinamico discreto:

$$\left[\begin{array}{c}{x_k}\\ {x_{k+1}}\end{array}\right]... ...y}{cc}{0}&{1}\\ {-1}&{2-\omega^2\,h^2}\end{array}\right]\;,$$


con equazione caratteristica:

$$det[A-\lambda\,I]= \lambda^2 -(2-\omega^2\,h^2)\,\lambda +1 =0$$


e autovalori

$$\lambda= 1-\frac{\omega^2\,h^2}2 \pm \frac 12\,\sqrt{\omega^4\,h^4-4\,\omega^2\,h^2}$$


(reali o complessi secondo il segno del discriminante $\omega^4\,h^4-4\,\omega^2\,h^2$). Quindi gli autovalori sono complessi coniugati e di modulo 1 per $\vert\omega\,h\vert<2$: in tal caso la soluzione del sistema dinamico discreto si può scrivere come

$$\left[\begin{array}{c}{x_k}\\ {x_{k+1}}\end{array}\right]... ... \left[\begin{array}{c}{x_0}\\ {x_1}\end{array}\right] \;,$$



dove la matrice $C$ ha per colonne la parte reale e la parte complessa di un autovettore di autovalore $\lambda$, ed $R_{k\theta}$ è una matrice di rotazione, con l'angolo di rotazione ad ogni passo dato da


$$\theta= arg(\lambda)= arcsin\left(\omega\,h\,\sqrt{1-\frac{\omega^2\,h^2}4}\right) =\omega\,h + O(\omega^3\,h^3)\;.$$



L'esempio qua sopra mostra che, nel caso di un sistema dinamico continuo dotato di due proprietà importanti come quella di essere conservativo e di avere un integrale primo, è talvolta possibile costruire una approssimazione discreta che mantiene entrambe le proprietà: infatti le rotazioni $R_{k\theta}$ che descrivono la soluzione del sistema discreto conservano l'area, e lasciano invariata la funzione energia. Queste proprietà però sono valide solo in un certo sistema di coordinate, definito dal cambiamento di coordinate $C^{-1}$. Inoltre il problema è capire se queste due proprietà possono essere conservate anche in casi non banali, nei quali la soluzione del sistema dinamico continuo non è nota esplicitamente.

Mappa standard: caso conservativo

Vogliamo generalizzare l'esempio precedente al caso di un qualsiasi sistema newtoniano ad un grado di libertà:

$$\frac{d^2{x}}{d{t}^2} =f(x) \Longleftrightarrow \left\{\beg... ... {\displaystyle=} &{\displaystyle f(x)} \end{array}\right.\;.$$

Se usiamo le differenze centrali otteniamo la discretizzazione

$$\Delta_0^2\, x_k= h^2\,f(x_k)\; .$$

Però, per essere in grado di usare le condizioni iniziali $x(0)=x_0, y(0)=y_0$ occorre trovare una discretizzazione anche per l'equazione $\dot x=y$. Se usiamo la differenza all'indietro

$$h\,y_k=\Delta_- x_k \Longleftrightarrow h\,y_{k+1}=\Delta_+ x_k\;,$$

si ottiene un sistema dinamico discreto decomponendo $\Delta_0^2=\Delta_+ -\Delta_-$

$$x_{k+1}-x_k=x_k-x_{k-1} +h^2\,f(x_k)\;,$$

quindi sostituendo $y_{k+1}, y_k$ si ha $h\,y_{k+1}=h\,y_k +h^2\,f(x_k)$; dividendo per $h$ questa equazione, e aggiungendoci $x_{k+1}=x_k+h\,y_{k+1}$ si trova il sistema dinamico discreto nonlineare

$$\left\{\begin{array}{lcl} {\displaystyle x_{k+1}} & {\disp... ...style=} &{\displaystyle y_k +h\,f(x_k)} \end{array}\right.\;,$$

che si chiama mappa standard del sistema newtoniano $\ddot x= f(x)$.

La mappa standard associata a $\ddot x= f(x)$, con $f$ di classe $C^1$ su un aperto $W\subset {\bf R}$, è conservativa su $W\times {\bf R}$. Dimostrazione:

Possiamo notare che, nelle equazioni che definiscono per ricorrenza la mappa standard, il calcolo della nuova $y_{k+1}$ deve essere eseguito prima del calcolo della nuova $x_{k+1}$. Perciò è conveniente spezzare la mappa che definisce il sistema discreto nella composizione di due mappe:

$$\left[\begin{array}{c}{x_k}\\ {y_k}\end{array}\right]\long... ...t[\begin{array}{c}{x_{k+1}}\\ {y_{k+1}}\end{array}\right]\;,$$


ciascuna delle quali è uno scorrimento, di una delle due forme

$$\left[\begin{array}{c}{x}\\ {y}\end{array}\right]\longmapsto \left[\begin{array}{c}{x}\\ {y+a(x)}\end{array}\right]$$


oppure

$$\left[\begin{array}{c}{x}\\ {y}\end{array}\right]\longmapsto \left[\begin{array}{c}{x+b(y)}\\ {y}\end{array}\right]\;,$$


con $a(x), b(y)$ funzioni $C^1$.

Ogni scorrimento è un sistema dinamico discreto conservativo: infatti una mappa conserva l'area se e solo se il determinante jacobiano della trasformazione vale ovunque $1$ (segue dal teorema del cambiamento di variabili negli integrali doppi). Nei due casi la matrice jacobiana è


$$\frac{\partial {(x,y+a(x))}}{\partial {(x,y)}}=\left[\begin{... ...[\begin{array}{cc}{1}&{b'(y)}\\ {0}&{1}\end{array}\right]\;,$$


con determinante comunque uguale ad $1$. Poiché la composizione di due trasformazioni che conservano l'area conserva l'area, la mappa standard è conservativa.

 C.D.D.

Figura: Mappa standard del pendolo con passo h=0.5; si noti che le soluzioni del sistema discreto non seguono le curve ad energia costante, ma sembrano seguire delle curve regolari distorte. In effetti alcune di queste curve sono biforcate in più ``isole'', come nel caso della condizione iniziale $(3,-2)$; inoltre esiste una traiettoria che sembra seguire entrambi i rami della separatrice. Si veda anche la Figura 6.14.

In conclusione la mappa standard costituisce una discretizzazione conservativa, che trasforma un sistema dinamico continuo conservativo in un sistema discreto pure conservativo. Però la risposta al quesito se la mappa standard abbia un integrale primo (analogo all'integrale dell'energia del sistema continuo) è negativa. Per mostrare questo è sufficiente un controesempio:

Esempio:

La mappa standard del pendolo si ottiene a partire dal sistema newtoniano $\ddot x= -\sin(x)$, cioè dal pendolo nonlineare conservativo:

$$\left\{\begin{array}{lcl} {\displaystyle x_{k+1}} & {\dis... ...style=} &{\displaystyle y_k-h\,\sin(x_k)} \end{array}\right.$$



Ricordarsi di calcolare prima la seconda equazione e poi la prima.

Figura 4.5: Mappa standard del pendolo con passo h=0.5; un ingrandimento della figura precedente, che mostra una sola orbita del sistema dinamico discreto, con condizione iniziale (3.141,0) molto vicina al punto di sella. È evidente da questa figura che non può esserci un integrale primo.



Figura 4.6: Mappa standard del pendolo con passo h=1; la regione caotica occupa gran parte del toro, una regione apparentemente regolare corrisponde alle oscillazioni a piccola ampiezza, e appaiono molte isole di risonanza.


La mappa è chiaramente periodica di periodo $2\pi$ nella variabile $x$, che può essere interpretata come una variabile angolo. Per semplificare la descrizione delle proprietà rispetto alla variabile $y$ conviene cambiare variabile, rimpiazzando $h\,y$ con $y$:

$$\left\{\begin{array}{lcl} {\displaystyle x_{k+1}} & {\dis... ...yle=} &{\displaystyle y_k-h^2\,\sin(x_k)} \end{array}\right.$$


in modo che la mappa sia periodica di periodo $2\pi$ anche rispetto ad $y$ (se $y$ aumenta di $2\pi$, l'incremento di $x$ aumenta di $2\pi$, il che è indifferente perché $x$ è una variabile angolo). In conclusione la mappa può essere interpretata come un diffeomorfismo del toro ${\bf T}^2$ su se stesso.

La mappa standard del pendolo dipende quindi da un solo parametro $\mu=h^2\geq 0$. Intuitivamente è chiaro che per $\mu$ molto piccolo il sistema discreto approssimerà il sistema continuo, e quindi la funzione energia $E(x,y)=y^2/2 - \cos(x)$ cambierà di poco ad ogni passo, e la successione $(x_k,y_k)$ si muoverà quasi lungo le curve di livello di $E$, cioè vicino alla soluzione del sistema continuo. Invece per $\mu$ grande il sistema discreto non avrà niente a che fare con il sistema continuo: per $\mu>4$ nell'origine diventa una mappa instabile, come risulta dall'analisi della discretizzazione del pendolo lineare.

Per osservare che cosa succede per valori intermedi di $\mu$, osserviamo le orbite della mappa standard per i valori $\mu=1/4$ (Figura 4.4) e $\mu=1$ (Figura 4.6).

Per valori relativamente piccoli di $\mu$, come nella Figura 4.4, sembra quasi che il sistema dinamico discreto ammetta un integrale primo, anche se con le curve di livello distorte rispetto a quelle della funzione energia del corrispondente sistema dinamico continuo. Ma se si guarda più in dettaglio, per esempio ingrandendo la figura, si nota che il comportamento qualitativo delle orbite del sistema dinamico discreto è molto diverso, ed in particolare alcune curve invarianti scompaiono: questi fenomeni di caos si verificano in modo più appariscente vicino alle separatrici che escono dal punto di sella (Figura 4.5), ma sono presenti anche vicino alle separatrici delle isole di risonanza.

Per valori più grandi di $\mu$, come in Figura 4.6, l'illusione che la mappa standard possa avere un integrale primo scompare del tutto. Per questo valore ($\mu=1$) l'ultima curva invariante corrispondente ad una circolazione è scomparsa, e le orbite caotiche possono muoversi liberamente lungo tutto l'asse $y$. Si veda la discussione nel Capitolo 6.

Problema Si consideri il problema dei due corpi, cioè il sistema newtoniano di dimensione 1

$$\frac{d^2{r}}{d{t}^2} = \frac {c^2}{r^3} -\frac{GM}{r^2}\;,$$

dove $GM>0$ e $c$ sono costanti. Si costruisca la mappa standard del problema dei due corpi, se ne trovino i punti fissi e se ne studi il comportamento al variare del passo $h$.

Suggerimento: Nel caso continuo, tutte le orbite limitate sono periodiche, con periodi che dipendono solo dall'integrale dell'energia; ponendo $p=\dot r$

$$E(r,p)=\frac 12 \,\left(p^2+\frac{c^2}{r^2}\right) -\frac{GM}{r}\;,$$

le orbite sono periodiche per $E<0$, ed il periodo è tanto più grande quanto l'energia si avvicina a 0 (per la terza legge di Keplero). Per $E>0$ le orbite vanno all'infinito. Quindi $E=0$ è una separatrice. Nel problema discretizzato la separatrice si trasforma in una regione caotica, tanto più ampia quanto più lungo è il passo $h$.

(Soluzione)

Mappa standard: caso dissipativo

L'impiego di una discretizzazione conservativa può risultare vantaggioso anche nel caso in cui il sistema dinamico continuo è dissipativo. Per esempio, applichiamo la mappa standard al sistema newtoniano dissipativo

$$\ddot x= f(x)-\gamma\, \dot x$$

usando le approssimazioni

$$D^2\simeq \frac{\Delta_0^2}{h^2} \hspace{5mm},\hspace{5mm}D \simeq \frac{\Delta_-}h$$

e otteniamo l'equazione alle differenze finite

$$\Delta_0^2 x_k= h^2 f(x_k)-\gamma h \, \Delta_- x_k \ .$$

Per ottenere un sistema dinamico usiamo, come nel caso della mappa standard del pendolo, $y_k=\Delta_-x_k$ e otteniamo il sistema dinamico discreto:

$$\left\{\begin{array}{lcl} {\displaystyle x_{k+1}} & {\displ... ...tyle h^2\, f(x_k) +(1-\gamma h)\, y_k} \end{array}\right. \ .$$

I punti fissi del sistema discreto soddisfano alle equazioni $x_{k+1}=x_k$ (da cui si deduce $y_k=0$) e $y_{k+1}=y_k$ (da cui $h^2\, f(x_k)=0$); ossia, i punti fissi sono della forma $(x_S, 0)$ con $f(x_S)=0$, esattamente i punti di equilibrio del sistema continuo. Le linearizzazioni in questi punti fissi hanno matrice

$$A=\left[\begin{array}{cc}{1+h^2\, f'(x_k)}&{1-\gamma h}\\ {h^2\, f'(x_k)}&{1-\gamma h}\end{array}\right]$$

con determinante $1-\gamma h$ e traccia $2+h^2\, f'(x_k)-\gamma h$. Perciò il polinomio caratteristico è

$$P(\lambda)=\lambda^2 -(2+h^2\, f'(x_k)-\gamma h)\, \lambda + 1-\gamma h\ .$$

Per determinare il carattere delle radici conviene cacolarne il valore in 0 e 1:

$$P(0)=1-\gamma h \hspace{5mm},\hspace{5mm}P(1)=-h^2\, f'(x_k) \ .$$

Possiamo assumere che $\gamma$ ed $h$ sono piccoli, per cui $\gamma h <1$, da cui $P(0)>0$ e $P(1)<0$ per $f'(x_k)>0$, cioè nel caso del massimo dell'energia potenziale. Allora le due radici sono una compresa tra 0 e 1 e una maggiore di 1, si ottiene un punto fisso iperbolico corrispondente alla sella del sistema continuo.

Se invece $f'(x_k)<0$, cioè nel caso del minimo dell'energia potenziale, $P(0)$ e $P(1)$ sono entrambe positivi. Se le radici sono complesse hanno modulo $\sqrt{1-\gamma h}<1$. Se invece sono reali non possono essere maggiori di $1$ perché $P'(1)=\gamma\,h -h^2\, f'(x_k)>0$ e, per $h$ abbastanza piccolo, $P'(0)<0$. Quindi le radici, se reali, devono essere entrambe in $[0,1]$, e il punto fisso è asintoticamente stabile e corrisponde al pozzo del caso continuo. In altre parole, per $h$ abbastanza piccolo il sistema discretizzato ha almeno localmente, nell'intorno dei punti fissi, le stesse proprietà qualitative del flusso integrale che sta approssimando.

Si noti che questa corrispondenza qualitativa è valida solo nell'intorno di un punto di equilibrio iperbolico, cioè nel piano nell'intorno di un pozzo, sella o sorgente. Per punti di equilibrio con esponenti di Lyapounov nulli le proprietà di stabilità non sono necessariamente conservate dalla discretizzazione.

Esercizio Dato il sistema newtoniano con dissipazione

$$\ddot x = x^2 -x^3 -\gamma \dot x$$

costruirne la discretizzazione conservativa mediante la la mappa standard. Determinare le proprietà di stabilità dei punti di equilibrio del sistema continuo e dei punti fissi del sistema discreto. Esiste un intervallo $(0,h_0)$ in cui sono le stesse? (Soluzione)