<Text-field layout="Heading 1" style="Heading 1"> Champ des \351l\351ments de contact</Text-field>

Mais, il y a plusieurs \351quations diff\351rentielles du premier ordre qui ne poss\350dent pas de solutions analytiques. Pour de telles \351quations, il y a plusieurs m\351thodes d'approximation num\351riques :

m\351thode d'Euler-Cauchy

m\351thode de Heun (ou m\351thode am\351lior\351e d'Euler-Cauchy)

m\351thode de Runge-Kutta

Dans un premier temps, nous alllons \351prouver la m\351thode d'Euler-Cauchy avec la solution analytique de l'\351quation diff\351rentielle du premier ordre NiMvKiYlI2R5RyIiIiUjZHhHISIiLCYqJiIiI0YmJSJ4R0YmRiglInlHRig=.

restart;

<Text-field layout="Heading 1" style="Heading 1"> Champ des \351l\351ments de contact</Text-field>

Une \351quation diff\351rentielle du premier ordre est une \351quation de la forme NiMvKiYlI2R5RyIiIiUjZHhHISIiLSUiZkc2JCUieEclInlH. En se rappelant l'interpr\351tation graphique de NiMqJiUjZHlHIiIiJSNkeEchIiI=, f(x,y) est donc une formule donnant la pente de la tangente \340 la courbe solution passant par le point (x,y). Ainsi, l'\351quation diff\351rentielle NiMvKiYlI2R5RyIiIiUjZHhHISIiLCYqJiIiI0YmJSJ4R0YmRiglInlHRig= peut \352tre visualis\351e par un graphique appel\351 champ des \351l\351ments de contact. Chaque \351l\351ment de ce champ est un petit segment de droite centr\351 au point (NiQmJSJ4RzYjIiIhJiUieUdGJQ==) d'orientation f(NiQmJSJ4RzYjIiIhJiUieUdGJQ==).

Ainsi, pour l'\351quation diff\351rentielle NiMvKiYlI2R5RyIiIiUjZHhHISIiLCYqJiIiI0YmJSJ4R0YmRiglInlHRig=, au point (-2,3), la pente de la tangente \340 la courbe solution y passant par ce point est NiMvJSN5J0csJiomIiIjIiIiLCRGJyEiIkYoRioiIiRGKg== = 1. Au point (6,2) la pente de la tangente \340 la courbe solution est NiMsJCIjOSEiIg==. Dans le premier cas, on illustre le r\351sultat par un petit segment de droite centr\351 en (-2,3) de pente 1 et, dans le second cas, par un autre segment de droite centr\351 en (6,2) de pente NiMsJCIjOSEiIg==. En r\351p\351tant ce pr\351cessus pour un certain nombre de couples (x,y), on obtient ce qu'on appelle un champ d'\351l\351ments de contact, parfois traduit de l'anglais par champ de pente.

La macro-commande dfieldplot de l'extension DEtools permet le trac\351 d'un champ d'\351l\351ments de contact. L'extension DEtools impose la formulation explicite de la variable d\351pendante. Dans le cas d'une fonction de deux variables x et y, si la variable y est d\351sign\351e comme d\351pendante, on doit l'\351noncer dans les requ\352tes avec la syntaxe y(x).

with(plots): with(DEtools):

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 \351galement repr\351senter dans un champ d'\351l\351ments de contact, des courbes appel\351es isoclines ( du grec iso qui signifie m\352me et du latin clinare qui signifie pencher). Les isoclines sont des courbes le long desquelles les \351l\351ments de contact ont une direction ( une inclinaison ) donn\351e.

Superposer, dans le champ d'\351l\351ments de contact pr\351c\351dent, les isoclines de pente 0 et \2611.

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'\351quations f(x,y)=c");

L'int\351r\352t du champ d'\351l\351ments de contact appara\356t clairement avec l'usage de l'ordinateur. Que l'\351quation diff\351rentielle poss\350de ou non une solution analytique, ce type de trac\351 permet de visualiser l'ensemble des courbes solution de l'\351quation diff\351rentielle. En effet, le champ d'\351l\351ments de contact pr\351c\351dent permet de visualiser ceci:

lorsque x < NiMsJComJSJ5RyIiIiIiIyEiIkYo, chaque courbe solution est croissante

lorsque x > NiMsJComJSJ5RyIiIiIiIyEiIkYo, chaque courbe solution est d\351croissante

lorsque x = NiMsJComJSJ5RyIiIiIiIyEiIkYo, chaque courbe solution admet un maximum relatif

lorsque x \256 NiMlKWluZmluaXR5Rw==, toutes les courbes solutions ont un comportement asymptotique

Puisque NiMvJSN5J0csJiUieEciIiIlInlHISIi est une \351quation diff\351rentielle du premier ordre et du premier degr\351, m\352me s'il est facile d'en obtenir directement la solution par un calcul \340 la main, r\351solvons quand m\352me cette \351quation avec Maple.

dsolve(ED);

La solution g\351n\351rale de cette \351quation est donc

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

Superposons au champ d'\351l\351ments de contact pr\351c\351dent, les solutions particuli\350res correspondant aux conditions initiales NiMvLSUieUc2IyIiISUiQ0c= pour C = NiMsJCIiIyEiIg==, NiMsJCIiIiEiIg==, 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\351 d'un champ d'\351l\351ments de contact avec ceux des solutions particuli\350res.

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]);

