Metodi di Runge-Kutta

In analisi numerica i metodi di Runge-Kutta [ˌʁʊŋəˈkʊta] sono un'importante famiglia di metodi iterativi impliciti ed espliciti per l'approssimazione delle soluzioni delle equazioni differenziali ordinarie. Queste tecniche furono sviluppate intorno al 1900 dai matematici tedeschi Carl Runge e Martin Wilhelm Kutta, dai quali prendono il nome.

I metodi di Runge-Kutta, spesso abbreviati con le iniziali "RK", sono una famiglia di metodi iterativi discreti utilizzati nell'approssimazione numerica di soluzioni di equazioni differenziali ordinarie (ODE), e più specificatamente per problemi ai valori iniziali. Fanno parte della famiglia più generale di metodi discreti per le equazioni differenziali ordinarie, ovvero di quella classe di metodi numerici che fornisce un'approssimazione della soluzione di un'equazione differenziale (o più precisamente di un problema di Cauchy) in un insieme discreto di punti.

Per trovare un'approssimazione della funzione ${\displaystyle y(t):\mathbb {R} \to \mathbb {R} ^{d}}$ che verifichi il generico problema di Cauchy:

${\displaystyle {\begin{cases}y'(t)=f(t,y(t))\\y(t_{0})=y_{0}\end{cases}}}$

in un insieme discreto di punti in cui si considera il problema (solitamente nell'intervallo ${\displaystyle \left[t_{0},t_{f}\right]}$), si considera un campionamento dell'intervallo ${\displaystyle \Delta }$ in un insieme di punti ${\displaystyle \{t_{i}\mid i=0\dots n\}}$, dove ${\displaystyle t_{i}=t_{0}+ih}$ e ${\displaystyle h=(t_{f}-t_{0})/n}$. Il metodo numerico fornisce allora l'approssimazione dei valori ${\displaystyle y(t_{i})}$, e per ottenere una ricostruzione abbastanza fedele della funzione il numero ${\displaystyle n}$ deve essere sufficientemente elevato.

L'idea che sta alla base dei metodi di Runge-Kutta è la trasposizione del problema dalla forma differenziale alla forma integrale, per la quale esistono metodi numerici (come le formule di quadratura) che permettono l'approssimazione della soluzione. In generale un metodo di Runge-Kutta è caratterizzato da tre parametri: un vettore ${\displaystyle b=(b_{i})_{i=0,\dots ,s}}$, una matrice ${\displaystyle a=(a_{i,j})_{i,j=0,\dots ,s}}$ e un vettore ${\displaystyle \theta =(\theta _{i})_{i=0,\dots ,s}}$. L'approssimazione è data dal sistema:

${\displaystyle {\begin{cases}y_{n+1}=y_{n}+h\sum _{i=1}^{s}b_{i}f(t_{n}+\theta _{i}h,Y_{i})\\Y_{i}=y_{n}+h\sum _{j=1}^{s}a_{ij}f(t_{n}+\theta _{j}h,Y_{j})\quad i=1,\dots ,s\end{cases}}}$

dove i valori ${\displaystyle y_{i}}$ approssimano il valore esatto ${\displaystyle y(t_{i})}$.

Considerando il generico problema di Cauchy:

${\displaystyle {\begin{cases}y'(t)=f(t,y(t))\\y(t_{0})=y_{0}\end{cases}}}$

Si può considerare una sua riformulazione in forma integrale:

${\displaystyle y(t)=y_{0}+\int _{t_{0}}^{t}y'(s)\,ds=y_{0}+\int _{t_{0}}^{t}f(s,y(s))\,ds}$

dove l'ultima sostituzione è legittima in quanto la funzione è soluzione dell'equazione differenziale. Da questa riformulazione segue in particolare che, per una suddivisione uniforme ${\displaystyle \Delta =\{t_{i}\}=\{t_{0}+ih\}_{\{i=0,\dots ,n\}}}$, il valore della soluzione nel punto ${\displaystyle t_{1}}$ è dato da:

${\displaystyle y(t_{1})=y_{0}+\int _{t_{0}}^{t_{0}+h}f(s,y(s))\,ds}$

Da questo risultato, attraverso la sostituzione ${\displaystyle s=t_{0}+\theta h}$ si può normalizzare l'intervallo di integrazione e ottenere:

${\displaystyle y(t_{1})=y_{0}+h\int _{0}^{1}f(t_{0}+\theta h,y(t_{0}+\theta h))\,d\theta }$

Da questa scrittura appare evidente che con opportune sostituzioni si può esprimere il valore della funzione in ogni punto intermedio tra ${\displaystyle t_{0}}$ e ${\displaystyle t_{1}}$.

Utilizzando un'approssimazione per via numerica dell'integrale attraverso delle formule di quadratura sui nodi ${\displaystyle \theta _{i}}$ e pesi ${\displaystyle b_{i}}$ (che dipendono dalla scelta della formula di quadratura) si ottiene una stima del valore ${\displaystyle y_{1}}$:

${\displaystyle y_{1}=y_{0}+h\sum _{i=0}^{\nu }b_{i}f(t_{0}+\theta _{i}h,Y_{i})}$

dove i valori ${\displaystyle Y_{i}}$ sono approssimazioni di ${\displaystyle y(t_{0}+\theta _{i}h)}$, che sono ignote a priori. Applicando nuovamente il ragionamento precedente si può scrivere che:

${\displaystyle y(t_{0}+\theta _{i}h)=y_{0}+h\int _{0}^{\theta _{i}}f(t_{0}+vh,y(t_{0}+vh))\,dv}$

Questi valori possono essere approssimati a loro volta utilizzando una formula di quadratura con pesi ${\displaystyle a_{ij}}$; si ottiene così che (per semplicità si considerano formule con i medesimi nodi di quadratura di quelle usate in precedenza):

${\displaystyle Y_{i}=y_{0}+h\sum _{j=1}^{\nu }a_{ij}f(t_{0}+\theta _{j}h,Y_{j})}$

Iterando tale costruzione si ottengono le approssimazioni nei punti ${\displaystyle t_{i}}$, e come conseguenza la formulazione generale per i metodi di Runge-Kutta.

Dato il problema ai valori iniziali:

${\displaystyle {\dot {y}}=f(t,y)\qquad y(t_{0})=y_{0}}$

dove i valori di ${\displaystyle t_{0}}$ e ${\displaystyle y_{0}}$ sono noti, si consideri un intervallo sufficientemente piccolo ${\displaystyle h>0}$ e si definiscano:

${\displaystyle {\begin{aligned}y_{n+1}&=y_{n}+{\tfrac {h}{6}}\left(k_{1}+2k_{2}+2k_{3}+k_{4}\right)\\t_{n+1}&=t_{n}+h\\\end{aligned}}}$

per ${\displaystyle n=1,2,3,\dots }$. In questo modo ${\displaystyle y(t_{n+1})}$ è approssimato con ${\displaystyle y_{n+1}}$, e ${\displaystyle y_{n+1}}$ è determinato da ${\displaystyle y_{n}}$ più la media pesata di quattro incrementi ${\displaystyle k_{1}}$, ${\displaystyle k_{2}}$, ${\displaystyle k_{3}}$ e ${\displaystyle k_{4}}$:

${\displaystyle {\begin{aligned}k_{1}&=f(t_{n},y_{n})\\k_{2}&=f(t_{n}+{\tfrac {h}{2}},y_{n}+{\tfrac {1}{2}}k_{1}h)\\k_{3}&=f(t_{n}+{\tfrac {h}{2}},y_{n}+{\tfrac {1}{2}}k_{2}h)\\k_{4}&=f(t_{n}+h,y_{n}+k_{3}h)\end{aligned}}}$

dove ogni incremento è il prodotto di ${\displaystyle h}$ e una stima della pendenza di ${\displaystyle f}$. Nello specifico:

• ${\displaystyle k_{1}}$ è l'incremento basato sulla pendenza all'inizio dell'intervallo, utilizzando ${\displaystyle {\dot {y}}_{n}}$ (metodo di Eulero)
• ${\displaystyle k_{2}}$ è l'incremento basato sulla pendenza alla metà dell'intervallo, utilizzando ${\displaystyle k_{1}}$
• ${\displaystyle k_{3}}$ è un altro incremento basato sulla pendenza alla metà dell'intervallo, utilizzando ${\displaystyle k_{2}}$
• ${\displaystyle k_{4}}$ è l'incremento basato sulla pendenza alla fine dell'intervallo, utilizzando ${\displaystyle k_{3}}$

Nel fare la media, gli incrementi valutati in un punto intermedio dell'intervallo hanno peso maggiore, ed i coefficienti sono scelti in modo che se ${\displaystyle f}$ è indipendente da ${\displaystyle y}$, sicché l'equazione dipende da un semplice integrale, allora il metodo RK coincide con la regola di Cavalieri-Simpson.

In generale, la famiglia di metodi espliciti di Runge-Kutta è data da:

${\displaystyle y_{n+1}=y_{n}+\sum _{i=1}^{s}hb_{i}k_{i}}$

dove:

${\displaystyle {\begin{aligned}k_{1}&=f(t_{n},y_{n})\\k_{2}&=f(t_{n}+c_{2}h,y_{n}+h(a_{21}k_{1}))\\k_{3}&=f(t_{n}+c_{3}h,y_{n}+h(a_{31}k_{1}+a_{32}k_{2}))\\&\vdots \\k_{s}&=f(t_{n}+c_{s}h,y_{n}+h(a_{s1}k_{1}+a_{s2}k_{2}+\cdots +a_{s,s-1}k_{s-1}))\\\end{aligned}}}$

Per specificare un particolare metodo si deve definire il numero ${\displaystyle s}$ ed i coefficienti ${\displaystyle a_{ij}}$ (per ${\displaystyle 1\leq j<i\leq s}$), ${\displaystyle b_{i}}$ (per ${\displaystyle i=1,2,\dots ,s}$) e ${\displaystyle c_{i}}$ (per ${\displaystyle i=2,3,\dots ,s}$). La matrice dei coefficienti ${\displaystyle a_{ij}}$ è detta matrice di Runge-Kutta. Si usa spesso rappresentare i coefficienti in una tabella:

	0
	${\displaystyle c_{2}}$	${\displaystyle a_{21}}$
	${\displaystyle c_{3}}$	${\displaystyle a_{31}}$	${\displaystyle a_{32}}$
	${\displaystyle \vdots }$	${\displaystyle \vdots }$		${\displaystyle \ddots }$
	${\displaystyle c_{s}}$	${\displaystyle a_{s1}}$	${\displaystyle a_{s2}}$	${\displaystyle \cdots }$	${\displaystyle a_{s,s-1}}$
		${\displaystyle b_{1}}$	${\displaystyle b_{2}}$	${\displaystyle \cdots }$	${\displaystyle b_{s-1}}$	${\displaystyle b_{s}}$

e si deve avere:

${\displaystyle \sum _{j=1}^{i-1}a_{ij}=c_{i}\qquad i=2,\ldots ,s}$

In generale, un metodo di Runge-Kutta di ordine ${\displaystyle s}$ può essere scritto come:

${\displaystyle y_{t+h}=y_{t}+h\cdot \sum _{i=1}^{s}a_{i}k_{i}+{\mathcal {O}}(h^{s+1})}$

dove:

${\displaystyle k_{i}=f\left(y_{t}+h\cdot \sum _{j=1}^{s}\beta _{ij}k_{j},t_{n}+\alpha _{i}h\right)}$

sono gli incrementi ottenuti valutando le derivate di ${\displaystyle y_{t}}$ all'${\displaystyle i}$-esimo ordine.

Un modo per derivare il metodo per ${\displaystyle s=4}$ utilizzando la formula generale è il seguente.^[1] Si sceglie:

${\displaystyle {\begin{array}{|l|l|}\hline \alpha _{i}&\beta _{ij}\\\hline \alpha _{1}=0&\beta _{21}={\frac {1}{2}}\\\alpha _{2}={\frac {1}{2}}&\beta _{32}={\frac {1}{2}}\\\alpha _{3}={\frac {1}{2}}&\beta _{43}=1\\\alpha _{4}=1&\\\hline \end{array}}}$

con ${\displaystyle \beta _{ij}=0}$ altrimenti. Si definiscono quindi le quantità:

${\displaystyle {\begin{aligned}y_{t+h}^{1}&=y_{t}+hf\left(y_{t},t\right)\\y_{t+h}^{2}&=y_{t}+hf\left(y_{t+h/2}^{1},t+{\frac {h}{2}}\right)\\y_{t+h}^{3}&=y_{t}+hf\left(y_{t+h/2}^{2},t+{\frac {h}{2}}\right)\end{aligned}}}$

dove:

${\displaystyle y_{t+h/2}^{1}={\dfrac {y_{t}+y_{t+h}^{1}}{2}}\qquad y_{t+h/2}^{2}={\dfrac {y_{t}+y_{t+h}^{2}}{2}}}$

Se si pone:

${\displaystyle {\begin{aligned}k_{1}&=f(y_{t},t)\\k_{2}&=f\left(y_{t+h/2}^{1},t+{\frac {h}{2}}\right)\\k_{3}&=f\left(y_{t+h/2}^{2},t+{\frac {h}{2}}\right)\\k_{4}&=f\left(y_{t+h}^{3},t+h\right)\end{aligned}}}$

sostituendo si ha a meno di ${\textstyle O(h^{2})}$:

${\displaystyle {\begin{aligned}k_{2}&=f\left(y_{t+h/2}^{1},t+{\frac {h}{2}}\right)=f\left(y_{t}+{\frac {h}{2}}k_{1},t+{\frac {h}{2}}\right)\\&=f\left(y_{t},t\right)+{\frac {h}{2}}{\frac {d}{dt}}f\left(y_{t},t\right)\\k_{3}&=f\left(y_{t+h/2}^{2},t+{\frac {h}{2}}\right)=f\left(y_{t}+{\frac {h}{2}}f\left(y_{t}+{\frac {h}{2}}k_{1},t+{\frac {h}{2}}\right),t+{\frac {h}{2}}\right)\\&=f\left(y_{t},t\right)+{\frac {h}{2}}{\frac {d}{dt}}\left[f\left(y_{t},t\right)+{\frac {h}{2}}{\frac {d}{dt}}f\left(y_{t},t\right)\right]\\k_{4}&=f\left(y_{t+h}^{3},t+h\right)=f\left(y_{t}+hf\left(y_{t}+{\frac {h}{2}}k_{2},t+{\frac {h}{2}}\right),t+h\right)\\&=f\left(y_{t}+hf\left(y_{t}+{\frac {h}{2}}f\left(y_{t}+{\frac {h}{2}}f\left(y_{t},t\right),t+{\frac {h}{2}}\right),t+{\frac {h}{2}}\right),t+h\right)\\&=f\left(y_{t},t\right)+h{\frac {d}{dt}}\left[f\left(y_{t},t\right)+{\frac {h}{2}}{\frac {d}{dt}}\left[f\left(y_{t},t\right)+{\frac {h}{2}}{\frac {d}{dt}}f\left(y_{t},t\right)\right]\right]\end{aligned}}}$

dove:

${\displaystyle {\frac {d}{dt}}f(y_{t},t)={\frac {\partial }{\partial y}}f(y_{t},t){\dot {y}}_{t}+{\frac {\partial }{\partial t}}f(y_{t},t)=f_{y}(y_{t},t){\dot {y}}+f_{t}(y_{t},t):={\ddot {y}}_{t}}$

è la derivata totale di ${\displaystyle f}$ rispetto a ${\displaystyle t}$. Scrivendo la formula generale con quanto ricavato:

${\displaystyle {\begin{aligned}y_{t+h}&=y_{t}+h\left\lbrace a\cdot f(y_{t},t)+b\cdot \left[f\left(y_{t},t\right)+{\frac {h}{2}}{\frac {d}{dt}}f\left(y_{t},t\right)\right]\right.+\\&{}+c\cdot \left[f\left(y_{t},t\right)+{\frac {h}{2}}{\frac {d}{dt}}\left[f\left(y_{t},t\right)+{\frac {h}{2}}{\frac {d}{dt}}f\left(y_{t},t\right)\right]\right]+\\&{}+d\cdot \left[f\left(y_{t},t\right)+h{\frac {d}{dt}}\left[f\left(y_{t},t\right)+{\frac {h}{2}}{\frac {d}{dt}}\left[f\left(y_{t},t\right)+\left.{\frac {h}{2}}{\frac {d}{dt}}f\left(y_{t},t\right)\right]\right]\right]\right\rbrace +{\mathcal {O}}(h^{5})\\&=y_{t}+a\cdot hf_{t}+b\cdot hf_{t}+b\cdot {\frac {h^{2}}{2}}{\frac {df_{t}}{dt}}+c\cdot hf_{t}+c\cdot {\frac {h^{2}}{2}}{\frac {df_{t}}{dt}}+\\&{}+c\cdot {\frac {h^{3}}{4}}{\frac {d^{2}f_{t}}{dt^{2}}}+d\cdot hf_{t}+d\cdot h^{2}{\frac {df_{t}}{dt}}+d\cdot {\frac {h^{3}}{2}}{\frac {d^{2}f_{t}}{dt^{2}}}+d\cdot {\frac {h^{4}}{4}}{\frac {d^{3}f_{t}}{dt^{3}}}+{\mathcal {O}}(h^{5})\end{aligned}}}$

e confrontando con l'espansione di serie di Taylor di ${\displaystyle y_{t+h}}$ intorno a ${\displaystyle t}$:

${\displaystyle {\begin{aligned}y_{t+h}&=y_{t}+h{\dot {y}}_{t}+{\frac {h^{2}}{2}}{\ddot {y}}_{t}+{\frac {h^{3}}{6}}y_{t}^{(3)}+{\frac {h^{4}}{24}}y_{t}^{(4)}+{\mathcal {O}}(h^{5})=\\&=y_{t}+hf(y_{t},t)+{\frac {h^{2}}{2}}{\frac {d}{dt}}f(y_{t},t)+{\frac {h^{3}}{6}}{\frac {d^{2}}{dt^{2}}}f(y_{t},t)+{\frac {h^{4}}{24}}{\frac {d^{3}}{dt^{3}}}f(y_{t},t)\end{aligned}}}$

si ottiene un sistema di equazioni per i coefficienti:

${\displaystyle {\begin{cases}&a+b+c+d=1\\[6pt]&{\frac {1}{2}}b+{\frac {1}{2}}c+d={\frac {1}{2}}\\[6pt]&{\frac {1}{4}}c+{\frac {1}{2}}d={\frac {1}{6}}\\[6pt]&{\frac {1}{4}}d={\frac {1}{24}}\end{cases}}}$

che fornisce:

${\displaystyle a={\frac {1}{6}}\quad b={\frac {1}{3}}\quad c={\frac {1}{3}}\quad d={\frac {1}{6}}}$

• (DE) C. Runge e H. Koenig Vorlesungen über numerisches Rechnen Springer - 1924 (capitolo 10)
• (EN) D. M. Yound e R. T. Gregory A survey of numerical methods (Dover, New York) ISBN 0486656926
• (EN) Kendall E. Atkinson. An Introduction to Numerical Analysis. John Wiley & Sons - 1989
• (EN) Ward Cheney e David Kincaid Numerical Mathematics and Computing Brooks/Cole - 1985 programmi Fortran
• (EN) David Kincaid e Ward Cheney Numerical Analysis Mathetatics of Scientific Computing Brooks/Cole -1991 programmi Fortran
• (EN) R.W. Brankin, I. Gladwell, e L.F. Shampine metodi di Runge Kutta in Fortran 77

• Metodi di soluzione numerica per equazioni differenziali ordinarie
• Metodo di Eulero
• Problema di Cauchy
• Integrazione numerica

• Wikimedia Commons contiene immagini o altri file su metodi di Runge-Kutta

• (EN) AA. VV., Runge-Kutta method, in Encyclopaedia of Mathematics, Springer e European Mathematical Society, 2002.
• (EN) Runge-Kutta 4th-order method textbook notes, PPT, Matlab Mathematica Maple Mathcad at Holistic Numerical Methods Institute