Equazioni differenziali al computer: Eulero, Runge-Kutta e problemi stiff
Metodo di Eulero e suo errore, Heun e Runge-Kutta del quarto ordine, passo adattivo, stabilità e problemi stiff, metodi impliciti in SPICE e metodi simplettici.
Le equazioni differenziali che si risolvono con carta e penna sono poche: quelle lineari a coefficienti costanti, quelle a variabili separabili. Il circuito con un diodo, il pendolo con oscillazioni ampie, un motore con attrito non lineare, l’orbita di un satellite attorno alla Terra schiacciata ai poli richiedono metodi numerici. L’idea è semplice: si parte dalla condizione iniziale e si avanza a piccoli passi nel tempo, seguendo il campo di direzioni. In questa lezione vediamo i metodi di Eulero e di Runge-Kutta, l’ordine di accuratezza e il problema della stabilità.
Il problema
![\[ y' = f(t, y),\qquad y(t_0) = y_0 \]](https://www.electronics-engineering.com/wp-content/ql-cache/quicklatex.com-22423eb556ce6796125e20f5420f0f4f_l3.png "Rendered by QuickLaTeX.com")
dove 
può essere anche un vettore: ogni equazione di ordine superiore si riscrive come sistema del primo ordine. Il pendolo 
diventa 
, 
. Si scelgono i tempi 
, con 
il passo, e si calcolano approssimazioni 
.
Il metodo di Eulero esplicito
Si segue la tangente per un passo:
![\[ y_{n+1} = y_n + h\,f(t_n, y_n) \]](https://www.electronics-engineering.com/wp-content/ql-cache/quicklatex.com-597ab5f9b3049db19a6b07cb48157665_l3.png "Rendered by QuickLaTeX.com")
È la definizione di derivata letta al contrario. A ogni passo si commette un errore locale proporzionale a 
(il termine trascurato dello sviluppo di Taylor); dopo 
passi l’errore globale è proporzionale a . Eulero è un metodo del primo ordine: dimezzando il passo l’errore si dimezza.
, 
, la soluzione esatta è 
(blu). Eulero con passo 
(rosso) e 
(arancione) resta sotto la curva vera, perché la tangente di una funzione convessa sta sotto la curva; dimezzando il passo l’errore si dimezza circa.Esempio. Con Eulero dà 
, cioè 
. Per arrivare a 
: con 
, 
invece di 
(errore 4,6%); con 
, 
(errore 0,5%). È la definizione di 
come limite di 
.
I metodi di Runge-Kutta
Per avere più precisione senza ridurre troppo il passo, si valuta la pendenza in più punti dentro il passo e se ne fa una media pesata. Il metodo di Heun (Runge-Kutta del secondo ordine) fa un passo di Eulero di prova, calcola la pendenza nel punto d’arrivo e usa la media delle due pendenze:
![\[ k_1 = f(t_n, y_n),\quad k_2 = f(t_n + h, y_n + hk_1),\quad y_{n+1} = y_n + \frac h2(k_1 + k_2) \]](https://www.electronics-engineering.com/wp-content/ql-cache/quicklatex.com-3a1be4c5287999c541e801612b90fa8a_l3.png "Rendered by QuickLaTeX.com")
Il metodo più usato è Runge-Kutta del quarto ordine (RK4):
![\[ \begin{aligned} k_1 &= f(t_n, y_n) & k_2 &= f\!\left(t_n + \tfrac h2, y_n + \tfrac h2k_1\right)\\ k_3 &= f\!\left(t_n + \tfrac h2, y_n + \tfrac h2k_2\right) & k_4 &= f(t_n + h, y_n + hk_3) \end{aligned} \]](https://www.electronics-engineering.com/wp-content/ql-cache/quicklatex.com-89e8f309a8a58016b73628c1a11db376_l3.png "Rendered by QuickLaTeX.com")
![\[ y_{n+1} = y_n + \frac h6\big(k_1 + 2k_2 + 2k_3 + k_4\big) \]](https://www.electronics-engineering.com/wp-content/ql-cache/quicklatex.com-0060e70b0bd6ec5cab8a421eeda941c8_l3.png "Rendered by QuickLaTeX.com")
Quattro valutazioni di 
per passo, errore globale proporzionale a 
: dimezzando il passo l’errore si riduce di sedici volte. Per con fino a , RK4 dà 
contro 
: errore di 2 milionesimi, contro il 4,6% di Eulero con lo stesso passo. La struttura dei pesi 
è la stessa della regola di Simpson: per che dipende solo da 
, RK4 è Simpson.
Passo adattivo. I programmi moderni (ode45 di MATLAB, solve_ivp di SciPy) calcolano a ogni passo due stime di ordine diverso, ne usano la differenza come stima dell’errore e regolano il passo: grande dove la soluzione è tranquilla, piccolo dove cambia in fretta.
La stabilità e i problemi stiff
Si applichi Eulero all’equazione 
con 
, la cui soluzione 
decresce. Eulero dà 
. Se 
, cioè se
![\[ h \gt \frac{2}{\lambda} \]](https://www.electronics-engineering.com/wp-content/ql-cache/quicklatex.com-003289e4c17ee3af9dbfae5b3ca01a48_l3.png "Rendered by QuickLaTeX.com")
i valori calcolati crescono e oscillano con segno alterno, mentre la soluzione vera va a zero. Il metodo è instabile. Non è una questione di precisione: anche se la soluzione dopo un po’ è quasi zero e non serve seguirla con dettaglio, il passo deve restare piccolo.
Un problema è stiff (rigido) quando contiene costanti di tempo molto diverse. Un circuito con una costante di tempo di 1 ns e un’altra di 1 ms, simulato per qualche millisecondo, costringerebbe un metodo esplicito a passi sotto i 2 ns per tutta la durata: milioni di passi, anche dopo che il transitorio veloce si è esaurito.
La soluzione sono i metodi impliciti, come Eulero implicito:
![\[ y_{n+1} = y_n + h\,f(t_{n+1}, y_{n+1}) \]](https://www.electronics-engineering.com/wp-content/ql-cache/quicklatex.com-572a6f79ba24d75c9f126e2cacb89bdb_l3.png "Rendered by QuickLaTeX.com")
dove la pendenza si valuta nel punto d’arrivo. Per dà 
, che decresce per qualsiasi passo. Il prezzo è che a ogni passo si deve risolvere un’equazione (o un sistema non lineare, con Newton). SPICE usa proprio metodi impliciti, la regola dei trapezi implicita e le formule di Gear, perché i circuiti elettronici sono quasi sempre stiff. Quando nelle opzioni di simulazione si cambia il metodo di integrazione da “trap” a “gear”, si sta scegliendo tra due metodi impliciti con proprietà di stabilità diverse.
La conservazione dell’energia
Per i sistemi oscillanti senza attrito, come un circuito LC ideale o un’orbita, Eulero esplicito fa crescere l’energia a ogni giro e la traiettoria si allarga a spirale; Eulero implicito la fa diminuire. Esistono metodi, detti simplettici (come il metodo di Verlet, usato nelle simulazioni di dinamica molecolare e nei motori fisici dei videogiochi), che conservano l’energia per tempi lunghissimi anche con passi grandi.
