Le double puits, en Python

1) L'adimensionnement. La force est $-E_{\mathrm{p}}^{\prime }(x)=-4a{x}^3+2bx$, d'où $m\ddot{x}=2bx-4a{x}^3$. Avec $x=x_{e}u$ et $t=\tau s$, on a $\ddot{x}=\frac{x_e}{{\tau }^2}{u}^{\prime \prime }$ (dériver deux fois par rapport à $t$, c'est dériver deux fois par rapport à $s$ et diviser par ${\tau }^2$). L'équation devient $$\frac{m}{{\tau }^2}{u}^{\prime \prime }=2bu-4a{x}_e^2{u}^3.$$ Or $\frac{m}{{\tau }^2}=2b$ et $4a{x}_e^2=4a\cdot \frac{b}{2a}=2b$ : en divisant par $2b$, $u^{\prime \prime }=u-u^3$. Les deux échelles ont été choisies pour cela : $x_e$ est la position d'un puits, $\tau$ fixe l'échelle de temps. Multiplier par $u^{\prime }$ donne $u^{\prime }{u}^{\prime \prime }-u{u}^{\prime }+u^3{u}^{\prime }=0$, c'est-à-dire $\frac{\mathrm{d}e}{\mathrm{d}s}=0$ : $e$ se conserve. Et l'on vérifie que $E_{\mathrm{m}}=\frac{1}{2}m{\dot{x}}^2+a{x}^4-b{x}^2=\frac{b^2}{a}e$ (chaque terme se calcule avec $m{x}_e^2/{\tau }^2=2b{x}_e^2=b^2/a$, $a{x}_e^4=b^2/(4a)$ et $b{x}_e^2=b^2/(2a)$) : $e$ est l'énergie mécanique mesurée dans l'unité naturelle $b^2/a$.

Pourquoi on adimensionne. Le problème dépendait de $m$, $a$, $b$ et d'une condition initiale ; il ne dépend plus que d'un seul nombre, $e$.

2) Les mouvements possibles. $e_p^{\prime }(u)=-u+u^3=u(u^2-1)$ : équilibres en $u=0$ (maximum local, $e_p=0$, instable) et en $u=\pm 1$ (minima, $e_p=-\frac{1}{4}$, stables). Allure : un W, deux puits de profondeur $\frac{1}{4}$ séparés par une bosse en $0$, remontant vers $+\infty$ aux deux bouts.

• Si $-\frac{1}{4}<e<0$ : le point est piégé dans le puits où il a été lancé, entre deux rebroussements.
• Si $e>0$ : il oscille entre deux rebroussements symétriques en passant par-dessus la barrière.
• Si $e=0$ : c'est la séparatrice ; le point met un temps infini à atteindre le sommet $u=0$.

Petites oscillations. $e_p^{\prime \prime }(u)=-1+3u^2$, donc $e_p^{\prime \prime }(\pm 1)=2$ : la pulsation réduite est $\sqrt{2}$, et la période réduite $2\pi /\sqrt{2}\approx 4,44$. En unités réelles, $\omega =\sqrt{2}/\tau =2\sqrt{b/m}$, valeur que l'on obtient aussi directement par $\sqrt{E_{\mathrm{p}}^{\prime \prime }(x_e)/m}$.

3) Le script. On pose $w=u^{\prime }$ : le système d'ordre 1 est $u^{\prime }=w$, $w^{\prime }=u-u^3$. Le point de départ. Partir sans vitesse de $u_0$ donne $e=-\frac{1}{2}{u}_0^2+\frac{1}{4}{u}_0^4$, soit $u_0^4-2u_0^2-4e=0$, du second degré en $u_0^2$ : $u_0^2=1\pm \sqrt{1+4e}$. La racine $1+\sqrt{1+4e}$ est positive pour tout $e>-\frac{1}{4}$ : $u_0=\sqrt{1+\sqrt{1+4e}}$ est le rebroussement extérieur du mouvement d'énergie $e$, du côté des $u>0$.

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def champ(Y, s):          # u'' = u - u^3 sous forme d'ordre 1
    u, w = Y              # w = u'
    return [w, u - u**3]

def depart(e):            # lacher au rebroussement exterieur
    return [np.sqrt(1 + np.sqrt(1 + 4*e)), 0.0]

s = np.linspace(0, 200, 400001)

