L'énergie potentielle d'un oscillateur harmonique est de la forme Ep(x) = ½.K.x². Pour ce type d'oscillateur, le mouvement est sinusoïdal, le portrait de phase v = f(x) est une ellipse et la période est indépendante de la valeur de l'énergie totale du système.
On considère ici un oscillateur (par exemple une masse m fixée à un ressort spécial) pour lequel l'expression de l'énergie potentielle est Ep(x) = ½.K.x² − 1/3.A.K.x³.
A est un coefficient sans dimension qui va caractériser l'anharmonicité de l'oscillateur.
A ce potentiel correspond la force F(x) = − dEp/dx = − K.x + A.K.x².
Cette force est nulle pour x = 0 et pour x = 1 / A qui sont des positions d'équilibre.
L'expression de la dérivée seconde de Ep est d²Ep / dx² = − K + 2.A.K.x
Pour x = 0 elle est positive et cette position est un équilibre stable. Par contre pour x = 1 / A, la dérivée seconde est négative et l'équilibre instable. La valeur correspondante de Ep est Ep(1/A) = Epmax = K / 6A².
Si l'on néglige toute forme d'amortissement l'équation du mouvement est : m.d²x / dt² + F.x = 0 (1)
Il n'existe pas de solution analytique au problème et il faut faire une intégration numérique. On utilise ici la méthode de Runge-Kutta à l'ordre 4.
Pour cette intégration, il faut préciser les conditions initiales. On prend à l'instant t = 0, v = 0.
L'énergie totale de m étant E = Ep + Ec = Ep(x) + ½.m.v², en t = 0, on a E = Ep. Pour trouver la valeur initiale de x, il faut résoudre E = Ep(x) par une méthode de zéro.
Il est aussi possible de prendre x = 0 pour t = 0. Comme Ep(0) = 0, on détermine v(0) à partir de E = Ec = ½.m.v².
La partie droite de l'animation montre la courbe Ep(x) = ½.K.x² − 1/3.A.K.x³ (en rouge), la courbe Ep(x) = ½.K.x² (en gris), le mouvement de la masse m. Le trait orange correspond à l'énergie potentielle et le trait cyan à l'énergie cinétique.
Le trait vert horizontal correspond à la valeur de l'énergie totale.
En haut à gauche on trouve la courbe v = f(x) (portrait de phase) et en bas à droite la courbe x = f(t).
Le programme affiche aussi la valeur de la période T₀ = 2π / (M / K)^½ de l'oscillateur harmonique.

Utilisation
Dans le programme, on pose m = 1.
Il est possible de modifier avec les curseurs la valeur de l'énergie totale, le rapport K / M et la valeur de A.
Le programme limite automatiquement la valeur de E à Epmax.
Pour A = 0, vérifier que v = f(x) est une ellipse, que x = f(t) est une sinusoïde et que la période est indépendante de la valeur de E.
Pour A ≠ 0, vérifier que si le puits de potentiel est étroit, (K grand) pour les faibles valeurs de E, l'oscillateur est pratiquement harmonique.
Tester la dépendance de la période avec E.