Transcript Introduction aux équations différentielles ordinaires (EDO) E. Grenier Exemple: dynamique des populations I.1.

Introduction aux équations
différentielles ordinaires (EDO)
E. Grenier
Exemple: dynamique des
populations
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).
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.
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.
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.
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
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
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)
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)
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.
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.
II.1. Populations en interaction: Lotka Volterra
II.1. Populations en interaction: Lotka Volterra
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.
II.2. Populations en interaction: Lotka Volterra modifié
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
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:
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.
II.4. Populations en interaction: mutualisme
Exemple: épidémies
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)
Epidémies: SIR
Epidémies: SIR
Petits systèmes d’EDOs
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)
Une EDO: un exemple
•
Exemple: comportements possibles pour
∂t u = u (1 – u)
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 !
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
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
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.
Deux EDO: cas linéaire: classification
Deux EDO: cas linéaire: classification
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
Trois ODE: Lorenz
Trois ODE: Lorenz
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 !
Simulations numériques
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
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
Schéma d’Euler (explicite)
•
Précision: proportionnelle à k
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.
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
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.
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 !
D’autres modèles
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.
Equations avec retard: exemple: mouches et moutons
•
•
Moutons australiens et mouches …
Oscillations avec une période 35 à 40 jours.
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.
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) ]
Equations avec retard: respiration de Cheyne Stoke
•
Simulation numérique
Equations avec retard: respiration de Cheyne Stoke
•
•
•
Changement de variables:
x = c/a, V* = V/Vmax, α = abVmax/p
T* = pT/a, t* = pt/a
ce qui donne
∂t x(t) = 1 – α x(t) x(t-T) / [1 + x(t-T)]
Etat stationnaire:
x indépendant du temps, égal à x*
α x*^2 = 1 + x*
Linéarisation
u = x – x* supposé très petit
∂t u(t) = - α V(x*) u(t) – α x* V’(x*) u(t-T)
on cherche des solutions
u(t) = u_0 exp(λ t)
ce qui donne
λ = - α V(x*) – α x* V’(x*) exp(-λ T)
équation en λ.
Equations avec retard: cellules sanguines
•
Dynamique du nombre de globules rouges ou blancs.
•
C(t) densité de cellules
•
Evolution:
–
–
•
Mortalité
Création par la moelle épinière, avec retard
Equation
∂t C(t) = f(C(t-T)) – g C(t)
où g est une constante et
f(x) = λ a^m x /(a^m + x^m).
Equations avec retard: cellules sanguines
Equations avec retard: cellules sanguines