Méthodes de Runge-Kutta

Les méthodes de Runge-Kutta sont des méthodes d'analyse numérique d'approximation de solutions d'équations différentielles. Elles ont été nommées ainsi en l'honneur des mathématiciens Carl Runge et Martin Wilhelm Kutta, lesquels élaborèrent la méthode en 1901.

Ces méthodes reposent sur le principe de l'itération, c'est-à-dire qu'une première estimation de la solution est utilisée pour calculer une seconde estimation, plus précise, et ainsi de suite.

Premiers exemples

Pour donner l'idée générale des méthodes de Runge-Kutta, on commence par donner une idée de la méthode sur quelques exemples de méthodes de Runge-Kutta avant de considérer le cas général. Soit le problème suivant :

${\displaystyle {dy \over dt}=f(t,y),\quad y(t_{0})=y_{0}}$

La fonction f est connue. y est une fonction de t inconnue. On cherche les valeurs de y lorsque t varie dans un ensemble discret t₀ < t₁ < ... < t_N.

En intégrant ce système, on montre que ${\displaystyle t\mapsto y(t)}$ vérifie les relations suivantes :

${\displaystyle {\begin{aligned}y(t_{n+1})&=y(t_{n})+\int _{t_{n}}^{t_{n}+h_{n}}\displaystyle f(t,y(t))\mathrm {d} t\\&=y(t_{n})+h_{n}\int _{0}^{1}f(t_{n}+h_{n}u,y(t_{n}+h_{n}u))\mathrm {d} u\end{aligned}}}$

où ${\displaystyle h_{n}=t_{n+1}-t_{n}>0}$ va désigner le pas entre deux approximations. Dans l'égalité précédente, l'intégrale n'est pas calculable car la fonction ${\displaystyle t\mapsto y(t)}$ n'est a priori pas connue. L'idée est de remplacer cette intégrale par un calcul approché.

En considérant une interpolation linéaire utilisant les points ${\displaystyle t_{n}}$ et ${\displaystyle t_{n+1}}$, on déduit grâce à la méthode des trapèzes :

${\displaystyle y(t_{n+1})\approx y(t_{n})+h_{n}{\frac {f(t_{n+1},y(t_{n+1}))+f(t_{n},y(t_{n}))}{2}}}$

Ainsi, en notant ${\displaystyle y_{i}}$ les valeurs approchant ${\displaystyle y(t_{i})}$, on déduit la relation suivante :

${\displaystyle y_{n+1}=y_{n}+{\frac {h_{n}}{2}}\left[f(t_{n+1},y_{n+1})+f(t_{n},y_{n})\right]}$.

En l'état, la méthode est appelée méthode de Crank-Nicolson. L'inconvénient est que pour calculer ${\displaystyle y_{n+1}}$ à partir de ${\displaystyle y_{n}}$, une équation non-linéaire doit être résolue à chaque fois : la méthode est implicite. Une alternative est d'utiliser un premier calcul approché pour estimer le terme ${\displaystyle y_{n+1}}$, noté ${\displaystyle {\tilde {y}}_{n+1}}$, dans le membre de droite puis d'avancer dans le calcul. C'est ce qui est fait dans le schéma de Heun :

${\displaystyle {\begin{aligned}{\tilde {y}}_{n+1}&=y_{n}+h_{n}f(t_{n},y_{n})\\y_{n+1}&=y_{n}+{\frac {h_{n}}{2}}\left(f(t_{n+1},{\tilde {y}}_{n+1})+f(t_{n},y_{n})\right)\end{aligned}}}$

Ainsi, ${\displaystyle y_{n+1}}$ est calculé immédiatement à partir de ${\displaystyle y_{n}}$, la méthode est explicite.

Dans les deux cas, il s'agit ici de l'idée des méthodes de Runge-Kutta : calculer ${\displaystyle y_{n+1}}$à partir de ${\displaystyle y_{n}}$ en passant par des approximations intermédiaires à des temps compris dans ${\displaystyle [t_{n};t_{n+1}]}$. La séquence générée ${\displaystyle (y_{i})_{0\leq i\leq N}}$ contient les approximations de ${\displaystyle (y(t_{i}))_{0\leq i\leq N}}$.

