Risoluzione numerica di un sistema di equazioni differenziali ordinarie

Lo scopo di questo progetto è la risoluzione numerica di un sistema lineare di ODE (ordinary differential equations) caratterizzato dall'avere un alto coefficiente di stiffness ed essere pertanto soggetto a generare una forte instabilità numerica se approcciato con un metodo numerico non consono. Il problema è stato risolto con vari algoritmi numerici e con vari tipi di mesh. Come mostrano i risultati, la A-stabilità ed, ancor più, la L-stabilità sono due proprietà fortemente desiderabili in un algoritmo risolutore, infatti, in loro assenza si rende necessario utilizzare un numero esorbitante di nodi per approssimare numericamente la soluzione. I metodi numerici che godono solo della zero-stabilità si sono rivelati assolutamente non in grado di approssimare efficacemente la soluzione, coerentemente con le previsioni teoriche.
I codici sono stati implementati in linguaggio MATLAB, tuttavia, essendo tutti gli algoritmi stati scritti from scratch (da zero) utilizzando unicamente la libreria standard di MATLAB, i codici risultano pienamente compatibili sia con le licenze base di MATLAB (senza la necessità di specifici toolbox a pagamento), sia con il software open-source GNU Octave.
Inoltre, una relazione molto più dettagliata di questo file README è stata inserita nel repository in formato PDF.
Sia i codici che la relazione sono stati sviluppati come progetto individuale.

Il problema

Si consideri il problema di Cauchy

$$ y \hspace{0.1cm}'(t)=M \cdot y(t), \qquad y(t_0)=y_0, $$

dove M è definita come la matrice di dimensione $15 \times 15$, avente come $j$-esimo elemento della diagonale principale il valore $-j^2$ e come elementi della sovradiagonale i valori $150$:

$$ M := \begin{pmatrix} -1 & 150 & 0 & \cdots & 0 \\ 0 & -4 & 150 & \cdots & 0 \\ 0 & 0 & -9 & & \vdots \\ \vdots & \vdots & & \ddots & 150 \\ 0 & 0 & \cdots & 0 &-225 \end{pmatrix}. $$

e dove

$$ \underline{y_0} = [ 1, 1, 1, ..., 1, 10 ]^T $$

da risolvere per $t \in [ t_0 , T ]$, con $T$ fissato e scelto nell'intervallo $[ 10 , 10000 ]$.
Per calcolare numericamente la soluzione, l'intervallo temporale $[t_0,T]$ viene suddiviso in una successione di nodi:

$$ t_0 < t_1 < \cdots < t_N=T , $$

chiamata mesh e su ogni nodo della mesh viene approssimato il valore della funzione (vettoriale).
La distanza $h_n=t_{n+1}-t_n$ tra due nodi consecutivi viene detta passo di discretizzazione, se il passo di discretizzazione è lo stesso per tutti i nodi della mesh, la mesh viene detta uniforme oppure omogenea.

Analisi teorica a priori

$M$ ha $15$ autovalori distinti ed ha spettro: $\sigma (M) = \{ −1,−4,…,−225 \}$, in particolare, tutti e $15$ i suoi autovalori sono reali negativi e pertanto il sistema di ODE in esame descrive una dinamica asintoticamente stabile per il teorema di Lyapunov.
Tuttavia, la scrittura della soluzione come combinazione lineare di termini della forma $e^{\lambda_i \cdot (t- t_0) } \cdot v_i$ mostra l'esistenza di transitori veloci e transitori lenti, in altre parole, il problema forza l'utilizzo di tantissimi nodi negli istanti iniziali, mentre negli istanti finali sono necessari molti meno nodi per approssimare efficacemente la soluzione. Conseguenza diretta di ciò è che utilizzare una mesh omogenea potrebbe non essere la scelta più vantaggiosa.
Come noto dall'analisi, i sistemi lineari di equazioni differenziali ordinarie ammettono come soluzione esatta $y(t) = \exp( M \cdot (t - t_0)) \cdot y_0$, dove:

$$ \exp(A) := \sum_{k=0}^{+ \infty} \frac{A^k}{k!} $$

La funzione $\exp(A)$ è detta esponenziale di matrice della matrice $A$, ed in MATLAB è calcolabile mediante la funzione di libreria standard expm.

Il coefficiente di stiffness del sistema vale:

