Série 9: Équations différentielles ordinaires

1. Exercices

Exercice 1

On considère le problème de Cauchy suivant:

\[\left\{\begin{array}{l} y^{\prime}(t)=-\arctan (y) \quad \text { pour } t>0 \\ y(0)=1 \end{array}\right.\]

1: Ecrire les schémas d’Euler progressif et rétrograde avec pas de temps $h$ pour l’approximation numérique de $y(t)$.
2: On peut montrer que la solution $y(t)$ reste positive pour $t>0$, qu’elle est décroissante, et que $y(t) \rightarrow 0$ lorsque $t \rightarrow \infty$. Donner une condition sur $h$ afin que la méthode d’Euler progressive soit absolument stable.
3: Soit $u_n$ la valeur de la solution numérique au temps $t=t_n=n h$; si on applique la méthode d’Euler rétrograde, montrer que la solution $u_{n+1}$ au pas $t^{n+1}$ vérifie une équation de la forme suivante:

\[F\left(u_{n+1}\right)=0 \text {. }\]

Ecrire la relation récursive qui définit les itérations $u_{n+1}^{(k)}$ de la méthode de Newton pour résoudre l’équation non-linéaire (4), ave c la donnée initiale $u_{n+1}^{(0)}=u_n$.

Solution

1: Le schéma d’Euler progressif s’écrit

\[u_{n+1}=u_n-h \arctan u_n .\]

Euler rétrograde:

\[v_{n+1}=v_n-h \arctan v_{n+1} .\]

Pour les deux schémas, on pose $u_0=v_0=y(0)=1$.

2: On a:

\[\frac{d}{d y}(-\arctan (y))=-\frac{1}{1+y^2} .\]

Soit $y(t)$ la solution du problème de Cauchy. On nous dit que $y(t)$ est positive. Puisque $y^{\prime}=-\arctan (y)<0$ pour $y>0, y(t)$ est décroissante. Donc $y(t) \in[0,1$(\operatorname{car} y(0)=1)], et si l’on pose $\lambda(t)=-\frac{1}{1+y(t)^2}$, on a

\[-\lambda_{\max } \leq \lambda(t) \leq-\lambda_{\min }\]

où $\lambda_{\min }=\frac{1}{1+1}=\frac{1}{2}>0, \lambda_{\max }=1<\infty$; alors, d’après la théorie, la condition de stabilité est la suivante:

\[h < h_0=\frac{2}{\lambda_{\max }}=2 .\]

3: On pose $F(x)=x-u_n+h \arctan (x)$, et on a alors $F^{\prime}(x)=1+h /\left(1+x^2\right)$; les itérations de Newton s’écrivent

\[\begin{aligned} x^{(k+1)} & =x^{(k)}-\frac{F\left(x^{(k)}\right)}{F^{\prime}\left(x^{(k)}\right)} \\ & =x^{(k)}-\frac{1+\left(x^{(k)}\right)^2}{1+h+\left(x^{(k)}\right)^2}\left(x^{(k)}-u_n+h \arctan \left(x^{(k)}\right)\right) . \end{aligned}\]

Exercice 2

On considère l’équation différentielle d’ordre 2 :

\[\left\{\begin{array}{l} y^{\prime \prime}(t)+c y^{\prime}(t)+k y(t)=0, \quad t>0 \\ y(0)=y_0, \\ y^{\prime}(0)=z_0, \end{array}\right.\]

où $c$ et $k$ sont des constantes réelles positives données.

1: On définit $w_1(t)=y(t)$ et $w_2(t)=y^{\prime}(t)$, puis le vecteur $\mathbf{w}(t)=\left(w_1(t), w_2(t)\right)^T$. Vérifier que l’équation (1) peut être écrite comme le système $2 \times 2$ d’équations différentielles du premier ordre qui suit

\[\mathbf{w}^{\prime}(t)=A \mathbf{w}(t),\]

$A$ étant une matrice carrée de taille 2 .

2: Ecrire les schémas d’Euler progressif et rétrograde pour résoudre le système (2).
3: On considère la diagonalisation de la matrice $A=V D V^{-1}$, où $D$ est la matrice diagonale des valeurs propres de $A$ et $V$ la matrice dont les colonnes sont les vecteurs propres correspondants. On introduit dans les deux schémas numériques trouvés à la question précédente la transformation $\mathbf{x}^n=V^{-1} \mathbf{w}^n$. Récrire les deux schémas précédents en fonction de la nouvelle inconnue $\mathbf{x}^n$.
4: On considère maintenant le cas où $c=5$ et $k=6$. Trouver les valeurs du pas de temps $h$ pour lesquelles la méthode d’Euler progressive est stable en utilisant les résultats pour le cas scalaire. (Sugg: On remarque que le système obtenu au point précédent est diagonal).

Solution

1: En utilisant la transformation proposée, le problème peut s’écrire sous la forme équivalente: On aboutit donc au système d’équations différentielles du premier ordre qui suit

\[\mathbf{w}^{\prime}(t)=A \mathbf{w}(t), \quad \mathbf{w}(0)=\left(y_0, z_0\right)^T,\]

$A$ étant la matrice $2 \times 2$ à coefficients constants

\[A=\left[\begin{array}{cc} 0 & 1 \\ -k & -c \end{array}\right] \text {. }\]

2: Les schémas d’Euler progressif et rétrograde s’écrivent respectivement: Euler progressif $\quad\left\{\begin{array}{l}\mathbf{w}^{n+1}=\mathbf{w}^n+h A \mathbf{w}^n=(I+h A) \mathbf{w}^n \\ \mathbf{w}^0=\mathbf{w}_0=\left(y_0, z_0\right)^T\end{array}\right.$ Euler rétrograde $\left\{\begin{array}{l}\mathbf{w}^{n+1}=\mathbf{w}^n+h A \mathbf{w}^{n+1} \\ \mathbf{w}^0=\mathbf{w}_0=\left(y_0, z_0\right)^T\end{array} \Longrightarrow\left\{\begin{array}{l}(I-h A) \mathbf{w}^{n+1}=\mathbf{w}^n \\ \mathbf{w}^0=\mathbf{w}_0=\left(y_0, z_0\right)^T\end{array}\right.\right.$ $h$ étant le pas d’intégration. On note que le deuxième schéma nécessite la résolution d’un système linéaire $2 \times 2$ à chaque itération.
3: Schéma d’Euler progressif. Si l’on introduit la transformation $\mathbf{w}^n=V \mathbf{x}^n$ dans [eq:3], on obtient

\[V \mathbf{x}^{n+1}=V \mathbf{x}^n+h A V \mathbf{x}^n,\]

En multipliant à gauche tous le termes $\operatorname{par} V^{-1}$, on obtient

\[\mathbf{x}^{n+1}=\mathbf{x}^n+h V^{-1} A V \mathbf{x}^n=(I+h D) \mathbf{x}^n,\]

avec la condition initiale $\mathbf{x}^0=\mathbf{x}_0=V^{-1} \mathbf{w}_0$. Schéma d’Euler rétrograde. En suivant le même raisonnement que pour le schéma explicite, on obtient aisément le système

\[(I-h D) \mathbf{x}^{n+1}=\mathbf{x}^n, \quad \mathbf{x}^0=V^{-1} \mathbf{w}_0\]

4: On s’aperçoit que le système [eq:5] est diagonal. Si on note par $\lambda_1$ et $\lambda_2$ les deux valeurs diagonales de la matrice $D$ (qui sont les valeurs propres de la matrice $A$ ), le système [eq:5] s’écrit

\[\left\{\begin{array}{l} x_1^{n+1}=\left(1+h \lambda_1\right) x_1^n \\ x_2^{n+1}=\left(1+h \lambda_2\right) x_2^n \end{array} .\right.\]

Notons que [eq:5] contient deux équations qui peuvent se résoudre de façon indépendante l’une de l’autre. Dans le cas où $c=5$ et $k=6$ on a

\[\operatorname{det}(\lambda I-A)=\operatorname{det}\left[\begin{array}{cc} \lambda & -1 \\ 6 & \lambda+5 \end{array}\right]=\lambda^2+5 \lambda+6 \quad \Longrightarrow \quad \lambda_1=-3, \quad \lambda_2=-2 .\]

Chaque équation du système [eq:6] correspond à un schéma d’Euler progressif scalaire et donc la condition de stabilité est $h<2 / 3$ pour la première équation et $h<1$ pour la deuxième. La condition à retenir c’est la plus restrictive, donc $h<2 / 3$. Puisque le schéma [eq:6] est équivalent au schéma [eq:3], cette dernière condition est aussi la condition de stabilité pour [eq:3].

Exercice 3

On considère la méthode de Crank-Nicolson

\[\left\{\begin{array}{l} \frac{u_{n+1}-u_n}{h}=\frac{1}{2}\left[f\left(t_n, u_n\right)+f\left(t_{n+1}, u_{n+1}\right)\right], \quad n=0,1, \ldots \\ u_0=y_0, \end{array}\right.\]

$h$ étant le pas de temps, pour approcher la solution du problème de Cauchy

\[\left\{\begin{array}{l} y^{\prime}(t)=f(t, y(t)), \quad t \in(0, T), \\ y(0)=y_0 . \end{array}\right.\]

On suppose que la fonction $f$ soit lipschitzienne par rapport à $y$, c’est-à-dire qu’il existe une constante $L>0$ telle que

\[\left|f\left(t, y_1\right)-f\left(t, y_2\right)\right| \leq L\left|y_1-y_2\right|\]

pour tout $t, y_1, y_2 \in \mathbb{R}$. On rappelle également que, si $y \in C^3(0, T)$, on $\mathrm{a}^1$

\[\frac{y\left(t_{n+1}\right)-y\left(t_n\right)}{h}=\frac{1}{2}\left[f\left(t_n, y\left(t_n\right)\right)+f\left(t_{n+1}, y\left(t_{n+1}\right)\right)\right]+\tau_{n+1}, \quad n=0,1, \ldots\]

où $\left|\tau_{n+1}\right| \leq \tau(h)=\frac{h^2}{12} \max _{t \in(0, T)}\left|y^{\prime \prime \prime}(t)\right|, \forall n=0,1, \ldots$

1: Trouver une équation de récurrence pour l’erreur $e_n=u_n-y_n$ en soustrayant les équations (7) et $(9)$. En supposant que la solution exacte $y(t)$ soit $C^3(0, T)$, prouver que pour $n=0,1, \ldots$ et $h<1 / L$ :

\[\left|e_n\right|=\left|u_n-y\left(t_n\right)\right| \leq C(T) \tau(h), \quad \text { où } C(T) \text { ne dépend que de } T \text {. }\]

2: En déduire l’ordre de convergence de la méthode de Crank-Nicolson. (Sugg. Utiliser (8) et (9) dans l’équation trouvée au point précédent).
3: On a appliqué les méthodes de Crank-Nicolson et d’Euler progressive pour résoudre le problème de Cauchy

\[\left\{\begin{array}{l} y^{\prime}(t)=-t y^2(t), \quad 0 < t < 1 \\ y(0)=2 \end{array}\right.\]

(dont la solution exacte est $y(t)=2 /\left(1+t^2\right)$ ). Le tableau suivant montre les erreurs commises par les deux méthodes à l’instant $t=5$, pour différentes valeurs du pas de temps $h$ :

h	méthode X	méthode Y
1	\(3.052276 \mathrm{e}-02\)	\(5.215832 \mathrm{e}-03\)
0.5	\(1.460741 \mathrm{e}-02\)	\(1.034686 \mathrm{e}-03\)
0.25	\(7.161485 \mathrm{e}-03\)	\(2.497104 \mathrm{e}-04\)
0.125	\(3.548112 \mathrm{e}-03\)	\(6.189158 \mathrm{e}-05\)
0.0625	\(1.765997 \mathrm{e}-03\)	\(1.543977 \mathrm{e}-05\)

méthode X

méthode Y

$3.052276 \mathrm{e}-02$

$5.215832 \mathrm{e}-03$

0.5

$1.460741 \mathrm{e}-02$

$1.034686 \mathrm{e}-03$

0.25

$7.161485 \mathrm{e}-03$

$2.497104 \mathrm{e}-04$

0.125

$3.548112 \mathrm{e}-03$

$6.189158 \mathrm{e}-05$

0.0625

$1.765997 \mathrm{e}-03$

$1.543977 \mathrm{e}-05$

Identifier quelle colonne entre les deux à été calculée en utilisant la méthode de Crank-Nicolson. Justifier la réponse.

4: Trouver la condition de stabilité sur $h$ pour la méthode d’Euler progressive appliquée au problème.
5: Implémenter la méthode de Crank-Nicolson pour résoudre le [eq11:problème] sur l’intervalle $t \in(0,1)$ avec $h=0.1$. Reporter le graphe de la solution numérique obtenue. Comparer les résultats avec la solution exacte.

Solution

1: En soustrayant les équations indiquées on trouve

\[\frac{e_{n+1}-e_n}{h}=\frac{1}{2}\left[f\left(t_n, u_n\right)-f\left(t_n, y\left(t_n\right)\right)\right]+\left[f\left(t_{n+1}, u_{n+1}\right)-f\left(t_{n+1}, y\left(t_{n+1}\right)\right)\right]-\tau_{n+1},\]

où on a utilisé le fait que $y_n-u_n=e_n$. Donc l’équation de récurrence pour l’erreur est:

\[e_{n+1}=e_n+h \frac{1}{2}\left[\left(f\left(t_n, u_n\right)-f\left(t_n, y\left(t_n\right)\right)\right)+\left(f\left(t_{n+1}, u_{n+1}\right)-f\left(t_{n+1}, y\left(t_{n+1}\right)\right)\right)\right]-h \tau_{n+1} .\]

Cette formule découle de l’estimation de l’erreur pour la formule d’intégration des trapèzes

2: A partir de (13), en sachant que la fonction $f$ est lipschitzienne, on trouve

\[\begin{aligned} \left|e_{n+1}\right| \leq\left|e_n\right|+h \frac{1}{2}\left[L\left|u_n-y\left(t_n\right)\right|+L\left|u_{n+1}-y\left(t_{n+1}\right)\right|\right]+ & h \tau(h) \\ & =\left|e_n\right|+\frac{L h}{2}\left[\left|e_n\right|+\left|e_{n+1}\right|\right]+h \tau(h), \end{aligned}\]

c’est-à-dire

\[\left(1-\frac{L h}{2}\right)\left|e_{n+1}\right| \leq\left(1+\frac{L h}{2}\right)\left|e_n\right|+h \tau(h) .\]

Si $h<1 / L$, on a $1-h L / 2>0$, donc on peut diviser par $1-h L / 2$ et écrire

\[\left|e_{n+1}\right| \leq a\left|e_n\right|+b, \quad n=1,2, \ldots\]

où

\[a=\frac{1+h L / 2}{1-h L / 2}=1+\frac{h L}{1-h L / 2}, \quad b=\frac{h}{1-h L / 2} \tau(h) .\]

Or, comme $e_0=y_0-u_0=0$, en appliquant $n$ fois la dernière inégalité on a:

\[\begin{aligned} \left|e_n\right| & \leq a\left|e_{n-1}\right|+b \\ & \leq a\left[a\left|e_{n-2}\right|+b\right]+b \\ & \vdots \\ & \leq\left[1+a+a^2+\ldots a^{n-1}\right] b \\ & =\frac{a^n-1}{a-1} b=\left[\left(1+\frac{h L}{1-h L / 2}\right)^n-1\right] \frac{1}{L} \tau(h) . \end{aligned}\]

On peut ensuite remarquer que, si $h<1 / L$, on a

\[\frac{h L}{1-h L / 2}<\frac{h L}{1-1 / 2}=2 h L\]

donc (par la même démarche que l’on a vue au cours pour la convérgence de la méthode d’Euler

\[\begin{aligned} \left|e_n\right| \leq\left[(1+2 h L)^n-1\right] \frac{1}{L} \tau(h) \leq\left(e^{2 h L n}-1\right) \frac{1}{L} \tau(h) & = \\ & \left(e^{2 L t_n}-1\right) \frac{1}{L} \tau(h) \leq\left(e^{2 L T}-1\right) \frac{1}{L} \tau(h), \end{aligned}\]

et finalement on obtient

\[\left|e_n\right| \leq C(T) \tau(h)\]

où

\[C(T)=\frac{e^{2 T L}-1}{L} .\]

La méthode de Cranck-Nicholson est donc d’ordre 2 ; en effet, on a $\tau(h)=\frac{h^2}{12} \max _{t \in[0, T$}\left|y^{\prime \prime \prime}(t)\right|] et donc l’erreur est majorée par une quantité proportionelle au carré du pas de temps:

\[\left|e_n\right| \leq\left[\frac{e^{2 T L}-1}{12 L} \max _{t \in[0, T]}\left|y^{\prime \prime \prime}(t)\right|\right] h^2\]

3: On a vu au point précédent que la méthode de Crank-Nicholson est d’ordre 2. Or, si l’on regarde le tableau, on remarque que, en divisant par deux le pas de temps, l’erreur dans la première colonne est (approximativement) divisée par 2, tandis que dans la seconde colonne elle est divisée par 4. On en déduit que la première colonne correspond à une méthode d’ordre 1 (par exemple Euler progressive ou Euler rétrograde) et la seconde à une méthode d’ordre 2 (donc celle de Crank-Nicholson).
4: On a:

\[\lambda(t)=\frac{\partial f}{\partial y}(t, y(t))=-2 t y(t)=-\frac{4 t}{1+t^2} \leq 0 \quad \text { pour } t \geq 0 .\]

On remarque que $\lambda^{\prime}(t)=4 \frac{1-t^2}{\left(1+t^2\right)^2}$ : les zéros de cette dérivé sont en $t= \pm 1$, donc

\[\max _t|\lambda(t)|=2 .\]

On peut donc utiliser la condition de stabilité pour la méthode d’Euler progressive que l’on a vue au cours: on a

\[h < \frac{2}{\max _t|\lambda(t)|}=1\]

Exercice 4

On considère l’équation différentielle linéaire suivante:

\[\left\{\begin{array}{l} y^{\prime}(t)=\lambda y(t), \quad t \in[0,20], \\ y(0)=1, \end{array}\right.\]

où $\lambda=-5$. On veut trouver la solution numérique de l’équation avec la méthode de Heun, avec un pas $h$ de discrétisation temporelle (et donc un nombre $N_h=\frac{20}{h}$ de sous-intervalles de $[0,20]$). Dans la suite, on désigne $u_n$ une approximation de $y\left(t_n\right)$ au temps $t_n=n h, n=0,1,2, \ldots$

La méthode d’Heun progressive n’est pas stable pour $h$ quelconque; pour assurer la stabilité, il faut que $h$ soit plus petit que $\bar{h}=\frac{2}{|\lambda|}$.

1: Calculer $\bar{h}$. Résoudre l’équation différentielle avec la méthode de Heun (commande heun) en choisissant $h=\bar{h}$ et $h=\bar{h} / 4$. Visualiser les graphes des solutions numériques en utilisant axis pour choisir les intervalles $[0,20]$ en abscisse et $[-0.5,1.5]$ en ordonnée. Commenter les résultats obtenus dans les deux cas.
2: Résoudre l’équation différentielle avec la méthode d’Euler progressive (commande feuler) et rétrograde (commande beuler) avec les mêmes valeurs de $N_h$ utilisées au point 1). Comparer les résultats et discuter la stabilité de chaque méthode.
3: Programmer une fonction meuler, avec la déclaration suivante function $[t, y]=\operatorname{meuler}(f, t s p a n, y, N h)$ dans le fichier meuler.m, qui implémente la méthode d’Euler modifiée:

\[u_{n+1}=u_n+h f\left(t_n+\frac{h}{2}, u_n+\frac{h}{2} f\left(t_n, u_n\right)\right), \quad n=0,1, \ldots, N_h-1,\]

pour l’approximation numérique du problème de Cauchy

\[\left\{\begin{array}{l} y^{\prime}(t)=f(t, y(t)), \quad t \in\left[t_0, t_1\right], \\ y(0)=y_0 . \end{array}\right.\]

Les variables d’entrée $\mathbf{f}$, tspan, y 0 et $\mathrm{Nh}$ de meuler sont respectivement la fonction inline $f(t, y)$, le vecteur $\left[t_0, t_1\right]$, la valeur initiale $y_0$, le nombre de pas temporals $N_h$. Suggestion: utilisez la même structure de heun.py pour écrire meuler.y. Utiliser votre code pour le problème (13), et calculer l’erreur

\[E_h=\max _n\left|y\left(t_n\right)-u_n\right|\]

pour $h=0.1$ et $h=0.01$. On rappelle que une méthode est dite d’ordre $p, p>0$, s’il existe $C>0$ tel que $E_h \leq C h^p$. Si l’on suppose $E_h \simeq C h^p$, d’après le calcul des erreurs, quelle est l’ordre de la méthode d’Euler modifiée?

Solution

1: Nous avons:

\[\bar{h}=\frac{2}{5}=0.4\]

Les valeurs de $N_h$ requises sont $N_{h_1}=20 / 0.4=50$ et $N_{h_2}=20 / 0.1=200$:

\[\mathrm{Nh} 1=50 ; \mathrm{Nh} 2=200 ;\]

Pour la méthode de Heun, utilisons la méthode de Heun pour calculer la solution numérique.

import numpy as np
import plotly.graph_objects as go
from scipy.integrate import solve_ivp
from tan.ode.heun import heun

# Définition de la fonction de l'EDO
def odefun(t, y):
    return -5 * y

# Paramètres
tspan = np.array([0, 20])
y0 = [1]
h1 = 0.4  # h_bar
Nh1 = int(20 / h1)
h2 = h1 / 4
Nh2 = int(20 / h2)

# Résolution par la méthode de Heun (à implémenter ou utiliser solve_ivp avec méthode 'Heun')
# Exemple d'utilisation de solve_ivp pour le contexte
t1,sol1 = heun(odefun, tspan, y0, Nh1)
t2,sol2 = heun(odefun, tspan, y0, Nh2)

# Plotting
fig = go.Figure()
fig.add_trace(go.Scatter(x=t1, y=sol1[:,0], mode='lines', name='Heun Nh=50'))
fig.add_trace(go.Scatter(x=t1, y=sol2[:,0], mode='lines', name='Heun Nh=200'))

fig.update_layout(title="Solution Numérique par la Méthode de Heun",
                  xaxis_title="Temps (t)",
                  yaxis_title="Solution (y)",
                  legend_title="Méthode")
fig.show()

Results

La Figure 1 affiche les résultats. Pour $h=\bar{h}$, la solution numérique est constante et ne tend pas vers zéro lorsque $t \rightarrow \infty$ : en effet, la condition de stabilité sur le pas de discrétisation $h$ n’est pas satisfaite. Pour $h=\bar{h} / 4$, c’est-à-dire $N_h=200$, la condition $h<\bar{h}$ est satisfaite et la solution numérique tend vers zéro.

2: Méthodes de Euler Progressive et Rétrograde. Appliquons la même démarche avec les méthodes Euler progressive et rétrograde. (Le code Python pour ces méthodes serait similaire, en ajustant la méthode de résolution dans solve_ivp ou en implémentant les méthodes spécifiquement.)

# Définition de la fonction de l'EDO
def odefun(t, y):
    return -5 *y
from tan.ode.feuler import feuler
from tan.ode.beuler import beuler

# Résolution par la méthode d'Euler progressive
t1,sol1 = feuler(odefun, tspan, y0, Nh1)
t2,sol2 = feuler(odefun, tspan, y0, Nh2)

# Résolution par la méthode d'Euler rétrograde
t3,sol3 = beuler(odefun, tspan, y0, Nh1)
t4,sol4 = beuler(odefun, tspan, y0, Nh2)

# Plotting
fig = go.Figure()
fig.add_trace(go.Scatter(x=t1, y=sol1[:,0], mode='lines', name='Euler Progressif Nh=50'))
fig.add_trace(go.Scatter(x=t2, y=sol2[:,0], mode='lines', name='Euler Progressif Nh=200'))
fig.add_trace(go.Scatter(x=t3, y=sol3[:,0], mode='lines', name='Euler Rétrograde Nh=50'))
fig.add_trace(go.Scatter(x=t4, y=sol4[:,0], mode='lines', name='Euler Rétrograde Nh=200'))

fig.update_layout(title="Solution Numérique par les Méthodes d'Euler",
                  xaxis_title="Temps (t)",
                  yaxis_title="Solution (y)",
                  legend_title="Méthode")
fig.show()

Results

3: Euler Modifié, voici le code

import numpy as np
import plotly.graph_objects as go

tspan = np.array([0, 20])
y0 = [1]
h1 = 0.4  # h_bar
Nh1 = int(20 / h1)
h2 = h1 / 4
Nh2 = int(20 / h2)

def meuler(f, tspan, y0, Nh):
    h = (tspan[1] - tspan[0]) / Nh  # le pas de temps
    t = np.linspace(tspan[0], tspan[1], Nh + 1)  # le vecteur des temps tn
    y = np.zeros((Nh + 1, len(y0)))
    y[0] = y0
    for i in range(Nh):
        y_mid = y[i] + h / 2 * f(t[i], y[i])  # Calcul du point milieu
        y[i + 1] = y[i] + h * f(t[i] + 0.5 * h, y_mid)  # Mise à jour avec la pente au point milieu
    return t, y

# Exemple d'utilisation avec une fonction d'EDO simple
def odefun(t, y):
    return -5*y

t1, y1 = meuler(odefun, tspan, y0, Nh1)
t2, y2 = meuler(odefun, tspan, y0, Nh2)

# Tracer avec Plotly
fig = go.Figure()
fig.add_trace(go.Scatter(x=t1, y=y1[:, 0], mode=f'lines+markers', name='Méthode Euler Modifiée h={h1}'))
fig.add_trace(go.Scatter(x=t2, y=y2[:, 0], mode=f'lines+markers', name='Méthode Euler Modifiée h={h2}'))
fig.update_layout(title='Solution Numérique par la Méthode d\'Euler Modifiée',
                  xaxis_title='Temps',
                  yaxis_title='Solution y(t)',
                  legend_title='Méthode')
fig.show()

Results

Soit $h_1=0.1, h_2=0.01 ;$ comme $h_1 / h_2=10$, on aura $p \simeq \log {10}\left(E{h_1} / E_{h_2}\right)$. Alors on calcule les solutions numériques et les erreurs dans les deux cas:

import numpy as np
yexact = lambda t: np.exp(-5 * t)
# Calcul des solutions numériques
t1, y1 = meuler(odefun, tspan, y0, Nh1)
t2, y2 = meuler(odefun, tspan, y0, Nh2)

# Calcul des erreurs
Eh1 = np.max(np.abs(y1[:,0] - yexact(t1)))
Eh2 = np.max(np.abs(y2[:,0] - yexact(t2)))
print(f"Erreur pour h1={h1} : {Eh1}")
print(f"Erreur pour h2={h2} : {Eh2}")
# Calcul de l'ordre de la méthode
p = np.log10(Eh1 / Eh2) / np.log10(h1 / h2)
print(f"ordre de Eurler Modifié: {p}")

Results

La méthode d’Euler modifiée est d’ordre 2.

Exercice 5

On considère l’équation différentielle suivante:

\[\left\{\begin{array}{l} y^{\prime}(t)=f(y(t))=-\operatorname{atan}(y(t)), \quad t \in(0,30) \\ y(0)=1 \end{array}\right.\]

On veut trouver la solution numérique de l’équation avec la méthode d’Euler progressive, avec un pas $h$ de discrétisation temporelle, aux temps $t_n=n h, n=1,2, \ldots, N_h$.

La méthode d’Euler progressive n’est pas stable pour toutes valeurs de $h$; pour assurer la stabilité, il faut que $h$ soit plus petit que $\bar{h}=\frac{2}{\max _{t \in[0,30]}|\lambda(t)|}$, où $\lambda(t)=\left|\frac{\partial f(y)}{\partial y}\right|_{y=y(t)}$.

1

Calculer $\bar{h}$, en utilisant l’approximation $y(30) \simeq 0$.

2

Appliquer la méthode d’Euler progressive (commande feuler) dans les deux cas suivants:

$\operatorname{avec} h=1.5 \cdot \bar{h}$;
avec $h=0.25 \cdot \bar{h}$. Visualiser la solution numérique dans les deux cas. Commenter les différences entre les résultats obtenus.

3

On veut maintenant résoudre l’équation différentielle avec la méthode de Heun (commande heun) pour $h=1.5 \cdot \bar{h}$ et $h=0.25 \cdot \bar{h}$. Visualiser les graphes des solutions numériques en utilisant axis pour choisir les intervalles $[0,30$] en abscisse et $[-1,1]$ en ordonnée. Commenter les résultats obtenus dans les deux cas.

4

On veut maintenant résoudre l’équation différentielle avec la méthode de Runge-Kutta d’ordre 4 (commande rk4) pour $h=1.5 \cdot \bar{h}$ et $h=0.25 \cdot \bar{h}$. Visualiser les graphes des solutions numériques en utilisant axis pour choisir les intervalles $[0,30$] en abscisse et $[-1,1]$ en ordonnée. Commenter les résultats obtenus dans les deux cas.

Solution

== Réponse à l’Exercice sur les Équations Différentielles Ordinaires

1: On a :

\[\frac{\partial f(y)}{\partial y} = - \frac{1}{1+y^2},\]

donc

\[\bar{h} = \frac{2}{\max_{t\in[0,30]} \frac{1}{1+y(t)^2} } \leq \frac{2}{1} = 2;\]

mais comme $y(30) \simeq 0$, on déduit que $\bar{h} \simeq 2$.

2: Soit $h = 1.5 \bar{h} = 3$; avec les commandes qui suivent on calcule la solution numérique par la méthode d’Euler progressive.

import numpy as np
import plotly.graph_objects as go

# Définition de la fonction dérivée
def dydt(t, y):
    return -1 / (1 + y**2)

# Paramètres
h_large = 3  # h = 1.5 bar{h}
h_petit = 0.5  # h = 0.25 bar{h}
tspan = np.array([0, 30])
y0 = [0]  # Condition initiale

# Solution numérique pour h_large
t_large = np.arange(tspan[0], tspan[1] + h_large, h_large)
t_large, y_large = feuler(dydt, tspan, y0, int(30 / h_large))


# Solution numérique pour h_petit
t_petit = np.arange(tspan[0], tspan[1] + h_petit, h_petit)
t_petit, y_petit = feuler(dydt, tspan, y0, int(30 / h_petit))

# Tracer les solutions numériques
fig = go.Figure()
fig.add_trace(go.Scatter(x=t_large, y=np.squeeze(y_large), mode='lines', name='h=1.5bar{h}'))
fig.add_trace(go.Scatter(x=t_petit, y=np.squeeze(y_petit), mode='lines', name='h=0.25bar{h}'))
fig.update_layout(title="Méthode d'Euler Progressive pour h=1.5bar{h} et h=0.25bar{h}",
                  xaxis_title="Temps",
                  yaxis_title="y(t)",
                  legend_title="Valeur de h")
fig.show()

Results

Les graphes que l’on obtient sont montrés. D’après le graphique, on peut vérifier que la méthode d’Euler progressive n’est pas stable pour tout choix de $h$: si on utilise la méthode avec $h=1.5\bar{h}$, on obtient des solutions avec des oscillations non-amorties. Par contre, avec $h = 0.25 \bar{h}$ on est dans la région de stabilité absolue de la méthode : les perturbations sont contrôlées, en particulier on n’a pas d’oscillations.

3: On veut maintenant résoudre l’équation différentielle avec la méthode de Heun (commande heun) pour $h=1.5 \cdot \bar{h}$ et $h=0.25 \cdot \bar{h}$.

from tan.ode.heun import heun

# Résolution par la méthode de Heun
t_large, y_large = heun(dydt, tspan, y0, int(30 / h_large))
t_petit, y_petit = heun(dydt, tspan, y0, int(30 / h_petit))

# Tracer les solutions numériques
fig = go.Figure()
fig.add_trace(go.Scatter(x=t_large, y=np.squeeze(y_large), mode='lines', name='h=1.5bar{h}'))
fig.add_trace(go.Scatter(x=t_petit, y=np.squeeze(y_petit), mode='lines', name='h=0.25bar{h}'))
fig.update_layout(title="Méthode de Heun pour h=1.5bar{h} et h=0.25bar{h}",
                  xaxis_title="Temps",
                  yaxis_title="y(t)",
                  legend_title="Valeur de h")
fig.show()

Results

4: On veut maintenant résoudre l’équation différentielle avec la méthode de Runge-Kutta d’ordre 4 (commande rk4) pour $h=1.5 \cdot \bar{h}$ et $h=0.25 \cdot \bar{h}$.

from tan.ode.rk4 import rk4

# Résolution par la méthode de Runge-Kutta d'ordre 4
t_large, y_large = rk4(dydt, tspan, y0, int(30 / h_large))
t_petit, y_petit = rk4(dydt, tspan, y0, int(30 / h_petit))

# Tracer les solutions numériques
fig = go.Figure()
fig.add_trace(go.Scatter(x=t_large, y=np.squeeze(y_large), mode='lines', name='h=1.5bar{h}'))
fig.add_trace(go.Scatter(x=t_petit, y=np.squeeze(y_petit), mode='lines', name='h=0.25bar{h}'))
fig.update_layout(title="Méthode de Runge-Kutta d'ordre 4 pour h=1.5bar{h} et h=0.25bar{h}",
                  xaxis_title="Temps",
                  yaxis_title="y(t)",
                  legend_title="Valeur de h")
fig.show()

Results

Exercice 6

L’évolution en temps d’une population $y$ de bactéries peut être obtenue en résolvant une équation différentielle du type

\[\left\{\begin{array}{l} y^{\prime}(t)=y(t) \cdot[y(t) \cdot(1-0.2 y(t))-z], \quad t>0 \\ y(0)=0.5 \end{array}\right.\]

où $z$ représente le nombre de phagocytes. On considère d’abord une population de phagocytes constante $z=0.1$ (toutes les unités sont normalisées).

1: Vérifier l’existence d’un état d’équilibre pour $t \rightarrow \infty$ en résolvant le problème (2) par la méthode d’Euler progressive (commande feuler) sur l’intervalle [0,50], avec les pas de temps $h=0.5$ et $h=0.1$. Au vu des graphes, y a-t-il une condition de stabilité à respecter sur le pas de temps $h$ ?
2: Si la population $z$ des phagocytes n’est pas constante, le modèle dynamique (2) peut être remplacé par le modèle suivant, dit de Odell:

\[\begin{cases}y^{\prime}(t)=y(t) \cdot[y(t) \cdot(1-0.2 y(t))-z], & t>0 \\ z^{\prime}(t)=z(t) \cdot(y(t)-1), & t > 0 \\ y(0)=0.5 & \\ z(0)=0.1 & \end{cases}\]

Résoudre ce modèle sur l’intervalle $[0,50$] avec $h=0.1$, puis reporter les graphes des populations $y(t)$ et $z(t)$. Est-ce que selon ce nouveau modèle les phagocytes arrivent à éliminer les bactéries?

3: On veut aller plus loin dans la modélidation de l’Interaction entre Bactéries et Phagocytes. Le système d’équations différentielles ordinaires (EDO) pour les bactéries et les phagocytes peut être utilisé pour modéliser les interactions entre ces deux populations dans un environnement. Un modèle simple pourrait prendre la forme suivante :
Bactéries ($y1$): La population de bactéries croît de manière logistique (croissance limitée par la capacité de l’environnement) et est réduite par la prédation des phagocytes.
Phagocytes ($y2$): La population de phagocytes croît en réponse à la présence de bactéries (leur source de nourriture) et décroît naturellement à un certain taux.

Ce modèle peut être représenté par le système d’EDO suivant :

\[\begin{align*} \frac{dy_1}{dt} &= y_1 \left( r_1 \left(1 - \frac{y_1}{K}\right) - a y_2 \right) \\ \frac{dy_2}{dt} &= r_2 y_2 \left(1 - \frac{y_2}{y_1}\right) \end{align*}\]

où : - $y1$ est la population de bactéries, - $y2$ est la population de phagocytes, - $r1$ est le taux de croissance intrinsèque des bactéries, - $K$ est la capacité de charge de l’environnement pour les bactéries, - $a$ est le taux auquel les phagocytes consomment les bactéries, - $r2$ est le taux de croissance des phagocytes en réponse à la disponibilité des bactéries.

Le modèle suppose que les bactéries ont une croissance logistique limitée par la capacité de l’environnement $K$ et que les phagocytes dépendent directement de la population de bactéries pour leur croissance. Le taux $a$ représente l’efficacité des phagocytes à réduire la population de bactéries.

Ce modèle peut être ajusté ou modifié pour inclure d’autres effets biologiques ou écologiques, tels que l’immunité spécifique, l’introduction de phages (virus attaquant les bactéries), ou d’autres mécanismes de régulation de la population.

Le choix des valeurs des paramètres pour un modèle impliquant des bactéries et des phagocytes peut avoir un impact significatif sur la dynamique observée dans le système. Voici quelques gammes ou valeurs typiques que vous pourriez envisager lors de la simulation du système. Gardez à l’esprit que ces valeurs devraient être ajustées en fonction du contexte spécifique de votre modèle, du réalisme biologique que vous recherchez et de toutes données empiriques que vous pourriez avoir :

Taux de croissance bactérien (r_1) : Il s’agit souvent d’une valeur positive reflétant la rapidité avec laquelle la population bactérienne peut croître dans des conditions idéales. Une gamme courante pour r_1 pourrait être de 0,1 à 2 par jour, selon le type de bactéries et les conditions environnementales.
Capacité porteuse (K) : Cela représente la taille maximale de la population de bactéries que l’environnement peut soutenir. Cela peut varier largement en fonction de la disponibilité des nutriments, de l’espace et d’autres facteurs environnementaux. Une gamme typique pourrait être de 10^3 à 10^9 cellules par mL ou par gramme de tissu, selon l’échelle et le contexte du modèle.
Taux de prédation (a) : Ce paramètre représente le taux auquel les phagocytes consomment les bactéries. Il pourrait varier de 0,001 à 0,1 par phagocyte par jour, reflétant l’efficacité avec laquelle les phagocytes peuvent réduire les populations bactériennes.
Taux de réponse à la croissance des phagocytes (r_2) : C’est le taux auquel la population de phagocytes peut augmenter en réponse à la présence de bactéries. Sa valeur pourrait varier de 0,05 à 1 par jour, indiquant la rapidité avec laquelle les phagocytes peuvent répondre à une augmentation de la population bactérienne.
Populations initiales : Les nombres de départ de bactéries (y_1) et de phagocytes (y_2) sont également cruciaux. Vous pourriez commencer avec un petit nombre de bactéries, disons 10^2 à 10^5 cellules, et un nombre relativement plus faible de phagocytes, puisque les phagocytes sont typiquement moins nombreux mais augmentent en réponse à la présence bactérienne.

Par exemple, vous pourriez commencer avec : - r_1 = 1.0 par jour (taux de croissance bactérien), - K = 10^6 cellules/mL (capacité porteuse), - a = 0.01 par phagocyte par jour (taux de prédation), - r_2 = 0.1 par jour (taux de réponse à la croissance des phagocytes), - Population bactérienne initiale y_1 = 10^5 cellules, - Population de phagocytes initiale y_2 = 10^2 cellules.

Implémenter le modèle et tester différentes valeurs des paramètres, Reporter votre exploration du modèle dans des graphes des populations $y(t)$ et $z(t)$. Commenter les résultats obtenus. NOTE: il sera profitable d’utiliser l’échelle semi-logarithmique pour visualiser les graphes du fait variation rapide des populations.

Solution

1: On a une équation différentielle décrivant la dynamique de la population :

import numpy as np
import plotly.graph_objects as go
from tan.ode.feuler import feuler

# Définition de la fonction dérivée
def odefun(t, y):
    return y * (y * (1 - y / 5) - 0.1)

y0 = np.array([1])
T = 10.

# Premier cas avec h = 0.5
h = 0.5
Nh = int(T / h)
t = np.linspace(0, T, Nh+1)
t, y = feuler(odefun, [0, T], y0, Nh)
print(y)

fig = go.Figure()
fig.add_trace(go.Scatter(x=t, y=y[:,0], mode='lines', name='h=0.5', line=dict(dash='dot')))

# Deuxième cas avec h = 0.1
h = 0.1
Nh = int(T / h)
t1 = np.linspace(0, T, Nh+1)
t1, y1 = feuler(odefun, [0, T], y0, Nh)

fig.add_trace(go.Scatter(x=t1, y=y1[:,0], mode='lines', name='h=0.1'))

fig.update_layout(title="Dynamique de la Population pour Différents Pas de Temps",
                  xaxis_title="Temps",
                  yaxis_title="Population",
                  legend_title="Valeur de h")
fig.show()

Results

La première remarque est qu’il existe une condition de stabilité à respecter : pour $h=0.5$, la méthode n’est pas stable tandis que pour $h=0.1$ elle l’est. Ensuite, la population atteint un équilibre à $y_\infty = 5$.

2: On considère maintenant un système de deux équations différentielles représentant les interactions entre les bactéries et les phagocytes.

def odefun_system(t, y):
    return np.array([y[0] * (y[0] * (1 - y[0] / 5) - y[1]), y[1] * (y[0] - 1)])

y0 = np.array([0.5, 0.1])
h = 0.1
Nh = int(T / h)
t = np.linspace(0, T, Nh+1)
t,y=feuler(odefun_system, [0, T], y0, Nh)

fig = go.Figure()
fig.add_trace(go.Scatter(x=t, y=y[:, 0], mode='lines', name='Bactéries'))
fig.add_trace(go.Scatter(x=t, y=y[:, 1], mode='lines', name='Phagocytes'))

fig.update_layout(title="Dynamique des Populations de Bactéries et Phagocytes",
                  xaxis_title="Temps",
                  yaxis_title="Population",
                  legend_title="Espèce")
fig.show()

Results

Selon le nouveau modèle, les phagocytes réussissent à éliminer les bactéries.

3: voici le code pour le modèle de bactéries et phagocytes, pour étudier les différentes valeurs des paramètres. Nous implémentons le modèle dans une fonction modele qui prend en entrée les paramètres et retourne les solutions numériques afin de pouvoir explorer les différentes valeurs des paramètres.

import numpy as np
import plotly.graph_objects as go
from tan.ode.feuler import feuler

def modele(h,T,r1, r2, K, alpha):
    def odefun_system_logistic(t, y, r1, r2, K, alpha):
        return np.array([y[0] * (r1*((1 - y[0] / K) - alpha*y[1])), r2*y[1] * ((1 - y[1] / y[0]))])

    T=100

    y0 = np.array([1e5, 1e2])
    h = 0.1
    Nh = int(T / h)
    t = np.linspace(0, T, Nh+1)

    t,y=feuler(odefun_system_logistic, [0, T], y0, Nh, r1, r2, K, alpha)

    return t,y

r1=1 # taux de croissance des bactéries en l'absence de phagocytes (1/jour) : entre 0.1 et 2
r2=0.1 # taux de croissance des phagocytes en l'absence de bactéries (1/jour) : entre 1e-3 et 1
K=1e6 # capacité de charge de l'environnement (nombre de bactéries) : entre 1e3 et 1e9
alpha=1e-2 # taux de prédation des phagocytes sur les bactéries (1/jour) : entre 1e-3 et 0.1

t,y=modele(0.1,100,r1, r2, K, alpha)

fig = go.Figure()
fig.add_trace(go.Scatter(x=t, y=y[:, 0], mode='lines', name='Bactéries'))
fig.add_trace(go.Scatter(x=t, y=y[:, 1], mode='lines', name='Phagocytes'))

fig.update_layout(title="Dynamique des Populations de Bactéries et Phagocytes avec Régulation Logistique",
                  xaxis_title="Temps",
                  yaxis_title="Population",
                  legend_title="Espèce")
# semi-logarithmic plot
fig.update_layout(yaxis_type="log")
fig.show()

Results

Nous pouvons maintenant explorer les différentes valeurs des paramètres pour le modèle de bactéries et phagocytes. Par exemple, pour $r_1=1$, $r_2=0.1$, $K=10^6$ et $\alpha=10^{-2}$, nous obtenons une dynamique des populations de bactéries et phagocytes comme montré ci-dessous

fig = go.Figure()
for alpha in [1e-3, 1e-2, 1e-1]:
    t,y=modele(0.05,100,np.array([1e5, 1e2]),r1, r2, K, alpha)

    fig.add_trace(go.Scatter(x=t, y=y[:, 0], mode='lines', name=f'Bactéries alpha={alpha}'))
    fig.add_trace(go.Scatter(x=t, y=y[:, 1], mode='lines', name=f'Phagocytes alpha={alpha}'))

fig.update_layout(title=f"Dynamique des Populations de Bactéries et Phagocytes avec alpha={alpha}",
                      xaxis_title="Temps",
                      yaxis_title="Population",
                      legend_title="Espèce")
    # semi-logarithmic plot
fig.update_layout(yaxis_type="log")
fig.show()

Results

Nous pouvons voir que pour des valeurs plus élevées de $\alpha$, la population de bactéries diminue plus rapidement, tandis que la population de phagocytes augmente plus rapidement. Cela suggère que les phagocytes sont plus efficaces pour éliminer les bactéries lorsque $\alpha$ est plus élevé.

Explorons maintenant les effets de la capacité de charge $K$ sur la dynamique des populations de bactéries et phagocytes. Nous démarrons avec des populations initiales de bactéries et de phagocytes de $10^3$ et $10^1$ respectivement, et nous explorons différentes valeurs de $K$ dans une échelle semi-logarithmique.

fig = go.Figure()
for K in [1e3, 1e6, 1e9]:
    t,y=modele(0.05,100,np.array([1e3, 1e1]),r1, r2, K, alpha)

    fig.add_trace(go.Scatter
    (x=t, y=y[:, 0], mode='lines', name=f'Bactéries K={K}'))
    fig.add_trace(go.Scatter(x=t, y=y[:, 1], mode='lines', name=f'Phagocytes K={K}'))

fig.update_layout(title=f"Dynamique des Populations de Bactéries et Phagocytes avec K={K}",
                      xaxis_title="Temps",
                      yaxis_title="Population",
                      legend_title="Espèce")
# semi-logarithmic plot
fig.update_layout(yaxis_type="log")
fig.show()

Results

Nous pouvons voir que pour des valeurs plus élevées de $K$, la population de bactéries atteint des niveaux plus élevés et diminue plus lentement, tandis que la population de phagocytes atteint des niveaux légèrement plus élevés et augmente plus rapidement. Cela suggère que les phagocytes sont plus efficaces pour éliminer les bactéries lorsque $K$ est plus élevé.

Dans les graphes ci-dessus, nous avons utilisé une échelle semi-logarithmique pour visualiser les populations de bactéries et de phagocytes, car les populations peuvent varier sur plusieurs ordres de grandeur. Cela nous permet de mieux visualiser les dynamiques des populations sur une échelle plus large.

en Master CSMI nous avons un cours de traitement des incertitudes et de l’optimisation, qui permet de faire des analyses de sensibilité et d’optimisation des paramètres. Ces analyses peuvent être utilisées pour explorer les effets des paramètres sur la dynamique du système et pour identifier les paramètres les plus influents sur les résultats du modèle et ceci de manière systématique.