Skip to content
Snippets Groups Projects
# author:
# - Orestis Malaspinas
title: Travail pratique sur les équations différentielles
autoSectionLabels: true
autoEqnLabels: false
eqnPrefix: 
  - "éq."
  - "éqs."
chapters: false
numberSections: false
chaptersDepth: 1
sectionsDepth: 3
lang: fr
documentclass: article
papersize: A4
cref: false
pandoc-numbering:
  - category: exercice
urlcolor: blue

\newcommand{\real}{\mathbb{R}}

Rappel théorique : approximation numérique d'équations différentielles ordinaires

En général il n'existe pas de solutions analytiques pour les équations différentielles d'un intérêt pratique en ingénierie, c'est pourquoi on est obligé d'utiliser des méthodes numériques pour approximer la solution. Nous présentons ici une famille de méthodes numériques : les méthodes de Runge-Kutta qui permettent d'approximer les solutions d'équations différentielles ordinaires

dydt=F(y,t),y(t0)=y0, \frac{\mathrm{d}y}{\mathrm{d}t} = F(y,t), \quad y(t_0)=y_0,
FF
est une fonction qui dépend de
yy
et de
tt
, et où
y0y_0\in\real
est la condition initiale de l'équation différentielle.

Par exemple, dans le cas de l'équation différentielle

y(t)=y(t)+8t2+1, y'(t)=y(t)+8t^2+1,
F(y,t)=y(t)+8t2+1F(y,t)=y(t)+8t^2+1
.

Afin d'avoir une solution unique, il est nécessaire de donner une condition initiale à une équation différentielle, de la forme

y(t0)=y0y(t_0)=y_0
. Ici on pourrait par exemple avoir
y(t=0)=0. y(t=0)=0.

Les méthodes de Runge-Kutta

Les méthodes de Runge--Kutta sont une famille de méthodes explicites pour résoudre des équations différentielles ordinaires. Elles ont la forme général suivante

