Méthodes de Runge-Kutta à formules emboîtées

1.a. Formules emboîtées

Il s'agit d'utiliser deux formules de Runge-Kutta simultanément pour chaque pas, utilisant les mêmes évaluations de dérivées mais avec des coefficients différents. L'une est à l'ordre p, l'autre à l'ordre p+1.

Considérons par exemple une méthode de Runge-Kutta d'ordre 5 avec les 6 évaluations suivantes :

$$\begin{array}{cc}\hfill & k_1=hf(t_n,y_n)\hfill \\ \hfill & k_2=hf(t_n+a_{2}h,y_n+b_{21}{k}_1)\hfill \\ \hfill & \cdots \hfill \\ \hfill & k_6=hf(t_n+a_{6}h,y_n+b_{61}{k}_1+b_{62}{k}_2+\cdots b_{65}{k}_5)\hfill \end{array}$$

Le vecteur des variables à l'instant $t_{n+1}$ est donné par :

$$y_{n+1}=y_n+c_1{k}_1+c_2{k}_2+c_3{k}_3+c_4{k}_4+c_5{k}_5+c_6{k}_6+O\left(h^6\right)$$

La seconde formule est d'ordre 4 :

$$y_{n+1}^{*}=y_n+c_1^{*}{k}_1+c_2^{*}{k}_2+c_3^{*}{k}_3+c_4^{*}{k}_4+c_5^{*}{k}_5+c_6^{*}{k}_6+O\left(h^5\right)$$

L'intérêt de ces deux formules emboîtées est de permettre une estimation de l'erreur locale pour la formule d'ordre 4 :

$$\Delta =y_{n+1}-y_{n+1}^{*}=\sum _{i=1}^6(c_i-c_i^{*})k_i$$

Par rapport à une formule d'ordre 4 classique, la méthode nécessite deux évaluations en plus. Cependant, l'évaluation de l'erreur est indispensable pour adapter le pas de temps h à la précision souhaitée.

1.b. Adaptation du pas

L'adaptation du pas consiste à modifier le pas h de manière à maintenir l'erreur inférieure à une limite fixée. L'erreur est évaluée à chaque pas selon la méthode données ci-dessus, puis le nouveau pas est calculé en fonction de celle-ci. Remarquons que le contrôle de l'erreur est local (à chaque pas). L'erreur cummulée après N pas ne peut être contrôlé directement.

Notons $y_n\left(i\right)$ la valeur obtenue pour la variable i à l'instant $t_n$ et $\Delta \left(i\right)$ l'estimation de l'erreur pour cette même variable.

On se fixe comme objectif de maintenir l'erreur inférieure à

$${\epsilon }_a+{\epsilon }_r\left|y\left(i\right)\right|$$

où ${\epsilon }_a$ est une tolérance absolue et ${\epsilon }_r$ une tolérance relative. Si l'on souhaite simplement une tolérance relative, on posera la première très faible, mais non nulle pour éviter une division par zéro lorque y_n(i) s'annule.

On calcule le vecteur suivant :

$$W\left(i\right)=\frac{\Delta \left(i\right)}{{\epsilon }_a+{\epsilon }_r\left|y_n\left(i\right)\right|}$$

Si des composantes de ce vecteur sont supérieures à 1, il faut diminuer le pas. Inversement, si toutes les composantes sont inférieures à 1, on peut augmenter le pas. Sachant que l'erreur est d'ordre h⁵, le nouveau pas est calculé par :

$$h_1=h(\frac{1}{\left|\right|W\left|\right|}{)}^{1/5}$$

où $\left|\right|W\left|\right|$ désigne une norme du vecteur W. La norme la plus simple est :

$$\left|\right|W\left|\right|=Max\left|W\left(i\right)\right|$$

qui revient à tenir compte uniquement de la variable dont l'erreur estimée est la plus forte (par rapport à la tolérance fixée).

On peut aussi utiliser la norme suivante :

$$\left|\right|W\left|\right|=(\sum _{i=1}^{N}W{\left(i\right)}^p{)}^{1/p}$$

Pour $p\to \infty$, on retrouve la norme précédente.

Si $h_1<h$, la valeur de $y_{n+1}$ doit être recalculée avec le nouveau pas. En pratique, on constate que l'exposant 1/5 conduit à plusieurs essais successifs avec réduction du pas. Un exposant 1/4 fonctionne mieux dans le cas où $\left|\right|W\left|\right|>1$

2.a. Mathematica

La résolution numérique d'un système différentiel du premier ordre est obtenue avec la fonction NDSolve.

Considérons comme exemple l'équation différentielle de l'oscillateur harmonique :

$$\frac{d^{2}u\left(t\right)}{d{t}^2}+u\left(t\right)=0$$

Il faut tout d'abord définir le système à deux variables du premier ordre correspondant :

$$\begin{array}{cc}\hfill & \frac{du}{dt}=v\hfill \\ \hfill & \frac{dv}{dt}=-u\hfill \end{array}$$

Dans Mathematica, il faut définir les équations dans une liste avec les conditions initiales nécessaires :

systeme={u'[t]==v[t],v'[t]==-u[t],u[0]==1,v[0]==0}

Effectuons la résolution numérique avec une méthode de Runge-Kutta à formules emboîtées :

tmax=100*Pi
solution=NDSolve[systeme,{u,v},{t,0,tmax},Method->{"ExplicitRungeKutta",DifferenceOrder->4,EmbeddedDifferenceOrder->5},
                        NormFunction->Infinity,AccuracyGoal->4,PrecisionGoal->10,MaxSteps->Infinity,StartingStepSize->0.01]

L'option NormFunction permet de définir la norme utilisée pour $\left|\right|W\left|\right|$ ou simplement l'entier p pour la norme définie plus haut. Ici, on choisit de prendre le maximum des valeurs absolues des composantes.

L'option AccuracyGoal définit la tolérance absolue, par exemple ici ${\epsilon }_a=10^{-4}$. De même, l'option PrecisionGoal définit la tolérance relative ${\epsilon }_r=10^{-10}$.

La solution de ce système étant connue $u\left(t\right)=cos\left(t\right)$, on trace l'erreur en fonction du temps :

Plot[u[t]-Cos[t]/.solution,{t,0,tmax}]

figA.pdf

2.b. Solveur C++

Ce solveur, écrit en C++, est décrit ici. le résultat du calcul est récupéré en Python.

On reprend l'exemple précédent. Instant final t=100π. Tolérance relative 10^-10, tolérance absolue 10^-4, pas initial 0.01.

from numpy import *
from pylab import *
                

erreur=y-cos(t)
plot(t,erreur)

figB.pdf

2.c. Scilab

Le système différentiel doit être défini par une fonction :

function [deriv]=systeme(t,y), 
    deriv(1)=y(2),
    deriv(2)=-y(1),
endfunction
                

Les instants à calculer sont définis dans un vecteur t puis on utilise la fonction ode qui renvoit une matrice dont la ligne n correspond aux valeurs de la variable n aux instants demandés.

t=[0:0.1:100*%pi];
tolA=1d-4;
tolR=1d-10;
y=ode('rkf',[1;0],0,t,tolR,tolA,systeme);
                

figC=scf();
plot2d(t,y(1,:)-cos(t))

figC.pdf