Solving second order ODEs

Aim:

To write a program in Matlab that solves 2nd order ODE of a Simple pendulum and to create a video of it.

Equation of Simple Pendulum:

Given:

Length of the string, L=1 meter

Mass of the bob, M=1 kg

Damping co-efficient=0.05

Acceleration due to gravity, g=9.81 m/sec^2.

Converting 2nd order ODE into 1st order ODE:

let us consider θ=θ1⇒1

diff. on both the sides

dθ/dt=dθ1/dt=θ2⇒2

therefore d^2θ/dt^2=d^2θ1/d^2t=d/dt(dθ1/dt)=dθ2/dt⇒3

sub-equation 1,2 and 3

dθ2/dt+b/m*θ2+g/L*sinθ1=0

Program:

clear all
close all
clc

% Equation of a simple pendulum. 
% inputs
l=1; % length of the string in meters
m=1; % Mass of the bob in kg
b=0.05; % Damping coefficient
g=9.81; % Accleration due to gravity in meter/second^2

%time span
t=linspace(0,40,100);

%displacement
theta=[0 3];

[t,theta]=ode45('second_order',t,theta); %[results]=ode45('fun',intial conditions)

% ploting velocity wrt time
plot(t,theta(:,1),'color','r','linewidth',2)
xlabel('time')
ylabel('Angular displacement/velocity')

hold  on 
% ploting acceleration wrt time
plot(t,theta(:,2),'color','c','linewidth',2)
xlabel('time')
ylabel('Angular displacement/velocity')
title('Angular displacement/velocity vs time')
grid on

% To create video file we are assigning

x0  = 0; % Fixed points of the string at x axis
y0 = 0;  % fixed point of the string at y axis
x1 = sin(theta(:,1)); %  bob position in x axis
y1 = -cos(theta(:,1)); %   bob position in y axis

Videofile = VideoWriter('pendulum.avi','Uncompressed AVI')
open(Videofile);
% loop for the time span 
for k=1:length(t)
    figure(1)
    plot(x0,y0,'Marker','x'); % fixed position
    hold on
    plot(x1(k),y1(k),'Marker','o'); % Bob position
    plot([x0 x1(k)],[y0 y1(k)]); % joining bob and the fixed point 
    axis([-2 2 -2 2]); % axis
    pause(0.02);
    hold off
    M(k) = getframe(gcf);
    writeVideo(Videofile,M(k));
  
  
end
close(Videofile);

function [thetadot]=second_order(t,theta)

thetadot=zeros(2,1);
b=0.05;
l=1;
m=1;
g=9.81;

%Rate of change of angular displacement
thetadot(1)=theta(2);

%Rate of angular accleration
thetadot(2)=(-b/m)*(theta(2))-(g/l)*sind(theta(1));
end

Graph

Result:

By using Ode45, We have solved the Equation of a simple pendulum and created a video of it.

Simple pendulum link