Pendulum, Solving second order ODEs

Aim:

The primary aim is to solve the given ordinary differential equation and calculate the displacement and velocity vector, while the secondary aim is to simulate the motion of simple pendulum.

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

Governing Equation:

According to newton's second law of motion, the equation of motion of oscillating body can be described as: `I(d^2theta)/ dt^2+b(d theta)/dt+theta=0`

Where: ‘I’ is moment of inertia, in the given case, `I = ml^2`

‘b’ is damping coefficient

`theta` is angular displacement

‘t’ is time time interval

Objective:

The objective was to solve and compute the displacement and velocity of the above differential equation in a time interval of 0 to 20 second, with initial displacement 0 and velocity 3rad/sec at t=0. The secondary objective was to generate a simulation of a simple pendulum with the above boundary conditions.

g = gravity i.e 9.81 m/s2,

l = length of the pendulum is 1 m,

m = mass of the ball is 1kg,

b=damping coefficient.

The code for function:

function [dtheta_dt]= ode_func(t,theta,b,g,l,m) % function dtheta/dt is defined with independent variable time't' and dependent variable 'theta' along with constants 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

The above code is described below:

The higher order differential equation:

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

Let: `theta_1= theta and theta_2= (d theta)/dt`

`(d theta_1)/dt=(d theta)/dt& (d theta_2)/dt= (d^2theta)/dt^2`

therefore: `theta_2=(d theta_1)/dt`

Converting to first order differential equation:

Now we can say that: `(d theta_2)/dt`= `-(b/m)theta_2 -(g/l)sin(theta_1)`

let:

`theta= [[theta_1],[theta_2]]`

and

`(d theta)/dt=[[(d theta_1)/dt],[(d theta_2)/dt]]`

First order vector valued ordinary differential equation:,

Hence,

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

Rest of the code:

clear all
close all
clc

%inputs
l=1; %length of pendulam
b=0.05; %damping coefficient
m=1; %mass of pendulam bob
g=9.81; %acceleration due to gravity

ct=1; %counter
%initial conditions
theta_0= [0;3]; %initial angual displacement is 0 and angular velocity is 3rad/sec

%time span
tlim= [0,20]; %time period of is of 20sec

[tsol,thetasol] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m), tlim, theta_0); %passing the command to solve ordinary differntial equation by defining the variables and function also giving the limiting conditions. 
%% Animation
O = [0 0]; % origin is defined
axis(gca, 'equal'); 
axis([-0.9 0.7 -1.5 .05]); %axis is definded as [x1 x2 y1 y2]
grid on 
%loop for animation
for  i= 1:length(tsol)%loop, in which 'i' extands from 1 to length of tsol
    P = l*[sin(thetasol(i,1)) -cos(thetasol(i,1))]; % the cordinates of bob is given by alloting the the different positions while in loop
    O_circle= viscircles([0 0],0.01); %hinge point for pendulum is created by alloting coordinates and radius
    pend = line([O(1) P(1)], [O(2) P(2)]); %link is created between bob and hinge
    ball= viscircles(P,0.05); %bob is created giving coordinate and radius as input
   
    pause(.0005); %a pause for frames
    M(ct)= getframe(gcf); %craeting a array to get the current frame of plot
    ct= ct+1; %incrementaion in counter after each for looop
   
   
    if i<length(tsol)
        delete(pend);
        delete(ball);
        delete(O_circle);
    end
end
movie(M)% creating the movie
videofile= VideoWriter('pendulam.avi','Uncompressed AVI');%writing the movie
open(videofile);%opening the video file
writeVideo(videofile,M);
close(videofile);% closing the video file.
hold on
plot(tsol,thetasol,'linewidth',3);

The link for video animation is provided below:https://www.youtube.com/watch?v=_1jxTaq9B4E

Conclusion:

1. Displacement and velocity vector is calculated at different instant of time.
2. With given damping coefficient the allotted time interval was not enough to bring the system to equilibrium condition.
3. Either by increasing the damping coefficient or the time interval we can reach upto the equilibrium state.
4. We can easily judge whether the system is underdamped, overdamped or critically damped by playing with different damping coefficients.