DOUBLE PENDULE

Le double pendule se caractérise par deux simples pendules: le second est suspendu au premier et le système est libre de se mouvoir dans un plan vertical.
Chacun des deux pendules est constitué d'une barre sans masse de longueur l et d'une masse m considérée comme ponctuelle attachée à son extrémité.

Pour l'intégration numérique, on considérera que: mgl=1 et mgl.l.=1

Dans un exercice ultérieur on démontrera que le mouvement d'un tel système est chaotique.
On demande:
- D'écrire les équations du mouvement
- D'écrire un code Matlab pour visualiser le mouvement du double pendule en temps réel.

Equation du mouvement

La masse au point A a comme coordonnée et vitesse: (l sin q, -l cos q) et (l dq/dt cos q, l dq/dt sin q)

De même pour la masse en B: (l.(sin q + sin f ), -l.(cos q + cos f))

et pour la vitesse:(l.(dq/dt cos q + df/dt cos f), l.( dq/dt sin q + df/dt sin f ) )

L'énergie cinétique s'écrit: 2T =

et l'énergie potentielle s'écrit: V= - mgl cos q - m g l (cos q + cos f )

La fonction de Lagrange L s'écrit: L = + mgl cos q + m g l (cos q + cos f )

Les équations du mouvement s'écrivent:

Ce système admet comme intégrale première T+V = E, ce qui peut servir à tester numériquement la précision de l'intégration par Runge-Kutta. Ceci est d'autant plus intéressant que la nature chaotique du mouvement peut entraîner des erreurs d'intégration conséquentes.

Code Matlab

Le Code Matlab donnant en temps réel le mouvement du double-pendule est donnée ci-après.

On prendra l=10 m

hold on;
echo off
options = odeset('RelTol',1e-8);
[t,teta] = ode45('doupen',[0:0.1:100],[1 1 10 0],options);

% Cette ligne procède au calcul du mouvement entre t=0 s et t=100 s par pas
% de 0.1 s.

% On a 1000 points d'intégration. A la fin de l'exécution de ce code, on
% obtient donc la position du pendule après 100 s. Il suffit de prendre cette
% position comme nouvelles conditions initiales et enserrer ce code dans une
% boucle plus grande pour ne pas spécifiquement limiter la simulation temps
% réel à 100 s. Si on opte pour cette stratégie, il y a intérêt à ne pas
% prendre 100 s mais quelques secondes car cette ligne demande tout de
% même une minute d'exécution, ce qui génèrerait des à-coups dans l'exécution
% du code, malvenus dans une simulation dite temps réel.

y=zeros(3,1);
z=zeros(3,1);
yold=zeros(3,1);
zold=zeros(3,1);
y(2)=10*sin(teta(1,1));
z(2)=-10*cos(teta(1,1));
y(3)=y(2)+10*sin(teta(1,3));
z(3)=z(2)-10*cos(teta(1,3));
axis([-30 30 -30 30])
line(y,z,'LineWidth',2);
line(yold,zold,'LineWidth',2);
p=plot(y,z,'EraseMode','none');
q=plot(yold,zold,'w','EraseMode','none');

% la boucle qui dessine la position du double pendule pour les 1000 points
% d'intégration calculés ci-dessus.

for i=1:1000
yold(2)=y(2);
zold(2)=z(2);
y(2)=10*sin(teta(i,1));
z(2)=-10*cos(teta(i,1));
yold(3)=y(3);
zold(3)=z(3);
y(3)=y(2)+10*sin(teta(i,3));
z(3)=z(2)-10*cos(teta(i,3));
tic

% tic enclenche une horloge interne

set(q,'XData',yold,'YData',zold,'LineWidth',2);
set(p,'XData',y,'YData',z,'LineWidth',2);

% Le graphe est considéré comme un object au sens programmation orientée
% objet. Les deux lignes qi précèdent remettent à jour les propriétés de
% l'objet qui ne sont rien d'autre que la position du double pendule.

while toc<0.1;end;

% cette ligne provoue une pause de 0.1 s. dans l'exécution du code.
% comme chaque point d'intégration est espacé de 0.1 s. cette ligne assure
% que chaque nouvelle position du double pendule est précisément affichée
% toutes les 0.1 s., ce qui nous garantit la simulation temps réel.

drawnow;
end
hold off
echo on

Le code doupen.m est donnée par:

function ydot = doupen(t,y)
ydot = zeros(4,1);
s = sin(y(1) - y(3));
c = cos(y(1) - y(3));
ydot(1) = y(2);
ydot(3) = y(4);
ydot(2) = (- y(4)*y(4)*s - 2*sin(y(1))- y(2)*y(2)*s*c + sin(y(3))*c)/(2 - c*c);
ydot(4) = (2*y(2)*y(2)*s - 2*sin(y(3))+ y(4)*y(4)*s*c + 2*sin(y(1))*c)/(2 - c*c);
end;

--o0o--