Metodi di Runge-Kutta

Da Wikipedia, l'enciclopedia libera.

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.

Introduzione[modifica | modifica wikitesto]

I metodi di Runge-Kutta (spesso abbreviati con "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_j), e per ottenere una ricostruzione abbastanza fedele della funzione il numero n deve essere sufficientemente elevato.

Formulazione[modifica | modifica wikitesto]

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) =  h\int_0^{\theta_i} f(t_0+v h,y(t_0 + v h))\,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[modifica | modifica wikitesto]

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{align}
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{align}

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{align}
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{align}

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} (metodo di Eulero)
  • k_2 è l'incremento basato sulla pendenza alla metà dell'intervallo, utilizzando  \dot{y} + \tfrac{h}{2}k_1
  • k_3 è un altro incremento basato sulla pendenza alla metà dell'intervallo, utilizzando  \dot{y} + \tfrac{h}{2}k_2
  • k_4 è l'incremento basato sulla pendenza alla fine dell'intervallo, utilizzando  \dot{y} + hk_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 h b_i k_i

dove:

 k_1 = f(t_n, y_n)
 k_2 = f(t_n+c_2h, y_n+h(a_{21}k_1))
 k_3 = f(t_n+c_3h, y_n+h(a_{31}k_1+a_{32}k_2))
 \vdots
 k_s = f(t_n+c_sh, y_n+h(a_{s1}k_1+a_{s2}k_2+\cdots+a_{s,s-1}k_{s-1}))

Per specificare un particolare metodo si deve definire il numero s ed i coefficienti a_{ij} (per i \le j < i \le 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[modifica | modifica wikitesto]

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.[1] Si sceglie:


\begin{array}{|l|l|}
\hline
\alpha_i & \beta_{ij} \\[8pt]
\hline
\alpha_1 = 0 & \beta_{21} = \frac{1}{2} \\[8pt]
\alpha_2 = \frac{1}{2} & \beta_{32} = \frac{1}{2} \\[8pt]
\alpha_3 = \frac{1}{2} & \beta_{43} = 1 \\[8pt]
\alpha_4 = 1 & \\[8pt]
\hline
\end{array}

con \beta_{ij} = 0 altrimenti. Si definiscono quindi le quantità:

\begin{align}
y^1_{t+h} &= y_t + hf\left(y_t, t\right) \\
y^2_{t+h} &= y_t + hf\left(y^1_{t+h/2}, t+\frac{h}{2}\right) \\
y^3_{t+h} &= y_t + hf\left(y^2_{t+h/2}, t+\frac{h}{2}\right)
\end{align}

dove:

y^1_{t+h/2} = \dfrac{y_t + y^1_{t+h}}{2} \qquad y^2_{t+h/2} = \dfrac{y_t + y^2_{t+h}}{2}

Se si pone:

\begin{align}
k_1 &= f(y_t, t) \\
k_2 &= f\left(y^1_{t+h/2}, t + \frac{h}{2}\right) \\
k_3 &= f\left(y^2_{t+h/2}, t + \frac{h}{2}\right) \\
k_4 &= f\left(y^3_{t+h}, t + h\right)
\end{align}

sostituendo si ha:

\begin{align}
k_2 &= f\left(y^1_{t+h/2}, 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^2_{t+h/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^3_{t+h}, t + h\right) = f\left(y_t + h f\left(y_t + \frac{h}{2} k_2, t + \frac{h}{2}\right), t + h\right) \\
&= f\left(y_t + h f\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{align}

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{align}
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 h f_t + b \cdot h f_t + b \cdot \frac{h^2}{2} \frac{df_t}{dt}   + c \cdot h f_t+ c \cdot \frac{h^2}{2} \frac{df_t}{dt} + \\
&{}+ c \cdot \frac{h^3}{4} \frac{d^2f_t}{dt^2} + d \cdot h f_t + d \cdot h^2 \frac{df_t}{dt} + d \cdot \frac{h^3}{2} \frac{d^2f_t}{dt^2} + d \cdot \frac{h^4}{4} \frac{d^3f_t}{dt^3} + \mathcal{O}(h^5)
\end{align}

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

\begin{align}
 y_{t+h} &= y_t + h \dot y_t + \frac{h^2}{2} \ddot y_t + \frac{h^3}{6} y^{(3)}_t + \frac{h^4}{24} y^{(4)}_t + \mathcal{O}(h^5) = \\
&= y_t + h f(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{align}

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}

Note[modifica | modifica wikitesto]

  1. ^ Ling-Hsiao Lyu - Higher-Order Numerical Integrations

Bibliografia[modifica | modifica wikitesto]

  • (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[modifica | modifica wikitesto]

Collegamenti esterni[modifica | modifica wikitesto]

matematica Portale Matematica: accedi alle voci di Wikipedia che trattano di matematica