Question

write a matlab script for double mass spring system with animation.

Answer 1

% call: >> double_pendulum([pi;0;pi;5;9.81;1;1;2;1],100,10,false)
% Or, simply call >> double_pendulum_init
%

function double_pendulum(ivp, duration, fps, movie)

% ---------------------------------------------------------------------
clear All; clf;
nframes=duration*fps;
sol=ode45(@double_pendulum_ODE,[0 duration], ivp);
t = linspace(0,duration,nframes);
y=deval(sol,t);
phi1=y(1,:)'; dtphi1=y(2,:)';
phi2=y(3,:)'; dtphi2=y(4,:)';
l1=ivp(8); l2=ivp(9);
% phi1=x(:,1); dtphi1=x(:,2);
% phi2=x(:,3); dtphi2=x(:,4);
% l1=ivp(8); l2=ivp(9);
h=plot(0,0,'MarkerSize',30,'Marker','.','LineWidth',2);
range=1.1*(l1+l2); axis([-range range -range range]); axis square;
set(gca,'nextplot','replacechildren');
for i=1:length(phi1)-1
if (ishandle(h)==1)
Xcoord=[0,l1*sin(phi1(i)),l1*sin(phi1(i))+l2*sin(phi2(i))];
Ycoord=[0,-l1*cos(phi1(i)),-l1*cos(phi1(i))-l2*cos(phi2(i))];
set(h,'XData',Xcoord,'YData',Ycoord);
drawnow;
F(i) = getframe;
if movie==false
pause(t(i+1)-t(i));
end
end
end
if movie==true
movie2avi(F,'doublePendulumAnimation.avi','compression','Cinepak','fps',fps)
end

% Simply run
% >> double_pendulum_init

%This script calls double_pendulum.
%
% ---------------------------------------------------------------------
phi1 = pi;
dtphi1 = 0;
phi2 = pi;
dtphi2 = 5;
g = 9.81;
m1 = 1;
m2 = 1;
l1 = 2;
l2 = 1;
duration = 100;
fps = 10;
movie = true;
clc; figure;
interval=[0, duration];
ivp=[phi1; dtphi1; phi2; dtphi2; g; m1; m2; l1; l2];
double_pendulum(ivp, duration, fps, movie);

% This function calls is called by double_pendulum.

% ---------------------------------------------------------------------

function xdot = double_pendulum_ODE(t,x)
g=x(5); m1=x(6); m2=x(7); l1=x(8); l2=x(9);
xdot=zeros(9,1);
xdot(1)=x(2);
xdot(2)=-((g*(2*m1+m2)sin(x(1))+m2(g*sin(x(1)-2*x(3))+2*(l2*x(4)^2+...
l1*x(2)^2*cos(x(1)-x(3)))*sin(x(1)-x(3))))/...
(2*l1*(m1+m2-m2*cos(x(1)-x(3))^2)));
xdot(3)=x(4);
xdot(4)=(((m1+m2)(l1*x(2)^2+g*cos(x(1)))+l2*m2*x(4)^2*cos(x(1)-x(3)))...
sin(x(1)-x(3)))/(l2*(m1+m2-m2*cos(x(1)-x(3))^2));