Introduction aux équations différentielles ordinaires (EDO)

1. Introduction aux équations différentielles ordinaires (EDO) E. Grenier

2. Exemple: dynamique des populations

3. I.1. Dynamique des populations • Population de N(t) individus au temps t • Décrire l’évolution de cette population: • Hypothèses sur la reproduction / prédation • Mise en équation • Simulations numériques • Discussion • Evolution de N ∂t N = naissances – décès • Le modèle le plus simple (Euler) • Le nombre de naissances est proportionnel à la population Naissances = a N • Le nombre de décès est proportionnel à la population Décès = b N • Bilan ∂t N = ( a – b ) N • Résolution N(t) = N(0) exp( (a-b) t).

4. I.1. Dynamique des populations: Euler • Résolution N(t) = N(0) exp( (a-b) t). • Discussion • a > b : natalité plus importante que la mortalité: croissance exponentielle de la population. • a < b: natalité plus faible que la mortalité: décroissance exponentielle de la population. • Discussion du modèle: • Simple à mettre en équations et à résoudre • La croissance exponentielle n’est pas réaliste: limitations dues au milieu ambiant • Pour être plus réaliste il faut faire dépendre a et b de N … • Changement d’hypothèses de mises en équations pour avoir un comportement plus réaliste. • Hypothèses à ajouter: limitation de la croissance dans un milieu fini.

5. I.1. Dynamique des populations: modèle logistique • Le modèle logistique (Verhulst, 1836): les nouvelles hypothèses • Hypothèse de milieu limité: le milieu peut nourrir K individus • Si N < K: il y a suffisamment de ressources: la population augmente car la natalité est supérieure à la mortalité. • Si N > K: il n’y a pas assez de ressources: des individus meurent de faim. La mortalité devient supérieure à la natalité. • Si N << K: cas d’Euler: croissance proportionnelle à N. • Mise en équations: ∂t N = f(N) où f(N) > 0 si N < K f(N) < 0 si N > K. f(N) ~ c N si N << K De plus, il n’y a pas de création spontanée d’individus: f(0) = 0. • Choix de f: le f le plus simple satisfaisant ces hypothèses est f(N) = r N (1 – N/K) où r est une constante.

6. I.1. Dynamique des populations: modèle logistique • Equation logistique: ∂t N = f(N) = r N (1 – N/K) • Discussion: • Parfois f(N) est donné par des mesures expérimentales. • Si on connaît bien la dynamique de reproduction / mort, on peut parfois en déduire f, mais il faut pour cela des hypothèses supplémentaires • La fonction proposée est la plus simple qui fonctionne • Signification des constantes • K: capacité du milieu: nombre d’individus qu’il peut nourrir. • r: • c’est l’inverse d’un temps. • Si N << K, f(N) ~ rN donc N(t) ~ N(0) exp(r t) • C’est la vitesse de croissance de la population quand N << K.

7. I.1. Dynamique des populations: modèle logistique • Equation logistique: ∂t N = f(N) = r N (1 – N/K) • Les solutions analytiques existent ….. • Allures des solutions: p2

8. I.1. Dynamique des populations: équilibre stable / instable • Equation logistique: ∂t N = f(N) = r N (1 – N/K) • Etat d’équilibre: N* tel que f(N*) = 0. • Notion de stabilité: • Équilibre stable: après une petite perturbation, le système revient à N* • Equilibre instable: une petite perturbation déstabilise le système. • Equation des petites perturbations: N(t) = N* + u(t) ∂t u = f(N*+u) = f ’(N*) u + O(u^2) soit en négligeant les termes de taille u^2 ∂t u = f ’(N*) u • Discussion: • f ’ (N*) > 0 : croissance de u : N* est un équilibre instable • f ’ (N*) < 0: décroissance exponentielle de u: N* est un équilibre stable. • Application à la logistique: deux états d’équilibre: 0 et K • f ’ (0) > 0 donc 0 est instable • f ’ (K) < 0 donc K est stable

9. I.2. Dynamique des populations: modèles avec prédation • Ajout d’un phénomène : la prédation. • Exemple: vers dans des arbres, mangés par des oiseaux. • Modélisation: • Modèle simple: pas d’équation sur le nombre de prédateurs • P(n) nombre d’individus morts par prédation par unité de temps ∂t N = f(N) – P(N) • Hypothèses sur la prédation: un premier exemple: • La prédation est proportionnelle au nombre de vers • Sauf s’il y a beaucoup de vers: effet de saturation: les oiseaux se gênent entre eux P(N) = BN / (A + N) • Système obtenu ∂t N = r N (1 – N/K) – BN / (A+N)

