% Parameters and initial conditions
% mass, link length, initial angles, simulation time
t = linspace(0, 10, 200);
% Solving ODE of a double pendulum
[T,Y] = ode45(@(t, x) double_pendulum(t, x, m, L), ...
t, [theta1, theta2, 0, 0]);
% Calculating joint coordinates for animation purposes
x = [L*sin(Y(:,1)), L*sin(Y(:,1))+L*sin(Y(:,2))];
y = [-L*cos(Y(:,1)), -L*cos(Y(:,1))-L*cos(Y(:,2))];
% Convert radians to degrees
tiledlayout(2,1) % Introduced in R2019b
nexttile % Introduced in R2019b
plot(T,ang,LineWidth=2) % Name=value syntax introduced in R2021a!!
hh1(1) = plot(T(1),ang(1,1),Marker=".",MarkerSize=20,Color="b");
hh1(2) = plot(T(1),ang(1,2),Marker=".",MarkerSize=20,Color="r");
hh2 = plot([0,x(1,1);x(1,1),x(1,2)],[0,y(1,1);y(1,1),y(1,2)], ...
Marker=".",MarkerSize=20,LineWidth=2);
axis([-2*L 2*L -2*L 2*L])
ht = title("Time: "+T(1)+" sec"); % String arrays introduced in R2016b
set(hh1(1),XData=T(id),YData=ang(id,1))
set(hh1(2),XData=T(id),YData=ang(id,2))
set(hh2(1),XData=[0,x(id,1)],YData=[0,y(id,1)])
set(hh2(2),XData=x(id,:),YData=y(id,:))
set(ht,String="Time: "+T(id)+" sec")
function dy = double_pendulum(t, y, m, L)
% System of ODEs for a double pendulum (mass m and link length L)
% See http://en.wikipedia.org/wiki/Double_pendulum for the differential
% Copyright 2010 The MathWorks, Inc.
dy(1) = 6/(m*L^2) * (2*p1-3*cos(theta1-theta2)*p2) / ...
(16-9*cos(theta1-theta2)^2);
dy(2) = 6/(m*L^2)*(8*p2-3*cos(theta1-theta2)*p1) / ...
(16-9*cos(theta1-theta2)^2);
dy(3) = -1/2*m*L^2*(dy(1)*dy(2)*sin(theta1-theta2)+3*g/L*sin(theta1));
dy(4) = -1/2*m*L^2*(-dy(1)*dy(2)*sin(theta1-theta2)+g/L*sin(theta2));
Comments
To leave a comment, please click here to sign in to your MathWorks Account or create a new one.