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 \[ \frac{\mathrm{d}y}{\mathrm{d}t} = F(y,t), \quad y(t_0)=y_0, \] où \(F\) est une fonction qui dépend de \(y\) et de \(t\), et où \(y_0\in\mathbb{R}\) est la condition initiale de l’équation différentielle.
Par exemple, dans le cas de l’équation différentielle \[ y'(t)=y(t)+8t^2+1, \] \(F(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(t_0)=y_0\). Ici on pourrait par exemple avoir \[ y(t=0)=0. \]
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 \[ y_{n+1}=y_0+\delta t\sum_{i=1}^s b_ik_i, \] où \(\delta t\) est le pas de temps, et \[ y_n\equiv y(t_n), \] avec \(t_n=t_0+n\delta t\). La valeur des \(k_i\) dépend de la précision de l’algorithme (et donc de la valeur du nombre “d’étages” \(s\)) et ils sont donnés de façon itérative par \[\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 \(a_{ij}\), \(b_i\) et \(c_i\) sont propres de nouveaux au nombre d’étages \(s\) de la méthode. Pour \(s=1\), on a simplement \(c_1=0\), \(b_1=1\) et donc il vient \[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=2\), on obtient \(c_1=0\), \(c_2=1/2\), \(a_{21}=1/2\), \(b_1=0\), et \(b_2=1\). On a donc \[ 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 \[ 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 est un système d’équations différentielles ordinaires agissant sur une variable vectorielle à trois composantes \(\vec y(t)=(y_x(t),y_y(t),y_z(t))\) \[\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 \(\beta=8/3\), \(\sigma=10\), et \(\rho=28\).
La solution de l’équation de Lorenz est dite chaotique. Cela signifie que si nous considérons deux conditions initiales \(\vec y^1\) et \(\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 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.
É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.
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).
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 \(t_{max}=1.5586522\): \(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 \(\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.
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é.
Afin de modéliser la propagation d’une épidémie, nous allons utiliser le modèle SEIR, qui est un modèle compartimental où \(S\) est la population susceptible d’être infectée, \(E\) la population exposée (pas encore infectieuse) qui correspond à la population où la population est en incubation, \(I\) la population infectieuse, et finalement \(R\) 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)=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ù \(\mathcal{R}_0\) est taux de reproduction de base, \(T_{inf}\) est le temps où un individu est infectieux, et \(T_{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.
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!
Me cassez pas les pieds avec vos évaluations!