<Text-field layout="Heading 1" style="Heading 1"> M\351thode d'Euler-Cauchy</Text-field>

Mettons en \351vidence la solution particuli\350re NiMvLSUieUc2IyIiISwkIiIiISIi pour x NiMlKGVwc2lsb25H [0; 1,6].

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

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

Ce trac\351 est exact car il a \351t\351 obtenu \340 l'aide de la solution g\351n\351rale. Obtenons maintenant une approximation de ce trac\351 en joignant les points que l'on obtient avec la m\351thode d'Euler-Cauchy. En fait, il s'agit d'approximer la courbe solution de l'\351quation diff\351rentielle NiMvKiYlI2R5RyIiIiUjZHhHISIiJSN5J0c= = NiMtJSJmRzYkJSJ4RyUieUc= avec la condition initiale NiMvLSUieUc2IyYlInhHNiMiIiEmRiVGKQ==.

Cette m\351thode repose sur le processus it\351ratif suivant. La condition initiale pr\351cise un premier point. Ce point est \351videmment exact. Soit (NiQmJSJ4RzYjIiIhJiUieUdGJQ==) ce premier point. En ce point, la pente de la tangente \340 la courbe solution est donn\351e par f(NiQmJSJ4RzYjIiIhJiUieUdGJQ==). Alors, l'\351quation de la tangente en ce point est

NiMvJSJ5RywmKiYtJSN5J0c2IyYlInhHNiMiIiEiIiIsJkYrRi5GKiEiIkYuRi4tRiRGKUYu

c'est-\340-dire

NiMvJSJ5RywmKiYtJSJmRzYkJiUieEc2IyIiISZGJEYsIiIiLCZGK0YvRiohIiJGL0YvRi5GLw==

L'expression (NiMsJiUieEciIiImRiQ2IyIiISEiIg==) sera vue comme un accroissement d'abscisse h qui sera \253assez\273 petit. \300 cet accroissement, correspondra une valeur d'ordonn\351e jusqu'\340 la tangente NiMvJiUieUc2IyIiIiwmKiYtJSJmRzYkJiUieEc2IyIiISZGJUYvRiclImhHRidGJ0YxRic=. On obtient ainsi un deuxi\350me point NiMvNiQmJSJ4RzYjIiIiJiUieUdGJzYkLCYmRiY2IyIiIUYoJSJoR0YoLCYqJi0lImZHNiRGLSZGKkYuRihGMEYoRihGNkYo.

De m\352me, un troisi\350me point de coordonn\351es NiMvNiQmJSJ4RzYjIiIjJiUieUdGJzYkLCYmRiY2IyIiIkYvJSJoR0YvLCYqJi0lImZHNiRGLSZGKkYuRi9GMEYvRi9GNkYv sera obtenu \340 l'aide de la direction de la tangente passant par le deuxi\350me point f(NiQmJSJ4RzYjIiIiJiUieUdGJQ==).

Et ainsi de suite: NiMvNiQmJSJ4RzYjJSJuRyYlInlHRic2JCwmJkYmNiMsJkYoIiIiRjAhIiJGMCUiaEdGMCwmKiYtJSJmRzYkRi0mRipGLkYwRjJGMEYwRjhGMA==.

REMARQUE: La solution d'une \351quation diff\351rentielle est une fonction. La m\351thode d'Euler-Cauchy ne donne pas une fonction: elle permet seulement d'obtenir une s\351quence de couple (x,y) qu'il suffit de relier de mani\350re appropri\351e pour approximer la courbe solution.

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

