Simple Pendulum Simulation by solving 2nd order ODE

Report : 

Simple Pendulum Simulation by solving second order ODE . 

.                                            

• The objective of this program is to simulate a simple pendulum by solving second order ODE into two first order ODE\'s. The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the \"ball\" w.r.t time. 
• The governing differential equation representing the motion of simple pendulum with damping is given as                                  

•              
•  where , `Theta` = Displacement 

         b    = damping Co-efficient 

         g     = gravity m/s^2

         L     = length of string in (m)

         m    = mass of ball (kg)

• Upon solving the given Second order ode , we convert it into two first order ode\'s , which are  :

              `(d(theta1))/dt = theta2`         &     `(d(theta2))/dt = - b/m* theta2 - g/l sin*theta1`

• Program : 

clear all
close all
clc

% input parameters
b = 0.05;
l = 1;
g = 9.81;
m = 1;

%displacement and velocity column wise
theta_0 = [0;3];

% time span 0-20 sec , 500 values in between
t_span = linspace(0,20,500);

% using ode45 command to solve second order ODE 
[t , results] = ode45(@(t , theta) pendufunc(t,theta,b,g,l,m),t_span,theta_0);

% plotting time vs displacement & time vs angular velocity
subplot(5,4, [13, 14, 17 ,18])
plot(t,results(:,1))

hold on 
plot(t,results(:,2))
hold off 

subplot(5,4,[15 16 19 20])
plot(results(:,1),results(:,2))
xlabel('Position (rad)')
ylabel('Angular velocity')


ct = 1;

% for loop to run the loop for simple pendulum angles for different values in time
for i = 1:length(results(:,1))
    x0 = 0;
    y0 = 0;
    x1 = -l*sin(results(i,1));
    y1 = -l*cos(results(i,1));
    subplot(5,4,[1 2 3 4 5 6 7 8 9 10 11 12]);
    
    % plotting the rigid base 
    plot([-1 1] ,[ 0 0] , 'linewidth' ,2,'color','g');
    hold on 
    % plotting string
    line([x0 x1] , [y0 y1], 'linewidth',1,'color','[0.5 0.5 0.5');
    hold on ;
     
    % plotting the mass
     plot(x1,y1,'o' ,'markers',20,'markerfacecolor','[ 0.5 1 1]');
    
     r = 0.25;
    rectangle('position', [(x1 -r/2), (y1 -r/2 ),r,r],'curvature',[1 1] ,'facecolor','g','linewidth',3)
    hold on
    axis equal
    
    grid on 
    axis ([-1.5 1.5 -1.5 0.5])
    pause(0.003);
    hold off
    M(ct) = getframe(gcf);
    ct = ct+1;
end

% movie making 
movie(M)
videofile = VideoWriter('simple_pendulum.avi','uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile) 

 Explanation of code: 

1. Input parameters for the Simple pendulum :
   1.  b = 0.05
   2. L  = 1 m 
   3. g = 9.81 m/s^2
   4. m = 1 kg

2.  

   1. Initial displacement = 0
   2. initial velocity = 3
3. Time span is 0 - 20 seconds for the motion of pendulum is selected for 500 intervals in between.
4. Now the ODE45 command is used to solve the second order ODE , which we resolved into 2 first order ode\'s , so we are solving for the later .
5. This command is executed using a function called \"pendufuc\" which we created seperately , which is :

1. Program :

2. function [dtheta_dt] = pendufunc(t,theta,b,g,l,m)
   
   theta1 = theta(1);
   theta2 = theta(2);
    
   % two first order ODE 
   dtheta1_dt = theta2;
   dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1);
   
   %clubbing them together
   dtheta_dt = [dtheta1_dt; dtheta2_dt];
   
   end
   

3. The function command is used to simplify the program , ODE is solved here and when we call it in the main program assigning its input parameters , it provides the output .
4. The two first order derivative solutions are written in the function code , and theta1 and theta2 values are assigned to theta(1) and theta(2) respectively to call them in the main function command using just \"theta\" .

1. Now in the main program , we plot time vs displacement and time vs angular velocity :
2.                         
3. Afterwards we see how the plot between Displacement and angular velocity looks like from the governing differential equation, the behavoir of these variables is :
4.                        

• For loop 

1. For loop is used to run the loop for simple pendulum at various angles of displacement for different values of time in the ODE . 
2. initial x0 y0 conditions are set and then x1 y1 conditions are set in minus for the direction of pendulum desired is below the x axis , vertical orientation .
3. Then the rigid base at the top is plotted from -1 to 1 and y = 0 , on x axis  , then string is plotted connecting from the fixed support to the mass of the ball at x2 y2.
4. Then marker command is used to plot the mass in the shape of a circle at position x2 y2 , a random color code is picked for it. 
5. All the figures are plotted using subplot command at first to orient all of them in a common main figure.
6.                            
7. Then the axis is specified and pause command is used . 
8. The frames of various positions of the pendulum is stored in ct = 1 from then on ct = ct + 1; which is called in getframe command to call all the plotted frames which can be displayed with the given value of pause making it appear like a movie. 

• Movie making

1. It is stored in variable "M" . 
2. videofile = VideoWriter('simple_pendulum.avi','uncompressed AVI') command is used to write the video and its name and type of extension needed is provided.
3.  Then its opened and the movie is written into it and then stored in the directory where this program is saved . The movie making code is self explanatory .

• Output : 

 END