Transcript Equations différentielles ordinaires

Equations différentielles ordinaires
Généralités
* EDO  𝐹 𝑡, 𝑡 , 𝑦 ′ 𝑡 , … , 𝑦 𝑛 𝑡 = 0 où F est connue, y(t) est la fonction inconnue
de l’unique variable indépendante t (ou x, …).
* L’ordre n de l’équation est celui de la dérivée de plus haut degré.
Exemple: (y-t) / (y+t) - y’ = 0 est une EDO du premier ordre.
On considérera les problèmes dits de Cauchy, ayant la forme :
𝑦 ′ 𝑡 = 𝑓(𝑡, 𝑦 𝑡 )
𝑦 𝑡0 = 𝑦0
Avec t  [t0, T] et où y0 est appelée ‘condition initiale’ (CI)
Si la fonction f(t, y) est continue par rapport à ses 2 variables et possède la propriété :
∃ 𝐿 ∈ ℝ+ 𝑡𝑒𝑙 𝑞𝑢𝑒 𝑓 𝑡, 𝑦1 − 𝑓 𝑡, 𝑦2 ≤ 𝐿. 𝑦1 − 𝑦2 , ∀𝑦1 , 𝑦2 ∈ ℝ
Alors la solution y = y(t) existe et est unique.
(Théorème de Cauchy-Lipschitz)
* Les EDO d’ordre n>1, peuvent se ramener à un système de n EDO d’ordre 1.
* Ces solutions sont rarement résolues analytiquement, on utilise donc des méthodes
numériques pour approcher ces solutions. Le domaine sur lequel la solution est à
déterminer doit être limité: t  [t0,T] et divisé en M intervalles de longueur h=(T-t0)/M
appelé pas de discrétisation. Pour chaque nœud ti, i=0, 1, 2, …, M on cherche la valeur
ui qui approche yi = y(ti). L’ensemble {u0 = y0, u1, ..., uM} est la solution numérique.
M=4 intervalles, mais
M+1=5 nœuds
t0
t1
t2
t3
t4=T
1
Equations différentielles ordinaires
Condition initiale (CI), unicité
Toutes ces courbes
ont la même dérivée
y’ = f(t, y)
Les CI assurent l’unicité de la solution
Si la CI est à l’intérieur, résoudre en 2 fois :
CIgauche, puis CIdroite
2
Equations différentielles ordinaires
Méthodes d’Euler
On peut obtenir la méthode d’Euler en intégrant l’EDO du
problème de Cauchy sur une étape (tn  tn+1) :
𝑦 ′ 𝑡 = 𝑓 𝑡, 𝑦 𝑡
On a :
𝑡𝑛+1
𝑑𝑦 = 𝑦 𝑡𝑛+1 − 𝑦 𝑡𝑛 = 𝑦𝑛+1 − 𝑦𝑛
𝑡𝑛
𝑡𝑛+1
𝑡𝑛+1
𝑑𝑦 =
𝑡𝑛
⟹ 𝑑𝑦 = 𝑓 𝑡, 𝑦 𝑡 . 𝑑𝑡
𝑓 𝑡, 𝑦 𝑡 . 𝑑𝑡 ≈ 𝑡𝑛+1 − 𝑡𝑛 . 𝑓 𝑡𝑛 , 𝑦𝑛 = ℎ. 𝑓 𝑡𝑛 , 𝑦𝑛
𝑡𝑛
En remplaçant « y » par la solution numérique « u » on a la méthode explicite d’Euler :
𝑢𝑛+1 − 𝑢𝑛 = ℎ. 𝑓 𝑡𝑛 , 𝑢𝑛 ⟹ 𝑢𝑛+1 = 𝑢𝑛 + ℎ. 𝑓𝑛
(E1)
L’approximation rectangulaire avec fn+1, donne la méthode implicite d’Euler :
𝑢𝑛+1 − 𝑢𝑛 = ℎ. 𝑓 𝑡𝑛+1 , 𝑢𝑛+1 ⟹ 𝑢𝑛+1 = 𝑢𝑛 + ℎ. 𝑓𝑛+1
(E2)
Cette dernière méthode ne permet pas d’obtenir un+1 directement, mais seulement après
la résolution de l’équation 𝑢𝑛+1 = 𝑢𝑛 + ℎ. 𝑓 𝑡𝑛+1 , 𝑢𝑛+1 où un+1 est l’inconnue. On peut
utiliser les fonctions fzero ou fsolve ou une méthode de point fixe à cette fin(voir chapitre
3
« Equations Non-Linéaires » ).
Equations différentielles ordinaires
Soit le problème de Cauchy ci-contre. On cherche à
le résoudre numériquement sur le domaine [0, 3]:
f=@(t, u) (u-t)./(u+4*t);
(E1)
Méthodes d’Euler: un exemple
𝑦′
𝑦−𝑡
𝑡 =
𝑦 + 4𝑡
𝑦 𝑡0 = 1
f=@(t, u) (u-t)./(u+4*t);
% u(n+1)=u(n) + h*f(t(n+1), u(n+1));
% ==>
% 0 = -u(n+1) + u(n) + h*f(t(n+1), u(n+1))
% u(n+1) est l'inconnue X de l'équation
% g(X) = 0 avec la fonction g:
% g=@(X) -X + u(n) + h*f(t(n+1), X);
h=1e-2;
t=0:h:3;
u(1)=1;
for n=1:length(t)-1
u(n+1)=u(n) + h*f(t(n), u(n));
end
plot(t, u)
h=1e-2;
t=0:h:3;
u(1)=1;
1.4
1.35
1.3
for n=1:length(t)-1
g=@(X) -X + u(n) + h*f(t(n+1), X);
u(n+1)= fzero(g, u(n));
end
plot(t, u )
1.25
1.2
1.15
1.1
1.05
1
0
0.5
1
1.5
2
2.5
3
(E2)
Equations différentielles ordinaires
Méthodes de Crank-Nicolson,
Runge-Kutta, multi-pas
On utilise maintenant une approximation trapézoïdale de l’aire sous la courbe f(t, y). On
obtient alors:
1
𝑢𝑛+1 − 𝑢𝑛 = ℎ. 𝑓 𝑡𝑛 , 𝑢𝑛 − ℎ. 𝑓 𝑡𝑛+1 , 𝑢𝑛+1 −𝑓 𝑡𝑛 , 𝑢𝑛
2
Ce qui constitue la méthode de Crank-Nicolson (méthode implicite).
1
𝑢𝑛+1 = 𝑢𝑛 + 2 ℎ. 𝑓𝑛+1 + 𝑓𝑛
(CN1)
En utilisant des approximations plus fines de quadrature (comme par exemple la formule
de quadrature de Simpson qui utilise le point milieu tn+1/2), on peut générer des schémas
plus sophistiqués donnant une plus grande précision dans la solution. Les méthodes de
Runge-Kutta sont basées sur ce principe et nécessitent l’évaluation de plusieurs valeurs de
la fonction f(t, y) sur chaque intervalle [tn, tn+1]. Toutes les méthodes vues ci-dessus sont
des méthodes dites « à un pas ».
Les méthodes multi-pas utilisent la même approche (quadrature) mais en dehors de
l’intervalle [tn,tn+1]. On a par exemple la formule d’Adams-Bashforth (explicite):
1
𝑢𝑛+1 = 𝑢𝑛 +
ℎ. 23𝑓𝑛 − 16𝑓𝑛−1 + 5𝑓𝑛−2
12
ou encore la formule implicite
1
à 3 pas d’Adams-Moulton:
𝑢𝑛+1 = 𝑢𝑛 + 24 ℎ. 9𝑓𝑛+1 + 19𝑓𝑛 − 5𝑓𝑛−1 + 𝑓𝑛−2
Equations différentielles ordinaires
Méthodes prédicteur-correcteur
Les méthodes implicites nécessitent pour chacune de leurs étapes la résolution d’une
équation où un+1 est l’inconnue.
On peut aussi rendre la méthode explicite en utilisant la méthode d’Euler (E1). Par
exemple, considérons la méthode de Crank-Nicolson, on a:
1
𝑢𝑛+1 = 𝑢𝑛 + ℎ. 𝑓 𝑡𝑛+1 , 𝑢𝑛+1 +𝑓 𝑡𝑛 , 𝑢𝑛
(𝐶𝑁1)
2
Utilisons (E1) pour produire une « prédiction » sur un+1,
𝑢∗ 𝑛+1 = 𝑢𝑛 + ℎ. 𝑓𝑛
(PC1)
que l’on injecte dans (CN1) pour « corriger » :
1
𝑢𝑛+1 = 𝑢𝑛 + ℎ. 𝑓 𝑡𝑛+1 , 𝑢∗ 𝑛+1 +𝑓 𝑡𝑛 , 𝑢𝑛
2
(𝑃𝐶2)
Les 2 phases (PC1) et (PC2) constitue une méthode dite prédicteur-correcteur (ici la
méthode est dite « de Heun »).
Equations différentielles ordinaires
Propagation de l’erreur,
consistance d’une méthode
Soit en+1 = yn+1 – un+1, l’erreur à l’étape n+1. On a d’une part en utilisant Taylor:
yn+1 = yn + h.y’n + (h2) = yn + h.f(tn, yn) + (h2)
et d’autre part (dans le cas de E1):
un+1 = un + h.f(tn, un)
soit donc:
en+1 = (yn - un )+ h.f(tn, yn) - h.f(tn, un) + (h2)
en+1 = en + h.{f(tn, yn) - f(tn, un)} + (h2)
en+1 = en + h.{f(tn, yn) - f(tn, yn- en)} + (h2)
en+1 = en + h.en.{f(tn, yn) - f(tn, yn- en)} / en + (h2)
en+1 = en + h. en.[f/y ]n + (h2)
en+1 = en { 1 + h.[f/y ]n } + (h2)
 = facteur d’amplification
