By: Tao Steven Zheng
Here I present the code I wrote for my computational physics project in 2013 (UBC PHYS 210).
Code for the Trajectory[]
function [] = compend(tmax, dt, theta1_0, theta2_0, omega1_0, omega2_0)
format long
%Constants of pendulum system:
% g -- gravity (m/s/s)
% L1 -- top rod length (m)
% L2 -- bottom rod length (m)
% m1 -- top mass (kg)
% m2 -- bottom mass (kg)
g = 9.8; L1= 1; L2 = 1; m1 = 1; m2 = 1;
%Initial parameters:
% theta1_0 -- top angle (degrees)
% theta2_0 -- bottom angle (degrees)
% omega1_0 -- top angular speed
% omega2_0 -- bottom angular speed
nt = tmax/dt;
t = linspace(0, tmax, nt-1);
t0 = linspace(0, tmax, nt);
%Initial conditions
omega1(1) = omega1_0;
omega2(1) = omega2_0;
theta1(1) = theta1_0;
theta2(1) = theta2_0;
%System of equations of the compound pendulum that is evolved to the final time
for i = 1:nt-1
k1(i) = -g*(2*m1+m2)*sin(theta1(i))-m2*g*sin(theta1(i)-2*theta2(i))-2*sin(theta1(i)-theta2(i))*m2*(L2*((omega2(i))^2)+L1*((omega1(i))^2)*cos(theta1(i)-theta2(i)));
k2(i) = 2*sin(theta1(i)-theta2(i))*(L1*((omega1(i))^2)*(m1+m2)+g*(m1+m2)*cos(theta1(i))+((omega2(i))^2)*L2*m2*cos(theta1(i)-theta2(i)));
z(i) = 2*m1+m2-m2*cos(2*theta1(i)-2*theta2(i));
alpha1(i) = k1(i)/(L1*(z(i)));
alpha2(i) = k2(i)/(L2*(z(i)));
omega1(i+1) = omega1(i) + dt * alpha1(i);
omega2(i+1) = omega2(i) + dt * alpha2(i);
theta1(i+1) = theta1(i) + dt * omega1(i);
theta2(i+1) = theta2(i) + dt * omega2(i);
%Coordinates of top mass and bottom mass
x1(i) = L1*sin(theta1(i));
y1(i) = -L1*cos(theta1(i));
x2(i) = x1(i)+L2*sin(theta2(i));
y2(i) = y1(i)-L2*cos(theta2(i));
%Energy equations that describe the Lagrangian equations of motion
T(i) = 0.5*m1*(L1^2)*((omega1(i))^2)+0.5*m2*[(L1^2)*(omega1(i))^2+(L2^2)*(omega2(i)^2)+2*L1*L2*(omega1(i))*(omega2(i))*cos(theta1(i)-theta2(i))];
V(i) = -(m1+m2)*g*L1*cos(theta1(i))-m2*g*L2*cos(theta2(i));
E(i) = T(i) + V(i);
end
end
Code for the Animation of the Compound Pendulum[]
function [] = animate(tmax, dt, theta1_0, theta2_0, omega1_0, omega2_0)
format long
%Constants of pendulum system:
% g -- gravity (m/s/s)
% L1 -- top rod length (m)
% L2 -- bottom rod length (m)
% m1 -- top mass (kg)
% m2 -- bottom mass (kg)
g = 9.8; L1= 1; L2 = 1; m1 = 1; m2 = 1;
%Initial parameters:
% theta1_0 -- top angle (degrees)
% theta2_0 -- bottom angle (degrees)
% omega1_0 -- top angular speed
% omega2_0 -- bottom angular speed
nt = tmax/dt;
%Initial conditions
omega1(1) = omega1_0;
omega2(1) = omega2_0;
theta1(1) = theta1_0;
theta2(1) = theta2_0;
%System of equations of that evolves the compound pendulum
for i = 1:nt-1
alpha1(i) = (-g*(2*m1+m2)*sin(theta1(i))-m2*g*sin(theta1(i)-2*theta2(i))-2*sin(theta1(i)-theta2(i))*m2*(L2*((omega2(i))^2)+L1*((omega1(i))^2)*cos(theta1(i)-theta2(i))))/(L1*(2*m1+m2-m2*cos(2*theta1(i)-2*theta2(i))));
alpha2(i) = (2*sin(theta1(i)-theta2(i))*(L1*((omega1(i))^2)*(m1+m2)+g*(m1+m2)*cos(theta1(i))+((omega2(i))^2)*L2*m2*cos(theta1(i)-theta2(i))))/(L2*(2*m1+m2-m2*cos(2*theta1(i)-2*theta2(i))));
omega1(i+1) = omega1(i) + dt * alpha1(i);
omega2(i+1) = omega2(i) + dt * alpha2(i);
theta1(i+1) = theta1(i) + dt * omega1(i);
theta2(i+1) = theta2(i) + dt * omega2(i);
%Coordinates of top mass and bottom mass
x1(i) = L1*sin(theta1(i));
y1(i) = -L1*cos(theta1(i));
x2(i) = x1(i)+L2*sin(theta2(i));
y2(i) = y1(i)-L2*cos(theta2(i));
%A plot of the compound pendulum with lines representing the rigid bars
clf;
plot(x1, y1, 'red', x2, y2, 'blue'); axis([-2 2 -2 2]);
hold on;
drawnow;
line([0 x1(i)], [0 y1(i)]);
line([x1(i) x2(i)], [y1(i) y2(i)]);
drawnow;
hold off;
end
end
Results[]
