Méthodes de Runge-Kutta

La courbe bleue est la solution exacte de l'équation différentielle. Les flèches rouges symbolisent les pentes utilisées impliquées dans une méthode de Runge-Kutta. La solution approchée est représentée en vert.

Les méthodes de Runge-Kutta sont des méthodes d'analyse numérique d'approximation de solutions de problèmes de Cauchy. 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.

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. On considère le problème de Cauchy suivant:

${\frac {\mathrm {d} y}{\mathrm {d} t}}=f(t,y),\quad y(t_{0})=y_{0}$

où la fonction f est connue, et y est une fonction de t inconnue. On suppose que la solution $y$ a été approchée à $t_{n}$ par $y_{n}$. Si $n=0$, ces valeurs sont simplement les valeurs initiales. On cherche à approcher la solution à $t_{n+1}=t_{n}+h$, où $h$ est un réel appelé le pas de temps.

Une première approche, celle historique, est d'utiliser un développement à l'ordre 1:

$$y(t_{n+1})\approx y(t_{n})+h{\frac {dy}{dt}}(t_{n})=y(t_{n})+hf(t_{n},y(t_{n}))$$

Ce qui mène à définir la méthode appelée méthode d'Euler explicite :

$$y_{n+1}=y_{n}+hf(t_{n},y_{n})$$

Cette méthode est le cas le plus simple de méthode de Runge-Kutta. Une version modifiée de ce schéma est la méthode d'Euler implicite :

$$y_{n+1}=y_{n}+hf(t_{n},y_{n+1})$$

La seule différence est que le terme de droite utilise $y_{n+1}$ au lieu de $y_{n}$. Cela implique qu'un système doit être résolu pour calculer $y_{n+1}$. C'est ce qui différencie les méthodes implicites qui ont besoin de résoudre un système d'équation des méthodes explicites qui ne consistent qu'en une séquence de calculs.

Plus généralement, une méthode de Runge-Kutta approche la dérivée de $y$ à plusieurs points, puis utilise ces approximations pour calculer $y_{n+1}$.

Un exemple plus complexe de méthode explicite est "RK4", la méthode de Runge-Kutta probablement la plus célèbre, qui est définie par :

$${\begin{aligned}k_{1}&=\ f(t_{n},y_{n}),\\k_{2}&=\ f\!\left(t_{n}+{\frac {h}{2}},y_{n}+h{\frac {k_{1}}{2}}\right),\\k_{3}&=\ f\!\left(t_{n}+{\frac {h}{2}},y_{n}+h{\frac {k_{2}}{2}}\right),\\k_{4}&=\ f\!\left(t_{n}+h,y_{n}+hk_{3}\right),\\y_{n+1}&=y_{n}+{\frac {h}{6}}\left(k_{1}+2k_{2}+2k_{3}+k_{4}\right).\\\end{aligned}}$$

Cette méthode est dites à 4 étapes car elle approche 4 fois la dérivée de $y$, ce qui est plus coûteux que la méthode d'Euler explicite, qui est à une étape. Bien que la méthode implicite d'Euler soit à une étape, le fait de devoir résoudre un système d'équation peut la rendre plus coûteuse à utiliser que RK4.

Principe général

On considère de nouveau un problème de Cauchy:

${\frac {\mathrm {d} y}{\mathrm {d} t}}=f(t,y),\quad y(t_{0})=y_{0}$

On suppose que la solution $y$ a précédemment été approchée à $t_{n}$ par une valeur $y_{n}$. Si $n=0$, ces valeurs sont simplement les valeurs initiales. On cherche à approcher la solution en $t_{n+1}=t_{n}+h$, où $h$ est un réel appelé le pas de temps.

En intégrant l'équation différentielle, on montre que la relation suivante est vérifiée :

$y(t_{n+1})=y(t_{n})+\int _{t_{n}}^{t_{n}+h}\displaystyle f(t,y(t))\,\mathrm {d} t\ {\stackrel {x\leftarrow {\frac {t-t_{n}}{h}}}{=}}\ y(t_{n})+h\int _{0}^{1}f(t_{n}+xh,y(t_{n}+xh))\,\mathrm {d} x$

Dans l'égalité précédente, l'intégrale n'est pas calculable car la fonction $t\mapsto y(t)$ n'est a priori pas connue, cette égalité est en réalité une équation. L'idée est de remplacer cette intégrale par une valeur approchée. On approche cette intégrale grâce à une méthode de quadrature à $s$ points qui a pour poids $b_{i}$ et pour nœuds $c_{i}$ :

$y(t_{n+1})\approx y_{n}+h\sum _{i=1}^{s}b_{i}f(t_{n}+c_{i}h,y(t_{n}+c_{i}h))$

