Equations différentielles ordinaires

Download Report

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]:
[email protected](t, u) (u-t)./(u+4*t);
(E1)
Méthodes dโ€™Euler: un exemple
๐‘ฆโ€ฒ
๐‘ฆโˆ’๐‘ก
๐‘ก =
๐‘ฆ + 4๐‘ก
๐‘ฆ ๐‘ก0 = 1
[email protected](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:
% [email protected](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
[email protected](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
[email protected](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:
[email protected](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;
>> [email protected](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 ».
[email protected](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