SIMULATION OF MOTION OF A SIMPLE PENDULUM USING OCTAVE

SIMULATION OF MOTION OF A SIMPLE PENDULUM USING OCTAVE

AIM:

To simulate the motion of a simple pendulum when the initial condition of angular displacement and angular velocity are provided.

INTRODUCTION:

A simple pendulum consists of a small mass suspended by a string from a pivot so that it can swing freely. When a pendulum is displaced from its equilibrium position, it is subject to a restoring force due to gravity such that it will accelerate back towards the equilibrium position.

When released, the restoring force acting on the pendulum mass causes it to accelerate back and forth about its mean position, thereby constituting a periodic motion. However, this is an ideal scenario. Generally factors such as air resistance, frictional losses etc. cause the amplitude of the simple pendulum to decrease continuously. For such a damped oscillation, equation of motion is given by the equation below:

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

Where theta is the angular displacement of the pendulum from its mean position.

g is the acceleration due to gravity (9.81 m/s)

m is the mass of the bob in kg

L is the length of the string in m

b is the damping coefficient in Ns/m

By solving the above ODE, the value of ‘theta’ and ‘d(theta)/(dt)’ can be obtained (i.e. instantaneous values of angular displacement/amplitude and angular velocity can be obtained)

ASSUMPTIONS:

• Motion of the pendulum is considered to be in one plane.
• No external force is acting on the pendulum once the oscillation begins.
• The string is assumed to be massless.

PROCESS:

The code is divided into 2 parts.

PART 1: Defining an ode function which takes the input values of mass of the bob, length of string, damping coefficient, acceleration due to gravity, angular displacement and time.

It creates and returns an angular velocity matrix, i.e. a second order differential equation is converted into 2 first order differential equations with the help of this function.

The given differential equation is:

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

Let `theta=theta1`

Then, `(d(theta))/(dt) =(d(theta1))/(dt) =theta2` ---------- eqn 1

Therefore, `(d^2theta)/(dt^2) = (d^2theta1)/(dt^2) =d/(dt)*(theta2) ` --------- eqn 2

On re-writing the eqns,

`(d(theta1))/(dt) =theta2` ------eqn 3

`(d^2theta1)/(dt^2) =(d(theta2))/(dt)`

`(d(theta2))/(dt)=-(b/m)*(theta2)=(g/L)*sin(theta1)` ------eqn 4

Hence, From eqn 3 and 4, we get

`d/(dt)[[theta1],[theta 2]]=[[theta2],[-(b/m)*theta2-(g/L)sin(theta1]]]` ---------eqn 5

The function generates this dtheta/dt matrix and returns it to the main program.

CODE: PART 1:

function [dtheta_dt] = odefunc(t,theta,b,m,g,l)
  theta1=theta(1)
  theta2=theta(2)
  dtheta1_dt=theta2
  dtheta2_dt=-(b/m)*theta2-(g/l)*sin(theta1)
  dtheta_dt=[dtheta1_dt dtheta2_dt]
endfunction

PART 2: This is the main program.

Here, the values of all inputs such as mass of the bob, length of string, damping coefficient, acceleration due to gravity, angular displacement and time are provided. Timespan is assumed to be 20 sec. The ode45 function is used to solve the differential equations to provide angular displacement and angular velocity matrix. The initial conditions (i.e. values at time t=0) of angular displacement and angular velocity are assumed to be 0 deg and 3rad/sec respectively. The code calculates and plots the values of angular displacement and velocity as a function of time.

In order to animate the motion of the pendulum, the instantaneous position of the bob is obtained by using the formula,

`x_theta=lsin(theta)`

`y_theta=-lcos(theta)`

CODE: PART 2:

% simple pendulum ode

clear all
close all
clc

%inputs for the program are g,l,b,m and time span
% g = acceleration due to gravity (m/s^2)
% m = mass of bob in kg
% b = damping co-efficient 
% l = length of pendulum in m
% t = time in sec

g=9.81
l=1
m=1
b=0.05

%initial conditions for displacement and velocity
theta_0=[0;3]

%time span 20 sec
timespan=linspace(0,20,200);

%ode solving using ode45 function by calling 'odefunc' function
[t ,results] = ode45(@(t,theta) odefunc(t,theta,b,m,g,l),timespan,theta_0)

%plotting the angular displacement and angular velocity plots w.r.t time
subplot(2,2,1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
xlabel('Time')
title('Plotting w.r.t time')
legend('\theta','\omega')
hold on

%plotting the angular displacement vs angular velocity plots 
subplot(2,2,2)
plot(results(:,1),results(:,2))
xlabel('Angular Displacement in rad')
ylabel('Angular velocity in rad/sec')
title('angular displacement vs angular velocity')
hold on

%Animation of the pendulum
leng=length(results(:,1))
theta_array=results(:,1)
ct=1


for i=1:leng
  x_e=l*sin(theta_array(i))
  y_e=-(l*cos(theta_array(i)))
  subplot(2,2,[3,4])
  %Plotting the string, bob and support
  plot([0 x_e],[0 y_e],'linewidth',3,'color','black',x_e,y_e,'.','markersize',40,'linewidth',3,'color','black',[-1 1],[0 0],'linewidth',1,'color' , 'r',[-1 1],[-.01 -0.01],'linewidth',1,'color' , 'r')  
  title('Animation')
  axis([-1.5 1.5 -1.5 0.5])
  pause(0.03)
  ani=sprintf("Pendulum%d.png",ct)
  saveas(1,ani)
  ct=ct+1
end

RESULTS:

The picture below is a snapshot of simple pendulum at an instant.

The URL for the animation is available below:

https://www.youtube.com/watch?v=uuczaKWmRdQ&feature=youtu.be