Étude du pendule amorti
Modélisation
On considère le système mécanique du pendule rigide amorti. On désigne par \(\theta\) l'angle entre la verticale passant par le point \(O\) de suspension et la direction de la barre rigide du pendule. La mécanique du point nous donne l'équation différentielle qui modélise l'évolution de l'angle \(\theta\) : \[\theta'' = -b\theta'-c\sin\theta\] où \(b\) et \(c\) sont des constantes strictement positives.
En augmentant la dimension du système, on se ramène à une équation différentielle du premier ordre: \[\begin{pmatrix}\theta\\\theta'\end{pmatrix}'=\begin{pmatrix}\theta'\\-b\theta'-c\sin\theta\end{pmatrix}=F(\theta,\theta'),\] où \(F\) est une fonction \(C^\infty\).
Nous considérerons ici le cas où \(b^2>4c\), c'est-à-dire le cas où les frottements sont assez petits. Dans nos exemples, on prendra les valeurs \(b=0.5\) et \(c=1\).
Les équilibres
Cherchons les équilibres du système, c'est-à-dire les zéros du champs de vecteur \(F\).
Comme \(\theta\) est un angle, nous allons nous restreindre à \(\theta\in]-\pi,\pi]\). Avec, ce domaine de définition, on a les deux équilibres \((0,0)\) et \((\pi,0)\).
Solution numérique et portrait de phase
Si on ne peut pas trouver de solution analytique à ce problème, on peut en calculer des solutions numériquement. Pour cela, on peut utiliser le langage Python, et ses bibliothèques numpy et scipy.
from scipy.integrate import odeint
import numpy as np
def pend(y, t, b, c):
theta, omega = y
dydt = [omega, -b*omega - c*np.sin(theta)]
return dydt
b=.5
c=1.0
t = np.linspace(0, 25, 201)
y0 = [np.pi-0.01, 0.0]
sol = odeint(pend, y0, t, args=(b, c))
np.savetxt("solution0.dat",sol, '%2.5f')
np.savetxt("solution0tps.dat",np.c_[t,sol[:,0]], '%2.5f')
On peut alors représenter la solution calculée.
Ici nous avons pris une condition initial proche d'un des deux équilibres. On observe que le système converge vers l'autre équilibre.
L'animation ci-dessus montre l'évolution du système physique (le pendule rigide amorti), le tracé de la fonction du temps \(\theta(\cdot)\), et enfin le portrait de phase de la solution dans le plan \((\theta,\theta')\).
Étude du linéarisé autour des équilibres
L'étude du système linéarisé autour des équilibres va nous permettre de comprendre ce qu'il se passe autour de ces équilibres.
Autour de \((0,0)\)
L'étude du système linéarisé autour d'un point d'équilibre \((\theta_0,\theta_0')\) est l'étude de l'équation différentielle linéaire autonome suivante : \[\begin{pmatrix}\theta\\\theta'\end{pmatrix}'=DF(\theta_0,\theta_0')\cdot\begin{pmatrix}\theta\\\theta'\end{pmatrix}.\] Puisque \(F\) est \(C^1\), le calcul de la différentielle peut se faire par le calcul des dérivées partielle. En \((0,0)\),cela donne : \[\begin{pmatrix}\theta\\\theta'\end{pmatrix}'=\begin{pmatrix}0&1\\-c&-b\end{pmatrix}\begin{pmatrix}\theta\\\theta'\end{pmatrix}=:A_0\begin{pmatrix}\theta\\\theta'\end{pmatrix}.\]
Lorsqu'on se place dans le cas où \(b^2-4c>0\), on a deux valeurs propres complexes conjuguées : \[\mu_\pm=\frac{-b\pm\mathrm{i}\sqrt{4c-b^2}}{2}.\]
En considérant \(w=a+\mathrm{i}b\in\mathbb{C}^2\) un vecteur propre non nul de \(A_0\) associé à \(\mu_+\), la base de \(\mathbb{R}^2\), \((a,b)\) et en constituant la matrice de passage \(P\) constituée des vecteurs \(a\) et \(b\) en colonne. Dans cette base, il est simple de montrer que \[P^{-1}\mathrm{e}^{tA_0}P=\begin{pmatrix}Re(\mathrm{e}^{t\mu_+})&Im(\mathrm{e}^{t\mu_+})\\-Im(\mathrm{e}^{t\mu_+})&Re(\mathrm{e}^{t\mu_+})\end{pmatrix}.\]
Ainsi, pour \(Y(t)=P^{-1}\begin{pmatrix}\theta\\\theta'\end{pmatrix}\), on obtient : \[Y(t)=\mathrm{e}^{-\frac{b}{2}t}\begin{pmatrix}\cos(\sqrt{4c-b^2}/2)&\sin(\sqrt{4c-b^2}/2)\\-\sin(\sqrt{4c-b^2}/2)&\cos(\sqrt{4c-b^2}/2)\end{pmatrix}Y(0).\]
À un changement de base près, le portrait de phase du système linéarisé autour du point d'équilibre \(0,0\) est donc :
Il s'agit bien là, d'un comportement similaire à celui autour de point d'équilibre dans le système non linéaire.
Autour de \((\pi,0)\)
Le système linéarisé autour du point d'équilibre \((\pi,0)\) est: \[\begin{pmatrix}\theta\\\theta'\end{pmatrix}'=\begin{pmatrix}0&1\\c&-b\end{pmatrix}\begin{pmatrix}\theta\\\theta'\end{pmatrix}=:A_\pi\begin{pmatrix}\theta\\\theta'\end{pmatrix}.\]
La matrice \(A_\pi\) a pour valeurs propres : \[\lambda_\pm=\frac{-b\pm\sqrt{\Delta}}{2},\] où \(\Delta = b^2+4c\) qui est strictement positif pour \(b\) et \(c\) strictement positif. On a que \(\lambda_+>0\) et \(\lambda_-<0\), on a donc un espace stable et un espace instable (généré par les vecteurs propres associés). Le portrait de phase est alors le suivant:
On voit donc qu'à part pour les conditions initiales sur l'espace stable, toutes les autres solutions diverge de l'équilibre. Ceci correspond bien à l'intuition donnée par le système physique dont l'équilibre «du haut» semble «très instable».
Notons tout de même que seul le côté des \(\theta\leq 0\) nous intéresse ici.
Conclusion
Grâce à l'étude des systèmes linéarisés autour des points d'équilibre, on comprend mieux le comportement du système non linéaire du départ.
Pour terminer, prenons une conditions initiale donc la vitesse initiale est positive pour que la solution se rapproche de l'équilibre \((\pi,0)\) pour ensuite converger vers l'équilibre stable \((0,0)\). Pour cela, nous prenons, par exemple, comme condition initiale \(\Theta = (\pi,0)+v_-\), où \(v_-\) est un vecteur propre instable de \(A_\pi\). Il faut prendre un \(v_-\) dont la coordonnée en \(\theta\) soit négative.
Si la vitesse initiale est trop élevée, alors le pendule dépasse la position verticale haute (\(\theta>\pi\)), pour converger vers \((2\pi,0)\) qui est le même point physique que \((0,0)\).
Bien sûr, le fait de prendre comme condition initiale \(\Theta = (\pi,0)+v_-\) ne garantit pas que la vitesse initiale ne soit pas suffisemment élevée pour que \(\theta(t)\) ne dépasse pas \(\pi\). Cela nous donne tout le même un ordre d'idée sur les conditions initiales proches du point d'équilibre qui vont permettre de se rapprocher (en la dépassant ou non) de la position d'équilibre verticale haute. D'ailleurs, ici, pour que le système reste dans l'espace \(\{\theta(t)\leq \pi\}\), il a fallu diminuer la valeur de la vitesse initiale donnée par le vecteur propre \(v_-\).
