TUTORIAL SUITE

Faculté des

Sciences

Appliquées

Mécanique Rationnelle

2e bachelier 2008-2009

TUTORIAL

Séminaire du 15.10.2008

On demande d’écrire l’équation du mouvement d’un double pendule (voir diagramme ci-après), de le modéliser par MATLAB et de simuler son mouvement pour différentes conditions initiales. On tracera la plan des phases et le section de Poincaré (la section de Poincaré consiste à ne prendre que des points bien choisis de ce plan des phases ; un exemple est de prendre des points espaces à intervalles réguliers).

Equation du mouvement

Modélisation par Matlab

Visualisation du mouvement pour E=-1.5

Visualisation du mouvement pour E=-2.5 (chaos)

Section Poincaré

Vers le TP2

Equation du mouvement : par Lagrange

¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

Les équations de Lagrange s’écrivent :

avec w²=g/l

Ce système est à énergie constante :

Le plafond où est fixé le double-pendule est la référence E=0.Souvent, la valeur de l’énergie est un critère pour l’apparition du chaos. De fait, c’est ce qu’on observe dans le cas du mouvement du double pendule. Pour E=-1.5, le chaos apparaît presque pour toutes les conditions initiales excepté quelques ilôts de stabilité. Pour E=-2.5, il n’y a pas de chaos. Si on se fixe l’énergie du systèmes, on ne peut plus fixer arbitrairement les 4 conditions initiales q (t=0) , f(t=0), dq/dt(t=0) et df/dt(t=0),. On doit utiliser la relation :

Pour la modélisation MATLAB, on pose y(1)=q, y(2)=dq/dt, y(3)=f, y(4)=df/dt

Retour haut de la page

Modélisation par MATLAB

¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

clear

%mov = avifile('c:\test12.avi','Compression','None','Quality',100,'fps',25)

figure('NumberTitle','on','Name','TUTORIAL Matlab 15/10/2008 - DOUBLE PENDULE','Renderer','OpenGL','Color','w','Position',[100 100 500 500])

% DEBUT CALCUL

E=-2.5

teta_dot0=1;

teta0=0;

phi0=teta0;

phi_dot0=-teta_dot0+sqrt(6*cos(teta0)-teta_dot0^2+2*E);

options = odeset('RelTol',1e-8);

[t,teta] = ode45('doupen',[0:0.05:500],[teta0 teta_dot0 phi0 phi_dot0],options);

% FIN CALCUL

% DEBUT DU DESSIN

y=zeros(3,1);z=zeros(3,1);

yold=zeros(3,1);zold=zeros(3,1);

for i=1:max(size(t))

tic

y(2)=10*sin(teta(i,1));z(2)=-10*cos(teta(i,1));

y(3)=y(2)+10*sin(teta(i,3));z(3)=z(2)-10*cos(teta(i,3));

subplot(1,1,1,'replace');grid on;box on;axis([-20 20 -20 20 ]);xlabel('m');ylabel('m');Title(['Valeur de l Energie =',num2str(E)]);

line([y(1) y(2)],[z(1) z(2)],'Color','r','LineWidth',3);

line(y(2),z(2),'Color','r','Marker','.','Markersize',50);

line([y(2) y(3)],[z(2) z(3)],'Color','b','LineWidth',3);

line(y(3),z(3),'Color','b','Marker','.','Markersize',50);

drawnow;

%F = getframe(gcf);

%mov = addframe(mov,F);

while toc<0.05;end;

end

mov=close(mov);

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;

Retour haut de la page

Simulation du mouvement pour E=-2.5 (aucun chaos - mouvement périodique)

¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

IMAGE LONGUE A CHARGER - PATIENCE
IMAGE LONGUE A CHARGER - PATIENCE
IMAGE LONGUE A CHARGER - PATIENCE
IMAGE LONGUE A CHARGER - PATIENCE
IMAGE LONGUE A CHARGER - PATIENCE
IMAGE LONGUE A CHARGER - PATIENCE
IMAGE LONGUE A CHARGER - PATIENCE
IMAGE LONGUE A CHARGER - PATIENCE
IMAGE LONGUE A CHARGER - PATIENCE
IMAGE LONGUE A CHARGER - PATIENCE

