Write a program to simulate the transient behaviour of a simple pendulum and to create an animation of it's motion.

Aim of the project:-- To solve 2nd order ODE of motion of Simple Pendulum by MATLAB & to simulate the transient behevior of the pendulum & to make a animation of the motion.

Code to solve the 2nd order ODE of motion of Simple Pendulum & to simulate the transient behevior of the pendulum :--

% Program to solve the differntial equation of motion of Simple Pendulum & 
% simulation of its motion.

clear all
clc

b=0.05;  %[damping co efficient,accelaration due to gravity,effective length 
g=9.81;  %of the simple pendulum & mass of the bob are declared here.]
l=1;
m=1;

% initial condition
theta_0=[0;5];     % initial angular diaplacement (0) & angular velocity(5)

%time points
t_span=linspace(0,20,500); % [Time range declared as 1 to 20 sec
                           % & within 20 seconds, for equaly divided 
                           % 500 moments data will be evluated]
% Solve ODE
[t,results]=ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);
                           % [ode45 solver is used to calculate angular displacement 
                           % & angular velocity]
                          
 plot(t,results(:,1))    %[after definite intregration by solver all the 
 hold on                 %angular displacement are stored in result(:,1)
                         % & angular velocity in result(:;2)as column matrix]

 plot(t,results(:,2))    %[The displacement & velocity are ploted w.r.t
                         % time]
                          

ylabel('Plot')
xlabel('time')
legend ('angular displacement','angular velocity')




%The code of function 'ode_func' which is called by solver ode45 from main program.

function[dtheta_dt]=ode_func(t,theta,b,g,l,m) % name of fuction 'ode_func' 
                                              % is declared with all input
                                              % varriables.]
 
theta1=theta(1);   % 1st value of the 'theta'column matrix is stored in 'theta1'                   
theta2=theta(2);   % 2nd value of the 'theta'column matrix is stored in 'theta2'   

dtheta1_dt=theta2;                           % 1st order ODE .....no-1
dtheta2_dt= -(b/m)*theta2-(g/l)*sin(theta1); %1st order ODE .... .no-2  

dtheta_dt=[dtheta1_dt;dtheta2_dt];  % [Differential equation in the form of column   
                                    %matrix'dtheta_dt'is declared which is required for 
                                    %solver ode45]
                                    
                                    
end





%The code of the program for  making an animation of the motion of simple pendulum using the out put of the solver ode45

b=0.05;
g=9.81;
l=1;
m=1;

% initial condition
theta_0=[0;3];

%time points
t_span=linspace(0,20,500);

% Solve ODE
[t,results]=ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);
d= results(:,1); % In the variables 'd' the elements of column matrix of 
                 % angular displacement is stored]



%Origin
O=[0 0];        % Co ordinate of Point of suspension is declared 

% loop for animation

ct=1;           %counter is declared

for i=1:length(t)
    
 %Bob centre
 c= l*[sin(d(i,1)) -cos(d(i,1))];  
 
               % [co ordinate of the centre of the bob is calculated for
               % different values of i ie through out the time span]         
   
    
 %pendulum 
 plot([O(1) c(1)],[O(2) c(2)],'linewidth',1,'color','b','marker','o','markersize',10,'markerfacecolor','r');
       
               % ploting the pendulum for through out the time span is done
 
  axis([-1 1 -1.5 0.2]);
  grid on        
  pause(0.01)
  M(ct)=getframe(gcf);
  
  ct=ct+1;    
    
end     
   movie(M);
videofile=VideoWriter('Damped SHM');
open(videofile);
writeVideo(videofile,M);    % video of damped motion of simple pendulum is made.
close(videofile); 


Link of animation of motion of Simple pendulum from Google Drive
https://drive.google.com/open?id=1ytbSKtoTuuaCAotjr94q_2wLK6iFnjzx

ODE:-- Ordinery Differential Equation (It may be of any order but must have only two varriables.)

Introduction:- The 2nd order ODE used here, is of motion of simple pendulum. The motion is damped simple harmonic motion. Damping coefficient (b) used here is 0.05. The restoring force mgSin(theta) is alwyas directed to the equlibrium point(mid point of the path)& is propotional to displacement of the bob at any moment, measured from equilibrium point. This restoring force is the cause of the motion of the pendulum. Damping force which is due to frictional force offered by surrounding air, or other resistance,alwayas opposes the motion. To maintain the oscillation against such resistive force, the pendulum has to do some work which cause the errosion of its energy. Due to this, amplitude of the oscillation diminishes gradualy & the pendulum comes to rest. This is called damped oscillation.

Process:-- Here from the 2nd order ODE , two 1st orer ODEs are made.

Here, theta= angular displacement (in radian); t=time in second

d(theta)/dt=angular velocity; d2(theta)/dt2 = angular accelaration

d2(theta)/dt2 = -(b/m)*d(theta)/dt - (g/l)*sin(theta)----- 2nd order ODE

Let, theta=theta1. Let, d(theta)/dt = d(theta1)/dt =theta2

d2(theta)/dt2=d2(theta1)/dt2=d(theta2)/dt

So we are getting two 1st order ODEs

d(theta1)/dt=theta2 ...........1

d(theta2)/dt = -(b/m)*theta2 - (g/l)*sin(theta1) ............2

Now in the function named 'ode_func', a column matrix 'theta' is made in which there are two elements .1st element, theta(1)= theta1,2nd element, theta(2)=theta2

Now from the function,differential equation in column-matrix-form, which is dtheta_dt =[dtheta1_dt ; dtheta_dt], is sent to solver ode45,to solve the ODE when the function is called by the solver in the main program.

As input, damping coeficient b=0.5, accelaration due to gravity g=9.81m/s2, effective length of pendulum l=1m, mass of bob of pendulum m=1kg.

Here the solver ODE45 actulay does definite intregation in the time range 0 to 20 & give 500 values of angular displacements & angular velocities of the pendulum respectively, in time range from 0 to 20sec. Here at initial condition (t=0), displacement & velocity are taken as 0 radian & 3 radian/sec respectively. As out put of integration, the displacement & velocity are stored in result1 & result2 as column matrix respectively. These are the solution of the ODE in the time range t=0 to 20

When we plot these angular diasplacements & angular velocities with respect to time, we get the graph as below which shows the simulation of transit behevior of simple pendulum.

The graph is showing how values of angular displacement & angular velocity change w.r.t time. When displacepent is zero, velocity is maximum & when displacement is maximum, velocity is zero. Displacement is represented by sine curve & velocity by cosine curve (at t=0, displacement=0, velocity is maximum).Besides this we can see in the graph, that both displacement & velocity of pendulum are diminishing due to damping effect.

For animation of motion of Simple pendulum, 500 nos of instanteneous angular displacements in the time range 0 to 20sec, which is one of the out put of the solution of 2nd order ODE (by solver ODE 45) of the motion, are used. For animation, other inputs is effective length l =1m.

Though the link of google drive for animation of the motion is given at end of coding, below is the screen shot of the amimation.

Here are some errors are given which have been faced in the course of making code.

Error1

Here parentheses after input parameter of function is not given.

Error2

Here in Error2 spelling of function name was wrong.

Error 3

Here in Error3 initial value of ct is wromgly taken as 0. It should be +ve integer.