Résolution numérique d'une équation différentielle d'ordre 1

Nous verrons ici quelques techniques de résolution numérique d'équations différentielles pour des équations d'ordre 1. Pour plus de renseignements, consultez l'aide de Maple en cherchant par exemple le terme suivant: dsolve[classical] . Toutes les méthodes classiques sont implantées dans Maple, mais on peut également définir soi-même la procédure qui nous intéresse (voir le fichier numeric-proc.mws).

Une des syntaxes possibles pour la commande est :

> dsolve(équation, numeric, method=classical[choix], vars, options )

où

équation représente notre équation différentielle (ou un système d'équation) avec condition(s) initiale(s)

numeric indique que l'on résoud par une méthode numérique

method=classical[choix] indique quelle méthode est utilisée :

foreuler pour Euler , heunform pour Euler améliorée et rk4 pour la méthode de Runge-Kutta d'ordre 4.

Nous résoudrons ici l'équation suivante

On cherche à déterminer y(1) ou à estimer une courbe solution quand x varie de 0 à 1.

La ligne ci-dessous permet de définir notre équation différentielle ( eqdif ) , la condition initiale ( init ) sous la forme y(a)=c ainsi que la valeur de x finale désirée ( b ) et le nombre d'itérations ( n ) ce qui permet de déterminer le pas ( h ) à chaque itération. On peut également fournir une condition initiale et un pas directement sans fournir une valeur finale (voir plus bas).

> eqdif :=diff(y(x),x)=x^2+y(x)^2; a:=0: init:=y(a)=0; b:=1: n:=10: h:=evalf((b-a)/n,3);

Méthode d'Euler

L'option " value=array( ) " dans la commande suivante permet de construire une matrice des valeurs obtenues. On utilise l'argument " foreuler " dans la commande classical[choix] .

Si vous travaillez avec un grand nombre d'itérations, vous devriez mettre un deux-points ( : ) au lieu du point-virgule ( ; ) à la fin de la commande dsolve( .....) pour éviter d'afficher le résultat, à savoir la matrice nommée repEuler .

> repEuler:=dsolve({eqdif,init},y(x),type=numeric,method=classical[foreuler],value=array([seq(a+i*h,i=0..n)]),stepsize=h);

Warning, the 'value' option is deprecated, use 'output' instead

Méthodes d'Euler amlioré et de Runge-Kutta

On utilise l'argument " heunform " pour la méthode d'Euler amélioré et l'argument " rk4 " pour Runge-Kutta d'ordre 4

> repHeun:=dsolve({eqdif,init},y(x),type=numeric,method=classical[heunform],value=array([seq(a+i*h,i=0..n)]),stepsize=h);

> repRK:=dsolve({eqdif,init},y(x),type=numeric,method=classical[rk4],value=array([seq(a+i*h,i=0..n)]),stepsize=h);

Warning, the 'value' option is deprecated, use 'output' instead

Warning, the 'value' option is deprecated, use 'output' instead

Comparaison graphique des 3 solutions numériques

La commande pour faire tracer une résolution numérique d'équation différentielle est odeplot( ) , qu'on peut utiliser après avoir chargé la librairie "plots".

Nous tracerons à l'aide de la commande display( ) les 3 solutions précédentes sur le même graphique.

> with(plots):

> A:=odeplot(repEuler,[x,y(x)],color=red):

> B:=odeplot(repHeun,[x,y(x)],color=blue):

> C:=odeplot(repRK,[x,y(x)],color=green):

> display({A,B,C});

On remarque que les solutions avec Euler amélioré et avec Runge-Kutta se confondent sur le graphique. Par contre, en considérant les tableaux de valeurs obtenues on peut voir que la valeur de y( 1 ) obtenue par ces deux méthodes diffère après la deuxième décimale.

Warning, the name changecoords has been redefined

Accès direct à une valeur y( b ) à partir d'une condition initiale y( a ) = c

On peut également associer cette commande dsolve( ) à un nom de variable et ainsi créer une procédure qu'on peut appeler au besoin en évaluant la valeur voulue, en appliquant la méthode numérique avec un pas " h " spécifié, sans créer une matrice de points intermédiaires.

> monEuler:=dsolve({eqdif,init},y(x),type=numeric,method=classical[foreuler],stepsize=0.1);

> monEuler(1);

> monRK:=dsolve({eqdif,init},y(x),type=numeric,method=classical[rk4],stepsize=0.01);

> monRK(1);

Ce dernier estimé pour y( 1 ) est obtenu par la méthode de Runge-Kutta avec un pas de h = 0.01 donc en 100 étapes. Cela permet d'apprécier la précision des résultats obtenus précédemment avec n = 10 étapes.