Question

Project on a double pendulum using MATLAB. I have simulated double pendulum chaos using the code I found in Chegg. Could you please guide in adding a code to determine the time it takes for the second arm of the pendulum to "flip" based on initial starting conditions?

Answer 1

clc %---------------------------------------------------------------------------------------------------%

clear % function:Animates the comparison of two double pendulum's with almost similar initial parameters %

% parameters: %

% th1(1) : the angle of the top mass of the first double pendulum %

% th2(1) : the angle of the bottom mass of the first double pendulum %

dt=0.01; % thd1(1): the speed of the top mass of the first double pendulum %

T=10; % thd2(1): the speed of the bottom mass of the first double pendulum %

N=int16(T/dt); % al1(1) : the angle of the top mass of the second double pendulum %

th1=zeros(N,1); % al2(1) : the angle of the bottom mass of the second double pendulum %

al1=zeros; % ald1(1): the speed of the top mass of the second double pendulum %

th2=zeros(N,1); % ald2(1): the speed of the bottom mass of the second double pendulum %

al2=zeros(N,1); % ___________________________________________________ %

thd1=zeros(N,1);

ald1=zeros(N,1);

thd2=zeros(N,1);

ald2=zeros(N,1);

thdd1=zeros(N,1);

aldd1=zeros(N,1);

thdd2=zeros(N,1);

aldd2=zeros(N,1);

X1=zeros(N,1);

Y1=zeros(N,1);

Z1=zeros(N,1);

D1=zeros(N,1);

X2=zeros(N,1);

Y2=zeros(N,1);

Z2=zeros(N,1);

D2=zeros(N,1);

figure;

hold off;

th1(1)=pi/2;

al1(1)=pi/2;

th2(1)=pi+0.02; %---------------------------------------------------------------------------------------------------%

al2(1)=pi;

thd1(1)=0;

ald1(1)=0;

thd2(1)=0;

ald2(1)=0;

g=10;

A=-2*g*sin(th1(1))-sin(th1(1)-th2(1))*thd2(1)^2;

J=-2*g*sin(al1(1))-sin(al1(1)-al2(1))*ald2(1)^2;

B=-g*sin(th2(1))+sin(th1(1)-th2(1))*thd1(1)^2;

K=-g*sin(al2(1))+sin(al1(1)-al2(1))*ald1(1)^2;

thdd1(1)=(A-B*cos(th1(1)-th2(1)))/(2-(cos(th1(1)-th2(1)))^2);

aldd1(1)=(A-B*cos(al1(1)-al2(1)))/(2-(cos(al1(1)-al2(1)))^2);

thdd2(1)=B-cos(th1(1)-th2(1))*thdd1(1);

aldd2(1)=B-cos(al1(1)-al2(1))*aldd1(1);

X1(1)=sin(th1(1));

Y1(1)=cos(th1(1));

D1(1)=sin(al1(1));

Z1(1)=cos(al1(1));

X2(1)=sin(th1(1))+sin(th2(1));

Y2(1)=cos(th1(1))+cos(th2(1));

D2(1)=sin(al1(1))+sin(al2(1));

Z2(1)=cos(al1(1))+cos(al2(1));

for i=2:N

thd1(i)=thd1(i-1)+dt*thdd1(i-1);

ald1(i)=ald1(i-1)+dt*aldd1(i-1);

thd2(i)=thd2(i-1)+dt*thdd2(i-1);

ald2(i)=ald2(i-1)+dt*aldd2(i-1);

th1(i)=th1(i-1)+dt*thd1(i);

al1(i)=al1(i-1)+dt*ald1(i);

th2(i)=th2(i-1)+dt*thd2(i);

al2(i)=al2(i-1)+dt*ald2(i);

A=-2*g*sin(th1(i))-sin(th1(i)-th2(i))*thd2(i)^2;

J=-2*g*sin(al1(i))-sin(al1(i)-al2(i))*ald2(i)^2;

B=-g*sin(th2(i))+sin(th1(i)-th2(i))*thd1(i)^2;

K=-g*sin(al2(i))+sin(al1(i)-al2(i))*ald1(i)^2;

thdd1(i)=(A-B*cos(th1(i)-th2(i)))/(2-(cos(th1(i)-th2(i)))^2);

aldd1(i)=(J-K*cos(al1(i)-al2(i)))/(2-(cos(al1(i)-al2(i)))^2);

thdd2(i)=B-cos(th1(i)-th2(i))*thdd1(i);

aldd2(i)=K-cos(al1(i)-al2(i))*aldd1(i);

X1(i)=sin(th1(i));

D1(i)=sin(al1(i));

Y1(i)=cos(th1(i));

Z1(i)=cos(al1(i));

X2(i)=sin(th1(i))+sin(th2(i));

D2(i)=sin(al1(i))+sin(al2(i));

Y2(i)=cos(th1(i))+cos(th2(i));

Z2(i)=cos(al1(i))+cos(al2(i));

plot([0, X1(i), X2(i)], [0, -Y1(i), -Y2(i)],'-o');

hold on

plot([0, D1(i), D2(i)], [0, -Z1(i), -Z2(i)],'-o');

axis([-2 2 -2 2]);

title(['t = ', num2str(double(i)*dt, '% 5.3f'), ' s']);

hold on

plot(X2(1:i), -Y2(1:i), 'r');

plot(D2(1:i), -Z2(1:i), 'g');

hold off

drawnow;

end