SOLVING SECOND ORDER ODE : SIMPLE PENDULUM

AIM:

A PROGRAM IN OCTAVE TO SOLVE SECOND ORDER ODE FOR STIMULATING TRANSIENT BEHAVIOUR OF SIMPLE PENDULUM AND CREATE ANIMATION OF IT’S MOTION

INTRODUCTION:

ODE (Ordinary Differential Equation) helps us to describe the transient behavior of a system. Simple Pendulum is the best example to get second order ODE i.e. from Newton’s second law equation which describes the position of bob with respect to time.

GOVERNING EQUATIONS:

Fig.1. Simple Pendulum

The ordinary differential equation representing the motion of the pendulum with damping -

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

Where,

g = acceleration due to gravity (m/s^2)

l = length of pendulum (m)

m = mass of the ball/ bob (kg)

b = damping coefficient

For solving ODE in octave, the function of ‘ode45/lsode’ is used.

To solve second order ODEs in OCTAVE / MATLAB, we have to convert it into two first order ODEs.

OBJECTIVES:

1) TO STUDY THE TRANSIENT BEHAVIOUR OF SIMPLE PENDULUM

2) TO CREATE ANIMATION GIF OF PENDULUM MOTION

Let us assume the values of pendulum,

l = 1m

g = 9.81 m/s^2

b = 0.05

m = 1 Kg

Common commands used in starting of every program-

Clear all command for clearing workspace

Close all command for closing plots and figures in background

Clc command for clearing the command window

Figure1.1. Flowchart for animation

PROGRAM:

OCTAVE PROGRAM1- FUNCTION PROGRAM

% Defining function
function [dtheta_dt] = ode_func(t,theta,b,g,l,m)
  
  theta1 = theta(1);
  theta2 = theta(2);
  
% Second order ODE is converted into two single ODE
  dtheta1_dt = theta2;
  dtheta2_dt = -(b/m)*theta2-(g/l)*sin(theta1);
  
% Array with two single ODE
  dtheta_dt = [dtheta1_dt;dtheta2_dt];
endfunction

Explanation:

Line1- Function is created that define the input variables (t,theta,b,g,l,m)

Line2&3- The solution variable theta(1)&(2)

Line4- Conversion into first single ODE

Line5- Conversion into second single ODE

Line6- Making array with two single ODEs

Line7- Closing the function

OCTAVE PROGRAM2- MAIN FUNCTION

clear all
close all
clc

% Input parameter variables
b = 0.05;
g = 9.81;
l = 1;
m = 1;

% Specifying initial conditions
theta = [0;3];

% Specifying Time domain
t_span = linspace(0,50,500);

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

% Plotting
figure(1);
plot(t,results(:,1));
hold on;
plot(t,results(:,2));
grid on; 
grid minor;

% Titles and labels
title('TRANSIENT BEHAVIOUR')
xlabel('TIME')
ylabel('DISPLACEMENT & VELOCITY')
legend('Displacement','Velocity')

% Animating the motion
ct = 1;
disp = results(:,1);
% Iteration
for i = 1:length(disp)
  RESULTS = results(i);
  % Coordinates of start point
  x0 = 0;
  y0 = 0;
  % Coordinates of end point
  x1 = l*sin(RESULTS);
  y1 = -l*cos(RESULTS);
  % Plotting
   figure(2);
  % String
  plot([x0 x1],[y0 y1],'linewidth',5,'color','b');
  hold on;
  % Ceiling
  plot([-1.5 1.5],[0 0],'linewidth',10,'color','g');
  % Bob
  plot(x1,y1,'marker','o','markersize',15,'markerfacecolor','r');
  axis([-1.5 1.5 -1.5 1]);
  title('SIMPLE PENDULUM MOTION');
  pause(0.05);
  
  % Creating gif animation 
  t_span = 0.1 
  printf('creating GIF-Progress...%d/%d\n',i,length(t_span)); 
  img = print ('-RGBImage'); 
  imwrite(img,'simple pendulum.gif','DelayTime',0.005,'Compression','bzip','Writemode','Append'); 
  ct = ct + 1;   
endfor

Explanation:

Line1- Defining damping coefficient value

Line2- Value defined for acceleration due to gravity

Line3- Defining length of the pendulum string

Line4- Defining mass of bob

Line5- Defining initial value for angular displacement and angular velocity

Line6- Time array for time period

Line7- Function to solve 2nd order ODE, we use ‘ode45’.

Line8- The figure(1) will represent the transient behaviour of pendulum i.e. representing the solved equation.

Line9- Plotting displacement vs time first column

Line10- Retaining the data from first plot

Line11- Plotting displacement vs time second column

Line12- Showing grids on the graph

Line13- For making the grid lines neat and small

Line14- Putting the title on the first plot showing the transient behaviour

Line15- Labelling x-axis

Line16- Labelling y-axis

Line17- Defining the lines on the graph

Line18- Counter with starting value assigned as 1

Line19- Assigning the first column result vector

Line20- Running the for-loop for all values of ‘disp’

Line21- Assigning all the displacement values in the first column

Line22- Coordinates of initial positions, x-axis=0

Line23- Coordinates of initial positions, y-axis=0

Line24- Coordinates from trigonometry for x-axis

Line25- Coordinates from trigonometry for y-axis

Line26- Opening new window for showing motion of pendulum

Line27- Plot representing pendulum’s string

Line28- Retaining the data in second plot

Line29- Plot representing the base or ceiling from which the pendulum in hanged

Line30- Plot representing the bob

Line31- Giving title to second plot

Line32- Pause command is used to cease the animation from jumping between the figures.

Line32- A new plot is obtained for each value of t_span, i.e. 0.1

Line33- The statement helps us to keep the track of GIF creation progress

Here, printf is used to print particular string in command window

%d is format specifier

\n creates a new line every time it gets executed

i is a loop variable

length() returns total number of elements inside an array

Line34- The print command saves the most recent plot generated

Line35- Here, imwrite()- This is used to write images in various file format

img- The image to be written into the gif file is stored inside img

pendulum.gif- Name of the output gif file

DelayTime- This is used to specify the delay between each frame

Compression- The raw image is a big file hence we need to compress to bring each image into minimum size as it is attached to gif file.

Writemode- Specifies if img is appended to the file

Append- Used to attach images one by one into gif file

Line36- Loop counter value increased by +1

Line37- Closing the for-loop

Figure1.2. Damped motion of pendulum

Reducing the damping coefficient factor by 10 i.e. b= 0.5

Figure1.3. Undamped motion of pendulum

Figure1.4. GIF creation progress

ANIMATION LINK:https://drive.google.com/file/d/1fKFG55cs4ajP4N-rvFREdiFpHNaPWk4H/view?usp=sharing

CONCLUSION: SECOND OR HIGHER ORDER DIFFERENTIAL EQUATION ARE EASILY SOLVED ON PROGRAMMING PLATFORMS LIKE MATLAB GIVING US THE CLEAR UNDERSTANDING ABOUT THE HARMONIC MOTION OF SIMPLE PENDULUM.