Hi,
I have edited this code to give timeArray, theta1 and theta2 as outputs.
First save edited version, then call the function with its outputs like below:
(You will see the t,theta1 and theta2 as outputs in your workspace)
[t, theta1, theta2] = double_pendulum_simulation
Then to plot do this:
timeArrayNew = linspace (t(1), t(end),length(theta1));
figure;
plot(timeArrayNew,theta1);
hold on;
plot(timeArrayNew,theta2);
legend('theta1','theta2');
Edited function:
function [timeArray, theta1Array, theta2Array] = double_pendulum_simulation(theta1_0,theta2_0,L1,L2,m1,m2,g,tail)
%DOUBLE_PENDULUM_SIMULATION Plots double pendulum dynamics until stopped
%
% Simulates and plots the dynamics behavior of the double simple
% pendulum until the user presses a key or clicks in the figure.
% If any of the input parameters is not given, it will be assigned a
% default value.
%
% Example uses:
% DOUBLE_PENDULUM_SIMULATION
% DOUBLE_PENDULUM_SIMULATION(THETA1,THETA2,L1,L1,M1,M2,G,TAIL)
% THETA1_0 and THETA2_0 initial angles (default 3*pi/4 and pi/2)
% L1 and L2 rod lengths (default 3 and 2)
% M1 and M2 bob masses (default 2 and 3)
% G acceleration of gravity (default 9.81)
% TAIL lengh of visualization tail (default 100)
%
% Author: vand@dtu.dk, 2014
if nargin<8 || isempty(tail)
tail = 20;
end
if nargin<7 || isempty(g)
g = 9.81;
end
if nargin<6 || isempty(m2)
m2 = 10;
end
if nargin<5 || isempty(m1)
m1 = 2;
end
if nargin<4 || isempty(L2)
L2 = 1;
end
if nargin<3 || isempty(L1)
L1 = 1;
end
if nargin<2 || isempty(theta2_0)
theta2_0 = pi/2;
end
if nargin<1 || isempty(theta1_0)
theta1_0 = 0;
end
% preparing for loop until user either keypresses or clicks
global USER_RESPONDED
USER_RESPONDED = 0;
figure
set(gcf,'WindowKeyPressFcn',@userRespondFcn,'WindowButtonDownFcn',...
@userRespondFcn,'DeleteFcn',@userRespondFcn)
% initial state
x = [mod(theta1_0,2*pi),mod(theta2_0,2*pi),0,0];
t = 0;
i_end = 0.02; % initial estimate of a average iteration duration
tail = repmat(L1*[sin(x(1)),-cos(x(1))],[tail,2]) + ...
repmat([0 0 L2*[sin(x(2)),-cos(x(2))]],[tail,1]); % initial tail
% ode function
double_penudlum = @(t,x)double_pendulum_system(x,L1,L2,m1,m2,g);
% preparing for loop
axis xy equal, box on, hold on
axis(1.1*[-1 1 -1 1]*(L1+L2))
[r1,r2] = bob_drawing_scale(m1,m2,L1,L2); % scale for drawing bobs
iter = 0;
tic
theta1Array = [];
theta2Array = [];
timeArray = [];
while ~USER_RESPONDED
[t,x] = ode45(double_penudlum,t(end)*[1 0.5 0] + i_end*[0 0.5 1] ,x(end,:)');
[tail, theta1, theta2] = patch_double_pendulum(t,x(end,:),L1,L2,r1,r2,tail);
theta1Array = [theta1Array;theta1];
theta2Array = [theta2Array;theta2];
timeArray = [timeArray;t];
iter = iter+1;
i_end = max(toc*(1+1/iter),t(end)+2*eps); % end time for next iteration
end
end
function userRespondFcn(~,~)
% callback function which execute when user responds
global USER_RESPONDED
USER_RESPONDED = 1;
end
function dx = double_pendulum_system(x,L1,L2,m1,m2,g)
% a system of differential equations defining a double pendulum
% from http://www.myphysicslab.com/dbl_pendulum.html
theta1 = x(1);
theta2 = x(2);
omega1 = x(3);
omega2 = x(4);
dtheta1 = omega1;
dtheta2 = omega2;
domega1 = (-g*(2*m1+m2)*sin(theta1)-m2*g*sin(theta1-2*theta2)-...
2*sin(theta1-theta2)*m2*(omega2^2*L2+omega1^2*L1*cos(theta1-theta2)))...
/(L1*(2*m1+m2-m2*cos(2*theta1-2*theta2)));
domega2 = (2*sin(theta1-theta2)*(omega1^2*L1*(m1+m2)+...
g*(m1+m2)*cos(theta1)+omega2^2*L2*m2*cos(theta1-theta2)))...
/(L2*(2*m1+m2-m2*cos(2*theta1-2*theta2)));
dx = [dtheta1;dtheta2;domega1;domega2];
end
function [tail, theta1, theta2] = patch_double_pendulum(t,x,L1,L2,r1,r2,tail)
% using patch to plot pendulum state
cla
% plotting tail
alpha = linspace(0,1,size(tail,1)+1)';
patch([tail(:,1);NaN],[tail(:,2);NaN],0,'EdgeColor','r','FaceColor',...
'none','FaceVertexAlphaData',alpha,'EdgeAlpha','interp','LineWidth',2);
patch([tail(:,3);NaN],[tail(:,4);NaN],0,'EdgeColor','r','FaceColor',...
'none','FaceVertexAlphaData',alpha,'EdgeAlpha','interp','LineWidth',2);
% plotting rods
theta1 = x(1);
theta2 = x(2);
xm1 = L1*sin(theta1);
ym1 = -L1*cos(theta1);
xm2 = xm1 + L2*sin(theta2);
ym2 = ym1 - L2*cos(theta2);
patch([0, xm1, xm2, NaN],[0, ym1, ym2, NaN],0,'EdgeColor','b',...
'FaceColor','none','LineWidth',2)
% plotting bobs
p = linspace(0,2*pi,17);
sint = sin(p);
cost = cos(p);
patch(xm1+r1*cost,ym1+r1*sint,0,'EdgeColor','b','FaceColor','b')
patch(xm2+r2*cost,ym2+r2*sint,0,'EdgeColor','b','FaceColor','b')
title(sprintf('time = %0.1f',t(end)))
drawnow
tail = [tail(2:end,:);xm1,ym1,xm2,ym2];
end
function [r1,r2] = bob_drawing_scale(m1,m2,L1,L2)
% finding a scale for plotting the size of the bobs: the bigger bob has
% a radius which is a given fraction of the length of the shorter rod.
r_max = max(m1^(1/3),m2^(1/3));
l_min = min(L1,L2);
scale = 0.1*l_min/r_max;
r1 = scale*m1^(1/3); % mass is proportional with the volume
r2 = scale*m2^(1/3);
end
