Simple Pendulum Simulation

Aim: Simulation of a Simple Pendulum using Matlab

Objectives: 

1. Determining the angular displacement, velocity and acceleration of a simple pendulum by solving second order ordinary differential equation.
2. Animate the solution while displaying instantaneous velocity and acceleration.

Background:

Simple pendulum motion is governed by the following equation:

`(d^2theta)/(dt^2)+b/m(d theta)/(dt)+g/lsin(theta) = 0`

Procedure:

Differential Equation Solution:

The second order differential equation is written as system of linear equations:

`theta1 = theta`

`v=(d theta1)/dt = theta2`

`a=(d^2 theta)/(dt^2) = -(b/m*theta2+g/l*sin(theta1))`

Using the equations following ODE function was defined:

%ODE function as first derivatives
function [Dtheta] = SimplePendulumODE(time,theta,l,m,b,g)
theta1 = theta(1);
theta2 = theta(2);
dtheta1_dt = theta2;
dtheta2_dt = - (b/m)*theta2 - (g/l)*sind(theta1);
Dtheta = [dtheta1_dt;dtheta2_dt];
end

Assumptions:

`l = 1m`

`m = 1kg`

`b = 0.05`

`g = 9.81m/s^2`

`theta circ = 0`

`v circ= (d theta circ)/dt = 3 m/s`

`t = 0 to 20 s`

Animation:

The Animation was using the following pieces:

1. Pendulum length

The length was plotted by using the angular displacement(theta) of the pendulum and its length(l)

Plotted at line 50 of the main program

x2 = -l*sind(kinematics(:,1));
y2 = -l*cosd(kinematics(:,1));

2. Mass of the pendulum with radius r at the free end (x2,y2) (called in line 51 of main program)

A solid circle was plotted by generation lines of length r around x2,y2 at angles t = 0 to 360 in intervals of 4 degrees.

%solid sphere at x2,y2 with a radius r
function pendulumMass(x2,y2,r)
hold on
for t=0:4:360
plot([x2,x2+r*cosd(t)],[y2,y2+r*sind(t)],'color',[0.5 0 0])
end
end

3. Velocity and Acceleration vectors starting from mass center with mangnitude and direction obtained from the ODE solver. Orientation of the vectors is pependicular to the pendulum.

Arrow Head:

The definition of arrow head depends on the direction of the vector, if-condition was used for this purpose. For each case two lines are needed to plot an arrow head. The orientation of the arrow was defined using axis transformation as it changes according to the orientation of the vector for each instance.

`[(X),(Y)] = [[cos(theta),sin(theta)],[-sin(theta),cos(theta)]]*[(x),(y)]`

function A = axistransform(a1,b1,t)
A=[cosd(t) sind(t);-sind(t) cosd(t)]*[a1;b1];
end

The function was called in the velocity and acceleration functions while defining the Arrows

Velocity Vector: 

One end of the vector being the free end of pendulum, the other end is defined by using magnitude and direction of velocity obtained from ODE solver and is plotted perpendicular to pendulum length.

Arrays x3,y3 denote the other end of the vector. 

x3 =  x2-0.2*(kinematics(:,2)).*cosd(-kinematics(:,1));%using angular velocity function to define velocity vector, magnitude reduced to 1/5th
y3 =  y2-0.2*(kinematics(:,2)).*sind(-kinematics(:,1));%using angular velocity function to define velocity vector, magnitude reduced to 1/5th

Velcoity function was called (line 56 of main.m) to plot the vector with arrow head.

function velocity(x1,y1,x2,y2,V)
t= V(1);
v= V(2);
plot([x1 x2],[y1 y2],'color',[0.5 0 0.5])%vector magnitude and orientation
hold on
%Arrow Heads using axis transformation
if v<0
T = 150;
a1=0.05*cosd(T);
b1=0.05*sind(T);
c = axistransform(a1,b1,t);
plot([x2,x2+c(1)],[y2,(y2+c(2))],'color',[0.5 0 0.5])
T = -150;
a1=0.05*cosd(T);
b1=0.05*sind(T);
c = axistransform(a1,b1,t);
plot([x2,x2+c(1)],[y2,(y2+c(2))],'color',[0.5 0 0.5])
elseif v>0   
T = 30;
a1=0.05*cosd(T);
b1=0.05*sind(T);
c = axistransform(a1,b1,t);
plot([x2,x2+c(1)],[y2,(y2+c(2))],'color',[0.5 0 0.5])
T = -30;
a1=0.05*cosd(T);
b1=0.05*sind(T);
c = axistransform(a1,b1,t);
plot([x2,x2+c(1)],[y2,(y2+c(2))],'color',[0.5 0 0.5])
end
end

Acceleration Vector:

One end of the vector being the free end of pendulum, the other end is defined by using magnitude and direction of velocity obtained from ODE solver and is plotted perpendicular to pendulum length.

Arrays x4,y4 denote the other end of the vector. 

x4 =  x2-0.5*(kinematics(:,3)).*cosd(-kinematics(:,1));%using angular acceleration function to define acceleration vector, magnitude reduced to 1/2th
y4 =  y2-0.5*(kinematics(:,3)).*sind(-kinematics(:,1));%using angular acceleration function to define acceleration vector, magnitude reduced to 1/2th

