MENUCours D'introduction à l'analyse numérique

Les techniques de Runge-Kutta sont des schémas numériques à un pas qui permettent de résoudre les équations différentielles ordinaires. Elles font parties des méthodes les plus populaires de part leur facilité de mise en œuvre et leur précision. C'est Carle Runge et Martin Kutta qui, au début du XX^e siècle, ont inventé ces méthodes.
Nous décrivons ici deux algorithmes assez utilisés : celles de Runge-Kutta d'ordre 2 et 4.

Introduction

Dans de nombreux cas, les systèmes d'équations différentielles que l'on rencontre en science peuvent se mettre sous la forme d'une équation différentielle ordinaire du premier ordre du type : \[ \left\{\begin{array}{lccc} \dfrac{\mathrm{d}y}{\mathrm{d}t} & = & f(t,y(t)), & 0\leq t \leq T \\[1mm] y(0) & = & y_{0} & \\ \end{array}\right. \] où \(y(t)\) est la fonction que l'on recherche, \(y_0\) sa valeur initiale et \(f\) une fonction connue suffisamment régulière pour que l'existence et l'unicité de la solution ne pose pas de problème. Notez que \(y(t)\) peut être un scalaire ou un vecteur ( cf. exemple d'une équa-diff d'ordre deux).

À l'instar de la méthode d'Euler, celles de Runge-Kutta sont des schémas numériques à un pas basés sur la discrétisationde la variable \(t\). On note \(h\) ce pas et \(y_{n}\) la valeur approchée de \(y(t_{n})\) pour les différents instants \(t_{n}=nh\).

En intégrant l'équation différentielle entre \(t_n\) et \(t_{n+1}\) on a la relation \begin{equation} y(t_{n+1})-y(t_n)=\int_{t_n}^{t_{n}+h} f(t,y(t))\, \mathrm{d}t \end{equation} L'idée consiste à approcher cette intégrale de façon plus précise que ne le fait la méthode d'Euler. Mais avant de voir comment, revenons sur la méthode d'Euler.

Retour sur Euler

Euler explicite

L'intégrale (1) peut s'approcher par la méthode du rectangle à gauche : \[ \int_{t_{n}}^{t_{n+1}}f(t,y(t))\, \mathrm{d}t \simeq h\times f(t_{n},y(t_{n})) \] D'où le schéma itératif suivant \[ \left\{ \begin{array}{ccc} y_{n+1} & = & y_{n}+hf(t_{n},y_{n})\\[2mm] y_{0} & = &y(0) \\ \end{array} \right. \] Il s'agit de la fameuse méthode d'Euler explicite. Graphiquement, on voit immédiatement que cette méthode sous-estime(sur-estime) l'aire quand la fonction \(f\) croît(décroît) au cours du temps.

L'erreur produite correspond à l'aire grisée de forme quasi triangulaire et de dimension \(h\times ph\) où \(p\) est la pente de \(f\) à l'instant \(t_n\). L'erreur vaut donc à peu près \begin{equation} e_E\simeq \frac12 p\,h^2 \end{equation} Après \(N\) itérations, on commet une erreur globale de l'ordre de \(N \frac12 p\,h^2=\frac12 Tp\,h\) où \(T\) est la durée totale. Aussi, pour une durée donnée, l'erreur globale augmente linéairement avec le pas \(h\) : on dit que la méthode d'Euler est d'ordre un.

Euler implicite

On aurait également pu approcher l'intégrale (1) par la méthode du rectangle à droite : \[ \int_{t_{n}}^{t_{n+1}}f(t,y(t))\, \mathrm{d}t \simeq h\times f(t_{n+1},y(t_{n+1})) \] D'où le schéma numérique, dit schéma d'Euler implicite : \begin{equation} \left\{ \begin{array}{ccc} y_{n+1} & = & y_{n}+hf(t_{n},y_{n+1})\\[2mm] y_{0} & = &y(0) \\[3mm] \end{array} \right. \end{equation}

Notez que dans ce cas, \(y_{n+1}\) est présent dans le terme de gauche et celui de droite. Contrairement à la méthode d'Euler explicite, la grandeur recherchée \(y_{n+1}\) est reliée à une fonction qui dépend de cette même grandeur. Autrement dit, \(y_{n+1}\) est défini implicitement d'où le nom de la méthode. Il faut donc résoudre l'équation implicite à chaque étape ce qui demande l'emploi d'un algorithme annexe (dichotomie, algorithme de Newton-Raphson, etc.).

Cette méthode, comme on le voit graphiquement, produit une erreur opposée à celle de la méthode d'Euler explicite ; elle est donc aussi d'ordre un. De surcroît, elle apporte une complication dans la résolution de l'équation (3). Cependant, cet inconvénient est compensé par le fait, qu'en général, les méthodes implicites sont plus stables que les méthodes explicites.

