Week 3 - Solving second order ODEs

SIMPLE PENDULUM

One of the classical example for a second order ordinary differential equation is the motion of the bob of a simple pendulum. From Newton's 2nd law, we can get the motion of the bob of a simple pendulum as a Second order ODE as below:

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

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.

The following is a Matlab program to calculate the and simulate the motion of a bob in a simple pendulum between 0 to 20 seconds with the following inputs and initial conditions

inputs : L=1m; m=1kg ; b=0.05 ; g=9.8 m/s^2

initial condition : at t=0, `theta=0 & omega=(dθ)/(dt)`=3 rad/sec

time span = 0-20 seconds

Main script:

clear all
close all
clc

%inputs
b=0.05;                                          %damping coefficient
m=1;                                             %mass of the ball in kg
l=1;                                             %length of the pendulum in m
g=9.8;                                           %acceleration due to gravity in m.s^-2

theta_0=[0,3];                                   %intitial conditions
t_span=linspace(0,20,200);                       %time interval

%calling ode_pendulum_func
[t,result]=ode45(@(t,theta) ode_pendulum_func(t,theta,b,g,m,l),t_span,theta_0);

%plotting
figure(1)
plot(t,result(:,1),'color','red')
hold on
plot(t,result(:,2),'color','blue')
xlabel('time(s)');
legend('angular displacement','angular velocity');
title('plot of the motion of a simple pendulum');

%Animation of the motion of the pendulum
a=result(:,1);
for i=1:length(a)
    angle=a(i);
    %fixed points
    x0=0; y0=0;
    %lose end ponts
    x1=l*sin(angle);
    y1=-l*cos(angle);
    figure(2)
    plot([x0,x1],[y0,y1],'r','linewidth',2);
    hold on
    plot([-2,2],[0,0],'color','green','linewidth',3);%ceiling line
    plot(x1,y1,'.','MarkerSize',50);
    axis auto
    M(i)=getframe(gcf);
end

movie(M)
video_file=VideoWriter('motion of a simple pendulum.avi','uncompressed AVI');
open(video_file);
writeVideo(video_file,M);
close(video_file);

Function Code:

function[Dtheta]=ode_pendulum_func(t,theta,b,g,m,l)
Dtheta_1=theta(2);
Dtheta_2=-(b/m)*theta(2)-(g/l)*sin(theta(1));
Dtheta=[Dtheta_1;Dtheta_2];

Results: