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
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
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:
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.
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
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\}.
$$
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}.
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
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
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,
| $\,$ |
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,
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
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,
| $\,$ | $\,$ |
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
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}.$$
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:
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.
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$).