# 1) portrait de phase et separatrice
for e in (-0.2, -0.01, 0.01, 0.5):
    Y = odeint(champ, depart(e), s[:40001], rtol=1e-10, atol=1e-12)
    plt.plot(Y[:, 0], Y[:, 1], label="e = %g" % e)
u = np.linspace(-np.sqrt(2), np.sqrt(2), 401)
r = np.sqrt(np.clip(1 - u**2/2, 0, None))
plt.plot(u, u*r, "k--"); plt.plot(u, -u*r, "k--")   # e = 0
plt.xlabel("u"); plt.ylabel("u'"); plt.legend(); plt.show()

# 2) periode = duree entre deux maxima successifs de u
def periode(e):
    u = odeint(champ, depart(e), s, rtol=1e-10, atol=1e-12)[:, 0]
    i = np.where((u[1:-1] > u[:-2]) & (u[1:-1] >= u[2:]))[0] + 1
    return s[i[1]] - s[i[0]]

for e in (-0.24, -0.1, -0.01, -0.001, 0.001, 0.5):
    print("e = %6g    periode reduite = %.2f" % (e, periode(e)))
print("2 pi / sqrt(2) = %.2f" % (2*np.pi/np.sqrt(2)))

La séparatrice s'obtient en écrivant $e=0$ : $u^{\prime 2}=u^2-\frac{1}{2}{u}^4$, soit $u^{\prime }=\pm u\sqrt{1-u^2/2}$ pour $|u|\le \sqrt{2}$ — une courbe en forme de huit couché qui passe par l'origine.

Le portrait de phase obtenu. Pour $e=-0,2$, une petite courbe fermée, presque une ellipse, autour de $(1,0)$ ; pour $e=-0,01$, une courbe qui épouse la boucle droite du huit ; pour $e=0,01$ et $0,5$, des courbes qui entourent les deux boucles. La séparatrice partage le plan de phase en trois régions : puits de droite, puits de gauche, au-dessus de la barrière.

4) Les périodes. L'exécution affiche :

e =  -0.24    periode reduite = 4.48
e =   -0.1    periode reduite = 5.24
e =  -0.01    periode reduite = 7.41
e = -0.001    periode reduite = 9.69
e =  0.001    periode reduite = 19.35
e =    0.5    periode reduite = 6.78
2 pi / sqrt(2) = 4.44

Le fond du puits. À $e=-0,24$, tout près du fond ($-\frac{1}{4}$), la période vaut $4,48$, à $1\%$ de la valeur harmonique $4,44$ : l'approximation du puits harmonique est confirmée, et sa limite de validité mesurée.

L'approche de la séparatrice. Quand $e\to 0^{-}$, la période croît sans limite : $5,24$, $7,41$, $9,69$. Elle augmente d'environ $2,3\approx \ln 10$ chaque fois que $|e|$ est divisé par dix : la divergence est logarithmique. Pourquoi. Le rebroussement intérieur s'approche du sommet $u=0$ de la barrière, où la force s'annule ; au voisinage de $u=0$, l'équation se réduit à $u^{\prime \prime }\simeq u$, dont les solutions s'écartent de $0$ comme ${\mathrm{e}}^s$ : pour s'éloigner d'une distance de l'ordre de $\sqrt{|e|}$ à une distance de l'ordre de $1$, il faut un temps de l'ordre de $\ln (1/\sqrt{|e|})$, qui diverge. Le point passe l'essentiel de sa période à « hésiter » au sommet de la bosse. Juste au-dessus ($e=0,001$), la période est environ double ($19,35$) : l'orbite fait le tour des deux puits et frôle la bosse deux fois par période. Loin au-dessus ($e=0,5$), la barrière ne compte plus guère et la période redescend ($6,78$).

L'effet du non linéaire. Pour un oscillateur harmonique, la période ne dépend pas de l'énergie (isochronisme). Ici elle en dépend, et elle peut même diverger : c'est le terme $-u^3$ qui crée le double puits, donc la bosse et la séparatrice. Aucune linéarisation ne peut rendre compte de ce comportement : c'est exactement pour cela qu'on résout numériquement.

Ce que l'exercice installe. Adimensionner avant de calculer : on réduit le nombre de paramètres et l'on fait apparaître les échelles naturelles. Mettre une équation d'ordre 2 sous forme d'un système d'ordre 1 pour la confier à un solveur. Lire un portrait de phase : courbes fermées (mouvements bornés), séparatrice. Et constater que, dans un potentiel non harmonique, la période dépend de l'amplitude — ce que le problème qui suit va quantifier pour le pendule.