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.

Introduzione

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 $y(t):\mathbb {R} \to \mathbb {R} ^{d}$ che verifichi il generico problema di Cauchy:

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

Formulazione

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 $b=(b_{i})_{i=0,\dots ,s}$, una matrice $a=(a_{i,j})_{i,j=0,\dots ,s}$ e un vettore $\theta =(\theta _{i})_{i=0,\dots ,s}$. L'approssimazione è data dal sistema:

${\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 $y_{i}$ approssimano il valore esatto $y(t_{i})$.

Considerando il generico problema di Cauchy:

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

Si può considerare una sua riformulazione in forma integrale:

$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 $\Delta =\{t_{i}\}=\{t_{0}+ih\}_{\{i=0,\dots ,n\}}$, il valore della soluzione nel punto $t_{1}$ è dato da:

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

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

$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 $t_{0}$ e $t_{1}$.

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

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

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

$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 $a_{ij}$; si ottiene così che (per semplicità si considerano formule con i medesimi nodi di quadratura di quelle usate in precedenza):

$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 $t_{i}$, e come conseguenza la formulazione generale per i metodi di Runge-Kutta.

Metodi espliciti

Dato il problema ai valori iniziali:

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

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

${\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 $n=1,2,3,\dots$. In questo modo $y(t_{n+1})$ è approssimato con $y_{n+1}$, e $y_{n+1}$ è determinato da $y_{n}$ più la media pesata di quattro incrementi $k_{1}$, $k_{2}$, $k_{3}$ e $k_{4}$:

${\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 $h$ e una stima della pendenza di $f$. Nello specifico:

• $k_{1}$ è l'incremento basato sulla pendenza all'inizio dell'intervallo, utilizzando ${\dot {y}}_{n}$ (metodo di Eulero)
• $k_{2}$ è l'incremento basato sulla pendenza alla metà dell'intervallo, utilizzando $k_{1}$
• $k_{3}$ è un altro incremento basato sulla pendenza alla metà dell'intervallo, utilizzando $k_{2}$
• $k_{4}$ è l'incremento basato sulla pendenza alla fine dell'intervallo, utilizzando $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 $f$ è indipendente da $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:

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

dove:

${\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 $s$ ed i coefficienti $a_{ij}$ (per $1\leq j<i\leq s$), $b_{i}$ (per $i=1,2,\dots ,s$) e $c_{i}$ (per $i=2,3,\dots ,s$). La matrice dei coefficienti $a_{ij}$ è detta matrice di Runge-Kutta. Si usa spesso rappresentare i coefficienti in una tabella:

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

e si deve avere:

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

Derivazione per il quarto ordine

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

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

dove:

$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 $y_{t}$ all'$i$-esimo ordine.

Un modo per derivare il metodo per $s=4$ utilizzando la formula generale è il seguente. Si sceglie:

${\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 $\beta _{ij}=0$ altrimenti. Si definiscono quindi le quantità:

${\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:

$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:

${\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})}$:

${\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:

${\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 $f$ rispetto a $t$. Scrivendo la formula generale con quanto ricavato:

${\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 $y_{t+h}$ intorno a $t$:

${\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:

${\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:

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

Bibliografia

• (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

Voci correlate

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

Altri progetti

Altri progetti

• Wikimedia Commons

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

Collegamenti esterni

• (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