A Simple pendulum motion by solving second order ODE (Ordinary Differential Equation).

1. Aim: To simulate the transient behaviour of a system (Simple Pendulum) using ODE.

2. Objectives:

• To solve the ODE representing the motion of a simple pendulum.
• To create an animation video of the simple pendulum motion.

3. Theory:

1. In mathematics, ODE is known as the Ordinary Differential Equation. This equation contains one or more functions of one inpendent variable and their derivatives of those functions. If there are more than one independent variables then the differential equation is defined as partial differential equation.
2. Let y be the independent variable and x be the dependent variable. So, the function is defined as y = f(x).
3. The motion of a simple pendulum is simulated with the help of second order ODE. This motion is generally based on the Newton's second of motion. The following second order ODE gives the different positions of the ball with respect to time:

`(ⅆ^2 θ)/(ⅆt^2 )+b/m(ⅆ θ)/(ⅆt) +g/l sinθ=0` ...equation 1

where θ = angular displacement

`(ⅆ θ)/(ⅆt)` = angular velocity

iv. The above ODE is break down to two first order ODE as first for the angular displacement (rad) and second for the angular velocity (rad/s). The initial conditions are:

• At t = 0, θ = 0
• At t = 0, `(ⅆ θ)/(ⅆt)` = 3

v. The following parameters are considered for the simple pendulum:

• Length of the pendulum = l = 1 metre
• Mass of the ball = m = 1 kg
• Damping coefficient = b = 0.05
• Gravitional acceleration = g = 9.81 m/s²

vi. The following state space method is used to solve the second order ODE:

Let θ = θ₁, so `(ⅆ θ)/(ⅆt)` = `(ⅆ θ1)/(ⅆt)` and `(ⅆ θ1)/(ⅆt)` = θ₂ so `(ⅆ^2 θ)/(ⅆt^2)` = `(ⅆ θ2)/(ⅆt)`

Then the equation 1 reduces to:

`(ⅆ θ2)/(ⅆt) + b/mθ2 +g/l sinθ1=0` ....equation 2

In matrix form,

`d/dt((θ1),(θ2)) = ((θ2),(-b/m θ2-g/l sinθ1))`

vii. The ode45 solver is used to solve ODE. This solver gives the angular displacement and angular velocity plots with respect to time. The following is the syntax of the ode45:

ode45@(variable1,variable2)

viii. The motion of the simple pendulum is simulated for 20 seconds.

ix. The position of the ball is defined by:

• On X - axis -> l*sinθ
• On Y - axis -> - l*cosθ (negative sign as per sign convention - `↑` as positive and `↓` as negative)

4. MATLAB function:

First ode_func named function is defined to solve the second order ODE as explained in the above method and the main program is descripted. The following is the MATLAB program of the ode_func function:

% creating an ode function to solve the second order differential equation
function [dtheta_dt] = ode_func(t, theta, b, g, l, m)
% theta(1) stores the different values of displacement with respect to
% time and storing these variables in the theta1 variable
theta1 = theta(1)

% theta(2) stores the different values of velocity with respect to
% time and storing these variables in the theta2 variable
theta2 = theta(2)
dtheta1_dt = theta2;

% considering the angle is in radians
dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1);

% dtheta_dt is the column matrix 
dtheta_dt = [dtheta1_dt;dtheta2_dt];

end

5. Main MATLAB Program:

The following is the MATLAB program which creates an animation video on the motion of the simple pendulum and plots the graphs of the angular displacement and angular velocity with respect to time. The program contains the code explanation and the corresponding comments.

% A program to solve the second order differential eqaution and creating
% an animation of a simple pendulum

% clearing the workspace 
clear all

% closing the plots
close all

% clearing the command window
clc

% inputs
% damping coefficient
b = 0.05;

% gravitational acceleration (m/s^2)
g = 9.81;

% length of the pendulum (m)
l = 1;

% mass (kg)
m = 1;

% initial conditions
% The system is of second order so making it into two first order ODEs - velocity and displacement
% initial condition - at t = 0 angular displacement is 0  and angular velocity is 3
theta_0 = [0;3];

% as integration is performed to solve the ODE, time array is defined in sec
% from 0 to 20 with 200 values in between
t_span = linspace(0,20,200);

% solve ODE
% t and theta are the main inputs as these are the dependent and
% independent variables of the ODE system
% @ is the symbol of the function 
[t, results] = ode45(@(t,theta) ode_func(t, theta, b, g, l, m), t_span, theta_0);

% plotting the time and displacement as results is a matrix of displacement
% and velocity
plot(t, results(:,1))

% holding the plot
hold on

% plotting the time and velocity
plot(t, results(:,2))

% labelling the curves
legend({'angular displacement', 'angular velocity'},'Location','Southeast')

% labelling X - axis
xlabel('time')

% labelling Y - axis
ylabel('plot')

% turning on the grid
grid on

% holding off the plot - figure 1
hold off

% opening the window of figure 2
figure(2)

% setting the origin 
origin = [0 0];

