23. Résolution approchée d'une équation différentielle

1Probleme de Cauchy et methode d'Euler

Le probleme de Cauchy

Etant donne \(f : [t_0, T]\times\mathbb R \to \mathbb R\) et \(y_0 \in \mathbb R\), trouver \(y\) derivable verifiant \[ \begin{cases} y'(t) = f\bigl(t,\, y(t)\bigr) & \text{pour } t \in [t_0, T],\\ y(t_0) = y_0. \end{cases} \] La donnee \(y(t_0) = y_0\) est la condition initiale.

Theoreme de Cauchy-Lipschitz

Si \(f\) est continue sur \([t_0,T]\times\mathbb R\) et lipschitzienne en sa seconde variable, uniformement en \(t\) : \[ \exists L > 0,\ \forall t,\ \forall y_1, y_2 : \quad |f(t,y_1) - f(t,y_2)| \le L|y_1 - y_2|, \] alors le probleme de Cauchy admet une unique solution sur \([t_0, T]\).

C'est ce theoreme qui donne un sens au calcul approche : sans unicite, on ne saurait pas ce que l'algorithme est cense approcher.

Quelques equations resolubles exactement

Equation	Solution
\(y' = ay\), \(y(0)=y_0\)	\(y(t) = y_0e^{at}\)
\(y' = ay + b\)	\(y(t) = \left(y_0+\frac ba\right)e^{at} - \frac ba\)
\(y' = t\), \(y(0)=0\)	\(y(t) = t^2/2\)
\(y' = -2ty\), \(y(0)=1\)	\(y(t) = e^{-t^2}\)
\(y' = y^2\), \(y(0)=1\)	\(y(t) = \frac{1}{1-t}\) : explose en \(t=1\)

Le dernier exemple montre qu'une solution peut n'exister que localement, meme si \(f\) est tres reguliere : \(f(t,y) = y^2\) n'est pas globalement lipschitzienne.

Methode d'Euler explicite

Principe

On subdivise \([t_0, T]\) en \(n\) pas de longueur \(h = \frac{T-t_0}{n}\), avec \(t_k = t_0 + kh\). Sur \([t_k, t_{k+1}]\), on approche la courbe par sa tangente : \[ y(t_{k+1}) \approx y(t_k) + h\,y'(t_k) = y(t_k) + h\,f\bigl(t_k, y(t_k)\bigr). \] D'ou le schema \[ \boxed{\;y_{k+1} = y_k + h\,f(t_k, y_k),\qquad y_0 \text{ donne}.\;} \]

def euler(f, y0, t0, T, n):
    """Methode d'Euler explicite. Renvoie (liste des t, liste des y)."""
    h = (T - t0) / n
    t = [t0 + k*h for k in range(n + 1)]
    y = [y0]
    for k in range(n):
        y.append(y[k] + h * f(t[k], y[k]))
    return t, y

from math import exp
# y' = -2y, y(0) = 1  ->  solution exacte y(t) = exp(-2t)
t, y = euler(lambda t, y: -2*y, 1.0, 0, 2, 20)
print(f"y(2) approche = {y[-1]:.6f}, exact = {exp(-4):.6f}")

Erreurs locale et globale

Erreur locale de consistance : en supposant \(y_k = y(t_k)\) exact, la formule de Taylor donne \[ y(t_{k+1}) - y_{k+1} = y(t_k) + hy'(t_k) + \frac{h^2}{2}y''(\xi) - \bigl(y(t_k) + hy'(t_k)\bigr) = \frac{h^2}{2}y''(\xi) = O(h^2). \]

Erreur globale : les \(n = (T-t_0)/h\) erreurs locales s'accumulent, mais pas simplement — chacune est amplifiee par la dynamique. On montre \[ |y(t_k) - y_k| \;\le\; \frac{h M}{2L}\left(e^{L(t_k - t_0)} - 1\right) \;=\; O(h), \qquad M = \max|y''|. \] La methode d'Euler est donc d'ordre 1 : diviser \(h\) par 2 divise l'erreur finale par 2 seulement.

Idee de la demonstration

Posons \(e_k = |y(t_k) - y_k|\). En comparant les deux relations de recurrence et en utilisant le caractere lipschitzien de \(f\) : \[ e_{k+1} \le e_k + hL e_k + \frac{h^2M}{2} = (1+hL)e_k + \frac{h^2M}{2}. \] Cette recurrence arithmetico-geometrique se resout en \(e_k \le \frac{hM}{2L}\bigl((1+hL)^k - 1\bigr)\), et l'on conclut avec \((1+hL)^k \le e^{khL} = e^{L(t_k-t_0)}\).

Consequence importante : le facteur \(e^{L(T-t_0)}\) croit exponentiellement avec la duree d'integration. Sur un intervalle long, ou pour une equation raide (\(L\) grand), l'erreur devient enorme meme avec un pas fin.

Autres schemas classiques

def euler_implicite(f, y0, t0, T, n, tol=1e-12, itmax=50):
    """Euler retrograde : y[k+1] = y[k] + h f(t[k+1], y[k+1]).

    Schema IMPLICITE : on resout a chaque pas par point fixe.
    """
    h = (T - t0)/n
    t = [t0 + k*h for k in range(n+1)]
    y = [y0]
    for k in range(n):
        z = y[k] + h*f(t[k], y[k])            # prediction par Euler explicite
        for _ in range(itmax):
            z_nouveau = y[k] + h*f(t[k+1], z)
            if abs(z_nouveau - z) < tol:
                z = z_nouveau; break
            z = z_nouveau
        y.append(z)
    return t, y

def point_milieu(f, y0, t0, T, n):
    """Runge-Kutta d'ordre 2 (methode du point milieu)."""
    h = (T - t0)/n
    t = [t0 + k*h for k in range(n+1)]
    y = [y0]
    for k in range(n):
        k1 = f(t[k], y[k])
        k2 = f(t[k] + h/2, y[k] + h/2*k1)
        y.append(y[k] + h*k2)
    return t, y

def heun(f, y0, t0, T, n):
    """Runge-Kutta d'ordre 2 (methode de Heun / trapezes explicite)."""
    h = (T - t0)/n
    t = [t0 + k*h for k in range(n+1)]
    y = [y0]
    for k in range(n):
        k1 = f(t[k], y[k])
        k2 = f(t[k+1], y[k] + h*k1)
        y.append(y[k] + h/2*(k1 + k2))
    return t, y

def rk4(f, y0, t0, T, n):
    """Runge-Kutta d'ordre 4 : la methode de reference."""
    h = (T - t0)/n
    t = [t0 + k*h for k in range(n+1)]
    y = [y0]
    for k in range(n):
        k1 = f(t[k],       y[k])
        k2 = f(t[k] + h/2, y[k] + h/2*k1)
        k3 = f(t[k] + h/2, y[k] + h/2*k2)
        k4 = f(t[k] + h,   y[k] + h*k3)
        y.append(y[k] + h/6*(k1 + 2*k2 + 2*k3 + k4))
    return t, y

Schema	Ordre	Evaluations de \(f\) par pas	Type
Euler explicite	1	1	explicite
Euler implicite	1	plusieurs (resolution)	implicite
Point milieu, Heun	2	2	explicite
RK4	4	4	explicite

RK4 offre le meilleur compromis : ordre 4 pour 4 evaluations. C'est la methode par defaut de la plupart des solveurs.

2Visualisation, schemas d'ordre superieur et stabilite

Methode d'Euler : la solution approchee suit les tangentes

n = 8 champ des tangentes

A observer

Sur \(y' = -2y\), Euler sous-estime systematiquement (la solution est convexe, les tangentes passent en dessous).
Sur \(y' = y\), Euler sous-estime aussi, et l'ecart croit exponentiellement : c'est le facteur \(e^{L(T-t_0)}\) de la majoration.
Avec \(n = 8\), RK4 est deja visuellement confondu avec la solution exacte, la ou Euler s'ecarte nettement.
Le champ des tangentes montre que toutes les solutions suivent les memes directions : la solution est la courbe « portee par le courant » depuis la condition initiale.

Stabilite

Probleme modele

On etudie \(y' = \lambda y\) avec \(\lambda < 0\), dont la solution \(y(t) = y_0e^{\lambda t}\) tend vers 0. Un bon schema doit reproduire cette decroissance.

Euler explicite donne \(y_{k+1} = (1 + h\lambda)y_k\), donc \(y_k = (1+h\lambda)^k y_0\). Il y a decroissance si et seulement si \[ |1 + h\lambda| < 1 \iff \boxed{\;0 < h < \frac{2}{|\lambda|}.\;} \]

def euler_simple(lam, y0, T, n):
    h = T/n
    y = y0
    for _ in range(n):
        y = y + h*lam*y
    return y

lam, T = -50, 1
print(f"{'n':>6}{'h':>10}{'|1+h*lam|':>12}{'y(1)':>16}{'exact':>12}")
for n in (20, 25, 40, 100):
    h = T/n
    print(f"{n:>6}{h:>10.4f}{abs(1+h*lam):>12.3f}"
          f"{euler_simple(lam, 1, T, n):>16.3e}{exp(-50):>12.2e}")

     n         h   |1+h*lam|            y(1)       exact
    20    0.0500       1.500       3.325e+03    1.93e-22
    25    0.0400       1.000       1.000e+00    1.93e-22
    40    0.0250       0.250       8.882e-25    1.93e-22
   100    0.0100       0.500       7.889e-31    1.93e-22

Pour \(\lambda = -50\), la condition de stabilite est \(h < 0{,}04\), soit \(n > 25\). Avec \(n = 20\), la solution numerique explose alors que la vraie solution tend vers 0 : c'est une instabilite numerique, sans rapport avec le probleme physique.

Equations raides

Une equation est raide (stiff) lorsque la condition de stabilite impose un pas beaucoup plus petit que celui qu'exigerait la seule precision. C'est le cas des systemes ou coexistent des echelles de temps tres differentes (cinetique chimique, circuits electroniques, combustion).

Euler implicite donne \(y_{k+1} = \frac{y_k}{1 - h\lambda}\), et \(\left|\frac{1}{1-h\lambda}\right| < 1\) pour tout \(h > 0\) quand \(\lambda < 0\) : il est inconditionnellement stable. C'est pourquoi on utilise des schemas implicites sur les problemes raides, malgre le cout de la resolution a chaque pas.

def euler_implicite_lineaire(lam, y0, T, n):
    """Pour y' = lam*y, la resolution est explicite : y = y/(1 - h*lam)."""
    h = T/n
    y = y0
    for _ in range(n):
        y = y / (1 - h*lam)
    return y

for n in (5, 10, 20):
    print(f"n = {n:>3} : implicite {euler_implicite_lineaire(-50, 1, 1, n):.3e}, "
          f"explicite {euler_simple(-50, 1, 1, n):.3e}")

n =   5 : implicite 9.765e-07, explicite -5.905e+04
n =  10 : implicite 1.024e-04, explicite 2.594e+05
n =  20 : implicite 4.437e-04, explicite 3.325e+03

L'implicite reste borne et decroissant meme avec 5 pas ; l'explicite oscille et diverge. L'implicite n'est pas precis pour autant (il reste d'ordre 1), mais il est qualitativement correct.

3Systemes differentiels et applications

Systemes differentiels et equations d'ordre superieur

Reduction a l'ordre 1

Toute equation d'ordre \(p\) se ramene a un systeme de \(p\) equations d'ordre 1. Pour \(y'' = g(t, y, y')\), on pose \(u = y\) et \(v = y'\) : \[ \begin{cases} u' = v,\\ v' = g(t, u, v),\end{cases} \qquad \text{soit } Y' = F(t, Y) \text{ avec } Y = \begin{pmatrix}u\\v\end{pmatrix}. \] Les schemas s'appliquent alors composante par composante, sans aucune modification.

import numpy as np

def rk4_systeme(F, Y0, t0, T, n):
    """RK4 pour un systeme. Y0 est un tableau NumPy. Renvoie (t, Y) avec Y de forme (n+1, d)."""
    h = (T - t0)/n
    t = np.linspace(t0, T, n+1)
    Y = np.zeros((n+1, len(Y0)))
    Y[0] = Y0
    for k in range(n):
        k1 = F(t[k],       Y[k])
        k2 = F(t[k] + h/2, Y[k] + h/2*k1)
        k3 = F(t[k] + h/2, Y[k] + h/2*k2)
        k4 = F(t[k] + h,   Y[k] + h*k3)
        Y[k+1] = Y[k] + h/6*(k1 + 2*k2 + 2*k3 + k4)
    return t, Y

# pendule simple : theta'' = -(g/L) sin(theta)
g, L = 9.81, 1.0
def pendule(t, Y):
    theta, omega = Y
    return np.array([omega, -(g/L)*np.sin(theta)])

t, Y = rk4_systeme(pendule, np.array([np.pi/3, 0.0]), 0, 10, 2000)

# verification : l'energie mecanique doit etre conservee
E = 0.5*L*L*Y[:,1]**2 + g*L*(1 - np.cos(Y[:,0]))
print(f"energie : min {E.min():.10f}, max {E.max():.10f}, "
      f"derive relative {(E.max()-E.min())/E[0]:.2e}")

energie : min 4.9050000000, max 4.9050000004, derive relative 8.2e-11

La conservation comme test

Verifier qu'une quantite conservee par le systeme physique (energie, moment cinetique, masse) reste constante dans la simulation est le meilleur test de validation d'un integrateur : il ne demande pas de connaitre la solution exacte, et detecte immediatement une erreur de schema.

Attention : RK4 conserve l'energie approximativement mais avec une derive lente seculaire. Pour des simulations sur des temps tres longs (mecanique celeste), on utilise des integrateurs symplectiques (Verlet, leapfrog) qui, bien que d'ordre 2 seulement, conservent exactement une energie modifiee et donc ne derivent pas.

def verlet(a, x0, v0, T, n):
    """Integrateur symplectique de Stormer-Verlet pour x'' = a(x). Ordre 2."""
    h = T/n
    x, v = np.zeros(n+1), np.zeros(n+1)
    x[0], v[0] = x0, v0
    acc = a(x0)
    for k in range(n):
        v_demi = v[k] + h/2*acc
        x[k+1] = x[k] + h*v_demi
        acc = a(x[k+1])
        v[k+1] = v_demi + h/2*acc
    return x, v

# oscillateur harmonique x'' = -x, energie = (x^2 + v^2)/2 = constante
x, v = verlet(lambda x: -x, 1.0, 0.0, 200, 4000)
E = 0.5*(x*x + v*v)
print(f"Verlet sur 200 periodes : derive d'energie {(E.max()-E.min())/E[0]:.2e}")

Systemes celebres

Proie-predateur (Lotka-Volterra)

\[ \begin{cases} x' = \alpha x - \beta xy & \text{(proies)},\\ y' = \delta xy - \gamma y & \text{(predateurs)}.\end{cases} \]

alpha, beta, delta, gamma = 1.0, 0.1, 0.075, 1.5

def lotka(t, Y):
    x, y = Y
    return np.array([alpha*x - beta*x*y, delta*x*y - gamma*y])

t, Y = rk4_systeme(lotka, np.array([10.0, 5.0]), 0, 40, 8000)

# invariant du systeme : V = delta*x - gamma*ln(x) + beta*y - alpha*ln(y)
V = delta*Y[:,0] - gamma*np.log(Y[:,0]) + beta*Y[:,1] - alpha*np.log(Y[:,1])
print(f"invariant : derive relative {(V.max()-V.min())/abs(V[0]):.2e}")
print(f"periode approximative : {40/ (np.diff(np.where(np.diff(np.sign(Y[:,0] - Y[:,0].mean())) > 0)[0]).size or 1):.2f}")

Les trajectoires sont fermees dans le plan de phase \((x,y)\) : les populations oscillent periodiquement. Le systeme admet un invariant \(V\), ce qui fournit un test de validation.

Attracteur de Lorenz

\[ x' = \sigma(y-x), \qquad y' = x(\rho - z) - y, \qquad z' = xy - \beta z. \]

sigma, rho, beta_l = 10.0, 28.0, 8/3

def lorenz(t, Y):
    x, y, z = Y
    return np.array([sigma*(y-x), x*(rho-z) - y, x*y - beta_l*z])

t, A = rk4_systeme(lorenz, np.array([1.0, 1.0, 1.0]), 0, 40, 40000)
t, B = rk4_systeme(lorenz, np.array([1.0 + 1e-9, 1.0, 1.0]), 0, 40, 40000)

ecart = np.linalg.norm(A - B, axis=1)
for k in (0, 5000, 10000, 20000, 39999):
    print(f"t = {t[k]:>5.1f} : ecart {ecart[k]:.3e}")

t =   0.0 : ecart 1.000e-09
t =   5.0 : ecart 1.161e-06
t =  10.0 : ecart 1.383e-03
t =  20.0 : ecart 1.878e+01
t =  40.0 : ecart 2.723e+01

Sensibilite aux conditions initiales. Une difference initiale de \(10^{-9}\) devient de l'ordre de la taille de l'attracteur (\(\approx 30\)) en \(t = 20\). L'ecart croit comme \(e^{\lambda t}\) avec \(\lambda \approx 0{,}9\) (exposant de Lyapunov) : chaque unite de temps multiplie l'erreur par 2,5 environ.

C'est l'effet papillon, decouvert par Lorenz en 1963 en modelisant la convection atmospherique. Consequence : au-dela d'un horizon de predictibilite, aucune amelioration du schema numerique ni de la precision machine ne permet de prevoir la trajectoire — seules les proprietes statistiques restent accessibles. C'est la raison fondamentale pour laquelle la meteorologie ne depasse pas une dizaine de jours de prevision fiable.

L'essentiel du chapitre

Schema d'Euler : \(y_{k+1} = y_k + h\,f(t_k, y_k)\), erreur locale \(O(h^2)\), erreur globale \(O(h)\) : ordre 1.
L'erreur globale contient un facteur \(e^{L(T-t_0)}\) : integrer longtemps coute cher en precision.
Heun et point milieu sont d'ordre 2 ; RK4 est d'ordre 4 pour 4 evaluations : c'est le bon compromis.
Euler explicite n'est stable que si \(h < 2/|\lambda|\) ; Euler implicite est inconditionnellement stable.
Toute equation d'ordre \(p\) se ramene a un systeme d'ordre 1 par \(u = y\), \(v = y'\), etc.
Verifier une quantite conservee est le meilleur test de validation.
Sur un systeme chaotique, la sensibilite aux conditions initiales borne l'horizon de predictibilite.