Principe général

On revient au cas général, à savoir la résolution du problème :

${\displaystyle {dy \over dt}=f(t,y),\quad y(t_{0})=y_{0}}$

que l'on va chercher à résoudre en un ensemble discret t₀ < t₁ < ... < t_N. Comme une solution explicite n'est pas toujours disponible, les méthodes de Runge-Kutta permettent de déterminer une solution approchée aux points ${\displaystyle (t_{i})_{0\leq i\leq N}}$. Dans la suite, on notera ${\displaystyle y_{i}}$ la solution approchée au temps ${\displaystyle t_{i}}$ et ${\displaystyle y(t_{i})}$ la valeur exacte.

Il faut donc introduire les points intermédiaires ${\displaystyle \{(t_{n,i},y_{n,i})\}_{1\leqslant i\leqslant q}}$ . Il s'agit de solutions approchées intermédiaires utilisées pour calculer par récurrence les valeurs (t_n , y_n) avec

${\displaystyle t_{n,i}=t_{n}+c_{i}\cdot h_{n}}$

où ${\displaystyle h_{n}=t_{n+1}-t_{n}}$ est le pas de temps et c_i est dans l'intervalle [0 ; 1].

Pour chaque point intermédiaire, on note la pente correspondante

${\displaystyle p_{n,i}=f(t_{n,i},y_{n,i}).}$

Ainsi, pour une solution exacte ${\displaystyle t\mapsto y(t)}$ du problème, on a

${\displaystyle y(t_{n,i})=y(t_{n})+\int _{t_{n}}^{t_{n,i}}f(t,y(t))\mathrm {d} t=y(t_{n})+h_{n}\int _{0}^{c_{i}}f(t_{n}+uh_{n},y(t_{n}+uh_{n}))\mathrm {d} u,\,\forall i=1,\dotsc ,q,}$

${\displaystyle y(t_{n+1})=y(t_{n})+\int _{t_{n}}^{t_{n+1}}f(t,y(t))\mathrm {d} t=y(t_{n})+h_{n}\int _{0}^{1}f(t_{n}+uh_{n},y(t_{n}+uh_{n}))\mathrm {d} u.}$

On calculera ces intégrales par une méthode de quadrature, qu'on peut choisir différentes pour deux valeurs distinctes de i:

${\displaystyle \int _{0}^{c_{i}}g(u)\mathrm {d} u\approx \sum _{k=1}^{i-1}a_{ik}g(c_{k}),\,\int _{0}^{1}g(u)\mathrm {d} u\approx \sum _{k=1}^{q}b_{k}g(c_{k}),}$

calculées ici pour g(u) = f(t_n + u h_n, y(t_n + u h_n)).

La méthode de Runge-Kutta à q étapes est donnée par :

${\displaystyle \forall i=1,\dotsc ,q,{\begin{cases}t_{n,i}&=t_{n}+c_{i}h_{n},\\y_{n,i}&=y_{n}+h_{n}\sum _{k=1}^{q}a_{ik}p_{n,k}\\p_{n,i}&=f(t_{n,i},y_{n,i})\end{cases}}}$${\displaystyle y_{n+1}=y_{n}+h_{n}\sum _{k=1}^{q}b_{k}p_{n,k}.}$

Tableau de Butcher