$$ C_{stiff} = \frac{ \max\limits_{i=1,\ldots,15} \left| Re( \lambda_i ) \right| }{ \min\limits_{i=1,\ldots,15} \left| Re( \lambda_i ) \right| } = 225 \gg 1 $$

Dal momento che il coefficiente di stiffness è così elevato si può catalogare questo sistema di ODE come problema stiff.

La soluzione esatta è la seguente:

Possiamo osservare che essa presenta un picco iniziale nei primissimi valori di 𝑡, dove raggiunge anche valori di $10^7$ nella prima componente e $10^5$ nella seconda componente, tuttavia dopo il picco iniziale essa decade molto velocemente, anche solo per $t=30$ il vettore soluzione ha norma euclidea dell’ordine di $10^{-5}$, mentre per $t=10,000$ tutte le componenti della soluzione esatta hanno valore indistinguibile da zero in aritmetica finita.
Dal momento che la funzione soluzione presenta valori così elevati, risulta significativo anche lo studio dell'errore relativo, oltre a quello dell'errore assoluto.

Il metodo di Radau IIA

Il metodo di Radau_IIA è un metodo numerico di Runge-Kutta basato sul metodo di collocazione Radau IIA a tre stadi, si tratta di un metodo di ordine $5$ ed L-stabile. Il presente metodo è stato implementato su mesh omogenea.
La L-stabilità permette di smorzare efficacemente le componenti associate agli autovalori grandi in modulo. Il prezzo da pagare è che il metodo è anche costoso da un punto di vista computazionale: per ogni nodo della mesh il metodo richiede di risolvere un sistema $45 \times 45$, tuttavia, in questo particolare caso dove il sistema di ODE è lineare e la mesh è omogenea, non solo i sistemi da risolvere numericamente sono lineari, ma la matrice $A$ che descrive il sistema $Ax_n = b_n$ è sempre la stessa per tutto il metodo, cambia solo il termine noto $b_n$, ciò consente di utilizzare la fattorizzazione $LU$ una volta per tutte all’inizio ed abbattere così il costo computazionale richiesto per risolvere i sistemi ad un costo quadratico di operazioni macchina (flops).

Seguono alcuni grafici per rappresentare la soluzione approssimata. I tre grafici seguenti mostrano la prima componente $y_1(t)$ del vettore soluzione $\underline{y}(t)$ e sono stati tutti e tre prodotti utilizzando come parametri $t_0 = 0$ e $T = 50$ , inoltre nella prima immagine $h = 1$, nella seconda $h = 0.1$, nella terza $h = 0.01$. È decisamente soddisfacente osservare come con anche solo $h = 0.1$ e dunque $N = 501$ nodi (un numero quasi sempre irrisorio quando si parla di cardinalità di una mesh) il metodo produca un errore relativo di $10^{-4}$, risultando quindi decisamente affidabile.

t0 = 0,   T = 50,   h = 1

t0 = 0,   T = 50,   h = 0.1

t0 = 0,   T = 50,   h = 0.01

Seguono inoltre le tabelle degli errori commessi:

T = 10 Numero di nodi Errore di discretizzazione globale Errore relativo 100 4.35928e+03 1.97120e-04 1,000 8.13346e-02 3.67199e-09 10,000 1.27448e-06 5.75384e-14 100,000 2.64496e-07 1.19411e-14

T = 10,000 Numero di nodi Errore di discretizzazione globale Errore relativo 100,000 4.35928e+03 1.97120e-04 1,000,000 8.13346e-02 3.67199e-09

Il metodo di Gauss-Legendre su mesh non omogenea

Come descritto sopra, i problemi stiff forzano ad utilizzare un passo molto più piccolo del necessario per approssimare la soluzione, tuttavia, un passo così piccolo è necessario per descrivere efficacemente la soluzione unicamente in un intervallo molto ristretto di valori, altrove si può tranquillamente utilizzare un passo molto più ampio. È qui che le mesh omogenee mostrano i loro limiti e si è di conseguenza deciso di implementare anche dei metodi su una mesh non omogenea. L'idea implementativa è stata quella di discretizzare con un passo $h_1$ l’intervallo $[t_0 , 20]$, con un passo $h_2$ l’intervallo $[20,100]$, e con un passo $h_3$ l’intervallo $[100,𝑇]$, con l’accortezza di scegliere $h_1$ estremamente piccolo ed $h_3$ estremamente grande.
L'idea è in questo caso vincente. Il prezzo da pagare è che essendo la mesh non omogenea non è più possibile utilizzare la fattorizzazione $LU$ per abbattere il costo computazionale, tuttavia, il costo computazionale da pagare per rinunciare alla fattorizzazione $LU$ è abbondantemente compensato dall’abbattere drasticamente il numero di sistemi lineari da risolvere.
Il metodo di Gauss-Legendre è anch'esso un metodo di Runge-Kutta basato su un metodo di collocazione, solo che ora viene utilizzata un'interpolazione di Gauss-Legendre. Il metodo implementato è un metodo di Gauss-Legendre a 3 stadi, e dunque di ordine 6. A differenza del metodo di Radau IIA, il metodo di Gauss-Legendre è solo A-stabile, e non L-stabile.

