30 Metodi di Runge-Kutta
Il metodo di Eulero1, sebbene utile in molti casi, ha limitazioni di accuratezza ed efficienza, soprattutto per problemi di grandi dimensioni, ad alta dimensionalità o con equazioni differenziali rigide. Per ottenere risultati più precisi e gestire tolleranze d’errore più stringenti, ci affidiamo ai metodi di Runge-Kutta2.
30.1 Formulazione Generale
Consideriamo un problema di Cauchy3 generico per un’equazione differenziale ordinaria (ODE):
\left\{ \begin{aligned} y' &= f(t,y), \\ y(a) &= y_a. \end{aligned} \right.
Un metodo di Runge-Kutta a un passo per approssimare la soluzione di questo problema può essere scritto come:
y_{k+1} = \Phi(t_k, y_k, h) = y_k + A h f(t_k, y_k) + B h f(t_k + C h, y_k + D h f(t_k, y_k)).
dove h è il passo di discretizzazione e i coefficienti A, B, C, D devono essere scelti per garantire l’accuratezza del metodo.
Richiamo: sviluppo di Taylor in due variabili
Tutto quello che segue poggia sullo sviluppo di f attorno a (t,y) in entrambi gli argomenti. Conviene enunciarlo con il resto quantitativo, perché è da lì che escono sia le condizioni d’ordine sia le costanti delle stime d’errore.
Proposizione 30.1 (Sviluppo di Taylor con resto) Sia f\in C^{1} in un intorno di (t,y). Allora
f(t+\alpha, y+\beta) = f(t,y) + \alpha\, f_t(t,y) + \beta\, f_y(t,y) + r_1(\alpha,\beta)
e, se inoltre il gradiente (f_t,f_y) è lipschitziano con costante M_1 in quell’intorno,
\boxed{\; \big| r_1(\alpha,\beta) \big| \leq \dfrac{M_1}{2}\big(\alpha^2+\beta^2\big) \;} \tag{30.1}
Se f\in C^{2} vale lo sviluppo al second’ordine
\begin{aligned} f(t+\alpha,y+\beta) &= f + \alpha f_t + \beta f_y \\ &\quad + \dfrac{1}{2}\Big( \alpha^{2} f_{tt} + 2\alpha\beta f_{ty} + \beta^{2} f_{yy} \Big) + r_2(\alpha,\beta) \end{aligned}
con |r_2| = \mathcal{O}\big((\alpha^2+\beta^2)^{3/2}\big), e |r_2|\leq \tfrac{M_2}{6}(\alpha^2+\beta^2)^{3/2} se l’hessiana è lipschitziana con costante M_2.
Nei metodi di Runge-Kutta si applica sempre con \alpha = c_i h e \beta = \mathcal{O}(h): ogni incremento è proporzionale ad h, quindi il resto 30.1 è \mathcal{O}(h^{2}) e, al second’ordine, \mathcal{O}(h^{3}). È questo che permette di troncare gli sviluppi all’ordine voluto e di raccogliere i termini in potenze di h.
Calcolo dell’Errore Locale di Troncamento
Per valutare l’accuratezza del metodo, calcoliamo l’errore locale di troncamento \tau(h), che rappresenta la differenza tra la soluzione esatta y(t+h) e la stima data dal metodo \Phi(t, y(t), h):
\begin{aligned} \tau(h) &= \color{blue}y(t+h)-\color{DarkGreen}\Phi(t,y(t),h)\\ & = {\color{blue}y(t+h)}- \color{DarkGreen} y(t) - A h f(t,y(t)) - B h f(t+Ch,y(t)+Dh f(t,y(t))) \end{aligned}
dove y(t) è la soluzione esatta per la quale possiamo scrivere la serie di Taylor4
y(t+h) = y(t) + y^\prime(t) h + y^{\prime\prime}(t)\dfrac{h^2}{2}+ \cdots + y^{(p)}(t)\dfrac{h^p}{p!}+\mathcal{O}(h^{p+1})
Dalla ODE possiamo scrivere i primi due termini della serie come segue
y(t+h) = y(t) + \underbrace{f(t,y(t))}_{y^\prime(t)} h + \underbrace{\left( \dfrac{\partial f(t,y(t))}{\partial t}+ \dfrac{\partial f(t,y(t))}{\partial y}f(t,y(t)) \right)}_{y^{\prime\prime}(t)}\dfrac{h^2}{2}+ \mathcal{O}(h^3)
Sviluppo dell’Errore Locale di Troncamento
Sostituiamo questa espansione in \tau(h):
\begin{aligned} \tau(h) & = {\color{blue} \cancel{y(t)} + f(t,y(t)) h+ \left(\dfrac{\partial f(t,y(t))}{\partial t}+ \dfrac{\partial f(t,y(t))}{\partial y}f(t,y(t)) \right) \dfrac{h^2}{2}+\mathcal{O}(h^3)} \\ & - \color{DarkGreen} \cancel{y(t)} - A h f(t,y(t))\\ & \color{DarkGreen} - B h \left( f(t,y(t))+ \dfrac{\partial f(t,y(t))}{\partial t}Ch+ \dfrac{\partial f(t,y(t))}{\partial y}Dh f(t,y(t))+ \mathcal{O}(h^2) \right) \\ & = h\cdot{\color{DarkRed}f(t,y(t))}(1-A-B) + \dfrac{h^2}{2}\cdot{\color{DarkRed}\dfrac{\partial f(t,y(t))}{\partial t}} \left(1-2BC\right)\\ &+ \dfrac{h^2}{2}\cdot{\color{DarkRed}\dfrac{\partial f(t,y(t))}{\partial y}f(t,y(t))} \left(1-2BD\right)+ \mathcal{O}(h^3) \end{aligned}
Per ottenere un metodo di secondo ordine, imponiamo le condizioni seguenti:
\left\{ \begin{aligned} 1-A-B & = 0 \\ 1-2BC & = 0 \\ 1-2BD & = 0 \end{aligned} \right.
Queste condizioni assicurano che l’errore di troncamento \tau(h) = \mathcal{O}(h^3), quindi il metodo globale sarà del secondo ordine.
Famiglia RK del secondo ordine
Risolvendo il sistema di equazioni sopra, possiamo esprimere i parametri A, C, e D in funzione di B come segue:
A = 1 - B, \qquad C = D = \dfrac{1}{2B}.
Pertanto, i metodi di Runge-Kutta del secondo ordine assumono la forma generale:
y_{k+1} = y_k + (1 - B) h f(t_k, y_k) + B h f\left(t_k + \dfrac{h}{2B}, y_k + \dfrac{h}{2B} f(t_k, y_k)\right).
Definiamo ora:
b_1 = 1 - B, \quad b_2 = B, \quad c_1 = 0, \quad c_2 = \dfrac{1}{2B}, \quad A_{21} = \dfrac{1}{2B}.
Possiamo riscrivere il metodo come segue:
\begin{aligned} y_{k+1} &= y_k + b_1 K_1 + b_2 K_2 \\[1em] K_1 &= h f(t_k + c_1 h, y_k) \\[1em] K_2 &= h f\left(t_k + c_2 h, y_k + A_{21} K_1\right) \end{aligned}
Questa è la formulazione classica dei metodi di Runge-Kutta, che permette di approssimare soluzioni di ODE fino a un errore del secondo ordine con combinazioni lineari delle valutazioni di f in punti intermedi.
Estensione della Procedura
Il metodo di Runge-Kutta può essere esteso includendo termini aggiuntivi per ottenere una maggiore precisione. Consideriamo una formulazione estesa del metodo:
\begin{aligned} y_{k+1} &= \Phi(t_k, y_k, h) \\ &= y_k + A h f(t_k + B h, y_k + C h f(t_k, y_k)) \\ &\qquad + D h f(t_k + E h, y_k + F h f(t_k, y_k)). \end{aligned}
Anche in questo caso, utilizzando l’espansione di Taylor e minimizzando l’errore locale di troncamento, otteniamo un sistema di condizioni non lineari sui coefficienti A, B, C, D, E e F che devono essere soddisfatte. Queste condizioni conducono a intere famiglie di metodi di Runge-Kutta con proprietà diverse in termini di accuratezza e stabilità.
Un metodo di Runge-Kutta generico è costruito a partire da una formula di avanzamento nella forma generale:
\begin{aligned} y_{k+1} &= y_k + \sum_{i=1}^{s} b_i K_i, \\ \text{dove} & \\ K_1 &= h f\Big(t_k + h c_1, y_k + \sum_{j=1}^{s} a_{1j} K_j\Big), \\[1em] K_2 &= h f\Big(t_k + h c_2, y_k + \sum_{j=1}^{s} a_{2j} K_j\Big), \\[1em] & \vdots \\[1em] K_s &= h f\Big(t_k + h c_s, y_k + \sum_{j=1}^{s} a_{sj} K_j\Big). \end{aligned}
Questa formulazione introduce un totale di 2s + s^2 coefficienti, dove:
- b_1, b_2, \ldots, b_s rappresentano i pesi per i termini K_i,
- c_1, c_2, \ldots, c_s rappresentano i coefficienti temporali,
- a_{ij} (con i, j = 1, 2, \ldots, s) rappresentano i coefficienti delle combinazioni lineari delle valutazioni di f nei punti intermedi.
Per determinare questi coefficienti, si applica la procedura della serie di Taylor in più variabili e si confrontano i termini della serie con l’espansione della soluzione esatta. Questo permette di ridurre l’errore locale di troncamento a un ordine \mathcal{O}(h^p) desiderato.
Poiché il calcolo di queste condizioni può essere complesso, esistono già tabelle che forniscono le condizioni di ordine per i coefficienti b_i, c_j, e a_{ij} fino a ordini molto elevati.
30.2 Condizioni di Ordine
Ordine 1 (Metodo di Eulero)
Un metodo di Runge-Kutta ha ordine 1 se soddisfa la seguenti condizioni:
\qquad\displaystyle\sum_{j=1}^{s} a_{ij} = c_i,\qquad i=1,2,\ldots,s
\qquad\displaystyle\sum_{i=1}^{s} b_i = 1
Ordine 2
Un metodo di Runge-Kutta ha ordine 2 se soddisfa in aggiunta la seguente condizione:
- \qquad\displaystyle\sum_{i=1}^{s} b_i c_i = \dfrac{1}{2}
Ordine 3
Un metodo di Runge-Kutta ha ordine 3 se soddisfa in aggiunta le seguenti condizioni:
\qquad\displaystyle\sum_{i=1}^{s} b_i c_i^2 = \dfrac{1}{3}
\qquad\displaystyle\sum_{i,j=1}^{s} b_i a_{ij} c_j = \dfrac{1}{6}
Ordine 4
Un metodo di Runge-Kutta ha ordine 4 se soddisfa in aggiunta le seguenti condizioni:
\qquad\displaystyle\sum_{i=1}^{s} b_i c_i^3 = \dfrac{1}{4}
\qquad\displaystyle\sum_{i,j=1}^{s} b_i c_i a_{ij} c_j = \dfrac{1}{8}
\qquad\displaystyle\sum_{i,j=1}^{s} b_i a_{ij} c_j^2 = \dfrac{1}{12}
\qquad\displaystyle\sum_{i,j,k=1}^{s} b_i a_{ij} a_{jk} c_k = \dfrac{1}{24}
Ordini Superiori
Le condizioni di ordine per p > 4 diventano sempre più complesse e richiedono un numero maggiore di condizioni sui coefficienti.
Ordine 5 richiede 17 condizioni.
Ordine 6 richiede 37 condizioni.
Ordine 7 richiede 85 condizioni.
Ordine 8 richiede 200 condizioni.
Ordine 9 richiede 486 condizioni.
Ordine 10 richiede 1205 condizioni. \large \color{magenta} \left. \begin{aligned} &\color{DarkGreen}\underbrace{ \color{DarkRed}\overbrace{\color{blue}\overbrace{ \underbrace{\color{blue}\bm{A}\bm{1}=\bm{c}}_{\text{consistenza}},\quad \color{blue}\bm{b}\cdot\bm{1}=1, }^{\text{ordine 1}}\quad \color{DarkRed}\bm{b}\cdot\bm{c} = \dfrac{1}{2}}^{\text{ordine 2}},\quad \color{DarkGreen}\bm{b}\cdot\bm{c}^2 = \dfrac{1}{3},\quad \bm{b}\cdot\bm{d} = \dfrac{1}{6} }_{\text{ordine 3}} \\ & \color{magenta}\bm{b}\cdot\bm{c}^3 = \dfrac{1}{4},\quad (\bm{b}\otimes\bm{c})\cdot\bm{d} = \dfrac{1}{8},\quad \bm{e}\cdot\bm{c}^2 = \dfrac{1}{12},\quad \bm{e}\cdot\bm{d} = \dfrac{1}{24} \quad \end{aligned} \right\} _{\color{magenta}\text{ordine 4}}
\bm{d} = \bm{A}\bm{c}, \qquad \bm{e} = \bm{A}^T\bm{b}, \qquad \bm{c}^2 = \bm{c}\otimes\bm{c},\qquad \bm{c}^3 = \bm{c}\otimes\bm{c}\otimes\bm{c}
Di seguito, riportiamo una lista dei principali metodi di Runge-Kutta comunemente utilizzati.
30.3 Metodi di Runge-Kutta Espliciti
La forma generale di un metodo di Runge-Kutta esplicito è un metodo a un passo che approssima la soluzione utilizzando più fasi intermedie all’interno di ciascun passo temporale. Il metodo è rappresentato come segue:
\begin{aligned} y_{k+1} &= y_k + \sum_{i=1}^{s} b_i K_i, \\ &\left\{\begin{aligned} K_1 &= h f\Big(t_k, y_k\Big), \\[1em] K_2 &= h f\Big(t_k + h c_2, y_k + a_{21} K_1\Big), \\[1em] K_3 &= h f\Big(t_k + h c_3, y_k + a_{31} K_1 + a_{32} K_2\Big), \\[1em] \vdots & \\[1em] K_s &= h f\Big(t_k + h c_s, y_k + \sum_{j=1}^{s-1} a_{sj} K_j\Big), \end{aligned} \right. \end{aligned}
dove K_1,K_2,\ldots,K_s sono valori intermedi chiamati stadi, e s è il numero di stadi. Le costanti a_{ij}, b_i e c_i definiscono il metodo di Runge-Kutta specifico e sono scelte per raggiungere un’accuratezza di ordine superiore. Il metodo calcola più somme pesate delle valutazioni della funzione (i valori K_i) per ottenere approssimazioni più accurate della soluzione.
I valori t_k+hc_i rappresentano punti di quadratura nell’intervallo temporale [t_k,t_{k+1}], e le quantità K_i/h approssimano i valori della funzione in questi punti.
Tableau di Butcher
Per semplificare la notazione dei metodi di Runge-Kutta, i coefficienti a_{ij}, b_i e c_i vengono spesso riassunti in una forma compatta chiamata tableau di Butcher5. Il tableau per il metodo di Runge-Kutta esplicito è dato da:
\begin{array} {c|ccccc} 0 & 0 & 0 & 0 & 0 & 0 \\ c_2 & a_{21} & 0 & 0 & 0 & 0 \\ c_3 & a_{31} & a_{32} & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & 0 & 0 \\ c_s & a_{s1} & a_{s2} & \ldots & a_{s,s-1} & 0 \\ \hline & b_1 & b_2 & \ldots & b_{s-1} & b_s \end{array}
Euler Esplicito
Il metodo di Euler esplicito rappresentato utilizzando un tableau di Butcher:
\def\arraystretch{1.5} \begin{array}{c|c} 0 & 0 \\ \hline & 1 \end{array} \tag{30.2}
Metodo del Punto Medio
Il Metodo del Punto Medio, noto anche come Metodo di Collatz, è una variante del metodo di Euler esplicito. Questo approccio è classificato come metodo del secondo ordine e consiste in due fasi, consentendo di migliorare l’accuratezza nell’integrazione numerica.
\def\arraystretch{1.5} \begin{array}{c|cc} 0 & 0 & 0 \\ 1/2 & 1/2 & 0 \\ \hline & 0 & 1 \end{array} \tag{30.3}
Metodo di Heun
Il Metodo di Heun, spesso noto come la regola del trapezio esplicita, rappresenta un miglioramento rispetto al metodo di Euler esplicito standard. Questa tecnica è classificata come un metodo del secondo ordine e utilizza due fasi per ottenere una maggiore accuratezza nell’integrazione numerica.
\def\arraystretch{1.5} \begin{array}{c|cc} 0 & 0 & 0 \\ 1 & 1 & 0 \\ \hline & 1/2 & 1/2 \end{array} \tag{30.4}
Metodo di Ralston
Il metodo di Ralston è una tecnica del secondo ordine che utilizza due fasi, progettata specificamente per ottenere un limite minimo all’errore locale.
\def\arraystretch{1.5} \begin{array}{c|cc} 0 & 0 & 0 \\ 2/3 & 2/3 & 0 \\ \hline & 1/4 & 3/4 \end{array} \tag{30.5}
Metodo generico del secondo ordine
\def\arraystretch{1.5} \begin{array}{c|cc} 0 & 0 & 0 \\ \alpha & \alpha & 0 \\ \hline\\[-1em] & 1-\dfrac{1}{2\alpha} & \dfrac{1}{2\alpha} \end{array} \qquad \alpha\in(0,1) \tag{30.6}
Runge-Kutta del terzo ordine (RK3)
\def\arraystretch{1.5} \begin{array} {c|ccc} 0 & 0 & 0 & 0 \\ 1/2 & 1/2 & 0 & 0 \\ 1 & -1 & 2 & 0 \\ \hline & 1/6 & 2/3 & 1/6 \end{array} \tag{30.7}
Metodo del terzo ordine di Heun
\def\arraystretch{1.5} \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ 1/3 & 1/3 & 0 & 0 \\ 2/3 & 0 & 2/3 & 0 \\ \hline & 1/4 & 0 & 3/4 \end{array} \tag{30.8}
Metodo del terzo ordine di Ralston
\def\arraystretch{1.5} \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ 1/2 & 1/2 & 0 & 0 \\ 3/4 & 0 & 3/4 & 0 \\ \hline & 2/9 & 1/3 & 4/9 \end{array} \tag{30.9}
Runge-Kutta SSP (SSPRK3)
\def\arraystretch{1.5} \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 1/2 & 1/4 & 1/4 & 0 \\ \hline & 1/6 & 1/6 & 2/3 \end{array} \tag{30.10}
Metodo generico del terzo ordine
Con \alpha \neq 0, 2/3, 1:
\def\arraystretch{1.5} \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ \alpha & \alpha & 0 & 0 \\ 1 & 1+\dfrac{1-\alpha}{\alpha(3\alpha-2)} & -\dfrac{1-\alpha}{\alpha(3\alpha-2)} & 0 \\[1em] \hline \\[-1em] & \dfrac{1}{2}-\dfrac{1}{6\alpha} & \dfrac{1}{6\alpha(1-\alpha)} & \dfrac{2-3\alpha}{6\alpha(1-\alpha)} \end{array} \tag{30.11}
Metodo di Runge-Kutta classico (RK4)
\def\arraystretch{1.5} \begin{array} {c|cccc} 0 & 0 & 0 & 0 & 0 \\ 1/2 & 1/2 & 0 & 0 & 0 \\ 1/2 & 0 & 1/2 & 0 & 0 \\ 1 & 0 & 0 & 1 & 0 \\ \hline & 1/6 & 2/6 & 2/6 & 1/6 \end{array} \tag{30.12}
Ralston metodo del quarto ordine
È un metodo del secondo ordine con due stadi, progettato per minimizzare l’errore di troncamento.
\def\arraystretch{1.5} \begin{array}{c|ccccc} 0 & 0 & 0 & 0 & 0 \\ 2/5 & 2/5 & 0 & 0 & 0 \\ c_3 & a_{3,1} & a_{3,2} & 0 & 0 \\ 1 & a_{4,1} & a_{4,2} & a_{4,3} & 0 \\ \hline\\ & b_1 & b_2 & b_3 & b_4 \end{array} \tag{30.13}
dove
\begin{aligned} b_1 & = \dfrac{263}{1812}+\dfrac{2}{151}\sqrt{5}\\ b_2 & = \dfrac{125}{3828}-\dfrac{250}{957}\sqrt{5} \\ b_3 & = \dfrac{3426304}{5924787}+\dfrac{553984}{1974929}\sqrt{5} \\ b_4 & = \dfrac{10}{41}-\dfrac{4}{123}\sqrt{5} \\ c_3 & = \dfrac{7}{8}-\dfrac{3}{16}\sqrt{5} \\ a_{3,1} & = \dfrac{357}{256}\sqrt{5}-\dfrac{2889}{1024} \\ a_{4,1} & = \dfrac{1047}{3020}\sqrt{5}-\dfrac{673}{1208} \\ a_{3,2} &= \dfrac{3785}{1024}-\dfrac{405}{256}\sqrt{5} \\ a_{4,2} &= -\dfrac{975}{2552}-\dfrac{1523}{1276}\sqrt{5} \\ a_{4,3} &= \dfrac{93408}{48169}+\dfrac{203968}{240845}\sqrt{5} \end{aligned}
Runge-Kutta del quarto ordine: regola 3/8
\def\arraystretch{1.5} \begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0 \\ 1/3 & 1/3 & 0 & 0 & 0 \\ 2/3 & -1/3 & 1 & 0 & 0 \\ 1 & 1 & -1 & 1 & 0 \\ \hline & 1/8 & 3/8 & 3/8 & 1/8 \end{array}
30.4 Controllo adattivo del passo
Finora il passo h è stato tenuto costante. Nei codici reali non lo è quasi mai: dove la soluzione varia poco conviene un passo grande, dove varia rapidamente un passo piccolo. Per farlo serve una stima dell’errore locale che costi poco.
Coppie annidate
L’idea delle coppie annidate (embedded pairs) è di combinare gli stessi stadi K_i con due diversi vettori di pesi, \bm{b} e \widehat{\bm{b}}, ottenendo due approssimazioni di ordine diverso:
y_{k+1} = y_k + \sum_{i=1}^{s} b_i K_i, \qquad \widehat{y}_{k+1} = y_k + \sum_{i=1}^{s} \widehat{b}_i K_i
La loro differenza
\boxed{\; e^{\mathrm{est}}_{k+1} = y_{k+1} - \widehat{y}_{k+1} = \sum_{i=1}^{s} \big( b_i - \widehat{b}_i \big) K_i \;} \tag{30.14}
è una stima dell’errore locale che non costa nessuna valutazione aggiuntiva di f: è questo il punto. Il tableau si scrive con due righe sotto la linea:
\def\arraystretch{1.5} \begin{array}{c|c} \bm{c} & \bm{A} \\ \hline & \bm{b}^T \\ \hline & \widehat{\bm{b}}^T \end{array}
Accettazione del passo e controllore
Per un sistema la stima si normalizza componente per componente con una tolleranza assoluta e una relativa:
\mathrm{err} = \left[ \dfrac{1}{d}\sum_{j=1}^{d} \left( \dfrac{e^{\mathrm{est}}_j} {\mathrm{atol}_j + \mathrm{rtol}_j \max\big(|y_{k,j}|,|y_{k+1,j}|\big)} \right)^{2} \right]^{1/2}
Il passo è accettato se \mathrm{err}\leq 1, altrimenti è rifiutato e ripetuto. In entrambi i casi il passo successivo è
\boxed{\; h_{\mathrm{new}} = \eta\, h\, \mathrm{err}^{-1/q} \;} \tag{30.15}
dove q è l’ordine della stima 30.14 e \eta\in(0,1) è un fattore di sicurezza (tipicamente 0.8–0.9). L’esponente -1/q segue dall’ipotesi e^{\mathrm{est}}\approx Ch^{q}: per avere \mathrm{err}=1 occorre scalare h di \mathrm{err}^{-1/q}. In pratica il fattore si limita a un intervallo, ad esempio [0.2,\,5], per evitare oscillazioni del passo.
Esempio: la coppia di Bogacki-Shampine 3(2)
È il metodo usato da ode23 di MATLAB e da BS3 in molte librerie.
\def\arraystretch{1.6} \begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0 \\ \dfrac{1}{2} & \dfrac{1}{2} & 0 & 0 & 0 \\ \dfrac{3}{4} & 0 & \dfrac{3}{4} & 0 & 0 \\ 1 & \dfrac{2}{9} & \dfrac{1}{3} & \dfrac{4}{9} & 0 \\ \hline & \dfrac{2}{9} & \dfrac{1}{3} & \dfrac{4}{9} & 0 \\ \hline & \dfrac{7}{24} & \dfrac{1}{4} & \dfrac{1}{3} & \dfrac{1}{8} \end{array} \tag{30.16}
Il primo vettore di pesi dà ordine 3, il secondo ordine 2, quindi e^{\mathrm{est}} = \mathcal{O}(h^{3}) e nella 30.15 si usa q=3.
Il controllo adattivo del passo è la ragione principale per cui i metodi espliciti moderni sono molto più efficienti di RK4 a passo fisso nei codici di uso generale: non perché siano più accurati a parità di passo, ma perché scelgono il passo da soli.
30.5 Metodi di Runge-Kutta Impliciti
Finora abbiamo discusso esclusivamente dei metodi di Runge-Kutta espliciti. Tuttavia, quando si affrontano equazioni rigide, è consigliabile utilizzare i metodi di Runge-Kutta impliciti.
Il metodo di Runge-Kutta implicito generico è un metodo a un passo caratterizzato dalla seguente rappresentazione:
\begin{aligned} y_{k+1} &= y_k + \sum_{i=1}^{s}b_iK_i \\ &\left\{\begin{aligned} K_1 &= hf\Big(t_k+hc_1,y_k + \displaystyle\sum_{j=1} ^{s}a_{1j}K_j\Big) \\[1em] K_2 &= hf\Big(t_k+hc_2,y_k + \displaystyle\sum_{j=1} ^{s}a_{2j}K_j\Big) \\[1em] \vdots \\[1em] K_s &= hf\Big(t_k+hc_s,y_k + \displaystyle\sum_{j=1} ^{s}a_{sj}K_j\Big) \end{aligned} \right. \end{aligned} \tag{30.17}
Il tableau per 30.17 sarà quindi della forma:
\def\arraystretch{1.5} \begin{array} {c|cccc} c_1 & a_{11} & a_{12} & \ldots & a_{1s} \\ c_2 & a_{21} & a_{22} & \ldots & a_{2s} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ c_s & a_{s1} & a_{s2} & \ldots & a_{ss} \\ \hline & b_1 & b_2 & \ldots & b_s \end{array}
Soluzione degli stadi impliciti
La 30.17 non è una formula ma un sistema non lineare: gli s stadi K_1,\ldots,K_s compaiono a destra e a sinistra. Per un problema di dimensione d si tratta di un sistema di s\cdot d equazioni da risolvere a ogni passo. È qui che i metodi dei capitoli precedenti tornano utili.
Raccogliendo gli stadi nel vettore \bm{K}=(K_1,\ldots,K_s), il sistema è
\bm{G}(\bm{K}) = \bm{0}, \qquad G_i(\bm{K}) = K_i - h f\Big( t_k + c_i h,\; y_k + \sum_{j=1}^{s} a_{ij}K_j \Big) \tag{30.18}
Iterazione di punto fisso
La scelta più semplice è riscrivere la 30.18 come iterazione di punto fisso (Capitolo 23):
K_i^{(m+1)} = h f\Big( t_k + c_i h,\; y_k + \sum_{j=1}^{s} a_{ij}K_j^{(m)} \Big), \qquad m=0,1,2,\ldots
partendo ad esempio da K_i^{(0)} = h f(t_k,y_k). Se f è lipschitziana con costante L rispetto a y, la mappa è una contrazione non appena
\boxed{\; h\, L \max_{i} \sum_{j=1}^{s} |a_{ij}| < 1 \;} \tag{30.19}
e in tal caso l’iterazione converge linearmente con fattore pari al membro sinistro della 30.19.
Metodo di Newton
Applicando il metodo di Newton (Capitolo 20) alla 30.18 si itera
\bm{J}_{\bm{G}}\big(\bm{K}^{(m)}\big)\, \Delta\bm{K}^{(m)} = -\bm{G}\big(\bm{K}^{(m)}\big), \qquad \bm{K}^{(m+1)} = \bm{K}^{(m)} + \Delta\bm{K}^{(m)}
dove il jacobiano ha struttura a blocchi s\times s, ciascun blocco di dimensione d\times d:
\boxed{\; \big[\bm{J}_{\bm{G}}\big]_{ij} = \delta_{ij}\,\bm{I} - h\, a_{ij}\, \bm{J}_{\bm{f}}\Big( t_k + c_j h,\; y_k + \sum_{l} a_{jl}K_l \Big) \;} \tag{30.20}
con \bm{J}_{\bm{f}} = \partial \bm{f}/\partial \bm{y} e \delta_{ij} la delta di Kronecker. Il sistema lineare ad ogni iterazione si risolve con la fattorizzazione LU (Capitolo 1) oppure, se sd è grande e la matrice è sparsa, con un metodo iterativo (Capitolo 26).
A differenza del punto fisso, la convergenza di Newton è quadratica e non richiede alcuna condizione del tipo 30.19: il passo h resta libero, che è esattamente il motivo per cui si è scelto un metodo implicito.
Metodo di Eulero Implicito
Il metodo di Eulero implicito rappresentato utilizzando un tableau:
\def\arraystretch{1.5} \begin{array}{c|c} 1 & 1 \\ \hline & 1 \end{array} \tag{30.21}
Metodo del Punto Medio Implicito
Il metodo del punto medio implicito è un’adattamento del metodo di Eulero esplicito, spesso chiamato metodo di Collatz. Si tratta di un metodo del secondo ordine caratterizzato da due stadi, che offre una maggiore precisione nella risoluzione delle equazioni differenziali ordinarie.
Metodo originale y_{k+1} = y_k + h f\left( t_k + \dfrac{1}{2}h, \dfrac{y_k+y_{k+1}}{2} \right)
formulazione come Runge-Kutta \begin{aligned} y_{k+1} & = y_k + K_1 \\ K_1 &= h f\left( t_k + \dfrac{1}{2}h, \dfrac{y_{k+1}+y_k}{2} \right) \\ &= h f\left( t_k + \dfrac{1}{2}h, y_k + \dfrac{y_{k+1}-y_k}{2} \right) \\ &= h f\left( t_k + \dfrac{1}{2}h, y_k + \dfrac{1}{2}K_1\right) \end{aligned}
per cui il tableau diventa \def\arraystretch{1.5} \begin{array}{c|c} 1/2 & 1/2 \\ \hline & 1 \end{array} \tag{30.22}
Metodo di Crank-Nicolson
Il metodo di Crank-Nicolson6 corrisponde alla regola trapezoidale implicita.
Metodo originale y_{k+1} = y_k + \dfrac{1}{2}h \left(f( t_k,y_k)+f( t_{k+1},y_{k+1})\right)
formulazione come Runge-Kutta \begin{aligned} y_{k+1} & = y_k + \dfrac{1}{2}K_1+ \dfrac{1}{2}K_2 \\ K_1 &= h f( t_k,y_k) \\ K_2 &= h f( t_k + h, y_{k+1}) \\ &= h f\left( t_k + h, y_k + \dfrac{1}{2}K_1+ \dfrac{1}{2}K_2\right) \\ \end{aligned}
per cui il tableau diventa
\def\arraystretch{1.5} \begin{array}{c|cc} 0 & 0 & 0 \\ 1 & 1/2 & 1/2 \\ \hline & 1/2 & 1/2 \end{array} \tag{30.23}
Metodi di collocazione
I tableau impliciti che seguono non sono stati tirati a sorte: nascono tutti da una sola idea, la collocazione, che mette insieme l’interpolazione (Capitolo 7) e la quadratura (Capitolo 10).
Si fissano s nodi distinti 0\leq c_1<\cdots<c_s\leq 1 e si cerca un polinomio u(t) di grado \leq s tale che
u(t_k) = y_k, \qquad u'\big(t_k + c_i h\big) = f\Big( t_k + c_i h,\; u(t_k+c_i h) \Big), \quad i=1,\ldots,s \tag{30.24}
cioè che soddisfi l’equazione differenziale esattamente negli s nodi. Si pone poi y_{k+1} = u(t_{k+1}).
Teorema 30.1 (Collocazione e Runge-Kutta) Il metodo di collocazione 30.24 coincide con il metodo di Runge-Kutta implicito di tableau
a_{ij} = \int_{0}^{c_i} \ell_j(\sigma)\,\mathrm{d}\sigma, \qquad b_j = \int_{0}^{1} \ell_j(\sigma)\,\mathrm{d}\sigma \tag{30.25}
dove \ell_j sono i polinomi elementari di Lagrange sui nodi c_1,\ldots,c_s.
Dimostrazione. Il polinomio u' ha grado \leq s-1 ed è noto negli s nodi, dove vale K_i/h con K_i = h f(t_k+c_ih, u(t_k+c_ih)). Per l’unicità del polinomio interpolante (Capitolo 6),
u'(t_k+\sigma h) = \dfrac{1}{h}\sum_{j=1}^{s} K_j\, \ell_j(\sigma)
Integrando da 0 a c_i e poi da 0 a 1 si ottengono rispettivamente
u(t_k+c_i h) = y_k + \sum_{j=1}^{s} K_j \int_{0}^{c_i}\ell_j(\sigma)\,\mathrm{d}\sigma, \qquad u(t_{k+1}) = y_k + \sum_{j=1}^{s} K_j \int_{0}^{1}\ell_j(\sigma)\,\mathrm{d}\sigma
che sono esattamente le equazioni degli stadi e del passo di un metodo di Runge-Kutta con i coefficienti 30.25. \square
La 30.25 spiega l’intera tabella dei metodi impliciti che segue, e spiega anche il loro ordine: poiché b_j sono i pesi della formula di quadratura interpolatoria sui nodi c_j, il metodo di collocazione eredita l’ordine della quadratura.
| Nodi c_i | Quadratura | Metodo RK | Stadi | Ordine |
|---|---|---|---|---|
| zeri di Legendre su (0,1) | Gauss-Legendre | Gauss-Legendre | s | 2s |
| estremi inclusi | Gauss-Lobatto | Lobatto IIIA | s | 2s-2 |
| un estremo incluso | Gauss-Radau | Radau IIA | s | 2s-1 |
| c_1=1 | — | Eulero implicito | 1 | 1 |
| c_1=\tfrac12 | punto medio | punto medio implicito | 1 | 2 |
Il salto da ordine s (quanto basterebbe a un metodo generico a s stadi) a ordine 2s è lo stesso fenomeno della quadratura gaussiana: la scelta ottimale dei nodi raddoppia l’ordine. Per questo i metodi di Gauss-Legendre sono i metodi impliciti di ordine massimo a parità di stadi.
Metodi di Gauss-Legendre
Questi metodi si basano sui punti della quadratura di Gauss-Legendre7.
Metodo di Gauss-Legendre del quarto ordine
\def\arraystretch{3} \begin{array}{c|cc} \dfrac{1}{2}-\dfrac{\sqrt3}{6} & \dfrac{1}{4} & \dfrac{1}{4}-\dfrac{\sqrt3}{6} \\ \dfrac{1}{2}+\dfrac{\sqrt3}{6} & \dfrac{1}{4}+\dfrac{\sqrt3}{6} & \dfrac{1}{4} \\ \hline & \dfrac{1}{2} & \dfrac{1}{2} \\ & \dfrac{1}{2}+\dfrac{\sqrt{3}}{2} & \dfrac{1}{2}-\dfrac{\sqrt{3}}{2} \end{array} \tag{30.26}
Metodo di Gauss-Legendre del sesto ordine
\def\arraystretch{3} \begin{array}{c|ccc} \dfrac{1}{2} - \dfrac{\sqrt{15}}{10} & \dfrac{5}{36} & \dfrac{2}{9}- \dfrac{\sqrt{15}}{15} & \dfrac{5}{36} - \dfrac{\sqrt{15}}{30} \\ \dfrac{1}{2} & \dfrac{5}{36} + \dfrac{\sqrt{15}}{24} & \dfrac{2}{9} & \dfrac{5}{36} - \dfrac{\sqrt{15}}{24} \\ \dfrac{1}{2} + \dfrac{\sqrt{15}}{10} & \dfrac{5}{36} + \dfrac{\sqrt{15}}{30} & \dfrac{2}{9} + \dfrac{\sqrt{15}}{15} & \dfrac{5}{36} \\[1em] \hline & \dfrac{5}{18} & \dfrac{4}{9} & \dfrac{5}{18} \\ & -\dfrac{5}{6} & \dfrac{8}{3} & -\dfrac{5}{6} \end{array} \tag{30.27}
Metodi di Lobatto
Metodo di Lobatto IIIA del quarto ordine
\def\arraystretch{1.5} \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ 1/2 & 5/24& 1/3 & -1/24 \\ 1 & 1/6 & 2/3 & 1/6 \\ \hline & 1/6 & 2/3 & 1/6 \end{array} \tag{30.28}
Lobatto IIIB del quarto ordine
\def\arraystretch{1.5} \begin{array}{c|ccc} 0 & 1/6 & -1/6 & 0 \\ 1/2 & 1/6 & 1/3 & 0 \\ 1 & 1/6 & 5/6 & 0 \\ \hline & 1/6 & 2/3 & 1/6 \end{array} \tag{30.29}
Lobatto IIIC del quarto ordine
\def\arraystretch{1.5} \begin{array}{c|ccc} 0 & 1/6 & -1/3 & 1/6 \\ 1/2 & 1/6 & 5/12 & -1/12 \\ 1 & 1/6 & 2/3 & 1/6 \\ \hline & 1/6 & 2/3 & 1/6 \end{array} \tag{30.30}
Lobatto IIIC^\star del quarto ordine
\def\arraystretch{1.5} \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ 1/2 & 1/4 & 1/4 & 0 \\ 1 & 0 & 1 & 0 \\ \hline & 1/6 & 2/3 & 1/6 \end{array} \tag{30.31}
Metodi di Radau
Metodo di Radau IA del quinto ordine
\def\arraystretch{3} \begin{array}{c|ccc} 0 & \dfrac{1}{9} & \dfrac{-1 - \sqrt{6}}{18} & \dfrac{-1 + \sqrt{6}}{18} \\ \dfrac{3}{5} - \dfrac{\sqrt{6}}{10} & \dfrac{1}{9} & \dfrac{11}{45} + \dfrac{7\sqrt{6}}{360} & \dfrac{11}{45} - \dfrac{43\sqrt{6}}{360}\\ \dfrac{3}{5} + \dfrac{\sqrt{6}}{10} & \dfrac{1}{9} & \dfrac{11}{45} + \dfrac{43\sqrt{6}}{360} & \dfrac{11}{45} - \dfrac{7\sqrt{6}}{360} \\ \hline & \dfrac{1}{9} & \dfrac{4}{9} + \dfrac{\sqrt{6}}{36} & \dfrac{4}{9} - \dfrac{\sqrt{6}}{36} \end{array} \tag{30.32}
Radau II del quinto ordine
\def\arraystretch{3} \begin{array}{c|ccc} \dfrac{2}{5} - \dfrac{\sqrt{6}}{10} & \dfrac{11}{45} - \dfrac{7\sqrt{6}}{360} & \dfrac{37}{225} - \dfrac{169\sqrt{6}}{1800} & -\dfrac{2}{225} + \dfrac{\sqrt{6}}{75} \\ \dfrac{2}{5} + \dfrac{\sqrt{6}}{10} & \dfrac{37}{225} + \dfrac{169\sqrt{6}}{1800} & \dfrac{11}{45} + \dfrac{7\sqrt{6}}{360} & -\dfrac{2}{225} - \dfrac{\sqrt{6}}{75}\\ 1 & \dfrac{4}{9} - \dfrac{\sqrt{6}}{36} & \dfrac{4}{9} + \dfrac{\sqrt{6}}{36} & \dfrac{1}{9} \\ \hline & \dfrac{4}{9} - \dfrac{\sqrt{6}}{36} & \dfrac{4}{9} + \dfrac{\sqrt{6}}{36} & \dfrac{1}{9} \end{array} \tag{30.33}
In sintesi
- I metodi di Runge-Kutta ottengono ordine alto valutando f in più stadi interni, senza mai derivare f; il tableau di Butcher (\bm{A},\bm{b},\bm{c}) li individua completamente.
- Le condizioni d’ordine si ricavano confrontando lo sviluppo di Taylor del metodo con quello della soluzione esatta; il loro numero esplode con p (17 condizioni per p=5, 1205 per p=10).
- Il numero di stadi necessari supera l’ordine da p\geq 5 in poi: non esiste un RK esplicito a 5 stadi di ordine 5.
- Le coppie annidate (Sezione 30.4) danno una stima dell’errore locale a costo nullo e permettono il controllo automatico del passo; la proprietà FSAL riduce ulteriormente le valutazioni.
- I metodi impliciti richiedono di risolvere un sistema non lineare di dimensione sd a ogni passo (Sezione 30.5.1): punto fisso sui problemi non stiff, Newton su quelli stiff.
- Gauss-Legendre, Lobatto e Radau non sono tableau arbitrari ma metodi di collocazione (Teorema 30.1): ereditano nodi, pesi e ordine dalle corrispondenti formule di quadratura.