Les améliorations de Runge-Kutta

Runge-Kutta d'ordre 2

On voit immédiatement que l'on peut améliorer l'estimation de l'intégrale en calculant l'aire d'un trapèze au lieu de celui d'un rectangle. La méthode du trapèze consiste en l'approximation suivante : \[ \int_a^b f(x)\, \mathrm{d}x\simeq \frac{b-a}{2}\left[f(a)+f(b)\right] \] Appliquée à l'intégrale (1), cela donne \[ \int_{t_{n}}^{t_{n+1}}f(t,y(t))\, \mathrm{d}t \simeq \frac{h}{2}\left[(f(t_{n},y(t_{n}))+f(t_{n+1},y(t_{n+1}))\right] \] Ici, l'intégrale dépend des valeurs de \(y_n\) et \(y_{n+1}\) ce qui, si on en restait là, donnerait lieu à une méthode implicite. Pour éviter ces complications, on utilise la méthode d'Euler (explicite) afin d'estimer la valeur \(y_{n+1}\) qui intervient dans \(f(t_{n+1},y(t_{n+1}))\). On obtient le schéma itératif suivant : \begin{equation} \fcolorbox{#FF9D00}{#FEF5EB}{\quad \(\displaystyle y_{n+1}=y_{n}+h\left(\frac{k_1}{2}+\frac{k_2}{2}\right) \quad\text{avec}\quad \left\lbrace \begin{array}{rcl} k_1 &=& f(t_{n},y_{n})\\ k_2 &=& f(t_n+h,y_n+hk_1)\\ y_{0} &=& y(0)\\ \end{array} \right. \notag \)\quad} \end{equation} Cette méthode présente l'avantage d'être précise et assez simple à programmer (voir algorithme ci-dessous).

Ce schéma numérique présente deux erreurs de troncature. Tout d'abord, l'approximation de l'intégrale par l'aire d'un trapèze produit une première erreur \(e_1\). On voit bien graphiquement que la méthode du trapèze neutralise l'erreur des méthodes d'Euler : l'erreur \(e_1\) est liée à la courbure de la fonction et non à sa pente c'est pourquoi \(e_1\) varie comme \(h^3\) (au lieu de \(h^2\) pour les méthodes d'Euler). Ensuite, on commet une deuxième erreur \(e_2\) lorsque l'on estime \(y_{n+1}\) par \(y_n+hf(t_n,y_n)\) dans le calcul de \(f(t_{n+1},y_{n+1})\). Dans ce cas, un calcul de propagation d'erreur mène au résultat suivant : \(e_2=\frac{h}{2}\times \frac{\partial f}{\partial y}\times e_E\) avec \(e_E\) l'erreur donnée par la formule (2). Pour finir, l'erreur produite à chaque étape s'écrit \(e_{RK2}=e_1+e_2=\mathrm{C^{te}}h^3\) de sorte que l'erreur globale, pour une durée \(T\) fixée, se comporte comme \(h^2\). Il s'agit donc bien d'une méthode d'ordre deux comme son nom l'indique.

Runge Kutta d'ordre 4

La méthode de Runge-Kutta d'ordre 4 est une étape supplémentaire dans le raffinement du calcul de l'intégrale (1). Au lieu d'utiliser la méthode des trapèzes, on utilise la méthode de Simpson. Celle-ci consiste à remplacer la fonction intégrée par une parabole passant par les points extrêmes et le point milieu. On a \[ \int_a^b f(x)\, \mathrm{d}x\simeq \frac{b-a}{6}\left[f(a)+4f(\frac{a+b}{2})+f(b)\right] \] Appliquée à l'intégrale (1), cela donne \[ \int_{t_{n}}^{t_{n+1}}f(t,y(t))\, \mathrm{d}t \simeq \frac{h}{6}\left[f(t_n,y(t_n))+4f(t_{n+1/2},y(t_{n+1/2}))+f(t_{n+1},y(t_{n+1}))\right] \] d'où la relation \[ y_{n+1}=y_n+\frac{h}{6}\left[f(t_n,y_n)+4f(t_{n+1/2},y_{n+1/2})+f(t_{n+1},y_{n+1})\right] \] Ici, une difficulté apparaît car l'équation présente deux inconnues : \(y_{n+1/2}\) et \(y_{n+1}\). Pour rendre le schéma explicite, il faut estimer \(4f(t_{n+1/2},y_{n+1/2})\) et \(f(t_{n+1},y_{n+1})\) à partir de \(y_n\), \(t_n\) et \(h\).

