Solving second order ODE: equation of motion of a simple pendulum using MATLAB

AIM

To solve the given ODE i.e., the equation of motion of a simple pendulum using MATLAB

OBJECTIVE

• To plot the variation of displacement and angular velocity with time
• To get the animation of the motion of the pendulum

THEORY

The equation of motion for a simple pendulum which is a second order ordinary differential equation (ODE) is given by,

`(d^2 θ)/(dt^2 )+b/m (dθ)/dt+g/L sin⁡θ=0`

Where,

`b=` Damping coefficient

`m=` Mass of the ball

`g=` Acceleration due to gravity

`L=` Length of the pendulum

Let consider,

`θ=θ_1` and

`(dθ)/dt=(dθ_1)/dt=θ_2` --- (1)

Therefore,

`(d^2 θ)/(dt^2 )=(d^2 θ_1)/(dt^2 )=(dθ_2)/dt`

Thus the equation of motion becomes,

`(dθ_2)/dt+b/m θ_2+g/L sin⁡θ_1=0`

The above equation is the representation of the equation of motion as first order ODE. Now,

`(dθ_2)/dt=-b/m θ_2-g/L sin⁡θ_1` --- (2)

From equation (1) and (2) we get the following matrix,

`d/dt [(θ_1),(θ_2 )]=[(θ_2),(-b/m θ_2-g/L sinθ_1 )]`

The above equation will give us the value of angular displacement and angular velocity.

Software used: - MATLAB R2014a

PROGRAM

• Program to define the function ‘ode_func’

  %Program for Function definition 
  
  function [dtheta_dt] = ode_func(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

• Program to create a plot the variation of displacement and angular velocity with time and animate the motion of the pendulum

% Program to create a plot the variation of Angular Displacement and Angular Velocity with time 
% and animate the motion of the pendulum

clear all 
close all
clc


% Inputs

b = 0.05;                             % Assigning the value of damping coefficient
g = 9.81;                             % Assigning the value of Acceleration due to gravity (Unit- metre/second^2)
L = 1;                                % Assigning the value of length of the pendulum (Unit- metre)
m = 1;                                % Assigning the value of mass of the bob (Unit- kilogram)

% Initial Condition
theta_0 = [0;3];                      % Assigning the value of initial angular displacement and angular velocity

% Time span
t_span = linspace(0,20,500);          % Assigning the value of time span

% Solving ODE

% Calling the function 'ode_func'
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,L,m),t_span,theta_0);
% [t,results] = Output arguments: time, angular displacement, angular velocity
% ode45 = Function keyword to solve second order ODE
% (@(t,theta) = The ODE is function of Time(t) and angular displacement(?)
% (t,theta,b,g,L,m) = Input arguments of 'ode_func'
 
% Plotting
plot(t,results(:,1),'color','r')              % To plot Time vs Angular Displacement
hold on
plot(t,results(:,2),'color','b')              % To plot Time vs Angular Velocity
xlabel('Time')                                % To label the x-axis as 'Time'
ylabel('Plot')                                % To label the y-axis as 'Plot'
legend('Angular Displacement','Angular Velocity')     % To give annotation to the respective curves
title('Variation of angular displacement and angular velocity with time')

% Iteration using for loop to create animation of simple pendulum 

ct = 1;                                       %Counter initiation     

for i =1:length(results(:,1))
    
    % Starting point of string at base
    x0 = 0;
    y0 = 0;
    
    % Ending point of string at bob
    x1 = L*sin(results(i,1));
    y1 = -L*cos(results(i,1));
    
    figure(2)
    plot([-1 1],[0 0], 'linewidth',4,'color','k')                % To create the base of the Pendulum system                
    hold on
    plot([x0 x1],[y0 y1],'linewidth',2,'color','b')              % To create the string of the Pendulum system                        
    hold on 
    plot(x1,y1,'o','markersize',20,'MarkerFaceColor','r')        % To create the bob of the Pendulum system  
    axis([-2 2 -2 1.5]);                                         % Assigning the range of axis to be shown in animation       
    pause(0.0003)                                                % Assigning command to pause each frame
    hold off
    
    M(ct) = getframe(gcf);                                       % Assigning command to get current figure into an array called 'M'
    ct = ct+1;                                                   % Assigning command to run the loop for each counter
    
end

movie(M)                                                         % Assigning command to generate a movie with all the frames
videofile = VideoWriter('Pendulum.avi','Uncompressed AVI');      % Assigning command to create a AVI file to store the movie
open(videofile)                                                  % Assigning command to open the video file
writeVideo(videofile,M)                                          % Assigning command to write the frames in the videofile
close(videofile)                                                 % Assigning command to close the videofile

RESULTS

Plot

Figure: Variation of Angular Displacement and Angular Velocity with time

Animation

Animation of the motion of simple pendulum

Link: https://youtu.be/fvdACj3fHq8

CONCLUSION

• At first the equation of motion of a simple pendulum which is a second order ODE is given. After that the second order ODE is broken down to first order ODE for MATLAB program analysis.
• A function ‘ode_func’ is created for solving the ODE.
• After derivation of Angular Displacement and Angular Velocity values, the variation of those parameters with time is plotted in same graph.
• At the end, the animation of motion of simple pendulum is created as AVI movie file.