On approche de manière similaire tous les $y(t_{n}+c_{i}h)$ par des valeurs ${\tilde {y}}_{i}$ définies par :

${\tilde {y}}_{i}=y_{n}+h\sum _{j=1}^{s}a_{ij}f(t_{n}+c_{j}h,{\tilde {y}}_{j})$

L'approximation d'une méthode de Runge-Kutta à $s$ étapes est donc donnée par : $$\forall i=1,\dotsc ,s,\,{\tilde {y}}_{i}=y_{n}+h\sum _{j=1}^{s}a_{ij}f(t_{n}+c_{j}h,{\tilde {y}}_{j})$$ $$y_{n+1}=y_{n}+h\sum _{j=1}^{s}b_{j}f(t_{n}+c_{j}h,{\tilde {y}}_{j})$$

Le fait que les ${\tilde {y}}_{i}$ utilisent les mêmes valeurs que l'approximation finale pour leur quadrature n'est pas une perte de généralité, n'importe quelle valeur utilisée pour l'une des quadratures peut être ajoutée aux autres quadratures en ajoutant un poids de valeur 0.

Tableau de Butcher

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

	$c_{1}$	$a_{1,1}$	$a_{1,2}$	$\cdots$	$a_{1,s}$
	$c_{2}$	$a_{2,1}$	$a_{2,2}$	$\cdots$	$a_{2,s}$
	$\vdots$	$\vdots$	$\vdots$	$\ddots$	$\vdots$
	$c_{s}$	$a_{s,1}$	$a_{s,2}$	$\cdots$	$a_{s,s}$
		$b_{1}$	$b_{2}$	$\cdots$	$b_{s}$

Une méthode est convergente si et seulement si $\sum _{j=1}^{s}b_{j}=1$. Dans ce cas, la précision de la méthode est améliorée en diminuant le pas.

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

• Si $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 $y_{n+1}$ ne fait intervenir que les approximations $y_{k,i}$ avec $i\leq k-1$ et $y_{n}$. Aucune équation n'est à résoudre pour calculer l'itérée suivante.
• Si $A$ n'est pas triangulaire inférieure stricte alors la méthode nécessite de résoudre des équations pour calculer $y_{n+1}$à partir de $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 désigne l'erreur réalisée par une unique itération de la méthode de Runge-Kutta. On la note $e_{n}$ et on a: $e_{n}=y(t_{n+1})-y_{n+1}$ avec $0\leq n\leq N$ et sous l'hypothèse $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

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

avec $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 $h=T/N$) est consistante d'ordre $k$ si $e_{n}={\mathcal {O}}(h^{k+1})$. Dans ce cas, on a

${\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 $y_{n+1}=y_{n}+hf\left(t_{n},y_{n}\right)$ avec un pas constant $h$, on a:

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

En notant que $y'(t_{n})=f(t_{n},y(t_{n}))$et grâce aux développement de Taylor, on montre que $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é

À cause de l'erreur effectuée par la méthode numérique, la valeur $y_{n}$est entachée d'une erreur notée $\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.

Supposons que f est k-lipschitzienne en y. Soit $\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

$\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 $n\in \mathbb {N}$, on a $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 :

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

La méthode RK1 utilise l'équation

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

	$0$	$0$
		$1$

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

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

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

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

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

	$0$	$0$	$0$
	$\alpha$	$\alpha$	$0$
		$1-{\tfrac {1}{2\alpha }}$	${\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-Nicolson

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

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

	$0$	$0$	$0$
	$1$	${\tfrac {1}{2}}$	${\tfrac {1}{2}}$
		${\tfrac {1}{2}}$	${\tfrac {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 :

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

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

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

où

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

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

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

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

	$0$	$0$	$0$	$0$	$0$
	${\tfrac {1}{2}}$	${\tfrac {1}{2}}$	$0$	$0$	$0$
	${\tfrac {1}{2}}$	$0$	${\tfrac {1}{2}}$	$0$	$0$
	$1$	$0$	$0$	$1$	$0$
		${\tfrac {1}{6}}$	${\tfrac {1}{3}}$	${\tfrac {1}{3}}$	${\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.

Solutions numériques de $y'=(\sin(t))^{2}y$ (avec $y(0)=1$) obtenue par différentes méthodes de Runge-Kutta et comparées à la solution exacte. Le pas est $h=1/3$. On note que l'ordre de précision de la méthode est lié à l'erreur : la méthode d'Euler Explicite a une erreur nettement plus importante que celle de RK4.

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 :

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

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

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

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

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

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

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

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

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

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