Résolution numérique

Il existe de très nombreuses méthodes numériques de résolution des équations différentielles : plus ou moins précises, plus ou moins coûteuses en temps de calcul, plus ou moins adaptées à tel ou tel type d'équation. Toutes ces méthodes ont été implémentées et les meilleures sont proposées dans les logiciels de calcul standard. Pour comprendre la problématique de la résolution numérique, nous détaillerons seulement la méthode la plus simple, qui est la méthode d'Euler.

Nous nous plaçons dans le cadre du théorème de Cauchy-Lipschitz. Le problème à résoudre est le problème de Cauchy (10) :

$$\left\{ \begin{array}{lcl} Y'(t) &=& G(Y(t))\\ Y(0) &=& y_0\;. \end{array}\right.$$

Nous supposons que les hypothèses du théorème 7 sont vérifiées : la fonction $G$ est continûment différentiable sur $\mathbb{R}^d$ et ses dérivées partielles sont bornées. Notre problème est de construire une approximation numérique de la solution. Il faut pour cela discrétiser le temps. Nous choisissons donc un pas de discrétisation $h>0$. Les multiples entiers de $h$ sont les instants de discrétisation.

Remarque : Dans certains cas, il est judicieux d'adapter la suite des instants de discrétisation aux valeurs de la solution : plus espacés là où la solution varie peu, ils seront plus rapprochés là où la solution varie rapidement. On parle dans ce cas de méthode à pas adaptatif.

Notre but est de calculer par récurrence dans $\mathbb{R}^d$ une suite de valeurs $z_n$ qui seront comprises comme des approximations de la fonction aux instants de discrétisation. En ces instants, la solution exacte vérifie :

$$Y((n+1)h)=Y(nh)+\int_{nh}^{(n+1)h} G(Y(s)) \mathrm{d}s\;.$$

Or si le pas $h$ est petit, la valeur de l'intégrale est proche du produit $hG(Y(nh))$ (la «surface du rectangle»  en dimension $1$). Il est naturel de définir la suite $z_n$ par $z_0=y_0$ et pour tout $n\geqslant 0$ :

$$z_{n+1} = z_n + hG(z_n)\;.$$

Exemple : Explicitons le calcul des vecteurs $z_n$ dans le cas particulier d'un système linéaire à coefficients constants :

$$\left\{ \begin{array}{lcl} Y'(t)&=&AY(t)\;,\\ Y(0)&=&y_0\;. \end{array}\right.$$

On a :

$$z_{n+1} = z_n + hAz_n = (I_d+hA)z_n\;.$$

Avec la condition initiale $z_0=y_0$, on obtient :

$$z_n=(I_d+hA)^ny_0\;.$$

La fonction censée approcher $Y(t)$ sera notée $Z_h(t)$. Elle est définie entre les instants de discrétisation par une interpolation linéaire :

$$\forall t\in [nh,(n+1)h]\;,\quad Z_h(t) = z_n + (t-nh)G(z_n)\;.$$

Théorème 8   Supposons que $G$ soit continûment différentiable, ses dérivées partielles étant bornées uniformément par $M$. Soit $Y$ la solution du problème de Cauchy (10). Fixons $T>0$. Sur l'intervalle compact $[0,T]$, la fonction continue $\Vert G(Y(t))\Vert$ atteint son maximum ; notons-le $K(T)$ :

$$K(T) = \sup_{t\in [0,T]} \Vert G(Y(t))\Vert\;.$$

Alors l'erreur maximale sur $[0,T]$ entre la solution approchée $Z_h$ et la solution exacte $Y$ est telle que :

$$\sup_{t\in [0,T]} \Vert Y(t)-Z_h(t)\Vert \leqslant \frac{K(T)}{2}\Big(\mathrm{e}^{MT}-1\Big) h\;.$$ (11)

En particulier $Z_h$ converge vers $Y$ uniformément sur $[0,T]$.

La majoration (11) montre que l'erreur commise en remplaçant la solution exacte par la solution approchée est de l'ordre de $h$, et tend donc vers 0 quand $h$ tend vers 0. La constante multiplicative montre que l'erreur sur $[0,T]$ peut dépendre fortement de $T$. On comprend intuitivement que, partant d'une condition initiale exacte, les erreurs commises à chaque pas puissent se cumuler, de sorte que l'erreur à l'instant $T$ est exponentielle en $T$.

Exemple : Considérons le cas de l'exponentielle dans $\mathbb{R}$ : $y'(t)=y(t)$, avec $y(0)=1$ (figure 12).

Figure 12: Approximation de l'exponentielle par la méthode d'Euler.

La suite approximante est :

$$z_n = (1+h)^n\;.$$

Pour $T=nh$ fixé, la valeur de $Z_h(T)$ est :

$$Z_h(T) = (1+h)^{T/h} = \mathrm{e}^T\left(1+ T\frac{h}{2} + o(h)\right)\;.$$

L'erreur commise est bien de l'ordre de $h$. Le tableau 1 donne les valeurs de l'erreur $\vert\mathrm{e}^{T}-Z_h(T)\vert$, pour $h=10^{-1}$, $10^{-2}$ et $10^{-3}$, et $T=nh$ allant de $1$ à $10$. Même pour $h=10^{-3}$, on constate que l'erreur en $T=10$ est très importante.

Tableau 1: Erreurs de la méthode d'Euler sur $\mathrm {e}^{nh}$.

La méthode d'Euler, très facile à programmer et peu coûteuse en temps de calcul, peut suffire pour certaines applications. De nombreuses autres méthodes ont été imaginées. Nous nous contenterons d'indiquer une voie de généralisation.

L'idée de base de la méthode d'Euler était d'approcher :

$$\int_{nh}^{(n+1)h} G(y(s)) \mathrm{d}s$$   par$$\quad hG(y(nh))\;.$$

Ceci est l'approximation la plus rudimentaire (méthode des rectangles à gauche) pour un calcul d'intégrale. On peut faire beaucoup mieux, par exemple par la méthode des trapèzes. Il s'agit alors d'approcher

$$\int_{nh}^{(n+1)h} G(y(s)) \mathrm{d}s$$   par$$\quad \frac{h}{2}\Big(G(y(nh))+G(y((n+1)h))\Big)\;.$$

Ceci conduit à définir la suite des valeurs approchées $z_n$ par $z_0=y_0$ et pour tout $n\geqslant 0$ :

$$z_{n+1} = z_n + \frac{h}{2}\Big(G(z_n)+G(z_{n+1})\Big)\;.$$

Sous cette forme, $z_{n+1}$ est définie comme la solution d'une équation, qu'il faut résoudre numériquement. On dit dans ce cas que la méthode est implicite. Pour la rendre explicite, on peut remplacer dans le membre de droite $z_{n+1}$ par l'approximation issue de la méthode d'Euler. Cela donne :

$$z_{n+1} = z_n + \frac{h}{2}\Big(G(z_n)+G(z_n+hG(z_n))\Big)\;.$$

Ce que l'on obtient est un cas particulier de la méthode de Runge-Kutta. Notons que les calculs seront nécessairement plus coûteux, car ils demandent à chaque pas deux évaluations de la fonction $G$, au lieu d'une. En revanche, on obtiendra un résultat plus précis : on démontre que l'erreur commise sur un intervalle de temps fixé $[0,T]$ est de l'ordre de $h^2$, contre $h$ pour la méthode d'Euler.

Exemple : Reprenons le cas de l'exponentielle dans $\mathbb{R}$ : $y'(t)=y(t)$, avec $y(0)=1$. La suite approximante est solution de l'équation de récurrence :

$$z_{n+1} = z_n+\frac{h}{2} (z_n+(z_n(1+h)) = z_n\Big(1+h+\frac{h^2}{2}\Big)\;,$$

soit :

$$z_n=\Big(1+h+\frac{h^2}{2}\Big)^n\;.$$

Le tableau 2 donne les valeurs de l'erreur $\vert\mathrm{e}^{nh}-z_n\vert$, pour $h=10^{-1}$, $10^{-2}$ et $10^{-3}$, et $T=nh$ allant de $1$ à $10$. On constate que les erreurs sont beaucoup plus faibles que celles de la méthode d'Euler.

Tableau 2: Erreurs de la méthode de Runge-Kutta sur $\mathrm {e}^{nh}$.