SECOND ORDER ODE - MOTION OF SIMPLE PENDULUM

OBJECTIVE:

To solve the equation of motion of simple pendulum which is a second order ODE and to simulate the motion of the simple pendulum with damping

INTRODUCTION:

Differential equations appear in a wide range of problems in physics and science. They are use to model many of the phenomenon that occur around us and we need to solve differential equations to analyse those phenomena. Some simple examples of modelling using differential equations include:

1. Motion of spring-mass systems and simple harmonic motion.
2. Mixture of substances and finding concentration.
3. Newton's law of cooling.
4. Population growth analysis.
5. Dynamics: finding acceleration, velocity and position of a body.
6. Analysis of current in capacitive resistive circuits and inductive-resistive circuits

Simple Pendulum:

A simple pendulum consists of a large mass (bob) suspended by a massless rod or a light string from a rigid support as shown in below figure. The length of the string is large compared to the dimensions of the bob. A stationary pendulum has its string in a vertical position with the bob hanging straight down. This is known as the pendulum's equilibrium position. Oscillations (to and fro motion) are set up in a pendulum by displacing the bob from its equilibrium positions and releasing it.

The point of pivot is taken as `(X_0, Y_0)` and Bob coordinates are taken as `(X_1, Y_1)`

INPUT DATA:

1. Length of Pendulum = `L` = 1 m

2. Mass of the Bob = `m` = 1 kg

3. Acceleration due to gravity = `g` = 9.81 `m/s^2`

4. Position of Pivot is at Origin.

`X_0` = 0

`Y_0` = 0

5. Different values of damping coefficient are taken in the problem to identify the effect of damping on simple pendulum

Damping Coefficient = `b` = 0.05, 0.1, 0.2, 0.3, 0.4, 0.5

6. Time span is taken as 20 seconds

7. Initial Angular Displacement = 0 radians

8. Initial Velocity of Pendulum = 3 m/s

CALCULATION FORMULAE:

i. Position of pendulum bob at a particular angle `theta` is given as

`X_1 = L Sintheta`

`Y_1 = - L Costheta`

ii. Equation of motion of simple pendulum with damping is given by

`(d^2 theta)/dt^2 + b/m (d theta)/dt + g/L sin theta = 0` ------(1)

iii. Conversion of Second order ODE in to system of First order ODEs

Let's consider

`theta = theta_1` ------(2)

`(d theta)/dt = (d theta_1)/dt = theta_2` ------(3)

which gives

`(d^2 theta)/dt^2 = d/dt((d theta)/dt) = d/dt(theta_2) = (d theta_2)/dt` ------(4)

Substituting above conditions in equation of simple pendulum with damping

`(d^2 theta_1)/dt^2 + b/m (d theta_1)/dt + g/l sin theta_1 = 0`

`(d theta_2)/dt + b/m theta_2 + g/l sin theta_1 = 0`

`(d theta_2)/dt = -b/m theta_2 - g/l sin theta_1` ------(5)

by combining ewuations (3) and (5) in to matrix form

`d/dt [(theta_1), (theta_2)] = [(theta_2),(- b/m theta_2 - g/l sin theta_1)]` ------(6)

CODE ALGORITHM or EXPLANATION:

1. Initial input parameters are assigned. Time span is divided in to particular number of divisions using linspace command.

2. Set of linear first order ODEs (Eq. 6) are solved using 'ode45' sover.

3. Results are plotted between Time and Angular Displacement and Time and Velocity using plot command.