Parametri utilizzati t0 = 0,   T = 10,000 h1 = 0.1,   h2 = 1,   h3 = 10 Un passo h1 troppo grande introduce delle oscillazioni nella componente u15 della soluzione approssimata, a causa della mancanza di L-stabilità. La A-stabilità del metodo provvede comunque a smorzare progressivamente tali oscillazioni.

Le oscillazioni compaiono solo sulla $15$-esima componente poiché è la dinamica con autovalore più alto in modulo, per osservare oscillazioni anche sulle altre componenti bisognerebbe prendere un passo $h_1$ ancora più grande, che però introdurrebbe un errore assoluto eccessivo, rendendo la soluzione approssimata trovata priva di qualunque significato. Non è un caso che oscillazioni sulle altre componenti siano difficili da osservare: la A-stabilità è una proprietà matematica decisamente forte e che comunque rende il metodo incondizionatamente assolutamente stabile.

L’immagine seguente mette in luce come appare la soluzione su una mesh non omogenea:

Parametri utilizzati t0 = 0,   T = 10,000 h1 = 0.01,   h2 = 0.1,   h3 = 10 La mesh è più densa negli istanti iniziali, dove la soluzione presenta i transitori più veloci, e diventa progressivamente più sparsa negli istanti finali.

Forniamo infine una tabella con gli errori fissati i parametri $t_0 = 0$ ed $T=10,000$.

Num. nodi	Err. Ass.	Err. Rel.	Scelta dei passi: h1 in [t0, 20], h2 in [20, 100], h3 in [100, T]
291	6.69894e+02	3.02915e-05	h1 = 0.1,   h2 = 1,   h3 = 1000
2,900	1.66442e-02	7.51429e-10	h1 = 0.01,   h2 = 0.1,   h3 = 100
20,180	5.34687e-07	2.41393e-14	h1 = 0.001,   h2 = 1,   h3 = 100
21,791	1.08034e-07	4.87737e-15	h1 = 0.001,   h2 = 0.1,   h3 = 10
100,000	6.69894e+02	3.02915e-05	h = 0.1, mesh omogenea
1,000,000	1.66442e-02	7.51429e-10	h = 0.01, mesh omogenea

È interessante osservare come, lavorare con $2,900$ nodi in mesh non omogenea, oppure lavorare con $1,000,000$ nodi in mesh omogenea, produca in realtà lo stesso errore. Il vantaggio di utilizzare una mesh non omogenea è qui inequivocabile.

Il metodo di Gauss-Legendre adattivo

Se non si vuole utilizzare una mesh omogenea, si può optare anche per una mesh adattiva.
I metodi adattivi sono metodi numerici in cui la mesh non è fissata a priori ma i suoi nodi sono calcolati in fase di esecuzione del programma, sulla base di una stima dell’errore. In ogni nodo si determina automaticamente quanto vale il passo di discretizzazione successivo. In questo particolare problema i metodi adattivi risultano particolarmente efficaci.
Il metodo di Gauss-Legendre 2(6) è un metodo di Runge Kutta implicito adattivo embedded, ciò significa che esistono due metodi di Gauss-Legendre, uno di ordine 2 ed uno di ordine 6, che condividono lo stesso Tableau di Butcher, eccezion fatta per il vettore $b$.
Il vantaggio è considerevole: risolvendo gli stessi sistemi lineari e cambiando solo i pesi $b$ delle combinazioni lineari (il costo computazionale di ciò è assolutamente irrisorio) abbiamo a disposizione due soluzioni una A-stabile di ordine 6, ed un’altra di ordine 2. Trattando la soluzione ottenuta col metodo di ordine superiore come se fosse la soluzione esatta e calcolando lo scarto tra le due soluzioni si ha a disposizione un valore che può essere considerato uno stimatore dell’errore assoluto.
Se lo scarto è inferiore ad una tolleranza fissata si accetta il valore trovato inserendolo nella mesh e poi si procede a calcolare il nuovo passo $h_n$, in caso alternativo, si ripete il calcolo riducendo il passo. Infine, si tiene memorizzata come soluzione numerica la soluzione fornita dal metodo di ordine 6.

