La méthode d'Euler : du schéma à la simulation

La stratégie. Avant de faire confiance à un schéma numérique, on le teste sur un cas dont on connaît la solution exacte. Le régime libre est ce cas : le schéma d'Euler s'y résout à la main, et l'on voit d'un coup sa précision et sa limite de stabilité. Ensuite seulement, on l'utilise là où il est utile : une excitation qui n'est pas un échelon.

1) Le schéma résolu. Avec $e=0$ : $u_{n+1}=u_n-\frac{h}{\tau }{u}_n=\left(1-\frac{h}{\tau }\right)u_n$. La suite est géométrique : $u_n=u_0{\left(1-\frac{h}{\tau }\right)}^n$. En $t=\tau$, $n=\tau /h$ :

Diviser le pas par $10$ divise l'erreur par $10$ : l'erreur de la méthode d'Euler est proportionnelle au pas. (On le voit sur le développement $\ln (1-x)=-x-\frac{x^2}{2}-\dots$ : ${\left(1-\frac{h}{\tau }\right)}^{\tau /h}=\exp \left(-1-\frac{h}{2\tau }-\dots \right)\approx e^{-1}\left(1-\frac{h}{2\tau }\right)$, d'où une erreur relative $\approx -h/(2\tau )$, soit $-10\%$, $-1\%$, $-0,1\%$.) Le schéma est toujours en dessous en décharge : il utilise la pente au début du pas, plus raide que la pente moyenne.

2) La stabilité. Le facteur $1-h/\tau$ commande tout.

• Pour $h<\tau$, il est entre $0$ et $1$ : la suite décroît vers $0$ sans changer de signe, comme la solution exacte.
• Pour $\tau <h<2\tau$, il est entre $-1$ et $0$ : la suite tend encore vers $0$, mais en changeant de signe à chaque pas. Avec $h=1,5\tau$ : $1$ ; $-0,5$ ; $0,25$ ; $-0,125$… Une décharge qui oscille : qualitativement faux.
• Pour $h>2\tau$, il est inférieur à $-1$ : la suite diverge en oscillant. Avec $h=2,5\tau$ : $1$ ; $-1,5$ ; $2,25$ ; $-3,4$ ; $5,1$…

Un pas trop grand ne donne pas un résultat imprécis, il donne un résultat absurde : une tension qui croît sans limite dans un circuit sans source. La règle pratique est $h\ll \tau$ ; et dans un circuit qui comporte plusieurs constantes de temps, c'est la plus petite qui impose le pas.

3) Le script. On choisit un pas, on avance de tangente en tangente, et l'on compare à la solution exacte. Pour $t>\theta$, la solution exacte repart de $u(\theta )$ vers $E$ avec la constante de temps $\tau$.

import numpy as np

E, tau, theta = 5.0, 1e-3, 5e-3    # V, s, s (theta = 5 tau)

def e(t):                           # rampe saturée
    return np.minimum(E * t / theta, E)

def euler(h, t_fin=10e-3):
    t = np.arange(0.0, t_fin + h / 2, h)
    u = np.zeros_like(t)            # u(0) = 0
    for n in range(len(t) - 1):
        u[n + 1] = u[n] + h * (e(t[n]) - u[n]) / tau
    return t, u

def exacte(t):                      # solution exacte
    u1 = E / theta * (t - tau + tau * np.exp(-t / tau))
    u_th = E / theta * (theta - tau + tau * np.exp(-theta / tau))
    u2 = E + (u_th - E) * np.exp(-(t - theta) / tau)
    return np.where(t <= theta, u1, u2)

for k in (10, 100, 1000):
    h = tau / k
    t, u = euler(h)
    n = int(round(theta / h))
    err = np.max(np.abs(u - exacte(t)))
    print(f"h = tau/{k:<4d}  u(theta) = {u[n]:.4f} V   "
          f"erreur max = {err:.5f} V")
print(f"exact        u(theta) = {exacte(np.array(theta)):.4f} V")
t, u = euler(tau / 1000)
n = int(round(theta / (tau / 1000)))
print(f"retard e - u en theta : {e(t[n]) - u[n]:.3f} V")

Résultat obtenu en l'exécutant :

h = tau/10    u(theta) = 4.0052 V   erreur max = 0.01920 V
h = tau/100   u(theta) = 4.0066 V   erreur max = 0.00185 V
h = tau/1000  u(theta) = 4.0067 V   erreur max = 0.00018 V
exact        u(theta) = 4.0067 V
retard e - u en theta : 0.993 V

On retrouve la proportionnalité de l'erreur au pas. La valeur exacte $u(\theta )=E\left(1-\frac{1}{5}+\frac{e^{-5}}{5}\right)=0,8013E=4,007\mathrm{V}$ est reproduite à $0,2\mathrm{mV}$ près avec $h=\tau /1000$.

4) Le retard. Pour $t\gg \tau$ pendant la rampe, le terme $\tau e^{-t/\tau }$ s'éteint et $$u(t)\approx \frac{E}{\theta }(t-\tau )=e(t-\tau ).$$ La sortie suit l'entrée avec un retard $\tau$ ; à un instant donné, elle est en dessous de l'entrée de $E\tau /\theta =1,0\mathrm{V}$ (le script donne $0,993\mathrm{V}$ en $\theta$, l'exponentielle n'étant pas encore tout à fait éteinte). Physiquement : pour que $u$ monte à la vitesse $E/\theta$, il faut un courant constant $i=CE/\theta$, donc une tension constante $Ri=E\tau /\theta$ aux bornes de la résistance. On retrouvera ce courant constant au problème : c'est la clé de la charge sans pertes.

Ce que l'exercice installe. Tester un schéma sur un cas exactement soluble avant de s'en servir ; savoir que l'erreur d'Euler est proportionnelle au pas et qu'un pas supérieur à $2\tau$ fait diverger le calcul ; utiliser la simulation pour une excitation qui n'est pas un échelon. Et une propriété physique : un circuit du premier ordre soumis à une rampe la suit avec le retard $\tau$.