Solving Second Order ODE

Objective:- To write MATLAB program to simulate transient behaviour of simple pendulum

Intoduction

A simple pendulum is an idealized mathematical model of a real pendulum. A weight or bob is attached at the end of a massless cord suspended from a pivot, without friction. There is no frictional energy loss. So, when given an initial displacement it will swing to and fro at a constant amplitude. The assuption considered for the mathematical modeling of simple pendulum are as follows:-

• The rod or cord on which the bob attached is massless, inextensible and always remains taut;
• The bob is a point mass;
• Motion occurs in two dimension i.e. the bob traces a straight line.
• The motion does not lose energy to friction or air resistance.
• The gravitational field is uniform.
• The support is fixed.

.

fig:- simple pendulum

The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "bob" w.r.t time. The differential equation which represents the motion of a simple pendulum is

`frac{d^2\theta}{dt^2}+frac{b}{m}frac{d\theta}{dt}+\frac{g}{\l} \sin\theta=0`

where `g=`accaeration due to gravity (`m/s^2`)

`l=`length of pendulum (`m`)

`b=`damping coefficient

`m=`mass of bob (`Kg`)

solving above diffential equation

let `theta=theta_1`

`frac{d\theta}{dt}=frac{d\theta_1}{dt}=theta_2`

`frac{d^2\theta}{dt^2}=frac{d^2\theta_1}{dt^2}=d/dt(frac{d\theta_1}{dt})=frac{d\theta_2}{dt}`

substituting the above value in differential equation

`frac{d^2\theta}{dt^2}=-theta_2(b/m)-(g/l)sintheta`

writing above equation in matrix form

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

Input parameter

`b=0.05`

`g=9.81 m/s^2`

`m=0.1 Kg`

`l=1 m`

MATLAB code for solving the differential equation

clear all;
close all;
clc;
%damping contants, acceleration due to gravity, length, mass
b=0.05;
g=9.81;
l=1;
m=0.1;
%initial condition (displacement,velocity)
theta_initial=[0;5];
%initial position one end of rod
x1=0;
y1=10;
%time span
t_span=linspace(0,10,500);
%solving ODE
[t,results]=ode45(@(t,theta)ode_function(t,theta,b,g,l,m),t_span,theta_initial);
%plot(t,results(:,1));
%hold on
%plot(t,results(:,2));
%xlabel('time');
%ylabel('plot');
%drawing a circular mass
ct=1
for i=1:length(results);
    x2=l*sin(results(i,1));
    y2=l*cos(results(i,1));
    %r=1;
    %xunit = r * cos(th)+x2;
    %yunit = r * sin(th)+y2;
    plot([x1,x2],[y1,y2],'linewidth',3,'color','r');
    hold on
    [x,y,z]=sphere;
    surf(x+x2,y+y2,z);
    shading interp;
    axis([-2,2,-2,10]);
    pause(0.03);
    M(ct)=getframe(gcf);
    ct=ct+1;
    hold off
end
%playing video file 
movie(M)
%creating uncompressed videofile and specified format
videofile=VideoWriter('pendulam.avi','Uncompressed AVI');
%open videofile for writing
open(videofile)
%write video in a videofile from movie array
writeVideo(videofile,M)
%close the videofile
close(videofile)

function for solving ODE

function[dtheta_dt]=ode_function(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

Result

We can simulate the motion of simple pendulum by solving the above differential equation. The above MATLAB code depicts the motion of simple pendulum. Animation video link is given below,

https://www.youtube.com/watch?v=nUgaistbm3Y