Simulation du mouvement pour E=-1.5 (mouvement chaotique - aucune périodicité

¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

Retour haut de la page

La sensibilité aux conditions initiales (condition nécessaire au chaos)

¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

Deux pendules doubles (bleu et rouge) dont las conditions initiales sont identiques à l’exception de dq/dt qui varie de 0.001 (dq/dt(t=0)=1 et dq/dt(t=0)=0.999). Le mouvement des deux double-pendules reste confondu.Le mouvement n'est d'ailleurs PAS chaotique

Deux pendules doubles (bleu et rouge) dont las conditions initiales sont identiques à l’exception de dq/dt qui varie de 0.001 (dq/dt(t=0)=1 et dq/dt(t=0)=0.999). Au bout d’un temps fini, le mouvement des deux double-pendules devient distinct. C’est la sensibilité aux conditions initiales (aussi connu comme l’effet papillon). Le mouvement des deux double-pendules est chaotique.

clear

figure('NumberTitle','on','Name','TUTORIAL Matlab semaine 15/10/2008 - DOUBLE PENDULE','Renderer','OpenGL','Color','w','Position',[100 100 500 500])

% DEBUT CALCUL

E=-1.5

teta_dot0_1=1;

teta_dot0_2=0.999;

teta0=0;

phi0=teta0;

phi_dot0=-teta_dot0_1+sqrt(6*cos(teta0)-teta_dot0_1^2+2*E);

phi_dot0_2=-teta_dot0_2+sqrt(6*cos(teta0)-teta_dot0_2^2+2*E);

options = odeset('RelTol',1e-8);

[t,teta] = ode45('doupen',[0:0.15:120],[teta0 teta_dot0_1 phi0 phi_dot0],options);

[t,teta2] = ode45('doupen',[0:0.15:120],[teta0 teta_dot0_2 phi0 phi_dot0_2],options);

% FIN CALCUL

% DEBUT DU DESSIN

y=zeros(3,1);z=zeros(3,1);

yy=zeros(3,1);zz=zeros(3,1);

for i=2:max(size(t))

tic

y(2)=10*sin(teta(i,1));z(2)=-10*cos(teta(i,1));

y(3)=y(2)+10*sin(teta(i,3));z(3)=z(2)-10*cos(teta(i,3));

yy(2)=10*sin(teta2(i,1));zz(2)=-10*cos(teta2(i,1));

yy(3)=yy(2)+10*sin(teta2(i,3));zz(3)=zz(2)-10*cos(teta2(i,3));

subplot(2,1,1,'replace');grid on;box on;axis([-20 20 -20 0 ]);xlabel('m');ylabel('m');Title(['Valeur de l Energie =',num2str(E)]);

line([y(1) y(2)],[z(1) z(2)],'Color','b','LineWidth',3);

line(y(2),z(2),'Color','b','Marker','.','Markersize',50);

line([y(2) y(3)],[z(2) z(3)],'Color','b','LineWidth',3);

line(y(3),z(3),'Color','b','Marker','.','Markersize',50);

line([yy(1) yy(2)],[zz(1) zz(2)],'Color','r','LineWidth',3);

line(yy(2),zz(2),'Color','r','Marker','.','Markersize',50);

line([yy(2) yy(3)],[zz(2) zz(3)],'Color','r','LineWidth',3);

line(yy(3),zz(3),'Color','r','Marker','.','Markersize',50);

drawnow;

subplot(2,1,2);grid on;box on;xlabel('t');ylabel('\theta et \phi');Title(['Valeur de l Energie =',num2str(E)]);

line([t(i-1) t(i)],[teta(i-1,1) teta(i,1)],'Color','r')

line([t(i-1) t(i)],[teta2(i-1,1) teta2(i,1)],'Color','b')

while toc<0.05;end;

end

Section de Poincaré

¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

La section de Poincaré permet de visualiser s’il y a chaos ou non et les zones de stabilité: le code MATLAB suivant n’est pas OPTIMISE. Il dessine la section de Poincaré sans faire appel à events pour prendre la valeur de la solution lorsque le critère de Poincaré est atteint (ceci a pour but d’illustrer les notations matricielles de MATLAB et leur concision).

Une belle section de Poincaré est donnée par

q - f =0

c’est-à-dire on trace dans le plan des phases les valeurs de q et dq/dt lorsque cette condition est respectée

Pour E=-1.5, les sections de Poincaré correspondant à chaque CI seront des courbes non continues (fractales) avec des ilôts de stabilité (courbes continues pour certaines CI ou Conditions Initiales) ; pour E=-2.5, les sections de Poincaré correspondant à chaque CI seront des courbes continues.

echo off

options = odeset('RelTol',1e-8);

figure('NumberTitle','on','Name','TUTORIAL Matlab 15/10/2008 - DOUBLE PENDULE','Renderer','OpenGL','Color','w','Position',[100 100 500 500])

hold on

grid on

box on

E=-1

xlabel('cos(\theta)');ylabel('d\theta/dt');Title(['Plan des Phases - Valeur de l Energie =',num2str(E)]);

i=0;

for teta_dot0=1:-0.1:-1

i=i+1;

teta0=0;

phi0=teta0;

phi_dot0=-teta_dot0+sqrt(6*cos(teta0)-teta_dot0^2+2*E);

[t,teta] = ode45('doupen',[0 3600],[teta0 teta_dot0 phi0 phi_dot0],options);

indices=find(abs(sin(teta(:,1)-teta(:,3)))<=1e-3);

plot(cos(teta(indices,1)),teta(indices,2),'+','MarkerSize',3)

drawnow;

F = getframe(gcf);

X = frame2im(F);

imwrite(X,['R:\pendule',num2str(i),'.png'],'png')

drawnow;

hold on;

end

E=-2.5 (non chaotique)

E=-1.5 (chaotique avec ilôts de stabilité selon les conditions initiales)

E=-1.5 (chaotique quasi aucun ilôt de stabilité selon les conditions initiales)

Retour haut de la page

Depuis le 15.10.2008, vous êtes le visiteur