La méthode d'Euler est-elle stable ?

La stratégie. Dans les deux situations, le schéma d'Euler engendre une suite arithmético-géométrique — réelle pour le condensateur, géométrique complexe pour l'oscillateur —, dont le comportement à long terme se lit sur la raison. C'est le cours sur la récurrence $u_{k+1}=a{u}_k+b$ et sur les suites géométriques complexes, qui convergent si et seulement si la raison est de module strictement inférieur à $1$, ou égale à $1$.

a) Le condensateur. L'équation s'écrit $y^{\prime }=F(y)$ avec $F(y)=\frac{E-y}{\tau }$. Le schéma d'Euler donne $$y_{k+1}=y_k+\frac{h}{\tau }(E-y_k)=\left(1-\frac{h}{\tau }\right)y_k+\frac{h}{\tau }E.$$ C'est une récurrence arithmético-géométrique $y_{k+1}=q{y}_k+b$, de coefficient $q=1-\frac{h}{\tau }$, différent de $1$ puisque $h>0$, et de terme constant $b=\frac{h}{\tau }E$. On suit la méthode du cours. La solution constante $c$ vérifie $c=qc+b$, soit $c=\frac{b}{1-q}=\frac{hE/\tau }{h/\tau }=E$ : c'est la tension finale, ce qui est rassurant. La suite $(y_k-E)$ est géométrique de raison $q$, donc $y_k-E=q^k(y_0-E)=-E{q}^k$, et $$y_k=E\left(1-{\left(1-\frac{h}{\tau }\right)}^k\right).$$ Par conséquent $y_k\to E$ si et seulement si $q^k\to 0$, c'est-à-dire, d'après les limites géométriques de référence du cours, si et seulement si $|q|<1$ : pour $q\le -1$, la suite $(q^k)$ n'a pas de limite. Or $$-1<1-\frac{h}{\tau }<1⟺0<\frac{h}{\tau }<2⟺0<h<2\tau .$$ Le schéma converge vers la bonne limite si et seulement si $h<2\tau$. Les cinq régimes se lisent sur le signe et la taille de $q$ :

• pour $0<h<\tau$ : $q\in \left]0,1\right[$, la suite $y_k=E\left(1-q^k\right)$ croît vers $E$, comme la solution exacte ;
• pour $h=\tau$ : $q=0$, et $y_k=E$ dès $k=1$ : la suite atteint en un pas une valeur que la solution exacte n'atteint jamais ;
• pour $\tau <h<2\tau$ : $q\in \left]-1,0\right[$, le signe de $y_k-E=-E{q}^k$ alterne : la suite converge en oscillant autour de $E$ ; pour $h=1,5\tau$, $\frac{y_k}{E}$ vaut $0$ ; $1,5$ ; $0,75$ ; $1,125$ ; $0,9375$ ;
• pour $h=2\tau$ : $q=-1$, et $y_k=E\left(1-(-1{)}^k\right)$ vaut alternativement $0$ et $2E$ : la suite diverge ;
• pour $h>2\tau$ : $q<-1$, et $|y_k-E|=E|q{|}^k\to +\infty$ : la suite diverge en oscillant, avec une amplitude qui explose.

Le point délicat. La solution exacte converge vers $E$ quel que soit $h$ : la divergence est fabriquée par le schéma. C'est la notion de stabilité numérique : un pas trop grand devant la constante de temps rend le calcul absurde — une tension qui dépasse $2E$ dans un circuit alimenté sous la tension $E$.

b) L'oscillateur, en complexe. La fonction $z=x^{\prime }+i\omega x$ est dérivable, de dérivée $z^{\prime }=x^{\prime \prime }+i\omega x^{\prime }$, avec $x^{\prime \prime }=-{\omega }^{2}x=i\omega \cdot i\omega x$ ; donc $$z^{\prime }=-{\omega }^{2}x+i\omega x^{\prime }=i\omega (x^{\prime }+i\omega x)=i\omega z.$$ Comme $x$ et $x^{\prime }$ sont réels, $|z{|}^2=x^{\prime 2}+{\omega }^2{x}^2$. Dérivons cette fonction réelle : $$(x^{\prime 2}+{\omega }^2{x}^2{)}^{\prime }=2x^{\prime }{x}^{\prime \prime }+2{\omega }^{2}x{x}^{\prime }=2x^{\prime }(x^{\prime \prime }+{\omega }^{2}x)=0.$$ Une fonction de dérivée nulle sur un intervalle est constante : $|z{|}^2$ est constant, c'est la conservation de l'énergie mécanique.