Acceleration function was called (line 57 of main.m) to plot the acceleration with arrow head.

function acceleration(x1,y1,x2,y2,V)
t= V(1);
v= V(3);
plot([x1 x2],[y1 y2],'color',[0 0.5 0.5])%vector magnitude and orientation
hold on
%Arrow Heads  using axis transformation
if v<0
T = 150;
a1=0.05*cosd(T);
b1=0.05*sind(T);
c = axistransform(a1,b1,t);
plot([x2,x2+c(1)],[y2,(y2+c(2))],'color',[0 0.5 0.5])
T = -150;
a1=0.05*cosd(T);
b1=0.05*sind(T);
c = axistransform(a1,b1,t);
plot([x2,x2+c(1)],[y2,(y2+c(2))],'color',[0 0.5 0.5])
elseif v>0   
T = 30;
a1=0.05*cosd(T);
b1=0.05*sind(T);
c = axistransform(a1,b1,t);
plot([x2,x2+c(1)],[y2,(y2+c(2))],'color',[0 0.5 0.5])
T = -30;
a1=0.05*cosd(T);
b1=0.05*sind(T);
c = axistransform(a1,b1,t);
plot([x2,x2+c(1)],[y2,(y2+c(2))],'color',[0 0.5 0.5])
end
end

Main:

1. Using ODE45 solver in matlab(line 19) to solve for time, theta with the assumed time span and initial conditions, kinematics array was obtained with angular displacement and velocity. Acceleration was added to the array using the governing equation(line 21) and all the variables were plotted(lines 24 through 29) with respect to time.

2. For loop was used to plot all the pieces of the animation for each instance of time i (lines 43 to 60)

3. Frames from each instance were saved to an array (line 58) and used to generate a .avi file using VideoWriter function(lines 61-66) 

clear all
close all
clc

%D2(theta)+(b/m)D(theta)+(g/l)sind(theta)=0 governing equation for simple
%pendulum, where D# is time derivative of order #

%input
l = 1 %m
m = 1 %kg
b = 0.05
g = 9.81 %m/s^2
%initial displacement = 0 initial velocity = 3 
%acceleration = - (b/m)* velocity - (g/l)sin(theta)
initial_condition = [0 3]
span = 0:0.1:20;

%ODE solver
[time,kinematics]=ode45(@(time,theta) SimplePendulumODE(time,theta,l,m,b,g),span,initial_condition);

kinematics(:,3) = - ((b/m)*kinematics(:,2)) - ((g/l)*sind(kinematics(:,1)));

%plot
plot(time,kinematics(:,1),'color',[0.7 0 0.9]) %pink
hold on 
plot(time,kinematics(:,2),'color',[0 0 0.7]) %blue
plot(time,kinematics(:,3),'color',[0 0.5 0.5]) %green
legend('Angular Displacement','Angular Velocity','Acceleration')
xlabel('time')

figure

%Penduum animation, pendulum linked to wall at 0,0

%animation of simple pendulum

x2 = -l*sind(kinematics(:,1));%using angular displacement to find the free end of pendulum
y2 = -l*cosd(kinematics(:,1));%using angular displacement to find the free end of pendulum 
x3 =  x2-0.2*(kinematics(:,2)).*cosd(-kinematics(:,1));%using angular velocity function to define velocity vector, magnitude reduced to 1/5th
y3 =  y2-0.2*(kinematics(:,2)).*sind(-kinematics(:,1));%using angular velocity function to define velocity vector, magnitude reduced to 1/5th
x4 =  x2-0.5*(kinematics(:,3)).*cosd(-kinematics(:,1));%using angular acceleration function to define acceleration vector, magnitude reduced to 1/2th
y4 =  y2-0.5*(kinematics(:,3)).*sind(-kinematics(:,1));%using angular acceleration function to define acceleration vector, magnitude reduced to 1/2th
for i = 1:length(kinematics(:,1))
    clf
    plot([-0.5,0.5],[0,0],'linewidth',6,'color',[0 0 0])%wall
    axis([-2 2 -1.6 0.4])
    hold on
    plot([-2,2],[0,0],[0,0],[0.4,-1.6],'linewidth',0.2,'color',[0 0 0])%axis reference lines
    plot(x2,y2,'color',[0.9 0 0])%path followed by center of mass of the spherical mass
    plot([0,x2(i)],[0,y2(i)],'color','b','linewidth',2)%instantaneous position of pendulum
    pendulumMass(x2(i),y2(i),0.1)%spherical mass representation 
    label = 'Acceleration vector(x0.5)';
    text(0.5,-1.3,label,'color',[0 0.5 0.5])
    label = 'Velocity vector(x0.2)';
    text(0.5,-1.4,label,'color',[0.5 0 0.5])
    velocity(x2(i),y2(i),x3(i),y3(i),kinematics(i,:))%velocity vector    
    acceleration(x2(i),y2(i),x4(i),y4(i),kinematics(i,:))%acceleration vector
    M(i)=getframe(gcf);%figure captured in array to make the animation
    pause(0.001)
end
movie(M);
f=VideoWriter('SimplePendulum.avi','uncompressed AVI');
v.FrameRate = 5;
open(f);
writeVideo(f,M);
close(f);

Results:

Solution Plot:

Animation:

SimplePendulum.avi