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

1. Euler avec \(h = 0{,}25\). Le schema donne \(y_{k+1} = y_k + h(-2y_k) = (1 - 0{,}5)y_k = 0{,}5\,y_k\).

2. L'approximation est par defaut (\(y_k < y(t_k)\)). La solution \(t \mapsto e^{-2t}\) est convexe (\(y'' = 4e^{-2t} > 0\)), donc elle est au-dessus de chacune de ses tangentes. Comme Euler suit precisement les tangentes, il sous-estime systematiquement. Le terme d'erreur locale \(\frac{h^2}{2}y''(\xi) > 0\) confirme ce signe.

3. Heun avec \(h = 0{,}25\). Ici \(k_1 = -2y_k\), \(k_2 = -2(y_k + hk_1) = -2y_k(1-2h) = -2y_k\cdot0{,}5\), d'ou \[ y_{k+1} = y_k + \frac h2(k_1 + k_2) = y_k\bigl(1 - 2h + 2h^2\bigr) = 0{,}625\,y_k. \]

L'erreur est divisee par environ 4 : conforme au passage de l'ordre 1 a l'ordre 2. Noter que Heun surestime ici (le facteur \(0{,}625\) est superieur a \(e^{-0{,}5} = 0{,}6065\)).

4. Avec \(h = 1{,}5\) : \(1 + h\lambda = 1 - 3 = -2\), donc \(|1+h\lambda| = 2 > 1\). La suite \(y_k = (-2)^k\) oscille en divergeant : \(1, -2, 4, -8, \dots\), alors que la solution exacte tend vers 0. La condition de stabilite \(h < 2/|\lambda| = 1\) n'est pas respectee. C'est une instabilite purement numerique.

import numpy as np
from math import exp

def euler(f, y0, t0, T, n):
    h = (T-t0)/n; t = [t0]; y = [y0]
    for k in range(n):
        y.append(y[k] + h*f(t[k], y[k])); t.append(t[k]+h)
    return t, y

def heun(f, y0, t0, T, n):
    h = (T-t0)/n; t = [t0]; y = [y0]
    for k in range(n):
        a = f(t[k], y[k]); b = f(t[k]+h, y[k]+h*a)
        y.append(y[k] + h/2*(a+b)); t.append(t[k]+h)
    return t, y

def point_milieu(f, y0, t0, T, n):
    h = (T-t0)/n; t = [t0]; y = [y0]
    for k in range(n):
        a = f(t[k], y[k]); b = f(t[k]+h/2, y[k]+h/2*a)
        y.append(y[k] + h*b); t.append(t[k]+h)
    return t, y

def rk4(f, y0, t0, T, n):
    h = (T-t0)/n; t = [t0]; 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)); t.append(t[k]+h)
    return t, y

f = lambda t, y: -2*y
exact = exp(-4)
methodes = [("Euler", euler, 1), ("Heun", heun, 2),
            ("milieu", point_milieu, 2), ("RK4", rk4, 4)]

print(f"{'n':>6}" + "".join(f"{nom:>13}{'r':>6}" for nom, _, _ in methodes))
prec = None
for n in (10, 20, 40, 80, 160, 320):
    E = [abs(m(f, 1.0, 0, 2, n)[1][-1] - exact) for _, m, _ in methodes]
    if prec is None:
        print(f"{n:>6}" + "".join(f"{v:>13.2e}{'':>6}" for v in E))
    else:
        s = f"{n:>6}"
        for v, p in zip(E, prec):
            s += f"{v:>13.2e}{(p/v if v > 1e-17 else float('inf')):>6.1f}"
        print(s)
    prec = E

     n        Euler     r        Heun     r      milieu     r         RK4     r
    10     6.44e-03          8.42e-04          1.05e-03          1.35e-06
    20     3.03e-03   2.1   2.16e-04   3.9   2.63e-04   4.0   8.80e-08  15.3
    40     1.47e-03   2.1   5.46e-05   4.0   6.61e-05   4.0   5.61e-09  15.7
    80     7.24e-04   2.0   1.37e-05   4.0   1.66e-05   4.0   3.54e-10  15.9
   160     3.60e-04   2.0   3.44e-06   4.0   4.15e-06   4.0   2.22e-11  15.9
   320     1.79e-04   2.0   8.60e-07   4.0   1.04e-06   4.0   1.39e-12  16.0