% defining a variable ct as 1 to capture each frame of the simple pendulum
ct = 1;

% using for loop to get the different frames varying with time 
for i = 1:length(t)
    
    % setting the length of the pendulum
    P = l*[sin(results(i,1)) -cos(results(i,1))];
   
    % plotting the hook on which pendulum is mounted, having red colour and
    % linewidth as 
    support = plot(origin(1), origin(2), 'O','color','r','linewidth',2);
    
    % holding the plot
    hold on
    
    % To construct the block as the support
    % for horizontal lines of the support
    x1 = [-1 1];
    y1 = [0 0];
    y2 = [0.1 0.1];
    y3 = [0.4 0.4];
    y4 = [0.7 0.7];
    y5 = [1 1];

    % for vertical lines of the support
    x2 = [-0.5 -0.5];
    x3 = [-0.3 -0.3];
    x4 = [0 0];
    x5 = [0.1 0.1];
    x6 = [0.3 0.3];
    x7 = [0.4 0.4];
    x8 = [0.6 0.6];
    x9 = [0.8 0.8];

    y6 = [0 0.1];
    y7 = [0.1 0.4];
    y8 = [0.4 0.7];
    y9 = [0.7 1];

    % for the central axis
    x10 = [0 0];
    y10 = [1.5 -1.5];
    
    % plotting the horizontal lines having black colour
    line (x1,y1,'color', 'black')
    line (x1,y2,'color', 'black')
    line (x1,y3,'color', 'black')
    line (x1,y4,'color', 'black')
    line (x1,y5,'color', 'black')

    % plotting the vertical lines having black colour
    line (x2,y6,'color', 'black')
    line (x4,y6,'color', 'black')
    line (x6,y6,'color', 'black')
    line (x8,y6,'color', 'black')

    line (x3,y7,'color', 'black')
    line (x5,y7,'color', 'black')
    line (x7,y7,'color', 'black')
    line (x9,y7,'color', 'black')

    line (x2,y8,'color', 'black')
    line (x4,y8,'color', 'black')
    line (x6,y8,'color', 'black')
    line (x8,y8,'color', 'black')

    line (x3,y9,'color', 'black')
    line (x5,y9,'color', 'black')
    line (x7,y9,'color', 'black')
    line (x9,y9,'color', 'black')

    % plotting the central line
    line (x10,y10,'LineStyle', '--','color', 'black')
    
    % plotting the link of the simple pendulum using black colour and
    % linewidth as 4 
    pend = line ([origin(1) P(1)], [origin(2) P(2)], 'color','black','linewidth',4);
    
    % plotting the ball either by using MarkerSize or increasing the
    % linewidth to 15 and having cyan colour
    % plot(P(1),P(2), '-O', 'MarkerSize', 10)
    plot(P(1),P(2), 'O','color','cyan', 'linewidth', 15);
    
    % plotting the moving central axis along with the motion of the
    % pendulum, having blue colour
    line ([-P(1) P(1)],[-P(2) P(2)],'LineStyle', '--','color', 'b');
    
    % holding off the plot
    hold off
    
    % setting the limits for the X - axis as -3 to 3 and that of Y - axis
    % as -2 to 2
    axis([-2 2 -2 2])
    
    % labelling X - axis
    xlabel('angular displacement on X - axis w.r.t time')

    % labelling Y - axis
    ylabel('angular displacement on Y - axis w.r.t time')
    
    % turning the grid on
    grid on
    
    % pausing the plot to see each frame
    pause(0.02)
  
   % to get the individual frame we create an array M using getframe
   % command
    M(ct)= getframe(gcf);
    ct = ct + 1;
end

%to create a movie using movie command
movie(M)

%to create a videofile using command VideoWriter with file name as
%simple_pend.avi and profile as Uncompressed AVI
videofile = VideoWriter('simple_pend.avi', 'Uncompressed AVI');

%to open the file
open(videofile)

%to wirte the video file of the frames M using writeVideo command
writeVideo(videofile,M)

%to close the videofile
close(videofile)

6. Results:

The outputs of the above MATLAB program are:

• Graph of angular displacement, angular velocity Vs time - It shows that both the angular displacement and angular velocity is damped (decreases) with time. More the damping coefficient, more will be the damping.

• Graph of the position of the ball with respect to time - The position of the ball is changed with respect to time and it comes to rest after few seconds, as damping is provided.

The animation video of a simple pendulum can be seen by the following video:

simple_pend.avi

7. Conclusion:

There were no errors encountered during the MATLAB program, just different trials were taken to watch the motion of the simple pendulum and also to build the block support. An animation video showing the motion of a simple pendulum was successfully completed.

8. References:

1. ode45 - https://www.mathworks.com/help/matlab/ref/ode45.html
2. getframe - https://www.mathworks.com/help/matlab/ref/getframe.html
3. movie - https://www.mathworks.com/help/matlab/ref/movie.html
4. legend - https://www.mathworks.com/help/matlab/ref/legend.html