Table des matières MathematicaScilabPython

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

1. Méthode

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 :

Le vecteur des variables à l'instant est donné par :

La seconde formule est d'ordre 4 :

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 :

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 la valeur obtenue pour la variable i à l'instant et l'estimation de l'erreur pour cette même variable.

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

est une tolérance absolue et 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 yn(i) s'annule.

On calcule le vecteur suivant :

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 h5, le nouveau pas est calculé par :

désigne une norme du vecteur W. La norme la plus simple est :

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 :

Pour , on retrouve la norme précédente.

Si , la valeur de 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ù

2. Exemples

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 :

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

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 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 . De même, l'option PrecisionGoal définit la tolérance relative .

La solution de ce système étant connue , on trace l'erreur en fonction du temps :

Plot[u[t]-Cos[t]/.solution,{t,0,tmax}]
figA.pngfigA.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)
figBfigB.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.pngfigC.pdf

3. Stabilité

La méthode de Runge-Kutta à formules emboîtées (ordres 4 et 5) semble être instable pour l'équation de l'oscillateur harmonique. Avec Mathematica, il est possible d'essayer un ordre supérieur :

tmax=10000*Pi
solution=NDSolve[systeme,{u,v},{t,0,tmax},Method->{"ExplicitRungeKutta",DifferenceOrder->6,EmbeddedDifferenceOrder->7},
                        NormFunction->Infinity,AccuracyGoal->4,PrecisionGoal->10,MaxSteps->Infinity,StartingStepSize->0.01]
Plot[u[t]-Cos[t]/.solution,{t,0,tmax}]
figD.pngfigD.pdf

L'erreur semble bornée dans un premier temps puis finit par augmenter. Cependant, l'augmentation semble linéaire (et non pas exponentielle).

Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.