10. I.2. Dynamique des populations: modèles avec prédation • Modélisation: • Hypothèses sur la prédation: un second exemple: • La prédation est proportionnelle au nombre de vers • Sauf s’il y a beaucoup de vers: effet de saturation: les oiseaux se gênent entre eux • Sauf s’il y a trop peu de vers: les oiseaux ne se déplacent pas: P(N) = BN^2 / (A^2 + N^2) • Système obtenu ∂t N = r N (1 – N/K) – BN^2 / (A^2+N^2)

11. I.2. Dynamique des populations: modèles avec prédation • Etats d’équilibre du second modèle: r N* (1 – N*/K) – BN*^2 / (A^2+N*^2)=0 • soit N* = 0 • Soit r (1 – N*/K) (A^2 + N*^2) – BN* = 0 qui est une équation polynomiale de degré 3 qui a • Soit trois racines réelles • Soit une racine réelle et deux racines complexes conjuguées.

12. II.1. Populations en interaction: Lotka Volterra • Deux populations: une de proies, une de prédateurs • N(t) nombre de proies, P(t) nombre de prédateurs • Le modèle de Lotka Volterra: • Naissances des proies proportionnelles à N • Morts par prédation: proportionnel à N et à P • Naissances de prédateurs: proportionne à P et à N • Mort de prédateur proportionnelle à P (mort naturelle). • Mise en équations ∂t N = a N – b N P ∂t P = c N P – d P • Remarque: pas de limitation par le milieu nourrisant les proies (herbivores). • Propriété remarquable: soit α = d / a. H = α c N / d + b P / a + log (N^α P) ne dépend pas du temps.

13. II.1. Populations en interaction: Lotka Volterra

14. II.1. Populations en interaction: Lotka Volterra

15. II.2. Populations en interaction: Lotka Volterra modifié • Deux populations: une de proies, une de prédateurs • N(t) nombre de proies, P(t) nombre de prédateurs • Le modèle de Lotka Volterra: • Modèle logistique + prédation pour les proies • Modèle logistique pour les prédateurs, la capacité du milieu étant proportionnelle au nombre de proies. • Mise en équations ∂t N = r N (1 – N/K) – B N P / (A + N) ∂t P = k (1 – h P/N) • Remarque: plus de quantité conservée comme H. • Notion de cycle limite: les trajectoires « s’enroulent » autour d’une solution périodique.

16. II.2. Populations en interaction: Lotka Volterra modifié

17. II.3. Populations en interaction: compétition • Deux populations en compétition pour la même ressource. • N(t) nombre d’individus de la première espèce, P(t) nombre d’individus de la seconde espèce. • Modélisation: • Les espèces se gênent • La capacité du milieu est partagée par les deux espèces • Mise en équations ∂t N = rn N (1 – N/ Kn – b P /Kn) ∂t P = rp P (1 – P/Kp – b’ N/Kp) • Signification des constantes: • Kn : nombre d’individus de la première espèce que peut nourrir le milieu • Kp : nombre d’individus de la seconde espèce que peut nourrir le milieu • bP: fraction des ressources du milieu utilisées par l’espèce 2 (y compris la gêne) • b’ N: fraction du milieu utilisé par l’espèce 1. • Etats d’équilibre: N + b P = Kn P + b’ N = Kp

18. II.3. Populations en interaction: compétition • Changement d’inconnues: u = N/Kn, v = P/Kp a = b Kp/Kn, a’ = b’ Kn/Kp • Discussion: principe d’exclusion compétitive:

19. II.4. Populations en interaction: mutualisme • Deux populations en symbiose se facilitent mutuellement l’accès à la ressource. • N(t) nombre d’individus de la première espèce, P(t) nombre d’individus de la seconde espèce. • Modélisation: • Les espèces en symbiose • La capacité du milieu est partagée par les deux espèces • Mise en équations ∂t N = rn N (1 – N/ Kn + b P /Kn) ∂t P = rp P (1 – P/Kp + b’ N/Kp) • Signification des constantes: • Kn : nombre d’individus de la première espèce que peut nourrir le milieu • Kp : nombre d’individus de la seconde espèce que peut nourrir le milieu • bP: fraction des ressources du milieu rendue utilisable par l’espèce 2 par symbiose. • b’ N: fraction du milieu rendue utilisable par l’espèce 1.

20. II.4. Populations en interaction: mutualisme

21. Exemple: épidémies

22. Epidémies: SIR • Maladie contagieuse. • Trois populations: • S(t): nombre d’individus sains • I(t): nombre d’individus malades • R(t): nombre d’individus morts, ou guéris et immunisés contre la maladie. • Modélisation: • Contamination proportionnelle au nombre de rencontres entre individus sains et malades. • Les malades ont une certaine probabilité de guérir par unité de temps. • Mise en équations ∂t S = - r S I ∂t I = r S I – a I ∂t R = a I • Remarque: S + I + R ne dépend pas du temps (conservation du nombre d’individus)

23. Epidémies: SIR

24. Epidémies: SIR

25. Petits systèmes d’EDOs

26. Une EDO • u un scalaire. • Comportements possibles • Explosion: réaction autocatalysée ∂t u = u^2 Solution en 1/ (T_0 – t) • Convergence vers un équilibre stable: ∂t u = f(u) u -> u* avec f(u*) = 0 convergence à vitesse exponentielle généralement • Solution constante, reste sur un équilibre instable (exceptionnel)

27. Une EDO: un exemple • Exemple: comportements possibles pour ∂t u = u (1 – u)

28. Une EDO: pas de comportement oscillant possible • Pas de solution périodique possible pour une seule EDO: • Preuve: ∂t u = f(u) On multiplie par ∂t u ce qui donne (∂t u)^2 = f(u) ∂t u Que l’on intègre entre t et t + T ∫ (∂t u)^2 = ∫ f(u) ∂t u Soit F définie par F’ = f Alors la dérivée de F(u(t)) vaut f(u) ∂t u donc la seconde intégrale vaut F(u(t+T)) – F(u(t)) = 0 Donc la première intégrale est nulle donc u est constante !

29. Deux EDO • u et v deux scalaires. ∂t u = f(u,v) ∂t v = g(u,v) • Comportements possibles • Explosion • Convergence vers un équilibre stable: u -> u* et v -> v* avec f(u*,v*) = g(u*,v*) = 0 convergence à vitesse exponentielle généralement • Solution constante, reste sur un équilibre instable (exceptionnel) • Convergence vers une solution périodique

30. Deux EDO • Exemple: solution périodique stable ∂t (u + i v) = i*(u+i v) + (u+i v)*(1 – u^2 – v^2) cycle limite stable u^2 + v^2 = 1

31. Deux EDO: cas linéaire • u et v deux scalaires. ∂t u = a u + b v ∂t v = c u + d v • Solution explicite: • Matrice M de coefficients a b c d • Recherche de vecteurs propres et valeurs propres M e_1 = λ_1 e_1 M e_2 = λ_2 e_2 (sauf cas particulier λ_1 = λ_2). • Solution est de la forme (u(t),v(t)) = a_1 e_1 exp(λ_1 t) + a_2 e_2 exp(λ_2 t). • Comportement asymptotique dépend des signes des parties réelles de λ_1 et λ_2 • 0 est stable si et seulement si Re(λ_1) < 0 et Re(λ_2) < 0.

32. Deux EDO: cas linéaire: classification

33. Deux EDO: cas linéaire: classification

34. Trois ODE • Trois scalaires u, v et w • Les comportements précédents ne sont pas les seuls possibles • Chaos possible: exemple le plus simple: le système de Lorenz: ∂t x = s (y-x) ∂t y = r x – y – xz ∂t z = x y – bz avec s = 10, r = 28, b = 8/3

35. Trois ODE: Lorenz

36. Trois ODE: Lorenz

37. Plusieurs ODE • Classification impossible … • Notion d’attracteur étrange • Grande variabilité en fonction des paramètres • Etude numérique est la seule possible en général • Sauf cas très rares, pas de solution explicite aux équations différentielles ordinaires !

38. Simulations numériques

39. Schéma d’Euler (explicite) • Le schéma le plus simple pour résoudre ∂t u = f(u). • Principe: • Calcul approché de u(t) pour t = 0, k, 2k, 3k, …. où k est le pas de temps. • On note u_n la valeur approchée de u au temps n k. • Erreur, d’autant plus petite que k est petit • Pour évaluer u au temps T il faut calculer u_(T/k) donc faire T/k calculs • Plus k est petit plus le calcul est précis et plus il est long (logique !) • Le schéma d’Euler explicite • Approche ∂t u au temps n k par ( u_(n+1) – u_n ) / k • Schéma u_(n+1) = u_n + k f(u_n) • u_0 donnée initiale • Calcul itératif de u_(n+1) en fonction de u_n