Les informations de la méthode sont souvent résumées par le tableau des différents poids de quadrature, appelé tableau de Butcher :

	${\displaystyle c_{1}}$	${\displaystyle a_{1,1}}$	${\displaystyle a_{1,2}}$
	${\displaystyle c_{2}}$	${\displaystyle a_{2,1}}$	${\displaystyle a_{2,2}}$	${\displaystyle \ddots }$
	${\displaystyle c_{3}}$	${\displaystyle a_{3,1}}$	${\displaystyle a_{3,2}}$	${\displaystyle \ddots }$
	${\displaystyle \vdots }$	${\displaystyle \vdots }$		${\displaystyle \ddots }$	${\displaystyle \ddots }$
	${\displaystyle c_{q}}$	${\displaystyle a_{q,1}}$	${\displaystyle a_{q,2}}$	${\displaystyle \cdots }$	${\displaystyle a_{q,q-1}}$	${\displaystyle a_{q,q}}$
		${\displaystyle b_{1}}$	${\displaystyle b_{2}}$	${\displaystyle \cdots }$	${\displaystyle b_{q-1}}$	${\displaystyle b_{q}}$

La consistance de la méthode est garantie si ${\displaystyle \forall i=2,\dotsc ,q,\sum _{k=1}^{q}a_{ik}=c_{i}}$. Dans ce cas, la précision de la méthode est améliorée en diminuant le pas.

Soit la matrice ${\displaystyle A=(a_{i,j})_{1\leq i,j\leq q}\in {\mathcal {M}}_{q}(\mathbb {R} )}$. On peut établir les remarques suivantes^[1] :

• Si ${\displaystyle A}$ est triangulaire inférieure stricte alors la méthode est dite explicite et les calculs sont explicites. En effet, dans ce cas, le calcul de ${\displaystyle y_{n+1}}$ ne fait intervenir que les approximations ${\displaystyle y_{k,i}}$ avec ${\displaystyle i\leq k-1}$ et ${\displaystyle y_{n}}$. Aucune équation n'est à résoudre pour calculer l'itérée suivante.
• Si ${\displaystyle A}$ n'est pas triangulaire inférieure stricte alors la méthode nécessite de résoudre des équations pour calculer ${\displaystyle y_{n+1}}$à partir de ${\displaystyle y_{n}}$. La méthode est implicite. Les méthodes implicites sont généralement plus lourdes en calculs.

Consistance

Définition

L'erreur de consistance relative^[2] désigne l'erreur réalisée par une unique itération de la méthode de Runge-Kutta. On la note ${\displaystyle e_{n}}$ et on a^[3]: ${\displaystyle e_{n}=y(t_{n+1})-y_{n+1}}$ avec ${\displaystyle 0\leq n\leq N}$ et sous l'hypothèse ${\displaystyle y_{n}=y(t_{n})}$.

L'erreur de consistance relative désigne l'erreur réalisée sur une seule itération. En pratique, on va répéter une succession d'itérations et les erreurs vont se cumuler. Pour cette raison, on dit qu'une méthode est consistante si

${\displaystyle \lim _{h_{\max }\rightarrow 0}\sum _{n=0}^{N}|e_{n}|=0}$

avec ${\displaystyle h_{\max }=\max _{0\leq n\leq N}h_{i}}$.

Plus précisément et comme une conséquence des définitions ci-dessus, on dit qu'une méthode (avec un pas constant ${\displaystyle h=T/N}$) est consistante d'ordre ${\displaystyle k}$ si ${\displaystyle e_{n}={\mathcal {O}}(h^{k+1})}$. Dans ce cas, on a

${\displaystyle {\begin{aligned}\lim _{h\rightarrow 0}\sum _{n=0}^{N}|e_{n}|&=\lim _{h\rightarrow 0}\sum _{n=0}^{N}|C_{n}h^{k+1}|&\\&=\lim _{h\rightarrow 0}\sum _{n=1}^{N}|C_{n}h^{k+1}|&{\text{ car }}e_{0}=0\\&\leq \lim _{h\rightarrow 0}N\max _{1\leq n\leq N}|C_{n}|h^{k+1}&\\&=\lim _{h\rightarrow 0}(T\max _{1\leq n\leq N}|C_{n}|)h^{k}&{\text{ car }}h=T/N\\&=\lim _{h\rightarrow 0}Ch^{k}{\text{ avec }}C=T\max _{1\leq n\leq N}|C_{n}|&\\&=0{\text{ si }}k>0.&\end{aligned}}}$

