Transient method of Simple Pendulum

//Code to generate ://
//i.	the angular displacement vs time plot//
//ii.	the angular velocity vs time plot//	
//iii.	the animation of pendulum motion//


clear all
close all
clc
 
b = 0.05;
g = 9.81;
L = 1;
m = 1;
i=1;
 
% Initial
theta_0 = [0;3];
 
% time limit
t_span = linspace(0,20,250);
 
%Solving
[t,results] = ode45(@(t,theta)ode_fun(t,theta,b,g,L,m),t_span,theta_0);
 
figure(1)
 
%angular velocity 
subplot(1,2,1)
plot(t,results(:,2),\'color\',\'g\') 
xlabel(\'time(s)\')
ylabel(\'angular velocity(rad/s)\')
title(\'velocity vs time\')
%angular displacement
subplot(1,2,2)
plot(t,results(:,1),\'color\',\'r\')  
xlabel(\'time(s)\')
ylabel(\'angular displacement\')
title(\'theta vs time\')
hold on
 
 
%Animation of the pendulum motion
 theta = results(:,1);
i =1;
for i= 1:length(t)
    THETA =  theta(i);
    
    %origin
    x0 = 0;
    y0 = 0;
    %End
    x1 = L*sin(THETA);
    y1 = -L*cos(THETA);
    
    
    figure(2)
    plot([x0 x1],[y0 y1],\'linewidth\',1.5,\'color\',\'black\')
    hold on
    plot(x1,y1,\'marker\',\'o\',\'markersize\',30,\'color\',\'black\',\'MarkerFaceColor\',\'black\')
    plot(x0,y0,\'-*\',\'linewidth\',2,\'color\',\'y\');
    title(\'Pendulum\')
    hold off
    axis([-3 3 -3 0])
    pause(0.01)
 
 
 
    % Taking frame
    M1(i)=getframe(gcf);
    i=i+1;
end
 
movie(M1);
videofile = VideoWriter(\'pendulum_animation.avi\',\'Uncompressed avi\');
open(videofile)
writeVideo(videofile,M1)
close(videofile)


//Function for generating the second degree equation which governs the pendulum motion//

function [dtheta_dt] = ode_fun(t,theta,b,g,L,m)
theta1= theta(1)
theta2= theta(2)
dtheta1_dt= theta2;
dtheta2_dt= -(b/m)*theta2 - (g/L)*sin(theta1);
dtheta_dt = [dtheta1_dt; dtheta2_dt];
end

STEPWISE EXPLANATION OF THE CODE

1. The pgrm starts with clear all close all and clc
2. Taking the inputs damping constant(b = 0.05), acceleration due to gravity in m/s^2(g = 9.81), length of the pendulum in m(L = 1) and mass of the pendulum bob in kg(m = 1). Along with this a counter (i=1) is given for animation.

• Initial value of angular displacement and velocity (boundary conditions) is taken to a column matrix theta_0 as theta_0 = [0;3]

Giving the time limits for the pendulum motion 0s to 10s and takes 50 values (t_span = linspace(0,20,250))

1. For the solving the differential equation using the function ode45, which takes the input t and theta, the function ode_fun (defined earlier, explained later) whose inputs are t,theta,b,g,L,m and the boundary conditions t_span and theta_0

[t,results] = ode45(@(t,theta)ode_fun(t,theta,b,g,L,m),t_span,theta_0)

1. Plotting the the angular displacement vs time graph and angular velocity vs time plot
2. Figure 1 is created figure(1)
3. The subplot command determines how many graphs needed inside the figure The form of the command subplot is subplot(no. of rows,no.of columns,selected subplot number)

For eg: subplot(1,2,1) creates space for 1x2=2 subplots and subplot 1 is selected here

1. Plot is done for the column 2 of the matrix (t,results) with color green.

plot(t,results(:,2),\'color\',\'g\')

1. xlabel, ylabel & name of graph is given

xlabel(\'time(s)\')

ylabel(\'angular velocity(rad/s)\')

title(\'velocity vs time\')

1. The same procedure is repeated for the theta v/s time plot

subplot(1,2,2)

plot(t,results(:,1),\'color\',\'r\')

xlabel(\'time(s)\')

ylabel(\'angular displacement\')

title(\'theta vs time\')

1. Plotting ends with using a hold on command which retains plots in the current axes so that new plots added to the axes do not delete existing plots.
2. Creating an animation of pendulum motion
3. A for loop is started which is governed by the value of a counter i

for i= 1:length(t)

1. The value theta(i) is read towards a variable THETA

THETA = theta(i)

1. Defining the start points and end points of the pendulum string

The origin is

x0 = 0

y0 = 0

The end point is

x1 = L*sin(THETA)

y1 = -L*cos(THETA)

1. Naming the pendulum figure as figure(2)
2. Plotting the line in black with width 1.5

plot([x0 x1],[y0 y1],\'linewidth\',1.5,\'color\',\'black\')

1. Using hold on command which retains plots in the current axes so that new plots added to the axes do not delete existing plots.
2. Plotting of the bob must be well specified. The plot is done at (x1,y1).

Different markers can be used using marker command.

For eg. s for square shape, p for star shape and o for circular which id needed in the case of a pendulum bob.

Then the size of the marker is specified using markersize command

And the color inside the marker is given using MarkerFaceColor commad

plot(x1,y1,\'marker\',\'o\',\'markersize\',30,\'color\',\'black\',\'MarkerFaceColor\',\'black\')

1. Plotting the fixed link at (x0,y0)
2. Naming the plot Pendulum
3. The using of hold off command erases the previous plot for the next plot for the next set of for loop operation.
4. Then axis and pause between frames is specified

axis([-3 3 -3 0])

pause(0.01)

1. The for loop ended by capturing the frame to a variable for movie creation using getframe command

M1(i)=getframe(gcf)

1. Finally the counter is incremented

i=i+1

1. The movie is created using the following commands

movie(M1);

videofile = VideoWriter(\'pendulum_animation.avi\',\'Uncompressed avi\');

open(videofile)

writeVideo(videofile,M1)

close(videofile)

1. ode_fun: The second ordr diff. eqn governing the pendulum motion is broke up into to first order diff eqn
2. theta(1) and theta(2) are taken into two variables theta1 and theta2 respectively.

theta1= theta(1)

theta2= theta(2)

1. theta(1) and theta(2) are differentiated wrt time.

dtheta1_dt= theta2

dtheta2_dt= -(b/m)*theta2 - (g/L)*sin(theta1)

1. created a column matrix to which the values of dtheta1_dt & dtheta2_dt are readed

dtheta_dt = [dtheta1_dt; dtheta2_dt]

1. The function[dtheta_dt] is ended which used in step

ERRORS AND SOLUTIONS

1. Each frame comes over other frame, the previous plot doesn’t disappears.

Solution: Usage of hold off command

hold off command erases the previous plot for the next plot for the next set of for loop operation.

1. Subplots are not created

Solution: Providing the subplots command in the format

subplot(no. of rows,no.of columns,selected subplot number) For eg: subplot(1,2,1)

NB: I have compressed the video using a third party app since the uploading took a lot of time due to the large size