40. Schéma d’Euler (explicite) • Implémentation informatique de • ∂t u = 2 u (1 – u) • u(0) = 0.1 • Programme Matlab ou R: u_0 = 0.1; % donnée initiale k = 0.01; % pas de temps Tmax = 5; % temps maximal de calcul u = zeros(Tmax/k,1); % vecteur qui va contenir la solution u(1) = u_0; for J=1:Tmax/k-1, % boucle de calcul u(J+1) = u(J) + k*2*u(J)*(1-u(J)); end plot(u); % affichage de la solution • Démonstration: programme eulerexplicite.m

41. Schéma d’Euler (explicite) • Précision: proportionnelle à k

42. Schéma d’Euler (explicite) • Limitations: • Erreur: proportionnelle à k … peut mieux faire -> Runge Kutta • Echoue sur les problèmes « raides ». Exemple ∂t u = N f(u) où N est très grand: réaction très rapide. • Avantages: • Très simple à mettre en œuvre • En particulier lorsque f est très complexe • Autres schémas: • Runge Kutta: ordre plus élevé: erreur en k^4 • Euler implicite: supporte mieux les problèmes raides.

43. Schéma d’Euler (implicite) • Le schéma d’Euler implicite • Approche ∂t u au temps n k par ( u_(n+1) – u_n ) / k • Schéma u_(n+1) = u_n + k f( u_(n+1) ) • u_0 donnée initiale • Implicite: il faut résoudre une équation pour obtenir u_(n+1) • Equation u_(n+1) - k f( u_(n+1) ) = u_n • Résolution de cette équation par une méthode de Newton

44. Runge Kutta • Schéma plus complexe u_(n+1) = u_n + k (k_1 + 2k_2 + 2k_3 + k_4) /6 avec k_1 = f (t_n, x_n) k_2 = f (t_n + k⁄2, x_n + k⁄2 k_1) k_3 = f (t_n + k⁄2, x_n + k⁄2 k_2) k_4 = f (t_n + k, x_n + k k_3) • Ordre 4: erreur en k^4 • Erreur beaucoup plus petite … mais schéma plus complexe ! • Très souvent implémenté dès que l’on veut plus de précision que pour Euler.

45. Autres méthodes • Schémas d’ordres plus élevés: précision en k^N avec N aussi grand que l’on veut … mais schéma plus complexes. • Schémas implicites • Schémas avec contrôle a posteriori d’erreur. • Bilan: • Débuter par Euler • Si nécessaire passer à Runge Kutta 4 • En cas d’échec … consulter un spécialiste !

46. D’autres modèles

47. Equations avec retard • Dynamique des populations: l’évolution de la population N au temps t dépend de N(t – T) où T est le temps de gestation pour les naissances, et de N(t) pour la mortalité • Epidémie: idem: T temps d’incubation • La dérivée de N(t) est une fonction de N(t – T) et de N(t) ∂t N(t) = f( N(t), N(t-T) ) • Exemple: une variante de l’équation logistique ∂t N(t) = r N(t) ( 1 – N(t-T) / K) avec K capacité du milieu • Solutions oscillantes possibles: exemple ∂t N(t) = π N(t-T) / 2T a pour solution N(t) = A cos (π t / 2 T) périodique de période 4T.

48. Equations avec retard: exemple: mouches et moutons • Moutons australiens et mouches … • Oscillations avec une période 35 à 40 jours.

49. Equations avec retard: respiration de Cheyne Stokes • Physiologie: • Respiration anormale • Périodes d’apnée • L’amplitude de la respiration augmente et diminue régulièrement, avec des périodes d’apnée.

50. Equations avec retard: respiration de Cheyne Stoke • C(t) niveau de CO2 dans les artères • La ventilation V(t) dépend de C(t) avec un retard T V = Vmax c(t-T) / [ a + c(t-T) ] • Vmax : ventilation maximale, T temps de retard • Evolution de c: ∂t C(t) = p – b V C(t) ce qui donne ∂t C(t) = p – b Vmax C(t) C(t-T) / [ a + c(t-T) ]