Equations différentielles ordinaires

1. Equations différentielles ordinaires Généralités * EDO  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. * 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. On considérera les problèmes dits de Cauchy, ayant la forme : 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é: Alors la solution y = y(t) existe et est unique. (Théorème de Cauchy-Lipschitz) M=4 intervalles, mais M+1=5 nœuds t0 t1 t2 t3 t4=T

2. 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

3. 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 : En remplaçant « y » par la solution numérique « u » on a la méthode explicite d’Euler : L’approximation rectangulaire avec fn+1, donne la méthode implicite d’Euler : Cette dernière méthode ne permet pas d’obtenir un+1 directement, mais seulement après la résolution de l’équation 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 « Equations Non-Linéaires »).

4. Equations différentielles ordinaires Méthodes d’Euler: un exemple 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); 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) (E1) 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 g=@(X) -X + u(n) + h*f(t(n+1), X); u(n+1)= fzero(g, u(n)); end plot(t, u ) (E2)

5. 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: Ce qui constitue la méthode de Crank-Nicolson (méthode implicite). 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): ou encore la formule implicite à 3 pas d’Adams-Moulton:

6. 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: Utilisons (E1) pour produire une « prédiction » sur un+1, que l’on injecte dans (CN1) pour « corriger » : Les 2 phases (PC1) et (PC2) constitue une méthode dite prédicteur-correcteur (ici la méthode est dite « de Heun »).

7. 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) * 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. Erreur de troncature ou erreur de consistance  = facteur d’amplification

8. Equations différentielles ordinaires Propagation de l’erreur, consistance d’une méthode .e1 e2 u1 (h2) e1 y1 u0 = y0

9. 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: génère 40 valeurs de la solution u, réparties entre t=0 et t=3 secondes. f=@(t, u) (u-t)./(u+4*t); tspan = [0, 3]; y0=1; [t, u] = ode45(f, tspan, y0); plot(t, u) Pour utiliser ODE45, la fonction f doit avoir la forme dy=f(t, y)

10. 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... … 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)

11. 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:

12. 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: Alternativement, on écrit la fonction dans un fichier: et on l’appelle avec ode45: >> 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‘) function dw=fpendule(t, w) g=9.81; L=1; dw=[w(2); -g/L*sin(w(1))]; end >> ode45(@fpendule, [0,10], [pi/4, 0]) Pour utiliser ODE45, la fonction f doit avoir la forme dy=f(t, y)

13. 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)

14. 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.

15. 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.

16. 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:

17. 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

18. Equations différentielles ordinaires Questions de cours 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.