4. In order to simulate the motion of simple pendulum, a 'for' loop is started to determine end positions of pendulum rod (`X_1, Y_1)` at each angle of inclination of pendulum.

`X_1 = L Sintheta`

`Y_1 = - L Costheta`

5. In order to get animated movie, first plot command is used to plot the position of pendulum rod, second plot command is used to plot pendulum bob at `(X_1, Y_1)`and last plot is used to plot the surface to which pendulum bob is suspended. This plot frame is stored in to array M using getframe command

% Plotting of Pendulum link
plot([x0 x1], [y0 y1], 'linewidth', 4, 'color', 'b');

% Plotting of Pendulum Bob
plot(x1, y1, 'o', 'MarkerSize', 10, 'MarkerfaceColor', 'r', 'MarkerEdgeColor', 'r', 'linewidth', 4);

% Plotting surface where Pendulum Bob is suspended
plot([-0.2 0.2], [0 0], 'linewidth', 5, 'color', 'k');

M(count) = getframe(gcf);

gcf returns the current frame handle and count varible is increased by one value

5. 'for' loop is ended and syntax for creating a movie file using VideoWriter command is used.

'ode45' solver syntax:

'ode45' is the most common ODE solver in MATLAB. It is a six-stage, fifth-order, Runge-Kutta method

[t,y] = ode45(odefun,tspan,y0),

where tspan = [`t_0` `t_f`], integrates the system of differential equations `y' = f(t,y)` from `t_0` to `t_f` with initial conditions `y_0`.

Each row in the solution array `y` corresponds to a value returned in column vector `t`.

Here we are solving two First Order ODEs simultaneously. Hence 'results' is an array of [Angular Displacement, Velocity]

[t, results] = ode45(@(t, theta) ode_func(t,theta,b,g,l,m), t_span, theta_0);

where

ode_func(t,theta,b,g,l,m) - system of first order ODE in Eq. 6

t_span - time span

theta_0 - Initial conditions of angular displacement and velocity

Upon solving ode_func for each value of time `t`, angular displacement and velocity are stored in 'results' array

Syntax for creating a movie using VideoWriter from an array M:

Command for VideoWriter

v = VideoWriter(filename,profile) or VideoWriter(filename)

open(v) % Open the video file 'v'
writeVideo(v,M) % Writes data from an array M to the video file 'v'
close(v) % Close the video file 'v'

Different profiles for VideoWriter

'Archival' - Motion JPEG 2000 file with lossless compression
'Motion JPEG AVI' - AVI file using Motion JPEG encoding
'Motion JPEG 2000' - Motion JPEG 2000 file
'MPEG-4' - MPEG-4 file with H.264 encoding
'Uncompressed AVI' - Uncompressed AVI file with RGB24 video
'Indexed AVI' - Uncompressed AVI file with indexed video
'Grayscale AVI' - Uncompressed AVI file with grayscale video

CODE:

% Second Order ODE (Simple Pendulum) Equation 
clear all
close all 
clc

% Input Parameters
b = 0.05;        % Damping Coefficient
l = 1;          % Length of the Pendulum in Meters
g = 9.81;       % Acceleration due to gravity in m/s^2
m = 1;          % Mass of the Pendulum in Kg

% Initial Conditions
theta_0 = [0; 3];   % [Angular Displacement, Velocity] & Angular displacement is taken in radians

% Time limit is divided in to particular number of divisions
t_span = linspace(0,20,500);

% Most commom ODE solver in MATLAB is 'ode45' solver
% 'ode45' is a six-stage, fifth-order, Runge-Kutta method
% [t,y] = ode45(odefun,tspan,y0), where tspan = [t0 tf], integrates the system of differential 
% equations y'=f(t,y) from t0 to tf with initial conditions y0. 
% Each row in the solution array y corresponds to a value returned in column vector t.

% Here we are solving two First Order ODEs simultaneously
% Hence 'results' is an array of [Angular Displacement, Velocity] 
[t, results] = ode45(@(t, theta) ode_func(t,theta,b,g,l,m), t_span, theta_0);

% Plotting
figure(1);
% Time vs Angular Displacement
plot(t,results(:,1), 'linewidth', 2, 'color', 'r')
hold on
% Time vs Velocity
plot(t,results(:,2), 'linewidth', 2, 'color', 'b')
axis([-0.5 20.5 -3.5 3.5])
xlabel('Time')
legend ('Angular Displacement', 'Velocity')
title('When damping coefficient = 0.05')

% Storing angular displacement values from Array 'results' in to seperate Array
theta = results(:,1);
% THETA = theta*(180/pi);
count = 1;

% Loop for creation of Animation
for i = 1:length(theta)
    
    % Position where Pendulum os suspended
    x0 = 0;
    y0 = 0;
    
    % Position of Pendulum Bob at angle theta(i)
    % theta(i) is in radians. Hence sin() and cos() are used instead of sind() and cosd() 
    x1 = l*sin(theta(i));   
    y1 = -l*cos(theta(i));
    
    % Plotting
    figure(2)
    % Plotting of Pendulum link
    plot([x0 x1], [y0 y1], 'linewidth', 4, 'color', 'b');
    hold on
    % Plotting of Pendulum Bob
    plot(x1, y1, 'o', 'MarkerSize', 10, 'MarkerfaceColor', 'r', 'MarkerEdgeColor', 'r', 'linewidth', 4);
    hold on
    % Plotting surface where Pendulum Bob is suspended
    plot([-0.2 0.2], [0 0], 'linewidth', 5, 'color', 'k');
    axis([-1.5 1.5 -1.5 0.2])
    title('Animation of Simple Pendulum (b = 0.05)');
    hold off
    
    % Storing the current frame into array M using getframe command
    % gcf returns current figure handle
    M(count) = getframe(gcf);
    count = count + 1;
end

% plays the movie frames in array M once & to play 'n' times use movie(M,n)
movie(M);

%%%%%%%%%%%%%%%%%%%%%%%%%% CREATING VIDEO FILE %%%%%%%%%%%%%%%%%%%%%%%%%

% Using VideoWriter object to create a video file from an array M
% Command for VideoWriter
% v = VideoWriter(filename,profile) or VideoWriter(filename)
% 
%     Different profiles for VideoWriter
%     'Archival' - Motion JPEG 2000 file with lossless compression
%     'Motion JPEG AVI' - AVI file using Motion JPEG encoding
%     'Motion JPEG 2000' - Motion JPEG 2000 file
%     'MPEG-4' - MPEG-4 file with H.264 encoding
%     'Uncompressed AVI' - Uncompressed AVI file with RGB24 video
%     'Indexed AVI' - Uncompressed AVI file with indexed video
%     'Grayscale AVI' - Uncompressed AVI file with grayscale video
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
v = VideoWriter('simple_pendulum.avi','Uncompressed AVI');
open(v) % Open the video file 'v'
writeVideo(v,M) % Writes data from an array M to the video file 'v'
close(v) % Close the video file 'v'

% ODE FUNCTION FOR SIMPLE PENDULUM
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

OUTPUT:

Time vs Angular Displacement and Velocity plots for various damping coefficients are shown below

It can be observed from above graphs that as the value of damping increases, the values of angular displacement and velocity values decreases with time. For higher damping values, simple pendulum comes to rest or equilibrium position quickly.

The final output animation of motion of the simple pendulum for different damping coefficients can be seen in the following links

1. When damping coefficient (`b`) = 0.05; https://youtu.be/UU5ev7lpC1g

2. When damping coefficient (`b`) = 0.1; https://youtu.be/d1Po5t3AuhQ

3. When damping coefficient (`b`) = 0.2; https://youtu.be/G83mfrFHukw

4. When damping coefficient (`b`) = 0.3; https://youtu.be/cCWtq-n484I

5. When damping coefficient (`b`) = 0.4; https://youtu.be/oxtycbPdNvU

6. When damping coefficient (`b`) = 0.5; https://youtu.be/-g2FJcS76_8

ERRORS:

No major errors has encountered in the program