Systèmes d’EDO

1. Extension au cas de systèmes

Considérons un système non homogène d’équations différentielles ordinaires linéaires à coefficients constants:

\[\label{sys-lin_edo} \begin{cases} \mathbf{y}'(t) = A\mathbf{y}(t) + \mathbf{b}(t) & t>0,\\ \mathbf{y}(0) = \mathbf{y}_0, \end{cases}\]

avec \(A \in \mathbb{R}^{p \times p}\) et \(\mathbf{b} (t) \in \mathbb{R}^p\), où l’on suppose que \(A\) possède \(p\) valeurs propres distinctes \(\lambda_j\), \(j=1,\ldots,p\).

Du point de vue numérique, les méthodes introduites dans le cas scalaire peuvent être étendues aux systèmes d’équations différentielles. Par exemple, le schéma d’Euler progressif ([ep]) devient:

\[\left\{ \begin{array}{l} \displaystyle{ \frac{\mathbf{u}_{n+1} - \mathbf{u}_n}{h} = A \mathbf{u}_{n} + \mathbf{b}_n} \qquad \text{pour} \; n=0,1,2,\ldots \\ \mathbf{u}_0 = \mathbf{y}_0 \; , \end{array} \right.\]

tandis que la méthode d’Euler rétrograde ([er]) devient

\[\left\{ \begin{array}{l} \displaystyle{ \frac{\mathbf{u}_{n+1} - \mathbf{u}_n}{h} = A \mathbf{u}_{n+1} + \mathbf{b}_{n+1}} \qquad \text{pour} \; n=0,1,2,\ldots \\ \mathbf{u}_0 = \mathbf{y}_0 \; , \end{array} \right.\]

En ce qui concerne la stabilité:

si \(\mathbf{b}\equiv \mathbf{0}\) et les valeurs propres \(\lambda_j\) (\(j=1,\ldots,p\)) de la matrice \(A\) sont strictement négatives: \(\lambda_j < 0\), \(j=1,\ldots,p\), alors \(\mathbf{y}(t) \to \mathbf{0}\) lorsque \(t \to + \infty\), et la méthode d’Euler progressive est stable (c-à-d \(\mathbf{u}_n \to \mathbf{0}\) si \(n \to + \infty\)) pourvu que

\[h < \frac{2}{\max_{j=1,\ldots,p} | \lambda_j |} = \frac{2}{\rho (A)} \; ,\]

où \(\rho (A)\) est le rayon spectral de \(A\), tandis que la méthode d’Euler rétrograde est inconditionnellement stable.

2. Système linéaire

Exemple: Système linéaire

Le système

\[\begin{aligned} y_1'(t) & = & -2 y_1(t) + y_2(t) + e^{-t} \nonumber\\ y_2'(t) & = & 3y_1(t) -4 y_2(t) \label{exemple-de-sys-lin-ode} \end{aligned}\]

avec les conditions initiales \(y_1(0)= y_{10}\), \(y_2(0)= y_{20}\), s’écrit sous la forme ([sys-lin_edo]), où

\[\mathbf{y}(t) = \begin{bmatrix} y_1(t) \\ y_2(t) \end{bmatrix}, \quad A = \begin{bmatrix} -2 & 1 \\ 3 & -4 \end{bmatrix}, \quad \mathbf{b}(t) = \begin{bmatrix} e^{-t} \\ 0 \end{bmatrix}, \quad \mathbf{y}_0 = \begin{bmatrix} y_{10} \\ y_{20} \end{bmatrix}.\]

Soit \(h>0\) le pas de temps; pour \(n\in\mathbb{N}\) on pose \(t_n = nh\), \(\mathbf{b}_n = \mathbf{b}(t_n)\), et on désigne \(\mathbf{u}_n\) une valeur approchée de la solution exacte \(\mathbf{y}(t_n)\) au temps \(t_n\).

Les schémas d’Euler progressif, rétrograde et de Crank-Nicolson pour approcher la solution \(\mathbf{y}(t)\) de ([exemple-de-sys-lin-ode]) s’écrivent respectivement:

\[\begin{aligned} &\text{Euler progressif } & &\begin{cases} \mathbf{u}_{n+1} = \mathbf{u}_n + hA\mathbf{u}_n +h\mathbf{b}_n = (I+hA)\mathbf{u}_n +h\mathbf{b}_n \end{cases} \\ &\text{Euler r\'etrograde } & &\begin{cases} (I-hA)\mathbf{u}_{n+1} = \mathbf{u}_n +h\mathbf{b}_{n+1} \\ \mathbf{u}_0 = \mathbf{y}_0 \end{cases} \\ &\text{Crank-Nicolson } & &\begin{cases} (I-\frac{h}{2}A)\mathbf{u}_{n+1} = (I+\frac{h}{2}A)\mathbf{u}_n + \frac{h}{2}\left( \mathbf{b}_n + \mathbf{b}_{n+1} \right)\\ \mathbf{u}_0 = \mathbf{y}_0 \end{cases} \end{aligned}\]

Il faut remarquer que à chaque étape des méthodes de ER et CN il faut résoudre un système linéaire de matrice \(I-hA\) et \(I-\frac{h}{2}A\) respectivement (il s’agit de méthodes implicites).

La méthode de EP est explicite (pas de système à résoudre) mais par contre elle est seulement conditionellement stable. Dans notre cas, les valeurs propres de \(A\) sont \(\lambda_1 = -1\), \(\lambda_2 = -5\); elles sont bien négatives, donc la condition ([stabilita]) sur \(h\) s’applique: comme \(\rho(A) = 5\), cette condition de stabilité est

\[h < \bar{h} = \frac{2}{5}.\]

Comportement du schéma de Euler progressif pour le système ([exemple-de-sys-lin-ode]) avec condition initiale \(\mathbf{y}_0 = [1,1\)^\top]: différentes valeur du pas de temps \(h\).

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

# Définir la matrice A et la fonction d'équation différentielle
A = np.array([[-2, 1], [3, -4]])

def dy(t, y, A):
    return np.dot(A, y) + np.array([np.exp(-t), 0])

# Calculer h_bar en fonction des valeurs propres de A
h_bar = 2 / np.max(np.abs(np.linalg.eigvals(A)))

# Résoudre l'EDO avec différents pas de temps
t_small, y_small_h = feuler(dy, [0, 10], np.array([1, 1]), int(10/(0.1*h_bar)), A)
t, y_h = feuler(dy, [0, 10], np.array([1, 1]), int(10/(h_bar)), A)
t_large, y_large_h = feuler(dy, [0, 10], np.array([1, 1]), int(10/(0.9*h_bar)), A)

# Créer le graphique avec Plotly
fig = go.Figure()

# Ajouter les trajectoires pour chaque pas de temps
fig.add_trace(go.Scatter(x=t_small, y=y_small_h[:,0], mode='lines', name='y_1 h = 0.1*h_bar'))
fig.add_trace(go.Scatter(x=t_small, y=y_small_h[:,1], mode='lines', name='y_2 h = 0.1*h_bar'))
fig.add_trace(go.Scatter(x=t, y=y_h[:,0], mode='lines+markers', name='y_1 h = h_bar', marker=dict(color='red')))
fig.add_trace(go.Scatter(x=t, y=y_h[:,1], mode='lines+markers', name='y_2 h = h_bar', marker=dict(color='blue')))
fig.add_trace(go.Scatter(x=t_large, y=y_large_h[:,0], mode='lines', name='y_1 h = 0.9*h_bar', line=dict(dash='dash', color='green')))
fig.add_trace(go.Scatter(x=t_large, y=y_large_h[:,1], mode='lines', name='y_2 h = 0.9*h_bar', line=dict(dash='dash', color='orange')))



# Mise en forme du graphique
fig.update_layout(title='Trajectoires des solutions avec différents pas de temps',
                  xaxis_title='y1',
                  yaxis_title='y2',
                  legend_title='Pas de temps')
fig.show()
Results

Tracons à présent \(y_1\) en fonction de \(y_2\) pour les différentes valeurs de \(h\).

# tracons y1 en fonction de y2

fig = go.Figure()
fig.add_trace(go.Scatter(x=y_small_h[:,0], y=y_small_h[:,1], mode='lines', name='h = 0.1*h_bar'))
fig.add_trace(go.Scatter(x=y_h[:,0], y=y_h[:,1], mode='lines+markers', name='h = h_bar', marker=dict(color='red')))
fig.add_trace(go.Scatter(x=y_large_h[:,0], y=y_large_h[:,1], mode='lines', name='h = 0.9*h_bar', line=dict(dash='dash', color='green')))
fig.update_layout(title='Trajectoires des solutions avec différents pas de temps',
                  xaxis_title='y1',
                  yaxis_title='y2',
                  legend_title='Pas de temps')
fig.show()
Results

3. Système non linéaire

On pourrait aussi considérer le cas d’un système non linéaire de la forme

\[\mathbf{y}'(t) = \mathbf{F} (t,\mathbf{y} (t))\]

(par exemple le système (#sys:ex-bio)). Si \(\frac{\partial \mathbf{F}}{\partial \mathbf{y}}\) est une matrice à valeurs propres réelles et négatives, alors la méthode d’Euler rétrograde est inconditionnellement stable, tandis que le schéma d’Euler progressif est stable sous la condition ([stabilita]), où maintenant \(A=\frac{\partial \mathbf{F}}{\partial \mathbf{y}}\).

Exemple: Système non-linéaire

Le système non-linéaire

\[\begin{aligned} y_1'(t) & = & -2 y_1(t) + \sin(y_2(t)) + e^{-t}\sin(t) \nonumber\\ y_2'(t) & = & \cos(y_1(t)) -4 y_2(t) \label{exemple-de-sys-nonlin-ode} \end{aligned}\]

avec les conditions initiales \(y_1(0)= y_{10}\), \(y_2(0)= y_{20}\), s’écrit sous la forme

\[\mathbf{y}'(t) = \mathbf{F} (t,\mathbf{y} (t)),\]

où

\[\mathbf{F} (t,\mathbf{y} (t)) = \begin{bmatrix} -2 y_1(t) + \sin(y_2(t)) + e^{-t}\sin(t) \\ \cos(y_1(t)) -4 y_2(t) \end{bmatrix}.\]

Soit \(h>0\) le pas de temps; pour \(n\in\mathbb{N}\) on pose \(t_n = nh\), et on désigne \(\mathbf{u}_n\) une valeur approchée de la solution exacte \(\mathbf{y}(t_n)\) au temps \(t_n\).

Les schémas d’Euler progressif, rétrograde et de Crank-Nicolson pour approcher la solution \(\mathbf{y}(t)\) de ([exemple-de-sys-nonlin-ode]) s’écrivent respectivement:

\[\begin{aligned} &\text{Euler progressif } & &\begin{cases} \mathbf{u}_{n+1} = \mathbf{u}_n + h\mathbf{F}(t_n,\mathbf{u}_n)\\ \mathbf{u}_0 = \mathbf{y}_0, \end{cases} \\ &\text{Euler r\'etrograde } & &\begin{cases} \mathbf{u}_{n+1} + h\mathbf{F}(t_{n+1},\mathbf{u}_{n+1}) = \mathbf{u}_n \\ \mathbf{u}_0 = \mathbf{y}_0 \end{cases} \\ &\text{Crank-Nicolson } & &\begin{cases} \mathbf{u}_{n+1} - \frac{h}{2} \mathbf{F}(t_{n+1},\mathbf{u}_{n+1}) = \mathbf{u}_n + \frac{h}{2}\mathbf{F}(t_n,\mathbf{u}_n) \\ \mathbf{u}_0 = \mathbf{y}_0 \end{cases} \end{aligned}\]

Il faut remarquer que à chaque étape des méthodes de ER et CN il faut résoudre un système non-linéaire.

La méthode de EP est explicite (pas de système à résoudre) mais par contre elle est seulement conditionellement stable. Dans notre cas, le jacobien de \(\mathbf{F}\) est donné par

\[J = \frac{\partial \mathbf{F}}{\partial \mathbf{y}} = \begin{bmatrix} -2 & \cos y_2 \\ -\sin y_1 & -4 \end{bmatrix}\]

et ses valeurs propres sont \(\lambda_{1,2} = -3 \pm \sqrt{1 - \sin y_1 \cos y_2};\) elles sont bien négatives, en particulier \(-3 - \sqrt{2} < \lambda_{1,2} < -3 + \sqrt{2} < 0\), et \(\rho(J)< 3 + \sqrt{2}\). Condition de stabilité:

\[h < \bar{h} = \frac{2}{\rho(J)}, \quad \text{satisfaite par example si} \quad h < \frac{2}{3+\sqrt{2}} \simeq 0.453.\]

Comportement du schéma de Euler progressif pour le système ([exemple-de-sys-nonlin-ode]) avec condition initiale \(\mathbf{y}_0 = [1,1\)^\top]: \(h=0.1\) (bleu) et \(h = 0.8\bar{h}\) (rouge). Si on prend \(h \geq \bar{h}\), on peut observer l’instabilité de la méthode.

image

On a utilisé les commandes suivantes:

import numpy as np
from numpy.linalg import eigvals
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from tan.ode.feuler import feuler

# Définir la fonction dy et la matrice Jacobienne J
def dy(t, y):
    return np.array([-2*y[0] + np.sin(y[1]) + np.exp(-t)*np.sin(t),
                     np.cos(y[0]) - 4*y[1]])

def J(y):
    return np.array([[-2, np.cos(y[1])],
                     [-np.sin(y[0]), -4]])


# Premier calcul avec h fixe
t, y = feuler(dy, [0, 100], [1, 1], int(100/0.1))

# Calculer rho pour chaque t
rho = [max(abs(eigvals(J(y[i,:])))) for i in range(len(t))]

# Calculer h_bar basé sur rho
h_bar = 2/max(rho)

# Deuxième calcul avec h basé sur h_bar
t_adjusted, y_adjusted = feuler(dy, [0, 100], [1, 1], int(100/(0.8*h_bar)))

# Visualisation avec Plotly
fig = make_subplots(rows=2, cols=1)

# Trajectoires des solutions
fig.add_trace(go.Scatter(x=y[:,0], y=y[:,1], mode='lines+markers', name='Original', marker=dict(color='blue')), row=1, col=1)
fig.add_trace(go.Scatter(x=y_adjusted[:,0], y=y_adjusted[:,1], mode='lines+markers', name='Adjusted', marker=dict(color='red')), row=1, col=1)

# Graphique de rho
fig.add_trace(go.Scatter(x=t, y=rho, mode='lines+markers', name='rho', marker=dict(color='green')), row=2, col=1)

# Mise à jour du layout
fig.update_layout(height=800, title='Système d\'EDO et Rho sur le temps', xaxis_title='t', yaxis_title='rho', showlegend=True)
fig.show()
Results