Cette page décrit une simulation en HTML5 Canvas modélisant le mouvement d'un pendule triple amorti. Nous obtiendrons les équations du mouvement grâce aux équations d'Euler-Lagrange puis implémenterons la méthode numérique
RK4 à pas variable (Runge Kutta ordre 4).
Le système différentiel est obtenu grâce aux équations d'Euler-Lagrange. En notant $L$ le Lagrangien du système, $q_{i}$ chacune des 3 coordonnées et $\dot{q_{i}}$ leur dérivée par rapport au temps, on peut écrire :
\begin{equation}
\dfrac{\text{d}}{\text{d}t}\bigg(\dfrac{\partial\,L}{\partial\,\dot{q_{i}}}\bigg)-\dfrac{\partial\,L}{\partial\,q_{i}}=0
\end{equation}
Cette équation est valable lorsque le système est conservatif (le travail des forces non-conservatives est nul). Quand il existe des forces non conservatives dont le travail est non nul, on les qualifie de forces d'amortissement. Il faut rajouter alors au membre de droite la force de frottement $\overrightarrow{F_{f}}=-\overrightarrow{\nabla_{v}}\,D_{r}$ avec $D_{r}$ la fonction de dissipation de Rayleigh. On obtient :
\begin{equation}
\dfrac{\text{d}}{\text{d}t}\bigg(\dfrac{\partial\,L}{\partial\,\dot{q_{i}}}\bigg)-\dfrac{\partial\,L}{\partial\,q_{i}}=-\dfrac{\partial\,D_{r}}{\partial\,\dot{q_{i}}}
\end{equation}
avec $v_{i}^{2}=\left\Vert\dfrac{\text{d}\overrightarrow{OM_{i}}}{\text{d}t}\right\Vert^{2}=\dfrac{\text{d}\overrightarrow{OM_{i}}}{\text{d}t}\cdot\dfrac{\text{d}\overrightarrow{OM_{i}}}{\text{d}t}$ et $h_{i}$ qui correspond à la hauteur par rapport au sol de la masse $i$.
Calcul de l'énergie cinétique - Les termes $v_{i}^{2}$ sont égaux à :
où $k_{i}$ correspond au coefficient d'amortissement selon l'angle $\theta_{i}$ et $v_{i}^{2}$ sont les vitesses calculées en \eqref{eq-vi2}. Le terme de droite dans l'équation d'Euler-Lagrange est alors donné par :
Le calcul du système d'équations peut se faire à la main mais il est relativement lourd. Une façon plus rapide consiste à utiliser la Symbolic Math Toolbox. Après avoir calculé l'expression symbolique du Lagrangien, on pose le système différentiel de manière automatique grâce au code Matlab suivant : compute_diff_system.m.
Nous décrivons brièvement ce code :
En partant des expressions de \eqref{eq-kinetic-potential} et de $v_{i}^{2}$, avec les variables $t_{i}=\theta_{i}$, $dt_{i}=\dot{\theta_{i}}$, on obtient le Lagrangien suivant :
% Kinetic Energy
T=1/2*l1*l1*dt1*dt1*(m1+m2+m3)...+1/2*l2*l2*dt2*dt2*(m2+m3)...+1/2*l3*l3*dt3*dt3*m3...+l1*l2*dt1*dt2*cos(t1-t2)*(m2+m3)...+l2*l3*dt2*dt3*cos(t2-t3)*m3...+l1*l3*dt1*dt3*cos(t1-t3)*m3;% Potential Energy
V=g*(l1+l2+l3)*(m1+m2+m3)...-g*l1*(m1+m2+m3)*cos(t1)...-g*l2*(m2+m3)*cos(t2)...-g*l3*m3*cos(t3);% Lagrangian
L=T-V;
Ecrivons maintenant l'expression de la fonction de dissipation :
% Rayleigh Dissipation function
Dr=1/2*l1*l1*dt1*dt1*(k1+k2+k3)...+1/2*l2*l2*dt2*dt2*(k2+k3)...+1/2*l3*l3*dt3*dt3*k3...+l1*l2*dt1*dt2*cos(t1-t2)*(k2+k3)...+l2*l3*dt2*dt3*cos(t2-t3)*k3...+l1*l3*dt1*dt3*cos(t1-t3)*k3;% Derivative of Rayleigh function
f1=diff(Dr,dt1);
f2=diff(Dr,dt2);
f3=diff(Dr,dt3);
Il faut maintenant établir le système différentiel à résoudre. Ceci se fait de la manière suivante :
% Get Lagrange equations and solve for expressing theta_second
Equations=Lagrange(L,[t1,dt1,ddt1,t2,dt2,ddt2,t3,dt3,ddt3]);[theta1_second theta2_second theta3_second]=solve(Equations(1)==-f1,...
Equations(2)==-f2,Equations(3)==-f3,ddt1,ddt2,ddt3);
On obtient finalement les 3 équations dynamiques du système.
Remarque :
Le langage Javascript ne reconnaît pas les exposants avec le caractère '^'. Pour pouvoir par exemple convertir l'expression 'm1^2' sous Matlab dans une syntaxe Javascript valide, c'est-à-dire 'm1*m1', nous pouvons utiliser la commande sed suivante :
avec les conditions initiales $t_{0}$ et $y_{0}$ connues. La formule itérative est la suivante :
\begin{equation}
y_{i+1}=y_{i}+\dfrac{(k_{1}+2k_{2}+2k_{3}+k_{4})}{6}
\end{equation}
avec
\begin{equation}
k_{1}=h\Psi(t_{i},y_{i})\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,k_{2}=h\Psi(t_{i}+\dfrac{h}{2},y_{i}+\dfrac{k_{1}}{2})\\
k_{3}=h\Psi(t_{i}+\dfrac{h}{2},y_{i}+\dfrac{k_{2}}{2})\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,k_{4}=h\Psi(t_{i}+h,y_{i}+k_{3})\\
\end{equation}
Pour un système d'équations faisant intervenir des dérivées seconde, on utilise une approche vectorielle; Dans notre cas, on a le système différentiel suivant :
En écrivant $u=\dfrac{\text{d}x}{\text{d}t}$, $v=\dfrac{\text{d}y}{\text{d}t}$, $w=\dfrac{\text{d}z}{\text{d}t}$, et en posant le vecteur $\overrightarrow{V}$ :
\begin{equation}
\overrightarrow{V}=\begin{bmatrix} x \\ u \\ y \\ v \\ z \\ w\end{bmatrix}
\end{equation}
On peut mettre le système différentiel sous une forme vectorielle équivalente à \eqref{eq-rk4-basic} :
\begin{equation}
\dfrac{\text{d}\overrightarrow{V}}{\text{d}t}=\overrightarrow{U}(\overrightarrow{V})=\begin{bmatrix} u(x,x',y,y',z,z') \\ F(x,x',y,y',z,z') \\ v(x,x',y,y',z,z') \\ G(x,x',y,y',z,z') \\ w(x,x',y,y',z,z') \\ H(x,x',y,y',z,z')\end{bmatrix}
\end{equation}
Les coefficients de Runge-Kutta se calculent alors de manière identique :
\begin{equation}
\overrightarrow{k_{1}}=h\overrightarrow{U}(\overrightarrow{V_{n}})\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\overrightarrow{k_{2}}=h\overrightarrow{U}(\overrightarrow{V_{n}}+\dfrac{\overrightarrow{k_{1}}}{2})\\
\overrightarrow{k_{3}}=h\overrightarrow{U}(\overrightarrow{V_{n}}+\dfrac{\overrightarrow{k_{2}}}{2})\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\overrightarrow{k_{4}}=h\overrightarrow{U}(\overrightarrow{V_{n}}+\overrightarrow{k_{3}})\\
\end{equation}
Voici un canvas permettant de simuler le mouvement d'un pendule triple. La longueur totale des axes est fixée à 3 mètres, la masse initiale de chaque masse prise à 1 kg et la vitesse de départ à $5\,\text{rad.s}^{-1}$. Vous pouvez faire un drag'n drop sur chaque masse du pendule pour leur affecter la position initiale désirée.
Le source du Canvas HTML5 est disponible sur ce lien triple_pendulum.js. On commence par initialiser l'interface
graphique puis on définit les fonctions correspondant au système différentiel. La méthode Runge-Kutta ordre 4 à pas variable est alors implémentée pour calculer les prochaines valeurs. Une fois ces valeurs obtenues, nous mettons à jour l'affichage du pendule.
Il se peut que selon les choix de la longueur des axes, ou une différence de rapports importante entre les 3 masses du pendule, les angles $\theta_{i}$ soient amener à varier très rapidement. Il faut dans ces cas là choisir un pas de temps plus petit que celui utilisé par défaut dans la simulation ($h=0.005$) afin d'éviter la divergence des valeurs calculées. On utilise alors la méthode RK4 à pas variable aussi appelée RK4 à pas adaptatif. Le principe est d'obtenir un pas de temps qui s'adapte en fonction de la variabilité des valeurs calculées : le pas de temps doit diminuer quand celles-ci changent brutalement et augmenter quand elles varient lentement.
Pour évaluer cette variabilité, on fixe une erreur minimale au delà de laquelle il faut diminuer le pas de temps. Pour cela, on compare les valeurs calculées sur le pas de temps courant $h$ avec les valeurs obtenues sur 2 itérations en prenant $h/2$. Si la différence est supérieure à l'erreur fixée, nous prenons un pas plus petit. Dans le cas contraire, nous l'augmentons.
