Week 3 - Solving second order ODEs

AIM

To solve the ODE which represents the motion of a simple pendulum

OBJECTIVE

To write a program that solves the following ODE and to create an animation of the oscillation of the pendulum

THEORY

The motion of a simple pendulum is example of simple harmonic motion,consisting a small bob and a massless string.movement of pendulum depends upon the newton's second law of motion.by this a second order derivative equation is formed with repect to time

EQUATION OF MOTION

In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

INPUTS

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

time span=0-20 sec

angular displacement=0

angular velocity=3 rad/sec @ time t=0

CONVERTING 2nd ORDER INTO 1st ORDER

thetha=thetha1

d_thetha1/d_t=thetha2

d_thetha^2/d_t^2=d_thetha2/d_thetha

MAIN PROGRAM/CODE WITH ANIMATION

clear all
close all
clc


%inputs
b=0.05
g=9.81
l=1
m=0.01

%initial conditions
thetha_0=[0;5];

%time points
t_span=linspace(0,10,200)

%solve ode

[t,results]=ode45(@(t,thetha) ode_func(t,thetha,b,g,l,m),t_span,thetha_0)

plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('plot')
xlabel('time')
displacement=results(:,1)
velocity=results(:,2)
hold off 

%Fixed coordinates

x0=1;
y0=1;
ct=1;

for i= 1:length(displacement)
%coordinates
x1=1-1*sin(displacement(i));
y1=1-1*cos(displacement(i));
figure(2)
%plot fixed link
plot([0 2],[1 1],'linewidth',15,'color','b');
hold on
plot([1 x1],[1 y1],[x1 x1], [y1 y1],'-O','markersize', 25,'markerfacecolor','r','linewidth',4);
plot([x0,y0]'-0','markersize',3,'markerfacecolor','b')
axis ([0 2 -0.5 2])

%movie


M(ct)=getframe(gcf);
ct=ct+1;
hold off
end
movie(M)
videofile=VideoWriter('simple_pendulum.avi','Uncompressed AVI')
open(videofile)
writeVideo(videofile, M)
close(videofile)

FUNCTION PROGRAM(ODE_FUNC)

function [dthetha_dt] = ode_func(t,thetha,b,g,l,m)
thetha1=thetha(1)
thetha2=thetha(2)
dthetha1_dt=thetha2
dthetha2_dt=-(b/m)*thetha2-(g/l)*sin(thetha1);
dthetha_dt=[dthetha1_dt;dthetha2_dt];
end