$$ \newcommand{\N}{\mathbb{N}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\Q}{\mathbb{Q}} \newcommand{\R}{\mathbb{R}} \newcommand{\C}{\mathbb{C}} \newcommand{\K}{\mathbb{K}} \newcommand{\Ns}{\mathbb{N}^*} % \newcommand{\ABA}{\mathscr{A}} \newcommand{\FF}{\mathscr{F}} \newcommand{\BB}{\mathscr{B}} \newcommand{\CC}{\mathscr{C}} \newcommand{\MM}{\mathscr{M}} \newcommand{\PP}{\mathscr{P}} \newcommand{\LL}{\mathscr{L}} \newcommand{\VV}{\mathscr{V}} \newcommand{\DD}{\mathscr{D}} \newcommand{\SB}{\mathscr{S}} \newcommand{\EE}{\mathscr{E}} \newcommand{\HH}{\mathscr{H}} \newcommand{\dx}{\mbox{d} x} \newcommand{\dt}{\mbox{d} t} %\newcommand{\dt}{\mathop{\mathrm{d}}\nolimits t} \newcommand{\du}{\mathop{\mathrm{d}}\nolimits u} \newcommand{\dg}{\mathop{\mathrm{deg}}\nolimits} \newcommand{\re}{\mathop{\mathrm{Re}}\nolimits} \newcommand{\im}{\mathop{\mathrm{Im}}\nolimits} \newcommand{\ii}{\mathop{\mathrm{\textbf{i}}}\nolimits} \newcommand{\ee}{\mathop{\mathrm{e}}\nolimits} \newcommand{\jj}{\mathop{\mathrm{\textbf{j}}}\nolimits} \renewcommand {\epsilon}{\varepsilon} \renewcommand {\le}{\leqslant} \renewcommand {\ge}{\geqslant} \renewcommand {\leq}{\leqslant} \renewcommand {\geq}{\geqslant} \newcommand{\transp}[1]{\sideset{^{t}}{}{\mathop{#1}}} \newcommand{\tendvers}[2]{\xrightarrow[#1 \rightarrow #2]{}} \newcommand{\tendversN}{\xrightarrow[n \rightarrow \infty]{}} \newcommand{\limiteX}[2]{\underset{#1\rightarrow #2}{\lim \,}} \newcommand{\limite}[1]{\underset{#1}{\lim \,}} \def\norme#1{\left\Vert#1\right\Vert} \def\normeinf#1{\left\Vert#1\right\Vert_\infty} \newcommand{\abs}[1]{\left| #1 \right|} \newcommand{\dsum}{\displaystyle \sum} \newcommand{\Norme}[1]{| \mkern-2mu | \mkern-2mu | {#1} | \mkern-2mu | \mkern-2mu |} \newcommand{\ud}{\,\mathrm{d}} \newcommand{\Rp}{\mathbb{R}_+} \newcommand{\enumsp}[3]{#1_{#2},\ldots,#1_{#3}} \newcommand{\inter}[1]{[\![#1 ]\!]} \newcommand{\fonct}[5]{#1: \left\{ \begin{array}{rcl} #2 & \longrightarrow & #3 \\ #4 & \longmapsto & #5 \end{array} \right.} \newcommand{\paren}[1]{\left( #1 \right)} $$

Résolution approchée des équations différentielles

Rappel

Les équations différentielles rencontrées dans les divers domaines des sciences (physique, mécanique...) ont très rarement des solutions explicites et ne s'intègrent donc pas de manière exacte. C'est pour cette raison qu'on fait appel à des méthodes numériques afin de trouver des solutions numériques approchées.
Soit $I=[t_0,t_0+T]$, $T>0$, un intervalle de $\R$, et soit $f$ une fonction continue de $I\times \R^m$ dans $\R^m$.

On appelle équation différentielle du premier ordre une équation de la forme: \begin{equation} \label{Cauchy1} y'(t)=f(t,y(t)),\;\;\forall t\in I, \end{equation}

Une condition initiale à ce problème est la donnée de la solution au point $t_0$: \begin{equation} \label{Cauchy2} y(0)=y_0,\;\;y_0\in \R^m. \end{equation} Le problème de Cauchy consite à déterminer une fonction $y$ de $I$ dans $\R^m$ continue et dérivable sur $I$; telle que $y$ vérifie \ref{Cauchy1} et \ref{Cauchy2}.
Selon la régularité de la fonction $f$, il y a existence et unicité d'une solution, mais il peut ne pas y avoir de solution. Dans ce TP on se limitera au cas o\`u la fonction $f$ est Lipschitzienne par rapport à sa deuxième variable. Dans ce cas, on a le résultat suivant:

Cauchy-Lipschitz

Supposons qu'il existe une constante $M\in \R_+^*$, telle que: \begin{equation} \label{Lip} \forall\, t\in I;\;\forall\,(x,y)\in \R^m,\;\;\|f(t,x)-f(t,y)\|\leq M\|x-y\|, \end{equation} alors le problème (\ref{Cauchy1},\,\ref{Cauchy2}) admet une et une seule solution $y$ définie sur I.

Schémas numériques, méthodes à un pas

Dans cette partie on étudie les méthodes à un pas.
Considérons le problème de Cauchy (\ref{Cauchy1},\ref{Cauchy2}). Dans la pratique, on ne peut pas résoudre explicitement l'équation différentielle, tout en sachant qu'elle admet une solution.
De ce fait la solution est construite numériquement, c'est-à-dire on calcule des valeurs approchées de la solution exacte en certains points $t\in I$. A partir de ces valeurs, on construit une solution approchée.
La méthode à un pas consiste à construire la suite $(y_n)_{0\leq n\leq N}$, $N\geq 1$, des valeurs approchées de la solution $y$. Pour cela on se donne une subdivision $(t_n)_n$ de $I$ telle que: $$ \forall i\in\{0,...,N\},\;\; t_i=t_0+h\times i, $$ $h=\frac{T}{N}$ est, par définition, le pas du maillage sur $I$. Dans cette méthode, le calcul de $y_{n+1}\simeq y(t_{n+1})$ se fait à partir de $y_n\simeq y(t_n)$ à l'aide de la relation de récurrence: \begin{equation} \label{unpas} \left\{ \begin{array}{l} y_0 \mbox{ est donné }=y(t_0)\\ y_{n+1}=y_{n}+h\Phi(t_n,y_n,h) \end{array} \right. \end{equation} où $\Phi$ est une fonction continue de $I\times \R^m \times [0,h^*]$ dans $\R^m$, et $h^*\geq h$.
Pour savoir la précision de cette méthode, on calcule une estimation de $e_n=y(t_n)-y_n$. On définit l'erreur globale de la solution numérique $E(y_0,h)$ par: $$ E(y_0,h)= \max\{\|e_n\|, 0\leq n\leq N\}. $$

La méthode à un pas est dite convergente si: $$ \forall y_0\in \R^m,\;\lim_{h\rightarrow 0}E(y_0,h)=0. $$

On définit aussi l'erreur locale de troncature par: $$ \varepsilon_n=y(t_{n+1})-y(t_n)-h\Phi(t_n,y_n,h), $$ $\varepsilon_n$ mesure avec quelle précision la solution exacte vérifie la schéma \ref{unpas}.

La méthode à un pas est dite consistante si pour toute solution $y$ du problème \ref{Cauchy1} on a: $$ \lim_{h\rightarrow 0}\sum_{n=0}^N\|\varepsilon_n\|=0. $$

On dit que la méthode à un pas est d'ordre $p>0$, si $p$ est le plus grand entier tel que: $$ \forall t\in I,\; \varepsilon(t)=y(t+h)-y(t)-h\Phi(t,y(t),h)=\mathcal{O}(h^{p+1})\;\;(h\rightarrow 0), $$ ceci pour tout $f\in \CC^p(I)$ et toute $y$ solution du problème du Cauchy.

Enfin, on dit que la méthode à un pas est stable s'il existe une constante $K\in \R_+^*$ telle que: pour tout suites $(y_n)_{0\leq n\leq N}$, $(z_n)_{0\leq n\leq N}$, et $(\varepsilon_n)_{0\leq n\leq N}$ vérifiant: $$ y_{n+1}=y_n+h\Phi(t_n,y_n,h),\quad z_{n+1}=z_n+h\Phi(t_n,z_n,h)+\varepsilon_n, $$ on a: $$ \max_{n}\|y_n-z_n\|\leq K\left(\|y_0-z_0\|+\sum_{j=0}^{N-1}\|\varepsilon_j\|\right). $$ Autrement dit, une petite perturbation sur les données du schéma n'entraîne qu'une petite perturbation sur la solution.

(Convergence) Si la méthode à un pas est stable et consistante alors elle est convergente.

Méthode d'Euler

Le schéma d'Euler est la forme la plus simple des méthodes numériques à un pas. Le schéma consiste à poser: $$ \forall (t,y,h)\in I\times \R^m\times [0,h^*]\;,\Phi(t,y,h)=f(t,y). $$ Ce qui revient à poser: $$ y_{n+1}=y_n+hf(t_n,y_n), $$ i.e. à la ($n+1$)\up{ième} étape on confond la courbe de $y$ au voisinage de $t_n$ avec ce que serait sa tangente si elle passait effectivement par le point $(t_n,y_n)$.
Le schéma d'Euler est consistant et stable (car $f$ est Lipschitzienne en $y$). Il est donc convergent et cette méthode est de l'ordre $1$.
Ce schéma est très simple à programmer, mais il est guère utilisé dans la pratique.

				def Euler(F, a, b, y_0, n):
h = (b - a)/n    # pas du maillage
X = [a]
x = a
Y = [y_0]
y = y_0
for k in range(0, n):
	y = y + h * F(x, y)
	x = x + h
	X.append(x)
	Y.append(y)
return X, Y
				
				
Exemple

On considère le problème suivant: $$\left\{\begin{array}{lcl} y'&=&y+1\\ y(t_0)&=&y_0. \end{array} \right. $$ Dans ce cas, on connaît la solution exacte, à savoir $y(t)=-1+(y_0+1)\ee^{t-t_0}$.
L'utilisation de la méthode d'Euler donne les résultat suivant,

$\,$


Exemple $(\star)$

On considère le problème suivant: $$\left\{\begin{array}{lcl} y'&=&\dfrac{3y}{t}-\dfrac{5}{t^3},\,\quad t>0\\ y(t_0)&=&y_0. \quad (t_0>0) \end{array} \right. $$ Dans ce cas, on connaît la solution exacte, à savoir $$\forall t>0,\quad y(t)=\dfrac{y_0t_0^2-1}{t_0^5}t^3-\dfrac{1}{t^2}.$$ L'utilisation de la méthode d'Euler donne les résultat suivant,


Méthode Euler-Cauchy

Cette méthode est une version améliorée de la méthode d'Euler.
Principe de la méthode:
Soit $y$ la solution du problème de Cauchy. Alors $y$ vérifie: $$ \forall (t,t+h)\in I^2,\;y(t+h)=y(t)+\int_t^{t+h}f(s,y(s))\ud s. $$ L'idée ici consiste à approcher l'intégrale ci-dessus par la méthode des trapèzes i.e. $$ \int_t^{t+h}f(s,y(s))\ud s\simeq \frac{h}{2}(f(t,y(t))+f(t+h,y(t+h))) $$ Or, on ne connaît pas la valeur de $y(t+h)$, donc on la remplace dans la relation précédente par son approximation au sens d'Euler (i.e., $y(t)+hf(t,y(t))$), soit: $$ y(t+h)\simeq y(t)+\frac{h}{2}(f(t,y(t))+f(t+h,y(t)+hf(t,y(t)))). $$ Alors $\Phi$ s'écrit de la manière suivante: $$ \forall (t,y,h)\in I\times \R^m\times [0,h^*]\;,\Phi(t,y,h)=(f(t,y)+f(t,y+hf(t,y)))/2. $$ Cette méthode donne de meilleurs résultats que celle d'Euler; de plus, elle est d'ordre 2.

Exemple $(\star)$

On considère le problème suivant: $$\left\{\begin{array}{lcl} y'&=&\dfrac{3y}{t}-\dfrac{5}{t^3},\,\quad t>0\\ y(t_0)&=&y_0. \quad (t_0>0) \end{array} \right. $$ Dans ce cas, on connaît la solution exacte, à savoir $$\forall t>0,\quad y(t)=\dfrac{y_0t_0^2-1}{t_0^5}t^3-\dfrac{1}{t^2}.$$ L'utilisation de la méthode d'Euler donne les résultat suivant,

$\,$
$\,$

Méthode de Rung-Kutta

Soit $q\in \N$. On considère $q$ nombres réels $\tau_1,...,\tau_q$ dans l'intervalle $[0,1]$ auxquels est associée la formule de quadrature numérique: \begin{equation} \label{qua_1} 1\leq i\leq q,\;\; \int_0^{\tau_i}\phi(\tau)d\tau\simeq \sum_{1\leq j\leq q}a_{ij}\phi(\tau_i) \end{equation} \begin{equation} \label{qua_2} \int_0^1\phi(\tau)d\tau\simeq \sum_{1\leq j\leq q}b_j\phi(\tau_j)\,. \end{equation} L'idée de cette méthode est d'approcher l'intégrale de $f$ entre $t_n$ et $t_n+h$ en utilisant ces formules.
Dans ce but, introduisons les instants intermédiaires: $$ t_{n,i}=t_n+\tau_i h,\;\;1\leq i\leq q. $$ Si $y$ est une solution du problème de Cauchy, on écrit en utilisant les relations \ref{qua_1} et \ref{qua_2}: $$ y(t_{n,i})=y(t_n)+\int_{t_n}^{t_{n,i}}f(s,y(s))ds\simeq y(t_n)+h\sum_{j=1}^{j=q}a_{ij}f(t_{n,j},y(t_{n,j})) $$ $$ y(t_{n+1})=y(t_n)+\int_{t_n}^{t_{n+1}}f(s,y(s))ds\simeq y(t_n)+h\sum_{j=1}^{j=q}b_jf(t_{n,j},y(t_{n,j}))\,. $$ Notons $y_{n,i}$ une approximation de $y(t_{n,i})$ et $K_{n,i}=f(t_{n,i},y_{n,i})$. La méthode de Rung-Kutta s'écrit alors sous la forme: $$ \left\{ \begin{array}{l} K_{n,i}=f\left(t_{n,i},h\sum_{1\leq j\leq q}a_{ij}K_{n,j}\right)\;\;(i\in\inter{1,q})\\ \\ y_{n+1}=y_n+h\sum_{j=1}^qb_jK_{n,j} \end{array} \right. $$ Dans ce TP, on choisit comme exemple la méthode de Runge-Kutta la plus classique (dite RK4), caractérisée par: $$(a_{ij})=\left( \begin{array}{cccc} 0&0&0&0\\ 1/2&0&0&0\\ 0&1/2&0&0\\ 0&0&1&0 \end{array} \right),\,(b_j)=\left( \begin{array}{cccc} 1/6\\ 1/3\\ 1/3\\ 1/6 \end{array} \right). $$ Et $\tau_1=0,\tau_2=1/2,\tau_3=1/2,\tau_4=1$.
Le schéma s'écrit alors: $$ \left\{ \begin{array}{l} K_{n,1}=f(t_n,y_n)\\ K_{n,2}=f(t_n+\frac{h}{2},y_n+\frac{h}{2}K_{n,1})\\ K_{n,3}=f(t_n+\frac{h}{2},y_n+\frac{h}{2}K_{n,2})\\ K_{n,4}=f(t_n+h,y_n+hK_{n,3})\\ \\ y_{n+1}=y_n+\frac{h}{6}\left(K_{n,1}+2K_{n,2}+2K_{n,3}+K_{n,4}\right) \end{array} \right. $$ On montre que cette méthode est d'ordre 4.

Exemple $(\star)$

On considère le problème suivant: $$\left\{\begin{array}{lcl} y'&=&\dfrac{3y}{t}-\dfrac{5}{t^3},\,\quad t>0\\ y(t_0)&=&y_0. \quad (t_0>0) \end{array} \right. $$ Dans ce cas, on connaît la solution exacte, à savoir $$\forall t>0,\quad y(t)=\dfrac{y_0t_0^2-1}{t_0^5}t^3-\dfrac{1}{t^2}.$$


Exemples concours

Mines 2016

Texte d'après le sujet de Mines 2015
L'étude de la propagation des épidémies joue un rôle important dans les politiques de santé publique. Les modèles mathématiques ont permis de comprendre pourquoi il a été possible d'éradiquer la variole à la fin des années 1970 et pourquoi il est plus difficile d'éradiquer l'apparition d'épidémies de grippe tous les hivers. Aujourd'hui, des modèles de plus en plus complexes et puissants sont développés pour prédire la propagation d'épidémies à l'échelle planétaire telles que le SRAS, le virus H5N1 ou le virus Ebola.
Ces prédictions sont utilisées par les organisations internationales pour établir des stratégies de prévention et d'intervention.

On s'intéresse ici à une première méthode de simulation numérique.
Les modèles compartimentaux sont des modèles déterministes où la population est divisée en plusieurs catégories selon leurs caractéristiques et leur état par rapport à la maladie. On considère dans cette partie un modèle à quatre compartiments disjoints:

  1. sains (S, "susceptibles"),
  2. infectés (I,"infected") ,
  3. rétablis (R,"recovered", ils sont immunisés)
  4. et décédés (D,"dead").
Le changement d'état des individus est gouverné par un système d'équations différentielles obtenues en supposant que le nombre d'individus nouvellement infectés (c'est-à-dire le nombre de ceux qui quittent le compartiment S) pendant un intervalle de temps donné est proportionnel au produit du nombre d'individus infectés avec le nombre d'individus sains.
En notant $S(t)$, $I(t)$, $R(t)$ et $D(t)$ la fraction de la population appartenant à chacune des quatre catégories à l'instant $t$, on obtient le système: \begin{equation}\label{eq1} \left. \begin{array}{rcl} \dfrac{\ud}{\ud t}S(t) & = & -rS(t)I(t) \\ \dfrac{\ud}{\ud t}I(t) & = & rS(t)I(t)-(a+b)I(t)\\ \dfrac{\ud}{\ud t}R(t) & = & a I(t) \\ \dfrac{\ud}{\ud t}D(t) & = & bI(t) \end{array} \right\} \end{equation} avec $r$ le taux de contagion, $a$ le taux de guérison et $b$ le taux de mortalité. On suppose qu'à l'instant initial $t=0$, on a $S(0)=0,95$, $I(0)=0,05$ et $R(0)=D(0)=0$.
On définit $$X(t)=\begin{pmatrix} S(t)\\I(t)\\R(t)\\D(t) \end{pmatrix},\,\fonct{f}{\R\times\R^4}{\R^4}{(t,\transp (x\,y\,z\,w))}{\begin{pmatrix} -rxy\\rxy-(a+b)y\\ay\\by\end{pmatrix}}$$ Ainsi, avec ces définitions, on obtient $$\forall t\geq 0,\, X'(t)=f(t,X(t)),\text{ et } X(0)=\begin{pmatrix}S(0)\\I(0)\\R(0)\\D(0)\end{pmatrix}.$$

Ci-après la solution approchée de ce système obtenue en utilisant la méthode de Rung-Kutta. Les paramètres intiales sont les mêmes que dans le sujet de concours, mais vous pouvez les modifier.



Centrale 2015

Modéliser les interactions physiques entre un grand nombre de constituants mène à l’écriture de systèmes différentiels pour lesquels, en dehors de quelques situations particulières, il n’existe aucune solution analytique. Les problèmes de dynamique gravitationnelle et de dynamique moléculaire en sont deux exemples. Afin d’analyser le comportement temporel de tels systèmes, l’informatique peut apporter une aide substantielle en permettant leur simulation numérique. L’objet de ce sujet, composé de quatre parties, est l’étude de solutions algorithmiques en vue de simuler une dynamique gravitationnelle afin, par exemple, de prédire une éclipse ou le passage d’une comète.

Soient 𝑦 une fonction de classe $\CC^2$ sur $\R$ et $𝑡_{min}$ et $𝑡_{max}$ deux réels tels que $𝑡_{min}< 𝑡_{min}$. On note $I$ l’intervalle $[𝑡_{min}, 𝑡_{max}]$. On s’intéresse à une équation différentielle du second ordre de la forme : \begin{equation} \label{EqC15} \forall t\in I,\quad y''(t) = f(y(t)). \end{equation} où $f$ est une fonction donnée, continue sur $\R$. De nombreux systèmes physiques peuvent être décrits par une équation de ce type. On suppose connues les valeurs $y_0 = y(t_{min})$ et $z_0 = y'(t_{min})$. On suppose également que le système physique étudié est conservatif. Ce qui entraine l’existence d’une quantité indépendante du temps (énergie, quantité de mouvement, …), notée $E$, qui vérifie l’équation (\ref{EqC15}) où $g'=-f$ . \begin{equation} \label{EqC15Bis} \forall t\in I,\quad \dfrac{1}{2}y'(t)^2+g(y(t))=E. \end{equation} Pour résoudre numériquement l’équation différentielle (\ref{EqC15}), on introduit la fonction $z:I\longmapsto\R$ définie par, $\forall t\in \R,\, z(t)=y'(t)$.
On pose alors $Y=\begin{pmatrix} y\\y'\end{pmatrix}$ et $\fonct{F}{\R^2}{\R}{\begin{pmatrix} x\\y\end{pmatrix}}{\begin{pmatrix}y\\f(x)\end{pmatrix}}$, alors $Y$ vérifie l'équation différentielle suivante: $$\forall t\geq 0,\,Y'(t)=F(Y(t)),\, \quad F(0)=\begin{pmatrix} y_0\\z_0 \end{pmatrix} $$ Soit $n\geq 2$, on pose $$h=\dfrac{t_{max}-t_{min}}{n},\,\,\forall k\in \inter{0,n},\, t_k=t_{min}+kh.$$ Avec ces notation, le schéma d'Euler est donné par les relations suivantes, $$ Y_0=\begin{pmatrix} y_0\\z_0\end{pmatrix},\,\forall k\in\inter{1,n-1},\, Y_{k+1}=Y_k+hF(Y_k)=\begin{pmatrix} y_i+hz_i\\ z_i+hf(y_i) \end{pmatrix} $$ Le physicien français Loup Verlet a proposé en 1967 un schéma numérique d’intégration d’une équation de la forme (\ref{EqC15}) dans lequel, en notant $f_i= f(y_i)$ et $f_{i+1} =f(y_{i+1})$, les relations de récurrence s’écrivent $$ y_{i+1}=y_i+hz_i+\dfrac{h^2}{2}f_i,\text{ et }z_{i+1}=z_i+\dfrac{h}{2}(f_i+f_{i+1}).$$ La simulation ci-après donne les résultats des schémas Euler et Verlet avec l'exemple donné dans le sujet ($f(y)=-\omega^2 y$).



Exercices