L'ordre de consistance est l'une des propriétés fondamentales d'un schéma de Runge-Kutta.

Exemple

Par exemple, on considère la méthode d'Euler ${\displaystyle y_{n+1}=y_{n}+hf\left(t_{n},y_{n}\right)}$ avec un pas constant ${\displaystyle h}$, on a:

${\displaystyle e_{n}=y(t_{n+1})-\left(y(t_{n})+hf(t_{n},y(t_{n}))\right)}$.

En notant que ${\displaystyle y'(t_{n})=f(t_{n},y(t_{n}))}$et grâce aux développement de Taylor, on montre que ${\displaystyle e_{n}={\frac {h^{2}}{2}}\max _{t_{n}\leq t\leq t_{n+1}}|f''(t,y(t))|={\mathcal {O}}(h^{2})}$. Donc la méthode est d'ordre 1 et consistante.

Stabilité

A cause de l'erreur effectuée par la méthode numérique, la valeur ${\displaystyle y_{n}}$est entachée d'une erreur notée ${\displaystyle \varepsilon _{n}}$ (erreurs d'arrondis, erreur de consistance, ...). La propagation de ces erreurs doit rester contrôlable. Cette propriété est liée à la stabilité de la méthode^[2].

Supposons que f est k-lipschitzienne en y. Soit ${\displaystyle \alpha =\max _{i=1,\dotsc ,q}\sum _{k=1}^{i}|a_{ik}|.}$ Les méthodes de Runge-Kutta sont stables sur [0,T], avec pour constante de stabilité e^ΛT, avec

${\displaystyle \Lambda =k\sum _{i=1}^{q}|b_{i}|\left(1+\alpha kh_{\max }+\ldots +(\alpha kh_{\max })^{i-1}\right).}$

Exemples

Pour les exemples ci-dessous, seul le cas d'un pas de discrétisation constant est considéré : pour tout ${\displaystyle n\in \mathbb {N} }$, on a ${\displaystyle h_{n}=h>0.}$

La méthode de Runge-Kutta explicite d'ordre 1 (RK1)

Cette méthode est équivalente à la méthode d'Euler Explicite, une méthode simple de résolution d'équations différentielles du premier degré.

Considérons le problème suivant :

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

La méthode RK1 utilise l'équation

${\displaystyle y_{n+1}=y_{n}+hf\left(t_{n},y_{n}\right)}$

où h est le pas de l'itération.

Le problème s'écrit sous forme du tableau de Butcher suivant :

	${\displaystyle 0}$	${\displaystyle 0}$
		${\displaystyle 1}$

Schéma prédicteur-correcteur explicite (RK2)

La méthode du point milieu est une composition de la méthode d'Euler :

${\displaystyle y_{n+1}=y_{n}+hf\left(t_{n}+{\frac {h}{2}},y_{n}+{\frac {h}{2}}f\left(t_{n},y_{n}\right)\right)}$

où h est le pas de l'itération.

Elle consiste à estimer la dérivée au milieu du pas d'intégration :

${\displaystyle {\tilde {y}}_{n+{\frac {1}{2}}}=y_{n}+{\frac {h}{2}}f\left(t_{n},y_{n}\right)}$

et à refaire le pas d'intégration complet à partir de cette estimation :

${\displaystyle y_{n+1}=y_{n}+hf\left(t_{n}+{\frac {h}{2}},{\tilde {y}}_{n+{\frac {1}{2}}}\right)}$

Ce schéma est couramment appelé schéma prédicteur-correcteur explicite.

C'est le cas particulier pour α = 1/2 de la méthode plus générale :

	${\displaystyle 0}$	${\displaystyle 0}$	${\displaystyle 0}$
	${\displaystyle \alpha }$	${\displaystyle \alpha }$	${\displaystyle 0}$
		${\displaystyle 1-{\tfrac {1}{2\alpha }}}$	${\displaystyle {\tfrac {1}{2\alpha }}}$

On reconnaît ainsi que la méthode de quadrature utilisée pour les temps intermédiaires est celle du point milieu.

C'est une méthode d'ordre 2 car l'erreur est de l'ordre de h³.

Une alternative à la méthode précédente est la méthode de Heun, correspondant au cas α = 1. La méthode de quadrature repose sur la méthode des trapèzes.

La méthode de Crank-Nicholson

La méthode de Crank-Nicholson est une méthode implicite déduite de la formule des trapèzes pour la quadrature. Elle s'énonce de la façon suivante :

${\displaystyle {\begin{aligned}y_{n+1}&=y_{n}+{\frac {h_{n}}{2}}\left((f(t_{n+1},y_{n+1})+f(t_{n},y_{n})\right)\end{aligned}}}$

Le tableau de Butcher associé est :

	${\displaystyle 0}$	${\displaystyle 0}$	${\displaystyle 0}$
	${\displaystyle 1}$	${\displaystyle {\frac {1}{2}}}$	${\displaystyle {\frac {1}{2}}}$
		${\displaystyle {\frac {1}{2}}}$	${\displaystyle {\frac {1}{2}}}$

C'est une méthode d'ordre 2: l'erreur est de l'ordre de h³.

La méthode de Runge-Kutta classique d'ordre quatre (RK4)

C'est un cas particulier d'usage très fréquent, noté RK4.

Considérons le problème suivant :

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

La méthode RK4 est donnée par l'équation :

${\displaystyle y_{n+1}=y_{n}+{\frac {h}{6}}\left(k_{1}+2k_{2}+2k_{3}+k_{4}\right)}$

où

${\displaystyle k_{1}=f\left(t_{n},y_{n}\right)}$

${\displaystyle k_{2}=f\left(t_{n}+{\frac {h}{2}},y_{n}+{\frac {h}{2}}k_{1}\right)}$

${\displaystyle k_{3}=f\left(t_{n}+{\frac {h}{2}},y_{n}+{\frac {h}{2}}k_{2}\right)}$

${\displaystyle k_{4}=f\left(t_{n}+h,y_{n}+hk_{3}\right)}$

L'idée est que la valeur suivante (y_n+1) est approchée par la somme de la valeur actuelle (y_n) et du produit de la taille de l'intervalle (h) par la pente estimée. La pente est obtenue par une moyenne pondérée de pentes :

• k₁ est la pente au début de l'intervalle ;
• k₂ est la pente au milieu de l'intervalle, en utilisant la pente k₁ pour calculer la valeur de y au point t_n + h/2 par le biais de la méthode d'Euler ;
• k₃ est de nouveau la pente au milieu de l'intervalle, mais obtenue cette fois en utilisant la pente k₂ pour calculer y ;
• k₄ est la pente à la fin de l'intervalle, avec la valeur de y calculée en utilisant k₃.

Dans la moyenne des quatre pentes, un poids plus grand est donné aux pentes au point milieu.

	${\displaystyle 0}$	${\displaystyle 0}$	${\displaystyle 0}$	${\displaystyle 0}$	${\displaystyle 0}$
	${\displaystyle {\tfrac {1}{2}}}$	${\displaystyle {\tfrac {1}{2}}}$	${\displaystyle 0}$	${\displaystyle 0}$	${\displaystyle 0}$
	${\displaystyle {\tfrac {1}{2}}}$	${\displaystyle 0}$	${\displaystyle {\tfrac {1}{2}}}$	${\displaystyle 0}$	${\displaystyle 0}$
	${\displaystyle 1}$	${\displaystyle 0}$	${\displaystyle 0}$	${\displaystyle 1}$	${\displaystyle 0}$
		${\displaystyle {\tfrac {1}{6}}}$	${\displaystyle {\tfrac {1}{3}}}$	${\displaystyle {\tfrac {1}{3}}}$	${\displaystyle {\tfrac {1}{6}}}$

La méthode RK4 est une méthode d'ordre 4, ce qui signifie que l'erreur commise à chaque étape est de l'ordre de h⁵, alors que l'erreur totale accumulée est de l'ordre de h⁴.

Ces formules sont aussi valables pour des fonctions à valeurs vectorielles.

La méthode de Runge-Kutta d'ordre quatre avec dérivée seconde

Dans le cas d'une équation différentielle ordinaire d'ordre 2 de la forme :

${\displaystyle {\begin{aligned}y''&=f(t,y,y'),\\y(t_{0})&=y_{0},\quad y'(t_{0})=y'_{0}\end{aligned}}}$

on peut décomposer le problème en un système d'équations :

${\displaystyle \left\{{\begin{matrix}y'=z\\z'=f(t,y,z)\end{matrix}}\right.}$

On se ramène alors à un système d'équations différentielles du premier ordre de la forme

${\displaystyle Y'=F(t,Y),\ \mathrm {avec} \ Y(t)={\begin{pmatrix}y(t)\\z(t)\end{pmatrix}},\ Y(t_{0})={\begin{pmatrix}y_{0}\\y'_{0}\end{pmatrix}}.}$

Le système peut être résolu numériquement en utilisant une méthode de Runge-Kutta. Par exemple, en appliquant la méthode RK4 à chacune de ces équations, puis en simplifiant on obtient^[4] :

${\displaystyle k_{1}=f\left(t_{n},y_{n},y'_{n}\right)}$

${\displaystyle k_{2}=f\left(t_{n}+{\frac {h}{2}},y_{n}+{\frac {h}{2}}y'_{n},y'_{n}+{\frac {h}{2}}k_{1}\right)}$

${\displaystyle k_{3}=f\left(t_{n}+{\frac {h}{2}},y_{n}+{\frac {h}{2}}y'_{n}+{\frac {h^{2}}{4}}k_{1},y'_{n}+{\frac {h}{2}}k_{2}\right)}$

${\displaystyle k_{4}=f\left(t_{n}+h,y_{n}+hy'_{n}+{\frac {h^{2}}{2}}k_{2},y'_{n}+hk_{3}\right)}$

On en déduit y_{n + 1} et y'_{n + 1} grâce à:

${\displaystyle y_{n+1}=y_{n}+hy'_{n}+{\frac {h^{2}}{6}}\left(k_{1}+k_{2}+k_{3}\right)}$

${\displaystyle y'_{n+1}=y'_{n}+{\frac {h}{6}}\left(k_{1}+2k_{2}+2k_{3}+k_{4}\right)}$

La méthode peut être étendue à des équations différentielles ordinaires d'ordre supérieur.

Références

Articles connexes

• Méthode leapfrog
• Résolution numérique des équations différentielles

Liens externes

• (en) Eric W. Weisstein, « Runge-Kutta Method », sur MathWorld
• (en) Eric W. Weisstein, « Predictor-Corrector Methods », sur MathWorld

v · m Analyse numérique
Recherche de zéro	Méthode de Josephy-Newton Méthode de la sécante Méthode de Newton Point fixe
Transformations de matrice	Matrice de Hessenberg Décomposition LU Factorisation de Cholesky
Résolutions de systèmes	Méthode de Gauss-Seidel Méthode de surrelaxation successive (SOR) Méthode de Jacobi Décomposition QR Décomposition LU
Intégration numérique	Méthode du point médian Méthode des trapèzes Méthode de Simpson Méthodes de quadrature de Gauss Formule de Newton-Cotes Méthode de Romberg Méthode de Monte-Carlo
Équations différentielles	Méthode d'Euler (et semi-implicite) Méthodes de Runge-Kutta Intégration de Verlet Leapfrog Méthode multi-pas linéaire (en) (Adams-Bashforth, backward differentiation formula (en))
Interpolation numérique	Spline Interpolation polynomiale Interpolation lagrangienne Interpolation d'Hermite Suite de polynômes orthogonaux

• Portail de l'analyse