c) Euler fabrique de l'énergie. Le schéma $z_{k+1}=(1+i\omega h)z_k$ définit une suite géométrique complexe de raison $1+i\omega h$, de module $$|1+i\omega h|=\sqrt{1+{\omega }^2{h}^2}>1.$$ Donc $|z_k|={\left(\sqrt{1+{\omega }^2{h}^2}\right)}^k|z_0|$, suite géométrique réelle de raison strictement supérieure à $1$, multipliée par $|z_0|>0$ : $|z_k|\to +\infty$. L'énergie calculée, proportionnelle à $|z_k{|}^2$, croît sans limite, quel que soit le pas : aucun pas n'est assez petit pour rendre ce schéma stable. (Ce n'est pas un artifice du passage au complexe : le schéma d'Euler appliqué au système réel $x^{\prime }=v$, $v^{\prime }=-{\omega }^{2}x$, soit $x_{k+1}=x_k+h{v}_k$ et $v_{k+1}=v_k-h{\omega }^2{x}_k$, donne exactement $z_{k+1}=(1+i\omega h)z_k$ pour $z_k=v_k+i\omega x_k$.)

Pour $\omega h=0,1$, une période $\frac{2\pi }{\omega }$ correspond à $\frac{2\pi }{\omega h}=20\pi \approx 62,8$ pas. L'énergie est multipliée à chaque pas par $|1+i\omega h{|}^2=1,01$, donc sur une période par $$1,01^{20\pi }\approx 1,87,$$ et l'amplitude par $\sqrt{1,87}\approx 1,37$ : au bout de dix périodes, l'énergie a été multipliée par plus de cinq cents.

Le schéma implicite. La relation s'écrit $(1-i\omega h)z_{k+1}=z_k$, et comme $1-i\omega h\ne 0$, $z_{k+1}=\frac{1}{1-i\omega h}{z}_k$ : suite géométrique de raison de module $\frac{1}{\sqrt{1+{\omega }^2{h}^2}}<1$, donc $z_k\to 0$. Le schéma implicite ne fabrique plus d'énergie, il en détruit : l'oscillateur calculé s'amortit — pour $\omega h=0,1$, l'amplitude est multipliée par environ $0,73$ à chaque période —, alors que le vrai ne s'amortit pas.

d) Le schéma des trapèzes. La relation s'écrit $\left(1-\frac{i\omega h}{2}\right)z_{k+1}=\left(1+\frac{i\omega h}{2}\right)z_k$, d'où $$z_{k+1}=\frac{1+\frac{i\omega h}{2}}{1-\frac{i\omega h}{2}}{z}_k.$$ Le numérateur et le dénominateur de la raison sont deux complexes conjugués, donc de même module : la raison est de module $1$, et $|z_{k+1}|=|z_k|$ pour tout $k$. Le schéma des trapèzes conserve exactement l'énergie — juste milieu entre l'explicite, qui prend la pente au début du pas, et l'implicite, qui la prend à la fin.

Ce que l'exercice installe. Le comportement à long terme d'une suite $u_{k+1}=a{u}_k+b$ se lit sur $|a|$ : c'est le critère de stabilité de tout schéma numérique linéaire. En physique, où le programme demande de mettre en œuvre la méthode d'Euler, on retiendra qu'un schéma peut diverger alors que la solution exacte converge, et que deux schémas également « raisonnables » peuvent créer ou détruire de l'énergie. On n'a étudié ici que la stabilité de la méthode ; la séance d'analyse asymptotique (TD 8) en mesurera l'erreur.