TD MAPLE n°2 Équations différentielles
I- Syntaxe Maple c. dsolve
Maple sait résoudre pas mal d’équations différentielles grâce à
Nous allons découvrir aujourd’hui comment Maple résout les équations dif-
dsolve({equa diff, conditions initiales}, fonction(variable)).
férentielles et comment se débrouiller quand même Maple ne parvient pas
Par exemple,
à déterminer une solution explicite : voilà qui devrait éveiller l’attention des
brillant(e)s physicien(ne)s que vous êtes. > dsolve({diff(y(t),t,t)+y(t)=0,D(y)(0)=1,y(0)=1},y(t));
est sans surprise. On peut même omettre les conditions initiales
a. Rappel : fonction - expression > dsolve({diff(y(t),t,t)+y(t)=0},y(t));
Nous aurons besoin de faire la différence entre une fonction et une expres- Au passage, notons que Maple permet de résoudre des équations différen-
sion, car Maple la fait... tielles dépendant de paramètres
On peut transformer une expression en fonction à l’aide de > dsolve({R*diff(q(t),t)+(q(t))/C=E,q(0)=0},q(t));
unapply(expression, variable)
Testez par exemple Maple sur les équations suivantes
1. (x + 1)y ′ − x y = 0 4. x y ′ + 2y = x/(1 + x 2 )
b. diff - D p
2. (1 + x 2 )y ′ + x y = 3x 3 + 3x 5. y ′ 1 − x 2 − y 2 − 1 = 0
Nous n’emploierons pas la même syntaxe selon que l’objet à dériver est
3. (1 + x 2 )y ′ + x y = 1 6. y ′′ − 3y ′ + 2y = xe x + sin(x)
une fonction ou une expression. Dans le cas d’une expression, on emploie
diff(expression, varaible1, variable2, ...). Par exemple,
d. plot - odeplot
> F:=sin(4*x^3+5*x); diff(F,x);
L’instruction de base pour représenter graphiquement une fonction est
On peut dériver plusieurs fois plot(expression(variable), variable=mini..maxi,options
> diff(F,x,x,x); diff(F,x$3); graphiques);
comme nous l’avons déjà vu.
pour dériver une fonction, on utilise l’opérateur D(fonction). Notez bien
Ce qui nous intéressera plus aujourd’hui, c’est de tracer des lignes brisées
que d(fonction) est une fonction. Par exemple :
> plot([[0,0],[1,3],[2,-4],[3,2]]);
f:=unapply(F,x); D(f); D(f)(x); D(f)(32); (D3)(f); (D3)
(f)(x); Notez bien l’utilisation des crochets : on a construit une liste de points diffé-
rente de
TD Maple n°2 2
> plot([[0,0],[2,-4],[1,3],[3,2]]); Nous ne nous pencherons pas ici sur le problème d’existence de solutions :
On peut comparer deux tracés sur un même graphique en utilisant la fonc- nous supposerons que f a les propriétés qui font des solutions existent là où
tion display qu’il faut aller chercher dans la bibliothèque plots il faut...
Nous allons par exemple « observer » u à intervalle de temps régulier : notons
> with(plots): h cette période d’échantillonnage.
> d1:=plot([[0,0],[1,3],[2,-4],[3,2]], olor=blue): Nous allons essayer d’obtenir une suite d’approximations u n de u(t n ) et
> d2:=plot(sin(2*x),x=0..3,style=point): d’évaluer la pertinence de cette approximation, le principe étant résumé sur
> display(d1,d2); la figure 1 de la présente page
Pour terminer ce petit tour d’horizon, parlons de la fonction odeplot qui, u
elle aussi, fait partie de la bibliothèque plots. Il faut d’abord définir une so- Solution numérique
u n+1
lution d’une équation différentielle, par exemple un
Solution exacte
> s:=dsolve({(x+1)*diff(y(x),x)-x*y(x)=0,y(0)=1},y(x), u1
u(t n+1 )
numeri ): u(t n )
L’argument numeric ajouté à la fin impose à Maple de résoudre l’équation u0
numériquement. Pour comprendre ce qui se passe, demandez
O t1 tn t n+1 t
> s(32); h
Maple renvoie les coordonnées d’un point de la courbe solution. Pour tracer
F IGURE 1 – Principe de discrétisation d’une EDO
la courbe correspondante, il suffit d’entrer
> odeplot(s,[x,y(x)],-0.5..2); Deux grandes possibilités s’offrent à nous pour « se débarrasser » de u ′ :
– soit remplacer une courbe par sa tangente, c’est à dire utiliser une approxi-
Le tracé semble manquer de précision. On peut rajouter comme option
mation affine de u. On remplacera alors u ′ (t ) par une fonction dépendant
numpoints=1000 sachant que la valeur par défaut est 50.
de u et ainsi obtenir une relation de récurrence définissant la suite (u n ).
Vous voici armé(e) pour entamer les hostilités
– soit intégrer l’EDO : nous n’aurons plus de u ′ ... mais il faudra calculer
une intégrale ! Nous pourrons alors utiliser des méthodes d’approximation
II - Résolution d’équations différentielles d’intégrales : méthodes des rectangles, des trapèzes, etc.
par discrétisation b. Méthode d’Euler
a. Principe Découverte sur un exemple
On considère une équation différentielle ordinaire (EDO pour les intimes) Cette méthode, que vous avez dû étudier au lycée, a pour idée directrice de
u ′ (t ) = f (t , u(t )). Lorsqu’on ne connait pas de solution exacte à cette EDO, on remplacer localement une courbe par sa tangente. Il y a des méthodes plus
essaye d’en avoir une bonne approximation par des méthodes numériques. efficaces a , mais pour l’instant hors de notre portée mathématique.
a. Notamment celle de Runge-Kutta qu’utilise d’ailleurs Maple
Guillaume Connan, PCSI 2 - Lycée Clemenceau, 2008-2009
TD Maple n°2 3
Soit donc u une fonction dérivable sur un intervalle I = [a ; b] vérifiant u(0) = S:=[x,y]; // on place le point de depart
1 et , pour tout t , u ′ (t ) = α · u(t ) (oui, oui, c’est bon, moi aussi je sais qui for k from ... to N do
c’est, on n’a pas besoin de Maple pour ça, mais mieux vaut découvrir sur un x:=x+h; // on avan e d’un pas
exemple simple b ...).
– Soit t n un élément de ]a ; b[ tel que, pour h un réel strictement positif suffi- ...
samment « voisin » de 0, t n + h soit encore dans I. On peut alors considérer
que le quotient P(k):=[x,y]: // on definit le point P(k)
S:=S,P(k): // on le rajoute a notre liste
u(t n + h) − u(t n ) od:
h CourbeEuler:=plot([S]):
réalise une approximation d’autant plus valable de u ′ (t n ) que h est petit. SolMaple:=dsolve(...
Que devient alors notre EDO « approchée » ? En d’autres termes, trouver CourbeMaple:=odeplot(SolMaple,...):
une relation de récurrence définissant (u n ). display(CourbeEuler,CourbeMaple);
Quelle est la nature de la suite (u n ) ? end:
Est-on sûr(e) que cette suite réalise une bonne approximation de la solu- Regardez par exemple
tion de l’EDO ?
– Si vous avez des doutes, regardons des exemples. Pour se simplifier la vie, > eu(-4,6,0,3);
nous réaliserons un maillage régulier de I : nous le découperons en N inter- > eu(-4,24,0,3);
valles de longueur (b − a)/N. Construisez alors une procédure dépendant et commentez...
des paramètres α et N qui compare c les valeurs de u(t n ) et u n pour diffé-
rentes valeurs de n. Méthode d’Euler : généralisation
– Pour rendre encore plus visible ce phénomène, nous allons faire apparaître
sur un même graphique des courbes obtenues par la méthode d’Euler et la Oui, c’est bien beau, mais l’intérêt de l’informatique c’est de régler des cas
solution exacte. les plus généraux possibles.
u(t + h) − u(t )
On construira pour cela une procédure eu:= proc(alpha,N,a,b) qui res- Donc réfléchissons un peu ... u ′ (t ) = f (u(t ), t ) et u ′ (t ) ≈ ...
h
semblera à
for k from 0 to N do
> with(plots); x:=x+h;
>eu:=pro (alpha,N,a,b,u0) y:=y+f(y,x)*h;
lo al ...: P(k):=[x,y]:
h:=(b-a)/N; // le pas regulier S:=S,P(k):
y:=u0; // y=u(0) au depart od:
x:=a; // x=a au depart
ln 2
b. On estime l’intensité des radiations émises par un corps radioactif en mesurant la concentration u (t ) d’un isotope instable. Soit τ sa demi-vie. Alors u ′ (t ) = α · u (t ) avec α = −
τ
c. C’est à dire qui calcule dn = u (tn ) − un
Guillaume Connan, PCSI 2 - Lycée Clemenceau, 2008-2009
TD Maple n°2 4
Estimation de l’erreur Quand vous aurez étudié les formules de Taylor, vous pourrez montrer que
l’approximation obtenue est meilleure puisque d’ordre 2, alors que la mé-
C’est un point extrêmement important mais les outils théoriques nous
thode d’Euler donne une approximation d’ordre 1.
manquent pour l’instant pour en discuter
Modifiez une des procédures précédentes pour utiliser cette méthode. Atten-
tion ! Il y a trois niveaux de récurrence...
Sans filet
Le but de notre nouveau jeu va être de discrétiser une équation différentielle x:=a;
du second ordre, de tracer le graphe correspondant et de comparer le résultat ymoins:=u0;
obtenu avec les solutions proposées par Maple grâce à dsolve et odeplot. y:=ymoins+h*f(ymoins,x);
Nous utiliserons le fait que pour une fonction y deux fois dérivable sur un for k from 1 to N do
intervalle I et h proche de 0, on a x:=a+k*h;
yplus:=ymoins+2*h*f(y,x);
y(x + h) ≈ y(x) + y ′ (x) × h P(k):=[x+h,yplus]:
et donc S:=S,P(k):
ymoins:=y;
y ′ (x + h) ≈ y ′ (x) + y ′′ (x) × h y:=yplus;
od:
Il s’agit maintenant de mettre au point une procédure EquaDiff:=
proc(f,xmin,xmax,y0,yp0,n) où y0 = y(0), y p0 = y ′ (0) et f est la fonction
Observez :
de 3 variables définie par
Mil((u,t)->-2*u+t,200,0,3,1);
y ′′ = f (y ′ , y, x)
Mil((u,t)->-2*u+t,600,0,3,1);
Le résultat doit faire apparaître sur le même graphique la courbe obtenue par Mil((u,t)->4*u+t,200,0,3,1);
la méthode d’Euler et celle obtenue à l’aide de odeplot.
Il ne restera plus qu’à tester votre procédure sur les nombreuses équations
MORALITÉ : l’approche intuitive de la notion de limite ( dx ≈ ∆x) permet de
différentielles du second ordre que vous avez résolues en cours avec Madame
modéliser facilement certaines situations, mais si on pousse un peu, on ne
Gornet.
peut se passer d’une théorisation rigoureuse pour expliquer certaines ano-
malies. Analyse numérique ne veut donc pas dire bidouillage... L’étude de
c. Méthode du point milieu l’estimation des erreurs commises doit être menée avec le plus grand soin
Les anglo-saxons l’appellent schéma Leapfrog d mais est encore hors de notre portée.
En fait, pour une fonction suffisamment régulière sur l’intervalle [a ; b], l’idée Du point de vue informatique, un autre facteur doit être pris en compte : la
est d’approcher la pente entre les points de coordonnées t − h, u(t − h) et vitesse de convergence, pour ne pas avoir à attendre quelques siècles le ré-
¡ ¢
t + h, u(t + h) par u ′ (t ). sultat de la compilation...
¡ ¢
d. Pourquoi ce nom ?...
Guillaume Connan, PCSI 2 - Lycée Clemenceau, 2008-2009
TD Maple n°2 5
u
d. À partir de l’intégration numérique
Premier exemple u((k + 1)h) = u k+1
u(kh) = u k
Soit (E) l’équation différentielle u ′ (t )+3u(t ) = 2, u(0) = 1. Nous savons parfai-
tement résoudre cette équation. Nous allons encore voir une méthode qui va
vous paraître beaucoup plus compliquée, mais qui s’avérera fort utile quand t
il s’agira de résoudre des équations différentielles plus sophistiquées où les kh (k + 1)h
méthodes habituelles de résolution ne nous permettront plus de conclure
mais que vous rencontrerez au hasard de votre carrière d’ingénieur. F IGURE 2 – Méthode des trapèzes
Nous allons « observer » la solution de l’équation à intervalle de temps ré-
Il semble raisonnable de penser que plus h sera petit, plus l’aire sombre le
gulier : la période d’échantillonnage qu’on notera encore h. On fixera par
sera également et donc meilleure sera l’approximation.
exemple h = 0, 2 s.
Maintenant, en vous souvenant de la formule A (trapèze) =
Nous allons intégrer l’équation (E) sur une période d’échantillonnage quel- (petite base + grande base) × hauteur
, donnez l’aire du trapèze bleu.
conque, c’est à dire sur un intervalle kh, (k + 1)h , avec k un entier naturel.
£ ¤
2
Cela donne, grâce à la linéarité de l’intégration :
Équation aux différences
Z(k+1)h Z(k+1)h Z(k+1)h Déduisez-en que la suite (u n ) vérifie l’équation aux différences (F)
u (t ) dt + 3
′
u(t ) dt = 2 dt
kh kh kh 13 7 2
u n+1 − u n =
10 10 5
On définit la suite (u k )k∈N par u k = u(kh). Vous reconnaissez une suite arithmético-géométrique. Exprimez u n en fonc-
tion de n.
Calculez le premier terme du premier membre ainsi que le second membre
de l’équation : pas de problème.
R (k+1)h Équation différentielle vs équation aux différences
Plus compliqué va être de calculer kh u(t )dt car nous ne connaissons pas
u, donc encore moins une de ses primitives. Nous allons malgré tout nous en Pour vérifiez si la solution de l’équation aux différences n’est pas très éloi-
tirer en déterminant une approximation du résultat. gnée de la solution théorique, nous allons observer la solution de (E) toutes
les 0,2 s : on pose y n = f (nh). Donnez l’expression de y n .
On s’attend à ce que u n et y n soient très proches mais pas égaux. En effet, le
passage de (E) à (F) s’opère à l’aide d’une approximation (celle des trapèzes)
Méthode des trapèzes et non d’une équivalence.
On pose d n = y n − u n . À l’aide de Maple, dressez un tableau où figureront les
Observez la figure de la présente page 10 premiers termes des suites (u n ), (y n ) et (d n ). Cela vous rassure-t-il ?
Guillaume Connan, PCSI 2 - Lycée Clemenceau, 2008-2009
TD Maple n°2 6
Montrez que (d n ) est en fait la différence entre deux suites géométriques. Uti- On pose encore d n = y n − u n . Pouvez-vous conclure comme dans l’exemple
lisez ce résultat pour confirmer ou infirmer votre conjecture. précédent ? Moralité ?
Épilogue : les maths, c’est pas Harry Potter Comparaison avec le méthode d’Euler
Cette méthode est-elle plus précise que la méthode d’Euler vue précédem-
Tout serait donc magique en maths, ou bien n’avez-vous été confronté(e)
ment ?
toute votre vie qu’à des problèmes truqués ?
Soit (E′ ) l’équation différentielle 2u ′ (t ) − 3u(t ) = t , u(0) = 2. Vous pouvez de-
mander à Maple de vérifier que la solution de cette équation est la fonction s e. Généralisation
définie par
1 2 20 En dehors des équations linéaires, la généralisation devient extrêmement
u(t ) = − t − + e3/2t difficile... Nous nous contenterons donc des équations du style
3 9 9
Vous pourrez peut-être un jour montrer que la suite solution de l’équation
aux différences associée à (E) est u ′ (t ) = αu(t ) + βt + γ
20 43 n 1
µ ¶
2 Écrivez donc une procédure ayant pour arguments α, β et γ, les bornes a et b
un = − n− de l’intervalle, le nombre N de « mailles du filet » et la valeur de u(0).
9 37 30 9
Guillaume Connan, PCSI 2 - Lycée Clemenceau, 2008-2009