Commençons par le terme \(4f(t_{n+1/2},y(t_{n+1/2})\) que l'on décompose en deux termes identiques \[ 2f(t_{n+1/2},y_{n+1/2})+2f(t_{n+1/2},y_{n+1/2}) \] Dans le premier, on remplace \(y_{n+1/2}\) par sa valeur déduite de la méthode d'Euler explicite ; à savoir \(y^a_{n+1/2}=y_n+h/2f(t_n,y_n)\). Dans le deuxième terme, on remplace \(y_{n+1/2}\) par sa valeur déduite de la méthode d'Euler implicite : \(y^b_{n+1/2}= y_n+h/2f(t_{n+1/2},y_{n+1/2})\) que l'on va approcher par \(y_n+h/2f(t_{n+1/2},y^a_{n+1/2})\). Les méthodes d'Euler implicite et explicite produisant des erreurs quasi opposées, on a ainsi l'espoir de minimiser l'erreur sur le calcul de \(4f(t_{n+1/2},y_{n+1/2})\). Pour résumer, on écrira \[ 4f(t_{n+1/2},y_{n+1/2})\simeq 2k_2+2k_3 \quad\text{avec}\quad \left\{\begin{array}{rcl} k_1&=&f(t_n,y_n)\\ k_2&=&f(t_n+\frac12 h,y_n+\frac 12 hk_1)\\ k_{3} & = & f(t_{n}+\frac12 h,y_{n}+\frac12 hk_{2}) \\ \end{array}\right. \]

Quant au terme \(f(t_{n+1},y_{n+1})\), on l'approche en estimant \(y_{n+1}\) par la méthode du point milieu, c'est-à-dire en appliquant la méthode du rectangle au milieu : \[ y_{n+1}\simeq y_n+hf(t_{n+1/2},y_{n+1/2})\simeq y_n+hf(t_{n+1/2},y^b_{n+1/2}) \] Finalement on obtient le schéma explicite, dit de Runge-Kutta d'ordre 4 : \begin{equation} \fcolorbox{#FF9D00}{#FEF5EB}{\quad \(\displaystyle y_{n+1}=y_n+h\left[\frac{k_1}{6}+\frac{k_2}{3}+\frac{k_3}{3}+\frac{k_4}{6}\right] \quad\text{avec}\quad \left\lbrace\begin{array}{lcl} k_{1} & = & f(t_{n},y_{n})\\[1mm] k_{2} & = & f(t_{n}+\frac12 h,y_{n}+\frac12 hk_{1}) \\[1mm] k_{3} & = & f(t_{n}+\frac12 h,y_{n}+\frac12 hk_{2}) \\[1mm] k_{4} & = & f(t_n+h,y_n+hk_{3}) \end{array}\right. \)\quad} \end{equation}

Par rapport à la méthode RK2, ce schéma numérique exige deux fois plus de calcul à chaque pas et donc un temps de calcul plus long, sans parler des erreurs d'arrondi qui s'accumulent plus vite. Cependant, ce défaut est compensé par un gain de précision. En effet, comme son nom l'indique, il s'agit d'une méthode d'ordre 4.

On peut le vérifier sur un exemple simple. Prenons l'équation \[ \frac{\mathrm{d}y}{\mathrm{d}t}=1-y \quad\text{avec}\quad y(0)=0 \] dont la solution est la fonction exponentielle \(y(t)=1-\mathrm{e}^{-t}\). Quelle valeur prévoit la méthode de Runge-Kutta-4 pour \(y_1=y(h)\) ? Sachant que l'on a ici \(f(t,y(t))=1-y(t)\), le schéma numérique (4) donne

d'où il vient \[ \begin{split} y_{1}=y_0+h\left[\frac16 k_1+\frac13 k_2+\frac13 k_3+\frac16 k_4\right]\\ =h-\frac12 h^2+\frac16 h^3-\frac{1}{24}h^4 \end{split} \]

qui n'est rien d'autre que le développement limité de la solution exacte \(y(h)\) à l'ordre 4. Ainsi, même avec un pas modeste, la précision reste bonne comme l'illustre la figure.

Terminons par l'algorithme complet :

Conclusion

Les techniques de Runge-Kutta, d'ordre 2 ou 4, ont l'avantage d'être simples à mettre en œuvre, précises et assez stables pour les fonctions courantes rencontrées en physique. C'est ce qui explique leur grande popularité. De nombreux logiciels de calcul utilisent par défaut la méthode RK4 dans sa version adaptative Le pas h est alors variable et s'adapte pour optimiser précision et temps de calcul..

Bien entendu ces méthodes ont aussi leurs défauts : elles sont assez gourmandes en temps de calcul et ne sont pas adaptés aux systèmes conservatifs aux temps longs (voir Méthode de Verlet).