Erreur de troncature ou erreur de consistance
* L’erreur se propage d’un pas au suivant, le facteur d’amplification doit être < 1 pour que
la solution numérique converge ce qui donne une condition sur h.
* On dit qu’une méthode est consistante si l’erreur de consistance est un infiniment petit
en h.
Equations différentielles ordinaires
Propagation de l’erreur,
consistance d’une méthode
.e1
u1
(h2)
e1
y1
u0 = y 0
e2
Equations différentielles ordinaires
ode45
La fonction ode45 de Matlab est un « solveur » d’EDO (ODE en anglais). Son algorithme
est basée sur une méthode du type Runge-Kutta explicite. Elle possède typiquement 2
sorties et 3 entrées:
[t, u] = ode45(f, tspan, y0)
avec:
u la solution numérique calculée aux temps contenus dans le vecteur t.
f est la poignée de la fonction f(t, y) du problème de Cauchy concerné.
tspan est un vecteur du type [tinitial, tfinal] indiquant à ode45 le domaine sur lequel
chercher la solution. Enfin y0 est la condition initiale. Une 4 ième entrée possible est
‘options’ que l’on utilise en conjonction avec la fonction ‘odeset’ afin par exemple de
régler la précision de la solution numérique.
Exemple:
1.4
1.35
1.3
f=@(t, u) (u-t)./(u+4*t);
tspan = [0, 3];
y0=1;
[t, u] = ode45(f, tspan, y0);
plot(t, u)
𝑦−𝑡
′
𝑦 𝑡 =
𝑦 + 4𝑡
𝑦 𝑡0 = 1
1.25
1.2
1.15
1.1
1.05
1
0
0.5
1
1.5
génère 40 valeurs de la solution u, réparties entre t=0 et t=3 secondes.
Pour utiliser ODE45, la fonction f doit avoir la forme dy=f(t, y)
2
2.5
3
Equations différentielles ordinaires
Système d’EDO
Exemple d’un système d’EDO où x et y sont des fonctions
d’une variable indépendante, par exemple du temps. Les
conditions initiales sont précisées.
On peut appliquer une méthode
du type (E1)...
... ou plus pratique du point de vue programmation,
écrire le système sous forme vectorielle...
5
… pour par exemple, utiliser ode45:
F=@(t,U) [U(1)-4*U(2); U(1)+U(2)];
[tsol, usol] = ode45(F, [0,1], [1;2]);
plot(tsol, usol(:,1), tsol, usol(:,2))
Pour utiliser ODE45, la fonction f doit avoir la forme dy=f(t, y)
0
-5
-10
-15
0
0.2
0.4
0.6
0.8
1
Equations différentielles ordinaires
EDO d’ordre > 1
Une EDO d’ordre m est une relation donnant la dérivée m-ienne en fonction du temps t,
de la variable inconnue y et de ses dérivées d’ordre inférieur. Une solution unique est
possible en connaissant m conditions initiales:
On peut la transformer en un système de m EDO du premier degré en posant:
et en dérivant ces nouvelles variables pour
obtenir une forme de Cauchy W’= F(t, W):
Les conditions initiales étant:
Equations différentielles ordinaires
EDO d’ordre > 1 : exemple
Soit un pendule simple de longueur L=1m lâché sans vitesse initiale
d’un angle y0=pi/4 compté par rapport à la verticale. Son mouvement
est régi par l’équation différentielle (g=9.81 m.s-2) du 2ième ordre:
On pose w1=y; w2=y’; on dérive  w’1= w2;
w’2= -g/L*sin(w1);
On programme:
>> g=9.81; L=1;
>> f=@(t, w) [w(2); -g/L*sin(w(1))];
>> [tsol, wsol]=ode45(f, [0,10], [pi/4, 0]);
>> plot(tsol, wsol(:,1), '-or', tsol, wsol(:,2), '-s‘)
2.5
2
1.5
1
0.5
Alternativement, on écrit la fonction dans un fichier:
function dw=fpendule(t, w)
g=9.81; L=1;
dw=[w(2); -g/L*sin(w(1))];
end
et on l’appelle avec ode45:
>> ode45(@fpendule, [0,10], [pi/4, 0])
Pour utiliser ODE45, la fonction f doit avoir la forme dy=f(t, y)
0
-0.5
-1
-1.5
-2
-2.5
0
1
2
3
4
5
6
7
8
9
10
Equations différentielles ordinaires
Exercice: convergence conditionnelle
Expérimenter sur la valeur du pas h de la méthode E1 ci-dessous et constater
qu’au dessus d’un pas critique, la solution numérique « diverge ».
f=@(t, u) (u-t)./(u+4*t);
h=1e-2;
t=0:h:3;
u(1)=1;
for n=1:length(t)-1
u(n+1)=u(n) + h*f(t(n), u(n));
end
plot(t, u)
Equations différentielles ordinaires
Exercice: stabilité absolue, conditionnelle
On appelle « problème modèle » l’EDO ci-contre dont la
solution exacte est:
Ce problème est utilisé pour définir la notion de stabilité
d’une méthode:
Une méthode est absolument stable si sa solution
numérique pour le problème modèle est telle que:
Il peut exister une condition sur le pas de discrétisation h pour que la méthode soit
stable. On parle alors de stabilité conditionnelle. Une méthode peut être instable.
1/ Montrer que E1 est stable à la condition :
2/ Montrer que E2 et CN1 possèdent une stabilité inconditionnelle.
Equations différentielles ordinaires
Exercice: méthode « mixte »
Que peut-on dire de la méthode:
lorsque =0 ? =1 ? =0.5 ?
Exercice: mouvement d’un mobile
Soit un mobile de coordonnées x(t), y(t) telles que x’=x-4y et y’=x+y. A l’instant
initial t=0, le mobile se trouve à la position (2, 3). Représenter sa trajectoire sur
l’intervalle t=[0, 3].
Exercice: pendule dans une gravitation évanescente
Résoudre le problème du pendule en supposant que l’accélération de la pesanteur
diminue après l’instant initial selon la loi g(t) = g0.2-t avec g0 = 9.81 m.s-2.
Même expérience si g(t) = g0.2t.
Indice: utiliser ode45, produire les graphiques y=y(t) où y est l’angle du pendule.
Equations différentielles ordinaires
Exercice: EDO 1ier ordre, 2ième ordre
1/ Pour t=[0, 10] résoudre numériquement l'équation
différentielle du 1er ordre suivante:
La solution exacte de ce système est :
Représenter graphiquement yExact et ynumérique
2/ Résoudre l’équation de van der Pol (prendre t=[0, 10], µ=1, y(0)=2, y’(0)=0) par
une méthode E1, puis comparer avec le résultat retourné par ode45:
Equations différentielles ordinaires
Exercice: trajectoire d’un satellite
On désire tracer la trajectoire d’un satellite artificiel de masse m autour de la Lune de
masse M. La lune est placée à l’origine d’un référentiel absolu où le satellite possède les
coordonnées x, y et on notera:
On applique le principe fondamental de la dynamique au satellite pour obtenir la force
F exercée par la Lune sur le satellite:
D’autre part, on a
Résoudre ce problème en expérimentant différentes conditions initiales. Produire les
graphiques de la trajectoires.
m=103 kg, M=7 1022 kg, K=6.67 10-11 N.m2.kg-2
17
Equations différentielles ordinaires
1/ Donner la forme du problème de Cauchy.
2/ Donner la méthode d’Euler explicite.
3/ Donner la méthode d’Euler implicite.
4/ Donner la méthode de Crank-Nicholson.
Questions de cours