Pourquoi Euler fait spiraler l'électron

La stratégie. Au lieu de constater la dérive sur un graphe, on calcule ce que fait un pas du schéma, exactement, en variable complexe : le défaut apparaîtra comme un facteur multiplicatif, qu'il suffira d'élever à la puissance du nombre de pas.

1) Un pas d'Euler. Avec les deux mises à jour, $${\underline{u}}_{n+1}=(v_x+\omega v_y\mathrm{d}t)+i(v_y-\omega v_x\mathrm{d}t)=(v_x+i{v}_y)-i\omega \mathrm{d}t(v_x+i{v}_y)=(1-i\omega \mathrm{d}t){\underline{u}}_n.$$ C'est la discrétisation de $\dot{\underline{u}}=-i\omega \underline{u}$ (exercice 5). Le module d'un produit est le produit des modules : $|{\underline{u}}_{n+1}|=|1-i\omega \mathrm{d}t||{\underline{u}}_n|=\sqrt{1+(\omega \mathrm{d}t{)}^2}|{\underline{u}}_n|$. Ce facteur est strictement supérieur à $1$ : à chaque pas, la vitesse grandit, quel que soit $\mathrm{d}t$. Le schéma exact multiplierait $\underline{u}$ par $e^{-i\omega \mathrm{d}t}$, de module $1$ : Euler remplace un point du cercle unité par un point de sa tangente, extérieur au cercle.

2) La spirale. Un tour compte $N=\frac{2\pi }{|\omega |\mathrm{d}t}$ pas. Le facteur sur un tour vaut $$(1+(\omega \mathrm{d}t{)}^2{)}^{N/2}=\exp \left(\frac{\pi }{|\omega |\mathrm{d}t}\ln \right(1+(\omega \mathrm{d}t{)}^2))\approx \exp (\pi |\omega |\mathrm{d}t),$$ puisque $\ln (1+\epsilon )\approx \epsilon$ pour $\epsilon$ petit. Avec les données, $|\omega |=\frac{1,6\times 10^{-19}\times 1,0\times 10^{-3}}{9,1\times 10^{-31}}\approx 1,76\times 10^8\mathrm{rad}{\mathrm{s}}^{-1}$, $|\omega |\mathrm{d}t\approx 1,76\times 10^{-3}$ (soit environ $3600$ pas par tour), et le facteur par tour vaut $e^{5,5\times 10^{-3}}$ : le rayon augmente de $0,55\%$ par tour. L'erreur est cumulative : au bout de cent tours, le facteur est $e^{0,55}\approx 1,74$. La trajectoire calculée n'est pas un cercle, c'est une spirale qui s'ouvre. Diviser le pas par dix divise la dérive par dix et multiplie le temps de calcul par dix : on recule le problème sans le supprimer.

3) Le schéma semi-implicite. Notons $(x,y)$ la vitesse avant le pas et $(x^{\prime },y^{\prime })$ après : $x^{\prime }=x+hy$, puis $y^{\prime }=y-h{x}^{\prime }$. On calcule $$x^{\prime 2}+y^{\prime 2}+h{x}^{\prime }{y}^{\prime }=x^{\prime 2}+y^2-2h{x}^{\prime }y+h^2{x}^{\prime 2}+h{x}^{\prime }y-h^2{x}^{\prime 2}=x^{\prime 2}-h{x}^{\prime }y+y^2=x^{\prime }(x^{\prime }-hy)+y^2.$$ Or $x^{\prime }-hy=x$, donc $x^{\prime }(x^{\prime }-hy)=x(x+hy)=x^2+hxy$, et finalement $$x^{\prime 2}+y^{\prime 2}+h{x}^{\prime }{y}^{\prime }=x^2+y^2+hxy:$$ la quantité $Q=v_x^2+v_y^2+h{v}_x{v}_y$ est exactement conservée, à chaque pas. Comme $|h{v}_x{v}_y|\le \frac{|h|}{2}(v_x^2+v_y^2)$, on a $$(1-\frac{|h|}{2})v^2\le Q\le (1+\frac{|h|}{2})v^2,$$ donc $v^2$ reste compris entre $\frac{Q}{1+|h|/2}$ et $\frac{Q}{1-|h|/2}$ : la norme de la vitesse oscille de moins de $|h|/2\approx 0,1\%$ autour de sa valeur, sans jamais dériver. Le schéma ne conserve pas exactement l'énergie cinétique, mais il conserve exactement une quantité qui lui est très voisine, et c'est suffisant pour que la trajectoire reste fermée.

4) Le script.

import numpy as np
q, m, B = -1.6e-19, 9.1e-31, 1.0e-3
w = q*B/m                          # pulsation algebrique
v0, dt = 1.0e5, 1.0e-11
N = int(round(2*np.pi/abs(w)/dt))  # pas par tour
def euler(vx, vy):                 # les deux composantes avec l'ancien v
    return vx + w*vy*dt, vy - w*vx*dt
def semi(vx, vy):                  # v_y calcule avec le NOUVEAU v_x
    vx = vx + w*vy*dt
    return vx, vy - w*vx*dt
for nom, pas in (("Euler explicite", euler), ("semi-implicite", semi)):
    vx, vy = v0, 0.0
    for tour in range(1, 101):
        for _ in range(N):
            vx, vy = pas(vx, vy)
        if tour in (1, 100):
            print(nom, tour, np.hypot(vx, vy)/v0)

Résultat obtenu en l'exécutant : avec Euler explicite, $|\overrightarrow{v}|/v_0=1,00554$ après un tour et $1,737$ après cent tours ; avec le schéma semi-implicite, $1,00000$ après un tour et $1,00007$ après cent tours. Les valeurs d'Euler sont exactement celles que prévoit le 2) ($e^{0,0055}$ et $e^{0,55}$) ; le schéma semi-implicite, pour le même coût, ne dérive pas.

Ce que l'exercice installe. L'analyse d'erreur d'un schéma numérique : on ne juge pas un schéma sur un graphe, on calcule ce qu'il fait en un pas, et l'on cherche s'il conserve — exactement ou presque — ce que la physique conserve. Un défaut d'un facteur $1+\epsilon$ par pas devient une exponentielle sur un grand nombre de pas : c'est la raison pour laquelle on n'intègre jamais une orbite par Euler explicite.