• Champ des lments de contact
• Mthode d'Euler-Cauchy
• Mthode d'Euler-Cauchy amliore ou mthode de Heun
• Mthode de Runge-Kutta

Mthodes numriques pour
                 les quations diffrentielles:

Mthode d'Euler-Cauchy

Mthode de Heun

Mthode de Runge-Kutta

Pierre Lantagne (avril 2001)

Collge de Maisonneuve

plantag@edu.cmaisonneuve.qc.ca

http://math.cmaisonneuve.qc.ca/plantagne

Nous connaissons plusieurs mthodes pour la rsolution analytique de certaines quations diffrentielles du premier ordre et du premier degr. Par exemples, pour les

quations diffrentielles variables sparables

quations diffrentielles homognes

quations diffrentielles linaires

quations diffrentielles de Bernouilli

quations diffrentielles exactes

quations diffrentielles non exactes

Mais, il y a plusieurs quations diffrentielles du premier ordre qui ne possdent pas de solutions analytiques. Pour de telles quations, il y a plusieurs mthodes d'approximation numriques :

mthode d'Euler-Cauchy

mthode  de Heun (ou mthode amliore d'Euler-Cauchy)

mthode de Runge-Kutta

Dans un  premier temps, nous alllons prouver la mthode d'Euler-Cauchy avec la solution analytique de l'quation diffrentielle du premier ordre .

Champ des lments de contact

Une quation diffrentielle du premier ordre est une quation de la forme .  En se rappelant l'interprtation graphique de ,  f(x,y) est donc une formule donnant la pente de la tangente la courbe solution passant par le point (x,y). Ainsi, l'quation diffrentielle peut tre visualise par un graphique appel champ des lments de contact. Chaque lment de ce champ est un petit segment de droite centr au point ( ) d'orientation f( ).

Ainsi, pour l'quation diffrentielle , au point (-2,3), la pente de la tangente la courbe solution y passant par ce point est   = 1.  Au point (6,2) la pente de la tangente la courbe solution est . Dans le premier cas, on illustre le rsultat par un petit segment de droite centr en (-2,3) de pente 1 et, dans le second cas, par un autre segment de droite centr en (6,2) de pente . En rptant ce prcessus pour un certain nombre de couples (x,y), on obtient ce qu'on appelle un champ d'lments de contact, parfois traduit de l'anglais par champ de pente.

La macro-commande dfieldplot de l'extension  DEtools permet le trac d'un champ d'lments de contact. L'extension DEtools impose la formulation explicite de la variable dpendante. Dans le cas d'une fonction de deux variables x et y, si la variable y est dsigne comme dpendante, on doit l'noncer dans les requtes avec la syntaxe y(x).

>	with(plots): with(DEtools):

Warning, the name changecoords has been redefined

>	ED:=diff(y(x),x)=-2*x-y(x);

>	Champ:=dfieldplot(ED,[y(x)], x=-2..2,y=-2..2,arrows=line,dirgrid=[25,25],scaling=constrained,color=orange): display(Champ);

>

On peut galement reprsenter dans un champ d'lments de contact, des courbes appeles isoclines ( du grec iso qui signifie mme et du latin clinare qui signifie pencher). Les isoclines sont des courbes le long desquelles les lments de contact ont une direction ( une inclinaison ) donne.

Superposer, dans le champ d'lments de contact prcdent, les isoclines de pente  0 et 1.

>	f:=(x,y)->-2*x-y;

>	j:=0: for i in [-1,0,1] do  iso||j:=implicitplot(f(x,y)=i,x=-2..2,y=-2..2,color=blue):  j:=j+1: od:

>	display({iso||(0..j-1),Champ},view=[-1..1,-1..1],title="Isoclines d'quations f(x,y)=c");

>

L'intrt du champ d'lments de contact apparat clairement avec l'usage de l'ordinateur. Que l'quation diffrentielle possde ou non une solution analytique, ce type de trac permet de visualiser l'ensemble des courbes solution de l'quation diffrentielle. En effet, le champ d'lments de contact prcdent permet de visualiser ceci:

lorsque x < , chaque courbe solution est croissante

lorsque x > , chaque courbe solution est dcroissante

lorsque x = , chaque courbe solution admet un maximum relatif

lorsque x , toutes les courbes solutions ont un comportement asymptotique

>

Puisque est une quation diffrentielle du premier ordre et du premier degr, mme s'il est facile d'en obtenir directement la solution par un calcul la main, rsolvons quand mme cette quation avec Maple.

>	dsolve(ED);

La solution gnrale de cette quation est donc

>	Sol:=y=-2*x+2+C*exp(-x);

Superposons au champ d'lments de contact prcdent, les solutions particulires correspondant aux conditions initiales pour C =   , , 0, 1 et 2.

>	Sol:=-2*x+2+C1*exp(-x): i:=0: for C in [-2,-1,0,1,2] do  C1:=C-2:  f||i:=plot([x,Sol,x=-2..2],color=navy,thickness=2):  i:=i+1 od: C1:='C1': C:='C':

>	display({f||(0..4),Champ},view=[-2..2,-2..2],title="dy/dx = -2x-y",titlefont=[TIMES,ROMAN,12]);

La macro-commande DEplot de l'extension DEtools permet directement la superposition du trac d'un champ d'lments de contact avec ceux des solutions particulires.

>	DEplot(ED,y(x),x=-2..2,[[y(0)=-2],[y(0)=-1],[y(0)=0],[y(0)=1],[y(0)=2]],y=-2..2,linecolor=navy,color=orange,scaling=constrained,arrows=line,dirgrid=[25,25]);

>

Mthode d'Euler-Cauchy

Mettons en vidence la solution particulire pour x [0; 1,6].

>	Sol_par:=dsolve({ED,y(0)=-1},y(x));

>	plot([x,rhs(Sol_par),x=0..1.6]);

>

Ce trac est exact car il a t obtenu l'aide de la solution gnrale. Obtenons maintenant une approximation de ce trac en joignant les points que l'on obtient avec la mthode d'Euler-Cauchy. En fait, il s'agit d'approximer la courbe solution de l'quation diffrentielle = avec la condition initiale .

Cette mthode repose sur le processus itratif suivant. La condition initiale prcise un premier point. Ce point est videmment exact. Soit ( ) ce premier point. En ce point, la pente de la tangente la courbe solution est donne par f( ). Alors, l'quation de la tangente en ce point est

c'est--dire

L'expression ( ) sera vue comme un accroissement d'abscisse h qui sera assez petit. cet accroissement, correspondra une valeur d'ordonne jusqu' la tangente . On obtient ainsi un deuxime point .

De mme, un troisime point de coordonnes sera obtenu l'aide de la direction de la tangente passant par le deuxime point f( ).

Et ainsi de suite: .

REMARQUE: La solution d'une quation diffrentielle est une fonction. La mthode d'Euler-Cauchy ne donne pas une fonction: elle permet seulement d'obtenir une squence de couple (x,y) qu'il suffit de relier de manire approprie pour approximer la courbe solution.

Appliquons maintenant ce raisonnement pour obtenir une srie de points que l'on reliera ensuite dans un trac qui sera superpos celui de la solution analytique.

Automatisons le calcul d'un nombre n de couples ordonnes (x,y). Rappelons-nous d'abord l'quation diffrentielle rsoudre.

>

ED;

La condition initiale impose le premier point.

>	X[0]:=0; Y[0]:=-1; Points_Euler:=[X[0],Y[0]];

Posons un accroissement .

>

h:=0.1;

Rpartissons le nombre supplmentaires de points par le calcul suivant: la largeur de l'intervalle d'approximation [0; 1,6] divise par la valeur de l'accroissement h.

>	N:=(1.6-0)/h;

>	for n from 0 to N-1 do X[n+1]:=X[n]+h: Y[n+1]:=Y[n]+h*f(X[n],Y[n]):  # f(x,y) est la formule de calcul de la pente Tampon:=[X[n+1],Y[n+1]]; Points_Euler:=Points_Euler,Tampon od:

>

Superposons  le trac de ces points (sans les relier) avec celui de la courbe solution.

>	Sol:=dsolve({ED,y(0)=-1},y(x)); Courbe_solution:=plot([x,rhs(Sol),x=0..1.6],color=navy,thickness=2): Points_approx:=plot([Points_Euler],style=point,symbol=circle,color=orange): display([Courbe_solution,Points_approx]);

Diminuons le valeur de l'accroissement pour une seconde approximation.

>	Points_Euler:=[X[0],Y[0]];

>

h:=0.05;

>	N:=(1.6-0)/h; for n from 0 to N-1 do X[n+1]:=X[n]+h: Y[n+1]:=Y[n]+h*f(X[n],Y[n]):  # f(x,y) est la formule de calcul de la pente Tampon:=[X[n+1],Y[n+1]]; Points_Euler:=Points_Euler,Tampon od:

>	Sol:=dsolve({ED,y(0)=-1},y(x)); Courbe_solution:=plot([x,rhs(Sol),x=0..1.6],color=navy,thickness=2): Points_approx:=plot([Points_Euler],style=point,symbol=circle,color=magenta): display([Courbe_solution,Points_approx]);

Relions finalement ces points avec la macro-commande pointplot de l'extension plots.

>	Points_approx:=pointplot([Points_Euler],connect=true,color=magenta): display([Courbe_solution,Points_approx]);

>

Mthode d'Euler-Cauchy amliore ou mthode de Heun

Cette mthode a comme caractristique que l'ordonne du prochain point est calcul non pas avec la pente de la tangente la courbe solution du point courant mais plutt avec la moyenne arithmtique des pentes des tangentes la courbe solution du point courant et du point suivant temporaire calcul dans la direction donne par la tangente au point suivant et ce, partir de l'estimation de l'ordonne du point courant.

Dans le cas de la mthode de Heun, on a:

  

  

  

Pour bien automatiser le calcul des points de la courbe solution estimer, automatisons d'abord la mthode d'Euler-Cauchy avec la procdure suivante.

>

Euler:=proc(f::procedure,y0::realcons,Intervalle::range,h::realcons) local n,N,Points_Euler,Tampon,X,Y; X[0]:=op(1,Intervalle); Y[0]:=y0; Points_Euler:=[X[0],Y[0]]:  # Point de la condition initiale N:=(op(2,Intervalle)-X[0])/h: for n from 0 to N-1 do X[n+1]:=X[n]+h: Y[n+1]:=Y[n]+h*f(X[n],Y[n]):  # f(x,y) est la formule de calcul de la pente Tampon:=[X[n+1],Y[n+1]]; Points_Euler:=Points_Euler,Tampon od end:

>

Maintenant, automatisons l'algorithme de la mthode d'Euler-Cauchy amliore, c'est--dire de la mthode de Heun.

>

Heun:=proc(f::procedure,y0::realcons,Intervalle::range,h::realcons) local n,N,Points_Heun,Tampon,X,Y,U; X[0]:=op(1,Intervalle); Y[0]:=y0; Points_Heun:=[X[0],Y[0]]:  # Point de la condition initiale N:=(op(2,Intervalle)-X[0])/h: for n from 0 to N-1 do X[n+1]:=X[n]+h: U[n+1]:=Y[n]+h*f(X[n],Y[n]):  # f(x,y) est la formule de calcul de la pente Y[n+1]:=Y[n]+h/2*(f(X[n],Y[n])+f(X[n+1],U[n+1])): Tampon:=[X[n+1],Y[n+1]]; Points_Heun:=Points_Heun,Tampon od end:

>

Comparons maintenant la mthode de Heun avec celle d'Euler-Cauchy avec l'quation diffrentielle du dbut avec comme condition initiale .

>	f:=(x,y)->-2*x-y;

>

Euler_Approx:=plot([Euler(f,-1,0..1.6,0.1)],style=point,symbol=circle,color=magenta): Heun_Approx:=plot([Heun(f,-1,0..1.6,0.1)],style=point,symbol=circle,color=green): Champ:=dfieldplot(ED,[y(x)], x=0..1.6,y=-2..-0.5,arrows=line,dirgrid=[15,15],scaling=constrained,color=orange): display([Heun_Approx,Euler_Approx,Courbe_solution,Champ]);

>

On constate effectivement que la mthode de Heun (points verts) amliore beaucoup l'estimation de la courbe solution obtenue avec la mthode d'Euler-Cauchy (points magenta).

>

Mthode de Runge-Kutta

La mthode de Runge-Kutta du quatrime ordre amliore davantage la mthode d'Euler amliore mais les calculs sont encore beaucoup plus labors. Cette mthode repose aussi en quelque sorte sur un calcul de moyenne arithmtique de pentes.  Prsentons, sans autres explications, la procdure automatisant l'algorithme de Runge-Kutta.

>

Runge_Kutta:=proc(f::procedure,y0::realcons,Intervalle::range,h::realcons) local n,N,p,Points_Runge_Kutta,q,r,s,Tampon,X,Y; X[0]:=op(1,Intervalle); Y[0]:=y0; Points_Runge_Kutta:=[X[0],Y[0]]:  # Point de la condition initiale N:=(op(2,Intervalle)-X[0])/h: for n from 0 to N-1 do X[n+1]:=X[n]+h: p[n]:=f(X[n],Y[n]): q[n]:=f(X[n]+h/2,Y[n]+h/2*p[n]): r[n]:=f(X[n]+h/2,Y[n]+h/2*q[n]): s[n]:=f(X[n]+h,Y[n]+h*r[n]): Y[n+1]:=Y[n]+h/6*(p[n]+2*q[n]+2*r[n]+s[n]): Tampon:=[X[n+1],Y[n+1]]; Points_Runge_Kutta:=Points_Runge_Kutta,Tampon od end:

Pour l'quation diffrentielle , la mthode de Heun montrait la section prcdente que les points calculs pousaient dj, de manire remarquable, la courbe solution. Le degr de prcision atteint avec cette mthode a pu tre apprci par une approche graphique l'gard de la mthode d'Euler. Par contre, avec la mthode du Runge-Kutta, le degr de prcision atteint l'gard de la mthode de Heun est difficilement apprciable graphiquement ( l'cran). En effet, constatons que les tracs des courbes en reliant les points calculs par les mthodes de Heun et de Runge-Kutta et la courbe solution se superposent presque parfaitement compte tenu de la rsolution cran. Il est ncessaire de procder au grosissant maximum de l'affichage (Ctrl-6), pour pouvoir   visualiser l'augmentation de la prcision.

>	A:=array(1..2): C1:=pointplot([Heun(f,-1,0..1.6,0.1)],connect=true,color=orange): C2:=pointplot([Runge_Kutta(f,-1,0..1.6,0.1)],connect=true,color=cyan): A[1]:=display([Courbe_solution,C1,C2]): A[2]:=display([Courbe_solution,C2,C1]): display(A);

>

Seul un calcul numrique pourra permettre clairement une comparaison entre les degrs de prcision atteint par chaque mthode l'gard de la courbe solution.

>