Discrétisation d'une équation différentielle

Pour contrôler la démarche, on prend une équation différentielle que l'on sait résoudre ! Soit par exemple : [Graphics:Images/discretiser_gr_1.gif]

[Graphics:Images/discretiser_gr_2.gif]
[Graphics:Images/discretiser_gr_3.gif]
vrai = Plot[solu, {x,0,2}, AxesOrigin->{0,0},
    Ticks->{{0,1,2},{0, 0.1, 0.2, 0.3}},
    PlotStyle->Hue[0],
    PlotRange->All];

[Graphics:Images/discretiser_gr_4.gif]

Le principe de la discrétisation est de remplacer un phénomène continu (ici une fonction y[x]) par un autre, discret (ici on va prendre une suite [Graphics:Images/discretiser_gr_5.gif]).
L'idée est que l'on choisit un pas δ— mettons 1/20 — et que [Graphics:Images/discretiser_gr_6.gif].
Il faut alors exprimer [Graphics:Images/discretiser_gr_7.gif] en termes de [Graphics:Images/discretiser_gr_8.gif].

La méthode d'Euler.

On prend [Graphics:Images/discretiser_gr_9.gif] ce qui donne la relation
[Graphics:Images/discretiser_gr_10.gif] ou encore [Graphics:Images/discretiser_gr_11.gif].

Implémentation

δ = 0.05;
Clear[u];
u[0] = 0.;
u[n_] := u[n] = (1-δ) u[n-1] + E^(-(n-1) δ) δ
euler = ListPlot[Array[u, {40}] ,
        PlotJoined->True];

[Graphics:Images/discretiser_gr_12.gif]

Show[GraphicsArray[{vrai, euler}] ];

[Graphics:Images/discretiser_gr_13.gif]

Ça a l'air nickel... mais il ne faut pas pousser trop loin:

solu/.x->20.
vrai1 = Plot[solu, {x,18,20},
    PlotStyle->Hue[0],
    AxesOrigin->{18,0},
    PlotRange->All];
[Graphics:Images/discretiser_gr_14.gif]

[Graphics:Images/discretiser_gr_15.gif]

u[200];
u[400]
euler1 = ListPlot[Table[u[n],
            {n, 360, 400}] ,
        PlotJoined->True];
[Graphics:Images/discretiser_gr_16.gif]

10% d'erreur donc…

[Graphics:Images/discretiser_gr_17.gif]

On peut en fait calculer exactement [Graphics:Images/discretiser_gr_18.gif]:

Needs["DiscreteMath`RSolve`"]
Clear[u, δ];
δ = 0.05;
u[n]/. First[ RSolve[{u[0] == 0.,
     u[n] == (1-δ) u[n-1] + δ E^(-(n-1) δ)}, u[n], n,
     
     Method->MethodEGF]//Simplify ]
     //Chop//Factor
[Graphics:Images/discretiser_gr_19.gif]

Au lieu d'une exponentielle polynôme on trouve la différence de deux exponentielles !

C'est plus ennuyeux avec une équation du second ordre, ou un système:

[Graphics:Images/discretiser_gr_20.gif]
[Graphics:Images/discretiser_gr_21.gif]
[Graphics:Images/discretiser_gr_22.gif]
[Graphics:Images/discretiser_gr_23.gif]
[Graphics:Images/discretiser_gr_24.gif]

Normalement on s'attend à des cercles:

soluces = {x[t], y[t]}/.
    First[DSolve[eq, {x[t], y[t]}, t] ]
[Graphics:Images/discretiser_gr_25.gif]
ParametricPlot[soluces, {t,0,20}, 
    AspectRatio->Automatic];

[Graphics:Images/discretiser_gr_26.gif]

Discrétisation par la méthode d'Euler

On a comme plus haut, [Graphics:Images/discretiser_gr_27.gif] et [Graphics:Images/discretiser_gr_28.gif] d'où

δ = 0.05;
Clear[u, v];
u[0] = 1; v[0] = 0;

u[n_] := u[n] = δ v[n-1]  +  u[n-1]
v[n_] := v[n] = - δ u[n-1] + v[n-1]

eulerSinus = ListPlot[Table[{u[n], v[n]},
     {n,500}] , PlotJoined->True];

[Graphics:Images/discretiser_gr_29.gif]

Gloups !!! !!!!!!&9&9&9&9&9&9&9

Heureusement il existe des méthodes de discrétisation plus efficaces, comme celles utilisées par la procédure Mathematica NDSolve:

solucesAméliorées = {x[t], y[t]}/.First[NDSolve[eq, {x[t], y[t]}, 
                            {t, 0, 20}] ]
[Graphics:Images/discretiser_gr_30.gif]
ParametricPlot[solucesAméliorées, {t,0,20},
                AspectRatio->Automatic];

[Graphics:Images/discretiser_gr_31.gif]


Converted by Mathematica      September 8, 2002