§ 1.2	Résolution numérique d'une équation différentielle ordinaire
	du premier ordre

Interprétation géométrique et résolution graphique

Dans le § 1.1, nous avons considéré une équation différentielle avec condition initiale de la forme

{({{v'=f(t,v)}, {v((Subscript[t, 0]))=(Subscript[v, 0])}})

Que faire lorsque la solution v ne nous est pas donnée ?

Pour distinguer la solution que nous allons construire de la solution donnée dans le § 1.1, utilisons un autre symbole et écrivons le système sous la forme

{({{y'=f(t,y)}, {y((Subscript[t, 0]))=(Subscript[y, 0])}})

où y(t) est la fonction cherchée.
A toute équation différentielle correspond un champ de directions associé à la fonction f. Pour l'exemple du § 1.1, voici la situation:

p=4.;r=1.;m0=8.;q=2.;Clear[f,t,y];
f[t_,y_]:=((p-r (y^2))/(m0-q t));
f[t,y]

tmin=0.;tmax=0.95(m0/q);ymin=0.;ymax=3.;

A tout point du plan, on peut associer le nombre réel

(t,y) ↦  f(t,y)

Chacun de ces nombres représente la pente que doit avoir la fonction  y(t)  au point  (t, y)

y'=f(t,y)

On représente la pente en chacun de ces points par un segment. Les vecteurs directeurs de ces segments sont

(({{1}, {f(t,y)}})).

L'ensemble de ces segments représente un champ de directions (ou champ de pentes) que l'on peut dessiner.

VectorPlot[
{1,f[t,y]},{t,tmin,tmax},{y,ymin,ymax},
VectorStyle->"Segment",VectorScale->{0.02, Automatic, None},
Frame->False,Axes->True,AxesLabel->{"t","y"},AspectRatio->Automatic
]

Le but est de construire une solution graphique de l'équation différentielle qui passe par le point ((Subscript[t, 0]),(Subscript[y, 0])). Il existe une et une seule courbe t ↦y(t) qui passe par le point donné et suit le champ de pentes, c'est-à-dire l'équation différentielle avec condition initiale possède une et une seule solution (voir § 1.1)

Pour montrer comment la solution dépend de la condition initiale, modifions-la.

{({{v'=f(t,y)}, {y(0)=3}})

On appelle solution générale l'ensemble des solutions de l'équation différentielle sans condition initiale

y'=f(t,y)

Off[Solve::"ifun"];
c[t_]={First[y[t]/.DSolve[{y'[t]==f[t,y[t]],y[0]==-1},y[t],t]],First[y[t]/.DSolve[{y'[t]==f[t,y[t]],y[0]==0},y[t],t]],First[y[t]/.DSolve[{y'[t]==f[t,y[t]],y[0]==1},y[t],t]],2,First[y[t]/.DSolve[{y'[t]==f[t,y[t]],y[0]==3},y[t],t]],First[y[t]/.DSolve[{y'[t]==f[t,y[t]],y[0]==4},y[t],t]]};
Plot[Evaluate[c[t]],{t,tmin,tmax},PlotRange->{0,3},AxesLabel->{"t\n[s]","y [m/s]"},ImageSize->{500,400}]
On[Solve::"ifun"];

La solution générale est une famille de fonctions. Poser une condition initiale équivaut à choisir une fonction de cette famille.

Méthode numérique d'Euler

Dans la méthode graphique, pour tracer la courbe, on part du point donné par la condition initiale, puis on trace la courbe de proche en proche en suivant les directions données par le champ de pentes.

Leonhard Euler est né à Bâle en 1707. Il reçut un enseignement mathématique dans la famille Bernoulli. Il mourut à Saint-Petersburg en 1783. Dans la méthode d'Euler, en partant d'un point et de la pente en ce point, on progresse un petit bout en ligne droite. On obtient ainsi un nouveau point à partir duquel on recommence l'opération. Finalement, on obtient une fonction affine par morceaux qui est une approximation de la solution cherchée.

Le point de départ est ((Subscript[t, 0]), (Subscript[y, 0])). En ce point, la pente est  y'((Subscript[t, 0]))=f((Subscript[t, 0]), (Subscript[y, 0])). On va effectuer "un petit pas" dans la direction (({{1}, {f((Subscript[t, 0]), (Subscript[y, 0]))}})).

On s'intéresse à la situation au temps (Subscript[t, 1])=(Subscript[t, 0])+h où h est un petit intervalle de temps donné. Effectuons le déplacement   h (({{1}, {f((Subscript[t, 0]),(Subscript[y, 0]))}})). La nouvelle position ((Subscript[t, 1]), (Subscript[y, 1])) sera telle que

(({{(Subscript[t, 1])-(Subscript[t, 0])}, {(Subscript[y, 1])-(Subscript[y, 0])}}))=h(({{1}, {f((Subscript[t, 0]),(Subscript[y, 0]))}}))

((Subscript[t, 1]),(Subscript[y, 1]))=((Subscript[t, 0])+h,(Subscript[y, 0])+h f((Subscript[t, 0]),(Subscript[y, 0])))

En répétant l'opération (voir Formulaires et tables)

((Subscript[t, 2]),(Subscript[y, 2]))=((Subscript[t, 1])+h,(Subscript[y, 1])+h f((Subscript[t, 1]),(Subscript[y, 1])))

((Subscript[t, k+1]),(Subscript[y, k+1]))=((Subscript[t, k])+h,(Subscript[y, k])+h f((Subscript[t, k]),(Subscript[y, k])))                    k ∈ ℕ

n=16;
h=((tmax-t0)/n);t0=tmin;y0=0;
euler[{t_,y_}]:={t+h,y+h f[t,y]};

Solution numérique d'Euler

solne=NestList[euler,{t0,y0},n]

Graphique de l'approximation d'Euler (fonction affine par morceaux)

SetOptions[ListPlot,AxesLabel->{"t","y"},PlotRange->{ymin,ymax},ImageSize->{500,300}];
ListPlot[solne]

ListPlot[solne,Joined->True]

Comparons la solution exacte (§ 1.1) et la solution numérique approchée, graphiquement et numériquement

Clear[v,t];
v[t_]:=((2 (-(t^2)+8 t))/((t^2)-8 t+32));
SetOptions[Plot,AxesLabel->{"t","y"},PlotRange->{ymin,ymax},ImageSize->{500,300}];
Plot[v[t],{t,tmin,tmax},Epilog->Point[solne]]

Liste des erreurs :                  (approximation d'Euler en (Subscript[t, i])) - (valeur exacte en (Subscript[t, i]))

{abscEuler,ordEuler}=Transpose[solne];
ordExact=Map[v,abscEuler];
ordErreur=ordEuler-ordExact;
Max[Abs[ordErreur]]

Si on prend une plus petite valeur de h, l'erreur tend à diminuer

n=128;
h=((tmax-t0)/n);
solne=NestList[euler,{t0,y0},n];
Plot[v[t],{t,tmin,tmax},Epilog->Point[solne]]

Liste des erreurs :                  (approximation d'Euler en (Subscript[t, i])) - (valeur exacte en (Subscript[t, i]))

{abscEuler,ordEuler}=Transpose[solne];
ordExact=Map[v,abscEuler];
ordErreur=ordEuler-ordExact;
Max[Abs[ordErreur]]

Lorsque h tend vers 0 c'est-à-dire pour n->∞, la méthode numérique d'Euler converge vers la solution exacte.

Autres méthodes numériques

Il existe des méthodes numériques plus performantes que la méthode d'Euler, mais leur étude n'est pas abordée ici. Sachez simplement que Mathematica fait appel à des méthodes plus élaborées qui convergent plus vite vers la solution et donnent des réponses plus précises.

Calcul numérique avec Mathematica

Soit à résoudre l'équation différentielle ordinaire du premier ordre avec condition initiale

?NDSolve

soln=NDSolve[{y'[t]==f[t,y[t]],y[t0]==y0},y,{t,tmin,tmax}]

Dans un premier temps, la méthode numérique produit une liste finie d'approximations de la forme {(Subscript[t, i]), (Subscript[y, i])}. Les valeurs intermédiaires sont ensuite calculées par interpolation.

Clear[sn];
sn=y/.soln[[1]]

Plot[sn[t],{t,tmin,tmax},PlotRange->{ymin,ymax},ImageSize->{500,300}]

Estimation numérique de l'erreur maximale:

Max[Table[Abs[sn[t]-v[t]],{t,tmin,tmax,((tmax-tmin)/100)}]]

Travaux dirigés du § 1.2

1.2 - TD 1

On donne l'équation différentielle avec condition initiale

{({{y'=1+(((y-t))^2)}, {y(0)=(1/2)}})

a) 	Résolution graphique
	Dans le rectangle  [0; 1]×[0; 2.5], calculez des valeurs numériques de la fonction

(t,y)↦ 1+(((y-t))^2)

avec Mathematica puis, à la main, représentez le champ de pentes.
	Tracez la solution  t↦y(t) qui vérifie la condition initiale.

b) 	Méthode d'Euler sans ordinateur
	Calculez une solution numérique approchée de l'équation différentielle par la méthode d'Euler
	sur l'intervalle [0; 1], avec le pas h=0.1

c) 	Méthode d'Euler avec ordinateur
	Calculez une solution numérique approchée de l'équation différentielle par la méthode d'Euler
	sur l'intervalle [0; 1], avec le pas h=0.01.

d) 	Résolution numérique avec Mathematica
	Calculez une solution numérique approchée de l'équation différentielle sur l'intervalle [0; 1]
	avec la méthode numérique de Mathematica.

e) 	Vérification
	Montrez que la fonction

y(t)=t+(1/(2-t))

est la solution de l'équation différentielle avec condition initiale.

f) 	Comparaison
	Comparez les cinq solutions qui apparaissent sous les points a, b, c, d et e.

g) 	Calcul de l'erreur d'approximation
	Pour les solutions c et d, calculez l'erreur d'approximation aux abscisses
	t=0, 0.05, 0.1, 0.15, ..., 0.95, 1.

1.2 - TD 2

a) 	Méthode de Heun		[avec Mathematica]
	Pour résoudre les équations différentielles, une autre méthode numérique est décrite dans les Formulaires et
	tables: la méthode de Heun. Utilisez cette méthode pour résoudre sur l'intervalle [0, 2] l'équation suivante avec
	condition initiale

{({{y'=((4 t)/(y+1))}, {y(0)=1}})

b) 	Méthode Runge-Kutta 4		[facultatif]
	Pour résoudre les équations différentielles, une autre méthode numérique est décrite dans les Formulaires et
	tables: la méthode Runge-Kutta 4. Utilisez cette méthode pour résoudre l'équation différentielle précédente
	sur l'intervalle [0, 2].

1.2 - TD 3           [facultatif]

Champ de pentes		[avec Mathematica]
	Dessinez le champ de pentes pour résoudre graphiquement l'équation suivante sur l'intervalle [0, 2]

{({{y'=((4 t)/(y+1))}, {y(0)=1}})