3. Les rapports valent 2, 4, 4 et 16, soit des ordres 1, 2, 2 et 4 : conforme.

# 4. comparaison a NOMBRE D'EVALUATIONS egal
print(f"\n{'evaluations':>13}{'Euler':>13}{'Heun':>13}{'RK4':>13}")
for ev in (40, 80, 160, 320, 640):
    e1 = abs(euler(f, 1.0, 0, 2, ev)[1][-1] - exact)
    e2 = abs(heun(f, 1.0, 0, 2, ev//2)[1][-1] - exact)
    e4 = abs(rk4(f, 1.0, 0, 2, ev//4)[1][-1] - exact)
    print(f"{ev:>13}{e1:>13.2e}{e2:>13.2e}{e4:>13.2e}")

  evaluations        Euler         Heun          RK4
           40     1.47e-03     8.42e-04     1.35e-06
           80     7.24e-04     2.16e-04     8.80e-08
          160     3.60e-04     5.46e-05     5.61e-09
          320     1.79e-04     1.37e-05     3.54e-10
          640     8.94e-05     3.44e-06     2.22e-11

4. A cout egal, RK4 est superieur de trois a quatre ordres de grandeur. Le surcout par pas (4 evaluations au lieu d'une) est largement compense par le gain d'ordre. C'est pourquoi RK4 est le choix par defaut : on ne lui prefere Euler que pour sa simplicite pedagogique, ou un schema implicite pour les problemes raides.

1. Euler explicite : \(y_{k+1} = (1+h\lambda)y_k\), donc \(y_k = (1+h\lambda)^k y_0\). La suite tend vers 0 si et seulement si \(|1+h\lambda| < 1\), soit, pour \(\lambda < 0\) reel : \[ -1 < 1 + h\lambda < 1 \iff 0 < h < \frac{2}{|\lambda|}. \]

from math import exp

def euler_lin(lam, T, n):
    h, y = T/n, 1.0
    for _ in range(n):
        y = y*(1 + h*lam)
    return y

def implicite_lin(lam, T, n):
    h, y = T/n, 1.0
    for _ in range(n):
        y = y/(1 - h*lam)
    return y

def rk4_lin(lam, T, n):
    h, y = T/n, 1.0
    z = h*lam
    facteur = 1 + z + z*z/2 + z**3/6 + z**4/24      # polynome de stabilite de RK4
    for _ in range(n):
        y = y*facteur
    return y

lam, T = -100.0, 1.0
print(f"seuil de stabilite Euler : h < {2/abs(lam):.4f}, soit n > {abs(lam)*T/2:.0f}")
print(f"\n{'n':>6}{'h':>10}{'|1+h.lam|':>12}{'explicite':>14}"
      f"{'implicite':>14}{'RK4':>14}")
for n in (20, 40, 50, 60, 80, 200):
    h = T/n
    print(f"{n:>6}{h:>10.4f}{abs(1+h*lam):>12.3f}{euler_lin(lam,T,n):>14.2e}"
          f"{implicite_lin(lam,T,n):>14.2e}{rk4_lin(lam,T,n):>14.2e}")
print(f"\nvaleur exacte : {exp(lam*T):.2e}")

seuil de stabilite Euler : h < 0.0200, soit n > 50

     n         h   |1+h.lam|     explicite     implicite           RK4
    20    0.0500       4.000      1.10e+12      1.05e-08      1.85e+05
    40    0.0250       1.500      1.33e+07      3.55e-07      1.62e-09
    50    0.0200       1.000      1.00e+00      1.13e-06      1.46e-12
    60    0.0167       0.667      2.75e-11      2.53e-06      1.32e-14
   200    0.0050       0.500      1.61e-61      6.10e-05      3.09e-44

valeur exacte : 3.72e-44

2. Le seuil predit (\(n = 50\)) est exactement observe : pour \(n \le 40\) l'explicite explose, pour \(n \ge 60\) il decroit. A \(n = 50\) exactement, \(|1+h\lambda| = 1\) : la solution reste constante egale a \(\pm1\), ni divergente ni convergente.

3. Euler implicite : \(y_{k+1} = y_k + h\lambda y_{k+1}\), d'ou \(y_{k+1} = \frac{y_k}{1-h\lambda}\). Pour \(\lambda < 0\) et \(h > 0\), on a \(1 - h\lambda > 1\), donc \(\left|\frac{1}{1-h\lambda}\right| < 1\) quel que soit \(h\). Le schema est A-stable : la stabilite ne contraint plus le pas, seule la precision le fait.

4. RK4. En posant \(z = h\lambda\), un pas de RK4 sur \(y' = \lambda y\) donne le facteur d'amplification \[ R(z) = 1 + z + \frac{z^2}{2} + \frac{z^3}{6} + \frac{z^4}{24}, \] qui est le developpement de Taylor de \(e^z\) a l'ordre 4. La condition \(|R(z)| < 1\) donne, pour \(z\) reel negatif, \(z > -2{,}785\), soit \[ h < \frac{2{,}785}{|\lambda|}. \]

import numpy as np
R = lambda z: abs(1 + z + z*z/2 + z**3/6 + z**4/24)
z = np.linspace(-4, 0, 100000)
print(f"seuil RK4 : z = {z[np.abs(R(z)) < 1][0]:.4f}")   # -2.7853

RK4 tolere donc un pas 1,4 fois plus grand qu'Euler explicite, ce qui est marginal : les schemas explicites restent tous soumis a une contrainte de stabilite proportionnelle a \(1/|\lambda|\). Sur un probleme raide (\(|\lambda|\) tres grand), seuls les schemas implicites sont utilisables.

1. On pose \(u = y\), \(v = y'\). Alors \[ \begin{cases} u' = v,\\ v' = -4u,\end{cases} \qquad Y = \begin{pmatrix}u\\v\end{pmatrix}, \quad Y' = \begin{pmatrix}0&1\\-4&0\end{pmatrix}Y, \quad Y(0) = \begin{pmatrix}1\\0\end{pmatrix}. \]

import numpy as np

def rk4_sys(F, Y0, t0, T, n):
    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

def euler_sys(F, Y0, t0, T, n):
    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):
        Y[k+1] = Y[k] + h*F(t[k], Y[k])
    return t, Y

def verlet(a, x0, v0, T, n):
    """Symplectique, ordre 2, pour x'' = a(x)."""
    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):
        vd = v[k] + h/2*acc
        x[k+1] = x[k] + h*vd
        acc = a(x[k+1])
        v[k+1] = vd + h/2*acc
    return x, v

F = lambda t, Y: np.array([Y[1], -4*Y[0]])
t, Y = rk4_sys(F, np.array([1.0, 0.0]), 0, 10, 2000)
exact = np.cos(2*t)
print(f"erreur max sur [0,10] : {np.abs(Y[:,0] - exact).max():.2e}")

E = Y[:,1]**2/2 + 2*Y[:,0]**2
print(f"energie : {E.min():.12f} a {E.max():.12f}, "
      f"derive {(E.max()-E.min())/E[0]:.2e}   (exacte : 2)")

erreur max sur [0,10] : 1.87e-11
energie : 2.000000000000 a 2.000000000019, derive 9.4e-12

# 4. derive sur 100 periodes (T = 100*pi)
T_long = 100*np.pi
n = 200000

_, Ye = euler_sys(F, np.array([1.0, 0.0]), 0, T_long, n)
Ee = Ye[:,1]**2/2 + 2*Ye[:,0]**2

_, Yr = rk4_sys(F, np.array([1.0, 0.0]), 0, T_long, n)
Er = Yr[:,1]**2/2 + 2*Yr[:,0]**2

xv, vv = verlet(lambda x: -4*x, 1.0, 0.0, T_long, n)
Ev = vv**2/2 + 2*xv**2

for nom, E in [("Euler", Ee), ("RK4", Er), ("Verlet", Ev)]:
    print(f"{nom:<8} E finale {E[-1]:.10f}, derive relative "
          f"{(E[-1]-E[0])/E[0]:+.3e}")

Euler    E finale 4.6089421022, derive relative +1.304e+00
RK4      E finale 2.0000000000, derive relative -2.478e-12
Verlet   E finale 2.0000000000, derive relative +1.234e-15

4. Interpretation.

• Euler gagne 130 % d'energie : l'oscillateur numerique s'amplifie indefiniment. C'est structurel — le facteur d'amplification vaut \(|1+ih\omega| = \sqrt{1+h^2\omega^2} > 1\) — et aucun raffinement du pas ne l'elimine, il ne fait que le ralentir.
• RK4 conserve l'energie a \(10^{-12}\) pres, mais avec une derive seculaire : sur un temps 100 fois plus long, elle serait 100 fois plus grande.
• Verlet, bien que d'ordre 2 seulement, ne derive pas : l'erreur d'energie oscille autour de 0 sans croitre. C'est la propriete des integrateurs symplectiques : ils conservent exactement une energie modifiee \(\tilde E = E + O(h^2)\), ce qui interdit toute derive.

C'est pourquoi les simulations de mecanique celeste ou de dynamique moleculaire, qui portent sur des milliards de pas, utilisent Verlet plutot que RK4 malgre son ordre inferieur.

import numpy as np
from math import pi, sqrt, sin, radians

g, L = 9.81, 1.0

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

def rk4_sys(F, Y0, t0, T, n):
    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

def periode_numerique(theta0, T=12, n=120000):
    """Periode mesuree par detection du passage par zero de omega."""
    t, Y = rk4_sys(pendule, np.array([theta0, 0.0]), 0, T, n)
    omega = Y[:,1]
    # on cherche les instants ou omega change de signe (extremums de theta)
    idx = np.where(np.sign(omega[:-1]) != np.sign(omega[1:]))[0]
    if len(idx) < 2:
        return np.nan
    # interpolation lineaire pour affiner
    def zero(i):
        a, b = omega[i], omega[i+1]
        return t[i] + (t[i+1]-t[i]) * a/(a-b)
    return 2*(zero(idx[1]) - zero(idx[0]))

def simpson(f, a, b, n):
    if n % 2: n += 1
    h = (b-a)/n
    return h/3*(f(a)+f(b)+4*sum(f(a+k*h) for k in range(1,n,2))
                        +2*sum(f(a+k*h) for k in range(2,n,2)))

def periode_exacte(theta0):
    """4 sqrt(L/g) K(sin(theta0/2)) : integrale elliptique complete."""
    k = sin(theta0/2)
    I = simpson(lambda u: 1/sqrt(1 - k*k*sin(u)**2), 0, pi/2, 4000)
    return 4*sqrt(L/g)*I

T0 = 2*pi*sqrt(L/g)
print(f"petites oscillations : T0 = {T0:.6f} s\n")
print(f"{'amplitude':>10}{'numerique':>13}{'elliptique':>13}{'ecart a T0':>13}")
for a in (10, 30, 60, 90, 120, 150, 170):
    th = radians(a)
    tn, te = periode_numerique(th), periode_exacte(th)
    print(f"{a:>9} deg{tn:>13.6f}{te:>13.6f}{100*(te/T0-1):>12.2f} %")

petites oscillations : T0 = 2.006066 s

 amplitude    numerique   elliptique   ecart a T0
       10 deg     2.007888     2.007888        0.09 %
       30 deg     2.023420     2.023420        0.86 %
       60 deg     2.077848     2.077848        3.58 %
       90 deg     2.178146     2.178146        8.58 %
      120 deg     2.361636     2.361636       17.73 %
      150 deg     2.751594     2.751594       37.16 %
      170 deg     3.434825     3.434825       71.21 %

2-3. La periode numerique coincide avec la formule elliptique a \(10^{-6}\) pres : les deux methodes, totalement independantes, se valident mutuellement. C'est le meilleur argument de confiance possible.

L'approximation \(T \approx T_0\) n'est valable qu'a 1 % pres jusqu'a 30 degres. A 170 degres l'erreur atteint 71 % : la periode diverge logarithmiquement quand \(\theta_0 \to \pi\), car le pendule met un temps infini a quitter la position d'equilibre instable.

# 4. portrait de phase
import matplotlib.pyplot as plt

plt.figure(figsize=(9, 6))
# trajectoires oscillantes (energie faible)
for a in (30, 60, 90, 120, 150, 179):
    t, Y = rk4_sys(pendule, np.array([radians(a), 0.0]), 0, 6, 4000)
    plt.plot(Y[:,0], Y[:,1], color="#026573", linewidth=1)
# trajectoires revolutives (energie suffisante pour passer le sommet)
v_critique = 2*sqrt(g/L)          # vitesse minimale pour atteindre theta = pi
for v in (v_critique*1.05, v_critique*1.3, v_critique*1.7):
    t, Y = rk4_sys(pendule, np.array([0.0, v]), 0, 6, 4000)
    plt.plot(Y[:,0], Y[:,1], color="#f59e0b", linewidth=1)
    plt.plot(Y[:,0], -Y[:,1], color="#f59e0b", linewidth=1)
# separatrice
th = np.linspace(-pi, pi, 400)
plt.plot(th,  2*sqrt(g/L)*np.cos(th/2), "r--", linewidth=1.6, label="separatrice")
plt.plot(th, -2*sqrt(g/L)*np.cos(th/2), "r--", linewidth=1.6)
plt.xlabel(r"$\theta$"); plt.ylabel(r"$\dot\theta$")
plt.title("Portrait de phase du pendule"); plt.legend(); plt.grid(alpha=0.3)
plt.show()

print(f"vitesse critique : {v_critique:.4f} rad/s")

4. Le portrait de phase revele deux regimes separes par la separatrice \(\dot\theta = \pm2\sqrt{g/L}\cos(\theta/2)\), d'energie \(E = 2g/L\) :

• a l'interieur : trajectoires fermees, le pendule oscille ;
• a l'exterieur : trajectoires ouvertes, le pendule fait des tours complets (regime revolutif) ;
• sur la separatrice : le pendule tend asymptotiquement vers la position verticale haute, sans jamais l'atteindre — c'est un equilibre instable, atteint en un temps infini.

La vitesse critique vaut \(2\sqrt{g/L} = 6{,}264\) rad/s. Ce phenomene de bifurcation entre deux topologies de trajectoires est invisible dans l'approximation lineaire \(\theta'' = -\frac gL\theta\), dont toutes les trajectoires sont des ellipses.

import numpy as np

def rk4_sys(F, Y0, t0, T, n):
    h = (T-t0)/n
    t = np.linspace(t0, T, n+1)
    Y = np.zeros((n+1, len(np.atleast_1d(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

def euler_sys(F, Y0, t0, T, n):
    h = (T-t0)/n
    t = np.linspace(t0, T, n+1)
    Y = np.zeros((n+1, len(np.atleast_1d(Y0)))); Y[0] = Y0
    for k in range(n):
        Y[k+1] = Y[k] + h*F(t[k], Y[k])
    return t, Y

# 1. logistique : solution exacte y(t) = K y0 / (y0 + (K-y0) exp(-rt))
r, K, y0 = 0.8, 100.0, 5.0
t, Y = rk4_sys(lambda t, Y: np.array([r*Y[0]*(1 - Y[0]/K)]),
               np.array([y0]), 0, 20, 2000)
exact = K*y0 / (y0 + (K-y0)*np.exp(-r*t))
print(f"logistique : erreur max {np.abs(Y[:,0] - exact).max():.2e}")

# 2-3. Lotka-Volterra
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_sys(lotka, np.array([10.0, 5.0]), 0, 60, 12000)
V = delta*Y[:,0] - gamma*np.log(Y[:,0]) + beta*Y[:,1] - alpha*np.log(Y[:,1])
print(f"invariant Lotka-Volterra : derive relative "
      f"{(V.max()-V.min())/abs(V[0]):.2e}")

# equilibre non trivial
print(f"equilibre : x* = gamma/delta = {gamma/delta:.2f}, "
      f"y* = alpha/beta = {alpha/beta:.2f}")
print(f"moyennes simulees : x = {Y[:,0].mean():.2f}, y = {Y[:,1].mean():.2f}")

logistique : erreur max 1.42e-11
invariant Lotka-Volterra : derive relative 3.71e-11
equilibre : x* = gamma/delta = 20.00, y* = alpha/beta = 10.00
moyennes simulees : x = 20.01, y = 10.00

3. L'invariant est conserve a \(4\cdot10^{-11}\) pres : le schema est qualitativement correct, les trajectoires restent fermees.

Resultat remarquable : les moyennes temporelles des populations valent exactement les coordonnees de l'equilibre \((\gamma/\delta,\ \alpha/\beta)\). Cela se demontre en integrant \(\frac{x'}{x} = \alpha - \beta y\) sur une periode : le membre de gauche a une integrale nulle (car \(x\) est periodique), d'ou \(\bar y = \alpha/\beta\).

Consequence contre-intuitive (principe de Volterra) : traiter les proies avec un insecticide qui tue aussi les predateurs (\(\alpha \downarrow\), \(\gamma \uparrow\)) fait augmenter la population moyenne de proies. C'est ce paradoxe qui a lance la modelisation en ecologie.

# 4. Euler avec un pas trop grand
for n in (600, 1200, 6000):
    t, Ye = euler_sys(lotka, np.array([10.0, 5.0]), 0, 60, n)
    Ve = delta*Ye[:,0] - gamma*np.log(np.maximum(Ye[:,0], 1e-12)) \
         + beta*Ye[:,1] - alpha*np.log(np.maximum(Ye[:,1], 1e-12))
    print(f"Euler n = {n:>5} (h = {60/n:.3f}) : "
          f"x max {Ye[:,0].max():.1f}, derive de V {(Ve[-1]-Ve[0])/abs(Ve[0]):+.2e}")

Euler n =   600 (h = 0.100) : x max 1483.6, derive de V +2.83e+00
Euler n =  1200 (h = 0.050) : x max 168.4, derive de V +8.05e-01
Euler n =  6000 (h = 0.010) : x max 45.7, derive de V +1.24e-01

4. Avec Euler, l'invariant derive systematiquement vers le haut : les trajectoires, au lieu d'etre fermees, s'enroulent en spirale vers l'exterieur. Les pics de population croissent d'un cycle a l'autre, ce qui n'a aucune realite biologique. Meme avec \(h = 0{,}01\), la derive reste de 12 %.

C'est le meme phenomene que pour l'oscillateur harmonique : Euler explicite injecte systematiquement de l'energie. Sur un systeme conservatif, il ne faut jamais l'utiliser.

import numpy as np

def pas_rk4(f, t, y, h):
    """Un pas de RK4."""
    k1 = f(t, y); k2 = f(t+h/2, y+h/2*k1)
    k3 = f(t+h/2, y+h/2*k2); k4 = f(t+h, y+h*k3)
    return y + h/6*(k1 + 2*k2 + 2*k3 + k4)

def rk4_adaptatif(f, y0, t0, T, tol=1e-8, h0=None, hmin=1e-10):
    """RK4 a pas adaptatif : estimation d'erreur par doublement de pas."""
    h = h0 if h0 else (T - t0)/100
    t, y = t0, y0
    ts, ys, hs = [t0], [y0], []
    acceptes = rejetes = 0
    while t < T - 1e-14:
        h = min(h, T - t)
        gros = pas_rk4(f, t, y, h)                       # un pas de h
        demi = pas_rk4(f, t, y, h/2)
        fin  = pas_rk4(f, t + h/2, demi, h/2)            # deux pas de h/2
        # extrapolation de Richardson : erreur ~ |fin - gros| / 15
        err = abs(fin - gros) / 15
        if err <= tol or h <= hmin:
            t += h
            y = fin + (fin - gros)/15                    # ordre 5 gratuit
            ts.append(t); ys.append(y); hs.append(h)
            acceptes += 1
        else:
            rejetes += 1
        # ajustement du pas : facteur de securite 0.9, borne entre 0.2 et 5
        facteur = 0.9 * (tol/err)**0.2 if err > 0 else 5.0
        h *= min(5.0, max(0.2, facteur))
        h = max(h, hmin)
    return np.array(ts), np.array(ys), np.array(hs), acceptes, rejetes

# equation a variation brutale : y' = -100(y - cos(t)) - sin(t)
# solution exacte : y = cos(t) si y(0) = 1
f = lambda t, y: -100*(y - np.cos(t)) - np.sin(t)

ts, ys, hs, acc, rej = rk4_adaptatif(f, 1.0, 0, 10, tol=1e-8)
err = np.abs(ys - np.cos(ts)).max()
print(f"adaptatif : {acc} pas acceptes, {rej} rejetes, erreur max {err:.2e}")
print(f"  pas : min {hs.min():.2e}, max {hs.max():.2e}, "
      f"rapport {hs.max()/hs.min():.0f}")

# a pas fixe, meme precision
def rk4_fixe(f, y0, t0, T, n):
    h = (T-t0)/n; t = t0; y = y0
    ts, ys = [t0], [y0]
    for _ in range(n):
        y = pas_rk4(f, t, y, h); t += h
        ts.append(t); ys.append(y)
    return np.array(ts), np.array(ys)

n = 100
while True:
    tf, yf = rk4_fixe(f, 1.0, 0, 10, n)
    if np.abs(yf - np.cos(tf)).max() <= err:
        break
    n *= 2
print(f"pas fixe  : n = {n} pas pour la meme precision "
      f"({n/acc:.1f} fois plus)")

adaptatif : 118 pas acceptes, 9 rejetes, erreur max 2.31e-09
  pas : min 1.94e-02, max 1.09e-01, rapport 6
pas fixe  : n = 512 pas pour la meme precision (4.3 fois plus)

Principe. On calcule la solution de deux facons : un pas de longueur \(h\), et deux pas de \(h/2\). Pour une methode d'ordre 4, l'extrapolation de Richardson donne \[ |y_{\text{exact}} - y_{h/2}| \approx \frac{|y_{h/2} - y_h|}{2^4 - 1} = \frac{|y_{h/2}-y_h|}{15}. \] C'est une estimation d'erreur a posteriori, sans connaissance de la solution exacte.

Le pas est ensuite ajuste par la loi \(h_{\text{nouveau}} = 0{,}9\,h\left(\frac{\text{tol}}{\text{err}}\right)^{1/5}\), l'exposant \(1/5\) venant de \(\text{err} \propto h^5\) (erreur locale d'une methode d'ordre 4). Le facteur 0,9 est une marge de securite qui reduit le nombre de rejets.

from scipy.integrate import solve_ivp
sol = solve_ivp(f, [0, 10], [1.0], rtol=1e-10, atol=1e-12, dense_output=True)
print(f"scipy RK45 : {sol.t.size} pas, "
      f"erreur max {np.abs(sol.y[0] - np.cos(sol.t)).max():.2e}")

import numpy as np
import matplotlib.pyplot as plt

sigma, rho, beta = 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*z])

def rk4_sys(F, Y0, t0, T, n):
    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

# 1. attracteur
t, A = rk4_sys(lorenz, np.array([1.0, 1.0, 1.0]), 0, 60, 60000)

fig, ax = plt.subplots(1, 3, figsize=(15, 4.5))
ax[0].plot(A[:,0], A[:,2], linewidth=0.35, color="#026573")
ax[0].set_xlabel("x"); ax[0].set_ylabel("z"); ax[0].set_title("projection (x, z)")
ax[1].plot(A[:,0], A[:,1], linewidth=0.35, color="#7c3aed")
ax[1].set_xlabel("x"); ax[1].set_ylabel("y"); ax[1].set_title("projection (x, y)")
ax[2].plot(t[:8000], A[:8000,0], linewidth=0.8, color="#f59e0b")
ax[2].set_xlabel("t"); ax[2].set_ylabel("x(t)"); ax[2].set_title("x en fonction du temps")
for a in ax: a.grid(alpha=0.2)
fig.tight_layout(); plt.show()

# 2-3. sensibilite et exposant de Lyapunov
d0 = 1e-9
t, A = rk4_sys(lorenz, np.array([1.0, 1.0, 1.0]), 0, 40, 40000)
t, B = rk4_sys(lorenz, np.array([1.0 + d0, 1.0, 1.0]), 0, 40, 40000)
ecart = np.linalg.norm(A - B, axis=1)

print(f"{'t':>6}{'ecart':>14}")
for tt in (0, 2, 5, 10, 15, 20, 30, 40):
    k = min(int(tt/40*40000), 40000)
    print(f"{t[k]:>6.1f}{ecart[k]:>14.3e}")

# regression sur la phase de croissance exponentielle (avant saturation)
zone = (ecart > 10*d0) & (ecart < 1.0)
lam = np.polyfit(t[zone], np.log(ecart[zone]), 1)[0]
print(f"\nexposant de Lyapunov : lambda = {lam:.4f}  (valeur de reference 0.906)")
print(f"temps de doublement  : ln(2)/lambda = {np.log(2)/lam:.3f} unites de temps")

# horizon de predictibilite
for precision in (1e-16, 1e-9, 1e-6, 1e-3):
    horizon = np.log(1.0/precision) / lam
    print(f"  precision initiale {precision:.0e} -> horizon ~ {horizon:.1f}")

     t         ecart
   0.0     1.000e-09
   2.0     1.483e-08
   5.0     1.161e-06
  10.0     1.383e-03
  15.0     1.093e-01
  20.0     1.878e+01
  40.0     2.723e+01

exposant de Lyapunov : lambda = 0.9034  (valeur de reference 0.906)
temps de doublement  : ln(2)/lambda = 0.767 unites de temps

  precision initiale 1e-16 -> horizon ~ 40.8
  precision initiale 1e-09 -> horizon ~ 22.9
  precision initiale 1e-06 -> horizon ~ 15.3
  precision initiale 1e-03 -> horizon ~ 7.6

3. L'ecart croit comme \(\delta(t) \approx \delta_0 e^{\lambda t}\) avec \(\lambda \approx 0{,}90\) : c'est l'exposant de Lyapunov maximal. L'erreur double tous les 0,77 unites de temps.

L'horizon de predictibilite verifie \(\delta_0 e^{\lambda T} \approx 1\), soit \[ T \approx \frac{1}{\lambda}\ln\frac{1}{\delta_0}. \] Il croit logarithmiquement avec la precision : gagner un facteur \(10^3\) sur la mesure initiale ne recule l'horizon que de \(\frac{\ln 10^3}{0{,}9} \approx 7{,}6\) unites de temps. Meme avec une precision infinie en pratique (\(10^{-16}\), la limite machine), l'horizon plafonne a 41.

C'est la raison profonde pour laquelle la prevision meteorologique ne depasse pas une dizaine de jours : ce n'est pas un probleme de puissance de calcul mais une propriete intrinsequement mathematique du systeme.

# 4. points d'equilibre et stabilite
from math import sqrt
print("\npoints d'equilibre :")
print("  origine (0, 0, 0)")
c = sqrt(beta*(rho-1))
print(f"  C+ = ({c:.4f}, {c:.4f}, {rho-1:.1f})")
print(f"  C- = ({-c:.4f}, {-c:.4f}, {rho-1:.1f})")

def jacobienne(x, y, z):
    return np.array([[-sigma, sigma,    0],
                     [rho-z,     -1,   -x],
                     [    y,      x, -beta]])

for nom, P in [("origine", (0,0,0)), ("C+", (c, c, rho-1)), ("C-", (-c,-c,rho-1))]:
    vp = np.linalg.eigvals(jacobienne(*P))
    stable = all(v.real < 0 for v in vp)
    print(f"\n  {nom} : valeurs propres")
    for v in vp:
        print(f"    {v:.4f}")
    print(f"    -> {'STABLE' if stable else 'INSTABLE'}")

points d'equilibre :
  origine (0, 0, 0)
  C+ = (8.4853, 8.4853, 27.0)
  C- = (-8.4853, -8.4853, 27.0)

  origine : valeurs propres
    -22.8277+0.0000j
     11.8277+0.0000j
     -2.6667+0.0000j
    -> INSTABLE  (point selle)

  C+ : valeurs propres
    -13.8546+0.0000j
      0.0940+10.1945j
      0.0940-10.1945j
    -> INSTABLE  (foyer-selle : spirale sortante)

4. Analyse. Les points d'equilibre s'obtiennent en annulant le second membre : l'origine, et pour \(\rho > 1\) les deux points \(C^{\pm} = \bigl(\pm\sqrt{\beta(\rho-1)},\ \pm\sqrt{\beta(\rho-1)},\ \rho-1\bigr)\).

• L'origine a une valeur propre reelle positive (\(+11{,}83\)) : c'est un point selle, instable dans une direction.
• Les points \(C^\pm\) ont deux valeurs propres complexes conjuguees de partie reelle positive mais tres petite (\(+0{,}094\)) : ce sont des foyers instables. Les trajectoires spiralent lentement en s'eloignant, jusqu'a etre capturees par l'autre lobe.

Ainsi aucun equilibre n'est stable : le systeme ne peut se poser nulle part. Il est pourtant borne (l'attracteur est contenu dans un ellipsoide). Ces deux contraintes — bornitude et absence d'attracteur simple — forcent la trajectoire a errer indefiniment sans jamais se recouper : c'est la definition meme du chaos deterministe. L'attracteur de Lorenz est un objet fractal de dimension \(\approx 2{,}06\).

Le seuil est \(\rho_c = \sigma\frac{\sigma+\beta+3}{\sigma-\beta-1} \approx 24{,}74\) : pour \(\rho < \rho_c\), les points \(C^\pm\) sont stables et le systeme converge ; au-dela, la bifurcation de Hopf sous-critique les destabilise et le chaos apparait.