Segue una tabella contenente i vari errori al variare della tolleranza, ($t_0 = 0$ , $T = 10,000$):

Tolleranza	Numero nodi	Err. Ass.	Err. Rel.	hmin	hmax
1e+5	38	1.32016e+02	5.97904e-06	0.001	7515.21744
1e+4	70	1.62036e+01	7.31537e-07	0.001	7552.67047
1e+3	138	3.22138e-01	1.45438e-08	0.001	7452.33204
1e+2	281	2.18029e-02	9.84334e-10	0.001	5030.59686
1e+1	589	7.42305e-04	3.35125e-11	0.001	8677.20076
1e+0	1,254	7.42405e-05	3.35171e-12	0.00096	8649.32396
1e-1	2,687	6.57684e-06	2.96922e-13	0.00045	7031.06471
1e-2	5,774	7.13678e-07	3.22201e-14	0.00021	4776.71363
1e-3	12,423	1.86265e-07	8.40921e-15	0.00010	7578.50364
1e-4	26,748	1.37839e-07	6.22294e-15	0.00004	7133.44752

È inoltre particolarmente interessante osservare come anche fissando una tolleranza enorme di $10^4$, il metodo numerico produca, con solo $70$ nodi, un errore assoluto dell’ordine delle decine per approssimare variabili che viaggiano su grandezze nell’ordine delle decine di milioni.
Il costo computazionale è assolutamente irrisorio, richiedendo al più qualche secondo per tolleranze inferiori a $10^{−4}$, e diventando pressoché istantaneo per tolleranze superiori a $10^{−2}$.

Segue una progressione di immagini che descrive con diversi livelli di zoom l’approssimazione di $y_1(t)$, con $t_0 = 0$ , $T = 10,000$, $toll = 10^4$, $Num nodi = 70$, le quattro immagini seguenti sono viste rispettivamente con zoom negli intervalli: $[0, 20]$, $[0, 50]$, $[0, 220]$ ed $[0, 10,000]$.

t ∈ [0, 20]	t ∈ [0, 50]
t ∈ [0, 220]	t ∈ [0, 10,000]

Altri metodi implementati

Il repository comprende inoltre i seguenti metodi numerici:

• BDF 2 e BDF 3 a passo costante;
• Crank-Nicolson su mesh non omogenea;
• Eulero implicito;
• Pareschi-Russo;
• Runge-Kutta classico di ordine 4;
• Runge-Kutta-Fehlberg 4(5) adattivo.

La trattazione, l'analisi ed i risultati relativi a questi metodi vengono discussi nella relazione completa.

Esecuzione

1. Aprire la cartella del repository in MATLAB o GNU Octave.

2. Aprire main.m e impostare la variabile metodo, per esempio:

   metodo = 'Radau_IIA_5';
   % metodo = 'Gauss_Legendre_6_mesh_non_omogenea';
   % metodo = 'Gauss_Legendre_2_6';

3. Modificare, se lo si desidera, i parametri t_0, T, il passo h uniforme, i tre passi della mesh non omogenea passo_1, passo_2 e passo_3, oppure la tolleranza toll del metodo adattivo.

4. Eseguire main.m.

Lo script genererà i dati iniziali, calcolerà la soluzione numerica e quella esatta, stamperà gli errori e produrrà i grafici delle componenti $y_1$, $y_2$ e $y_{15}$.

Limiti dell'implementazione

I solver numerici sono stati implementati ed ottimizzati per questo specifico esperimento: assumono un sistema autonomo lineare $y \hspace{0.1cm}'=My$ di dimensione $15$ e non costituiscono una libreria ODE per un generale sistema di equazioni differenziali.
Questa scelta permette di mettere in evidenza le proprietà numeriche dei metodi e di sfruttare direttamente la struttura lineare del problema per ottimizzare il costo computazionale richiesto dai metodi numerici.