yn+1=y0+δti=1sbiki, y_{n+1}=y_0+\delta t\sum_{i=1}^s b_ik_i,
δt\delta t
est le pas de temps, et
yny(tn), y_n\equiv y(t_n),
avec
tn=t0+nδtt_n=t_0+n\delta t
. La valeur des
kik_i
dépend de la précision de l'algorithme (et donc de la valeur du nombre "d'étages"
ss
) et ils sont donnés de façon itérative par
k1=F(yn,tn),k2=F(yn+a21δtk1,tn+c2δt),ks=F(yn+as1δtk1+as2δtk2++as,s1δtks1,tn+csδt).\begin{aligned} k_1&=F(y_n,t_n),\\ k_2&=F(y_n+a_{21}\delta t k_1,t_n+c_2\delta t),\\ &\cdots\\ k_s&=F(y_n+a_{s1}\delta t k_1+a_{s2}\delta t k_2+\cdots +a_{s,s-1}\delta t k_{s-1},t_n+c_s \delta t).\end{aligned}
et où les
aija_{ij}
,
bib_i
et
cic_i
sont propres de nouveaux au nombre d'étages
ss
de la méthode. Pour
s=1s=1
, on a simplement
c1=0c_1=0
,
b1=1b_1=1
et donc il vient
yn+1=yn+δtF(yn,tn).y_{n+1}=y_n+\delta t F(y_n,t_n).
Cette méthode s'appelle également méthode d'Euler explicite (vous l'avez peut-être reconnue). Pour
s=2s=2
, on obtient
c1=0c_1=0
,
c2=1/2c_2=1/2
,
a21=1/2a_{21}=1/2
,
b1=0b_1=0
, et
b2=1b_2=1
. On a donc
yn+1=yn+δtF(yn+δt2F(yn,tn),tn+δt2). y_{n+1}=y_n+\delta t F\left(y_n+\frac{\delta t}{2}F(y_n,t_n),t_n+\frac{\delta t}{2}\right).
Cette méthode s'appelle également méthode de point du milieu (elle est proche de la méthode du point du milieu pour les intégrales).

Finalement, une méthode bien connue pour les méthode de Runge--Kutta est celle d'ordre 4

yn+1=yn+16(k1+2(k2+k3)+k4), y_{n+1}=y_n+\frac{1}{6}(k_1+2(k_2+k_3)+k_4),
où \begin{align} k_1&=\delta tF(y_n, t_n),\ k_2&=\delta tF\left(y_n+\frac{k_1}{2}, t_n+\frac{\delta t}{2}\right),\ k_3&=\delta tF(y_n+\frac{k_2}{2}, t_n+\frac{\delta t}{2}),\ k_4&=\delta tF(y_n+k_3, t_n+\delta t). \end{align}

L'équation de Lorenz

L'équation de Lorenz est un système d'équations différentielles ordinaires agissant sur une variable vectorielle à trois composantes

y(t)=(yx(t),yy(t),yz(t))\vec y(t)=(y_x(t),y_y(t),y_z(t))
dyxdt=σ(yyyx),dyydt=yx(ρyz)yy,dyzdt=yyyxβyz,\begin{aligned} \frac{\mathrm{d}y_x}{\mathrm{d}t}&=\sigma(y_y-y_x),\\ \frac{\mathrm{d}y_y}{\mathrm{d}t}&=y_x(\rho-y_z)-y_y,\\ \frac{\mathrm{d}y_z}{\mathrm{d}t}&=y_yy_x-\beta y_z, \end{aligned}
où les constantes sont données par
β=8/3\beta=8/3
,
σ=10\sigma=10
, et
ρ=28\rho=28
.

Solution dans le plan x-z de l'équation de Lorenz pour une solution initiale donnée (image tirée de Wikipedia).{#fig_lorenz width="50%"}

La solution de l'équation de Lorenz est dite chaotique. Cela signifie que si nous considérons deux conditions initiales

y1\vec y^1
et
y2\vec y^2
qui sont très proches. Elles s'éloigneront l'une de l'autre exponentiellement vite au cours du temps.

Par ailleurs, comme on peut le constater sur la figure 1{reference-type="ref" reference="fig_lorenz"} la solution a une structure de "papillon". Cette structure est retrouvée indépendemment de la condition initiale (on aura toujours cette forme de solution) bien que la position à un temps donné ne soit pas la même pour deux conditions initiales différentes.

Travail à effectuer

Écriture des solveurs

Écrire un solveur pour la méthode d'Euler explicite (méthode de Runge--Kutta à un étage), un pour la méthode de Runge--Kutta à deux étages et un pour la méthode de Runge--Kutta d'ordre 4.

Équation de Lorenz

Afin de valider les solveurs, appliquer les trois solveurs à l'équation de Lorenz et reproduire une figure approchant celle de l'énoncé (faire un graphique éventuellement tri-dimmensionnel de la trajectoire obtenue pour une condition initiale quelconque).

Orbites périodiques

L'attracteur de Lorenz contient des trajectoires périodiques sur lesquelles le système revient au point initial après un certain temps. Ces trajectoires sont dites "instables", comme toute autre trajectoire de ce système : il suffit qu'on dévie un tout petit peu de la trajectoire, et on finit par s'en éloigner complètement. Le point suivant se trouve (dans la limite de la précision numérique fournie) sur une trajectoire périodique de période

tmax=1.5586522t_{max}=1.5586522
:
y0=(0.9101673912,1.922121396,18.18952097)y_0=(-0.9101673912,-1.922121396,18.18952097)
. Prendre ce point comme valeur initiale, et faites évoluer le système avec un pas de temps donné (par exemple
δt=0.001\delta t = 0.001
). Combien de temps arrivez-vous à rester sur l'orbite périodique avec chacun des solveurs? Tracez les deux trajectoires sur une figure 3D pour les comparer.

Modélisation d'épidémies

Une fois les solveurs validés vous vous attaquerez à la tâche d'utiliser vos solveurs pour résoudre un vrai problème d'actualité.

Le modèle SEIR

Afin de modéliser la propagation d'une épidémie, nous allons utiliser le modèle SEIR, qui est un modèle compartimental où

SS
est la population susceptible d'être infectée,
EE
la population exposée (pas encore infectieuse) qui correspond à la population où la population est en incubation,
II
la population infectieuse, et finalement
RR
la population rétablie.

La dynamique de ce système est caractérisée par la population initiale de chaque catégorie:

S(t=0)=S0,E(t=0)=E0,I(t=0)=I0,R(t=0)=R0, S(t=0)=S_0,\quad E(t=0)=E_0,\quad I(t=0)=I_0,\quad R(t=0)=R_0,
et par quatre équation différentielles ordinaires \begin{align} S'(t)&=-\frac{\mathcal{R}0}{T{inf}}I(t)\frac{S(t)}{N},\ E'(t)&=\frac{\mathcal{R}0}{T{inf}}I(t)\frac{S(t)}{N}-\frac{1}{T_{inc}}E(t),\ I'(t)&=\frac{1}{T_{inc}}E(t)-\frac{1}{T_{inf}}I(t),\ R'(t)&=\frac{1}{T_{inf}}I(t), \end{align} où
R0\mathcal{R}_0
est taux de reproduction de base,
TinfT_{inf}
est le temps où un individu est infectieux, et
TincT_{inc}
est le temps d'incubation de la maladie. La taux de reproduction de base peut être remplacé par le taux de reproduction effectif, \mathcal{R}_t=\mathcal{R}_0\frac{S(t)}{N} qui représente la susceptibilité effective au temps t.

Ce modèle peut être amélioré pour inclure d'autres termes (par exemple des termes reliés aux trajets de gens entrant ou sortant dans la population, ainsi que ternir compte des morts naturelles, etc). Pour le moment nous nous contentons de cette version simplifiée qui pourra être rendue encore plus réaliste par la suite.

Travail à effectuer

L'utilisation que vous ferez de votre modèle ne dépendra que de vous. Essayez de concevoir des scénarios différents de réaction des pouvoirs publics et simulez les!

Évaluation

Me cassez pas les pieds avec vos évaluations!