Résolution d'équations différentielles

Cours · Python (école d'ingénieur), chapitre 13 · école d'ingénieur, informatique · le début du chapitre, la suite sur Adloun

L'évolution temporelle de nombreux systèmes physiques, biologiques, ou chimiques est régie par des équations différentielles. S'il est parfois possible de trouver des solutions analytiques (comme la décharge simple d'un condensateur dans un circuit $R C$ ), la grande majorité des modèles réels (tels que le mouvement gravitationnel à trois corps ou le modèle épidémiologique SIR) ne possède aucune solution exacte sous forme de fonctions usuelles. Nous devons alors recourir aux schémas de résolution numérique. Ce chapitre présente les méthodes d'intégration à un pas pour les problèmes de Cauchy (méthodes d'Euler explicite et implicite, méthodes de Runge-Kutta d'ordres 2 et 4) et montre comment modéliser et simuler des systèmes d'équations différentielles d'ordres supérieurs en Python.

13.1 Problème de Cauchy et Champ de Vecteurs

On s'intéresse à la résolution du problème de Cauchy suivant sur un intervalle $[t_{0}, t_{0} + T]$ :

${x^{'} (t) = f (x (t), t) x (t_{0}) = x_{0}$

où $f$ est une fonction continue de deux variables et $x_{0}$ est la condition initiale.

13.1.1 Visualisation : Champ de Vecteurs (Direction Field)

Graphiquement, l'équation $x^{'} (t) = f (x, t)$ attribue à chaque point $(t, x)$ un vecteur de pente $f (x, t)$ . Une solution $x (t)$ est une courbe tangente en tout point à ce champ de vecteurs.

Exemple 13.1Champ de vecteurs avec quiver


import numpy as np
import matplotlib.pyplot as plt

def champ_vecteurs(f, t_range, x_range, norm=True):
    """Affiche le champ de vecteurs de l'équation x' = f(x, t)."""
    t = np.linspace(t_range[0], t_range[1], 20)
    x = np.linspace(x_range[0], x_range[1], 20)
    T, X = np.meshgrid(t, x)

    # Composantes du champ (dt=1, dx=f(x,t))
    U = np.ones_like(T)
    V = f(X, T)

    if norm:
        N = np.sqrt(U**2 + V**2)
        U, V = U/N, V/N

    plt.quiver(T, X, U, V, color="gray", alpha=0.7)
    plt.xlabel("t")
    plt.ylabel("x")
    plt.title("Champ de vecteurs")
    plt.grid(True, alpha=0.3)

# Exemple : x' = -x/tau (décroissance exponentielle)
tau = 2.0
f = lambda x, t: -x / tau

champ_vecteurs(f, (0, 10), (-2, 2))
plt.show()

13.2 Méthodes d'Euler

Pour approcher les valeurs de $x (t)$ sur une subdivision régulière de pas $h = \frac{T}{n}$ définie par $t_{k} = t_{0} + k h$ , nous pouvons intégrer l'équation entre $t_{k}$ et $t_{k + 1}$ :

$x (t_{k + 1}) - x (t_{k}) = \int_{t_{k}}^{t_{k + 1}} f (x (t), t) d t$

13.2.1 Méthode d'Euler Explicite

En approchant l'intégrale par le rectangle à gauche (valeur en $t_{k}$ ), on obtient le schéma explicite :

$x_{k + 1} = x_{k} + h f (x_{k}, t_{k})$

◆Théorème 13.2Consistance d'Euler

Lorsque la fonction $f$ est de classe $C^{1}$ , le schéma d'Euler explicite est consistant d'ordre 1.


def euler_explicite(f, x0, t0, tf, n):
    """Méthode d'Euler explicite."""
    h = (tf - t0) / n
    t = np.linspace(t0, tf, n+1)
    x = np.zeros(n+1)
    x[0] = x0

    for k in range(n):
        x[k+1] = x[k] + h * f(x[k], t[k])

    return t, x

13.2.2 Méthode d'Euler Implicite

En utilisant le rectangle à droite (valeur en $t_{k + 1}$ ), on obtient le schéma implicite :

$x_{k + 1} = x_{k} + h f (x_{k + 1}, t_{k + 1})$

Pour une équation $x^{'} = f (x, t)$ non linéaire, on doit résoudre une équation à chaque pas. Pour l'équation linéaire $x^{'} = - λ x$ , la solution implicite s'écrit :

$x_{k + 1} = \frac{x _{k}}{1 + λh}$


def euler_implicite_linear(f, x0, t0, tf, n, lambd):
    """Euler implicite pour x' = -lambda x."""
    h = (tf - t0) / n
    t = np.linspace(t0, tf, n+1)
    x = np.zeros(n+1)
    x[0] = x0

    for k in range(n):
        x[k+1] = x[k] / (1 + lambd * h)

    return t, x

13.2.3 Stabilité des méthodes d'Euler

Exemple 13.3Comparaison de la stabilité des méthodes d'Euler

Pour l'équation test $x^{'} = - λ x$ ( $λ > 0$ ), la solution exacte décroît exponentiellement vers 0.

Euler explicite : $x_{k + 1} = (1 - λh) x_{k}$ . Stabilité si $∣1 - λh ∣ < 1$ , soit $h < \frac{2}{λ}$ .
Euler implicite : $x_{k + 1} = \frac{x _{k}}{1 + λh}$ . Toujours stable ( $∣1/ (1 + λh) ∣ < 1$ pour tout $h > 0$ ).


def comparer_euler(lambd=5.0, t_max=2.0):
    """Compare Euler explicite et implicite pour x' = -lambda x."""
    x0 = 1.0
    t_exact = np.linspace(0, t_max, 200)
    x_exact = x0 * np.exp(-lambd * t_exact)

    plt.figure(figsize=(10, 6))
    plt.plot(t_exact, x_exact, 'k-', linewidth=2, label='Solution exacte')

    for h in [0.3, 0.5, 0.7]:
        n = int(t_max / h)

        # Euler explicite
        t_exp, x_exp = euler_explicite(lambda x, t: -lambd*x, x0, 0, t_max, n)
        plt.plot(t_exp, x_exp, '--', label=f'Euler expo (h={h})')

        # Euler implicite
        t_imp, x_imp = euler_implicite_linear(lambda x, t: -lambd*x, x0, 0, t_max, n, lambd)
        plt.plot(t_imp, x_imp, ':', label=f'Euler impli (h={h})')

    plt.xlabel('t')
    plt.ylabel('x(t)')
    plt.title('Comparaison Euler explicite vs implicite')
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.show()

13.3 Méthodes de Runge-Kutta

Les méthodes de Runge-Kutta (RK) évaluent la fonction $f$ en plusieurs points intermédiaires pour approcher la pente moyenne avec un ordre de précision supérieur.

13.3.1 Méthode du Point Milieu (RK2 - ordre 2)

$⎩ ⎨ ⎧ k_{1} = f (x_{k}, t_{k}) k_{2} = f (x_{k} + \frac{h}{2} k_{1}, t_{k} + \frac{h}{2}) x_{k + 1} = x_{k} + h k_{2}$

13.3.2 Méthode de Heun (RK2 alternative)

$⎩ ⎨ ⎧ k_{1} = f (x_{k}, t_{k}) k_{2} = f (x_{k} + h k_{1}, t_{k} + h) x_{k + 1} = x_{k} + \frac{h}{2} (k_{1} + k_{2})$

13.3.3 Méthode RK4 classique (ordre 4)

$⎩ ⎨ ⎧ k_{1} = f (x_{k}, t_{k}) k_{2} = f (x_{k} + \frac{h}{2} k_{1}, t_{k} + \frac{h}{2}) k_{3} = f (x_{k} + \frac{h}{2} k_{2}, t_{k} + \frac{h}{2}) k_{4} = f (x_{k} + h k_{3}, t_{k} + h) x_{k + 1} = x_{k} + \frac{h}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4})$


def rk4(f, x0, t0, tf, n):
    """Méthode de Runge-Kutta d'ordre 4."""
    h = (tf - t0) / n
    t = np.linspace(t0, tf, n+1)
    x = np.zeros(n+1)
    x[0] = x0

    for k in range(n):
        k1 = f(x[k], t[k])
        k2 = f(x[k] + 0.5 * h * k1, t[k] + 0.5 * h)
        k3 = f(x[k] + 0.5 * h * k2, t[k] + 0.5 * h)
        k4 = f(x[k] + h * k3, t[k] + h)
        x[k+1] = x[k] + (h / 6.0) * (k1 + 2.0*k2 + 2.0*k3 + k4)

    return t, x

13.4 Systèmes Différentiels et Équations d'Ordre Supérieur

Pour simuler des systèmes dynamiques complexes ou des équations différentielles d'ordre $p \geq 2$ , on les récrit sous la forme d'un système vectoriel de premier ordre.

13.4.1 Principe de Réduction

Soit l'équation du second ordre $y ¨ = g (y, y ˙, t)$ . En posant le vecteur d'état $X (t) = (y (t) y ˙ (t))$ , sa dérivée s'écrit :

$X^{'} (t) = (y ˙ (t) y ¨ (t)) = (X_{1} (t) g (X_{0} (t), X_{1} (t), t)) = F (X, t)$

Exemple 13.4Oscillateur harmonique amorti

L'équation $y ¨ + γ y ˙ + ω_{0}^{2} y = 0$ se récrit en système :

$X^{'} = (v - γ v - ω_{0}^{2} y)$


def oscillateur(X, t, gamma=0.5, omega0=2.0):
    """Système pour oscillateur harmonique amorti."""
    y, v = X
    return np.array([v, -gamma * v - omega0**2 * y])

def rk4_systeme(F, X0, t0, tf, n, args=()):
    """RK4 pour système vectoriel."""
    h = (tf - t0) / n
    t = np.linspace(t0, tf, n+1)
    X = np.zeros((n+1, len(X0)))
    X[0] = X0

    for k in range(n):
        k1 = F(X[k], t[k], *args)
        k2 = F(X[k] + 0.5 * h * k1, t[k] + 0.5 * h, *args)
        k3 = F(X[k] + 0.5 * h * k2, t[k] + 0.5 * h, *args)
        k4 = F(X[k] + h * k3, t[k] + h, *args)
        X[k+1] = X[k] + (h / 6.0) * (k1 + 2.0*k2 + 2.0*k3 + k4)

    return t, X

# Simulation
t, X = rk4_systeme(oscillateur, [1.0, 0.0], 0, 20, 1000, args=(0.5, 2.0))

plt.figure(figsize=(10, 6))
plt.plot(t, X[:, 0], 'b-', label='Position y(t)')
plt.plot(t, X[:, 1], 'r--', label='Vitesse v(t)')
plt.xlabel('t')
plt.ylabel('y, v')
plt.title('Oscillateur harmonique amorti (RK4)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

13.5 Modèle SIR (Épidémiologie)

Le modèle SIR (Susceptible-Infectious-Removed) est un classique de l'épidémiologie :

$⎩ ⎨ ⎧ \frac{d S}{d t} = - β S I \frac{d I}{d t} = β S I - γ I \frac{d R}{d t} = γ I$

où $S$ est la proportion de susceptibles, $I$ d'infectieux, $R$ de guéris.


def sir_model(X, t, beta, gamma):
    """Modèle SIR."""
    S, I, R = X
    dS = -beta * S * I
    dI = beta * S * I - gamma * I
    dR = gamma * I
    return np.array([dS, dI, dR])

# Paramètres
beta, gamma = 0.3, 0.1
S0, I0, R0 = 0.99, 0.01, 0.0
t, X = rk4_systeme(sir_model, [S0, I0, R0], 0, 100, 1000, args=(beta, gamma))

plt.figure(figsize=(10, 6))
plt.plot(t, X[:, 0], 'b-', label='Susceptibles S(t)')
plt.plot(t, X[:, 1], 'r-', label='Infectieux I(t)')
plt.plot(t, X[:, 2], 'g-', label='Guéris R(t)')
plt.xlabel('t')
plt.ylabel('Population')
plt.title('Modèle SIR')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

13.6 Utilisation de Solveurs Professionnels

La bibliothèque scipy.integrate fournit la fonction solve_ivp (Initial Value Problem) avec plusieurs méthodes intégrées.


from scipy.integrate import solve_ivp

def oscillateur_scipy(t, X, gamma, omega0):
    y, v = X
    return [v, -gamma * v - omega0**2 * y]

# Résolution avec solve_ivp
solution = solve_ivp(oscillateur_scipy, [0, 20], [1.0, 0.0],
                     args=(0.5, 2.0), method='RK45', dense_output=True)

t_scipy = np.linspace(0, 20, 500)
X_scipy = solution.sol(t_scipy)

plt.figure(figsize=(10, 6))
plt.plot(t_scipy, X_scipy[0], 'b-', label='solve_ivp (RK45)')
plt.xlabel('t')
plt.ylabel('y(t)')
plt.title('Oscillateur avec solve_ivp')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

13.7 Exercices résolus

Niveau 1 : Méthode d'Euler explicite

Exercice 1 - Décroissance exponentielle

Soit l'équation différentielle $x^{'} (t) = - 2 x (t)$ avec $x (0) = 1$ .

Donner la solution analytique.
Implémenter la méthode d'Euler explicite avec $h = 0.1$ sur $t \in [0, 3]$ .
Comparer avec la solution exacte.

Démonstration (Solution)

La solution analytique est $x (t) = e^{- 2 t}$ .
Implémentation et comparaison en Python :


import numpy as np
import matplotlib.pyplot as plt

# Solution analytique
def x_exact(t):
    return np.exp(-2 * t)

# Euler explicite
def euler_explicite(f, x0, t0, tf, n):
    h = (tf - t0) / n
    t = np.linspace(t0, tf, n+1)
    x = np.zeros(n+1)
    x[0] = x0

    for k in range(n):
        x[k+1] = x[k] + h * f(x[k], t[k])

    return t, x

# Paramètres
f = lambda x, t: -2 * x
x0 = 1.0
t0, tf = 0.0, 3.0
n = 30 # h = 0.1

t_euler, x_euler = euler_explicite(f, x0, t0, tf, n)
t_exact = np.linspace(t0, tf, 300)
x_exact_vals = x_exact(t_exact)

# Erreur
erreur = np.abs(x_euler[-1] - x_exact(t_euler[-1]))
print(f"Valeur exacte à t=3 : {x_exact(3):.6f}")
print(f"Euler explicite à t=3 : {x_euler[-1]:.6f}")
print(f"Erreur absolue : {erreur:.6f}")

# Visualisation
plt.figure(figsize=(10, 6))
plt.plot(t_exact, x_exact_vals, 'k-', linewidth=2, label='Solution exacte')
plt.plot(t_euler, x_euler, 'ro-', markersize=4, label='Euler explicite (h=0.1)')
plt.xlabel('t')
plt.ylabel('x(t)')
plt.title('Décroissance exponentielle')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

Figure : Décroissance exponentielle résolue par la méthode d'Euler explicite (

h = 0.1

Résultat dans la console :


Valeur exacte à t=3 : 0.002479
Euler explicite à t=3 : 0.002478
Erreur absolue : 0.000001

L'erreur diminue avec le pas $h$ (théoriquement en $O (h)$ ).

Exercice 2 - Étude de l'erreur d'Euler

Pour l'équation $x^{'} = - 2 x$ , étudier l'évolution de l'erreur en fonction du pas $h$ .

La suite du cours est sur Adloun

Le chapitre en entier, formules et figures comprises, avec le tuteur, les fiches, les QCM et les exercices.

Lire la suite du cours

En accès libre, sans inscription.

Exercices corrigés de ce chapitre

Lire la suite du cours →