Automatisons le calcul d'un nombre n de couples ordonn\351es (x,y). Rappelons-nous d'abord l'\351quation diff\351rentielle \340 r\351soudre.

ED;

La condition initiale impose le premier point.

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

Posons un accroissement NiQvJSJoRyIiISIiIg==.

h:=0.1;

R\351partissons le nombre suppl\351mentaires de points par le calcul suivant: la largeur de l'intervalle d'approximation [0; 1,6] divis\351e 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\351 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:

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]);

<Text-field layout="Heading 1" style="Heading 1"> M\351thode d'Euler-Cauchy am\351lior\351e ou m\351thode de Heun</Text-field>

Cette m\351thode a comme caract\351ristique que l'ordonn\351e du prochain point est calcul\351 non pas avec la pente de la tangente \340 la courbe solution du point courant mais plut\364t avec la moyenne arithm\351tique des pentes des tangentes \340 la courbe solution du point courant et du point suivant temporaire calcul\351 dans la direction donn\351e par la tangente au point suivant et ce, \340 partir de l'estimation de l'ordonn\351e du point courant.

Dans le cas de la m\351thode de Heun, on a:

NiMvJiUieEc2IywmJSJuRyIiIkYpRiksJiZGJTYjRihGKSUiaEdGKQ==

NiMvJiUidUc2IywmJSJuRyIiIkYpRiksJiYlInlHNiNGKEYpKiYlImhHRiktJSJmRzYkJiUieEdGLUYrRilGKQ==

NiMvJiUieUc2IywmJSJuRyIiIkYpRiksJiZGJTYjRihGKSomJSJoR0YpKiYsJi0lImZHNiQmJSJ4R0YsRitGKS1GMjYkJkY1NiMsJkY1RilGKUYpJiUidUdGJkYpRikiIiMhIiJGKUYp

Pour bien automatiser le calcul des points de la courbe solution \340 estimer, automatisons d'abord la m\351thode d'Euler-Cauchy avec la proc\351dure 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 m\351thode d'Euler-Cauchy am\351lior\351e, c'est-\340-dire de la m\351thode 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 m\351thode de Heun avec celle d'Euler-Cauchy avec l'\351quation diff\351rentielle du d\351but NiMvKiYlI2R5RyIiIiUjZHhHISIiLCYqJiIiI0YmJSJ4R0YmRiglInlHRig= avec comme condition initiale NiMvLSUieUc2IyIiISwkIiIiISIi.

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 m\351thode de Heun (points verts) am\351liore beaucoup l'estimation de la courbe solution obtenue avec la m\351thode d'Euler-Cauchy (points magenta).

<Text-field layout="Heading 1" style="Heading 1"> M\351thode de Runge-Kutta</Text-field>

La m\351thode de Runge-Kutta du quatri\350me ordre am\351liore davantage la m\351thode d'Euler am\351lior\351e mais les calculs sont encore beaucoup plus \351labor\351s. Cette m\351thode repose aussi en quelque sorte sur un calcul de moyenne arithm\351tique de pentes. Pr\351sentons, sans autres explications, la proc\351dure 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'\351quation diff\351rentielle NiMvKiYlI2R5RyIiIiUjZHhHISIiLCYqJiIiI0YmJSJ4R0YmRiglInlHRig=, la m\351thode de Heun montrait \340 la section pr\351c\351dente que les points calcul\351s \351pousaient d\351j\340, de mani\350re remarquable, la courbe solution. Le degr\351 de pr\351cision atteint avec cette m\351thode a pu \352tre appr\351ci\351 par une approche graphique \340 l'\351gard de la m\351thode d'Euler. Par contre, avec la m\351thode du Runge-Kutta, le degr\351 de pr\351cision atteint \340 l'\351gard de la m\351thode de Heun est difficilement appr\351ciable graphiquement (\340 l'\351cran). En effet, constatons que les trac\351s des courbes en reliant les points calcul\351s par les m\351thodes de Heun et de Runge-Kutta et la courbe solution se superposent presque parfaitement compte tenu de la r\351solution \351cran. Il est n\351cessaire de proc\351der au grosissant maximum de l'affichage (Ctrl-6), pour pouvoir \253 visualiser \273 l'augmentation de la pr\351cision.

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 num\351rique pourra permettre clairement une comparaison entre les degr\351s de pr\351cision atteint par chaque m\351thode \340 l'\351gard de la courbe solution.