Simulating a 2nd order ODE system of a simple pendulum

Aim: To write program and simulate the motion of simple pendulum by solving 2^nd order differential equation in MATLAB.

Simple Pendulum:

It can be defined as point mass suspended on mass less string or rod.

Now, when this point mass is applied with given acceleration, it is subjected to periodic motion of ‘to & fro’ from equilibrium point at the center. So, time period for each oscillation is given by,

`T=2π*√l/g 1)`

`where T=time taken to complete 1 oscillation ,l=length of pendulum (m),g=gravitation acceleration`

The acceleration of single pendulum is given by,

`a= -g*sin⁡θ 2)`

The way pendulum moves depends upon the Newton ‘s 2^nd law. When this law is written we get second order differential equation, which governs the position of ball w.r.t time.

This governing equation is,

`(d^2 θ)/(dt^2 )+b/m dθ/dt+g/l sinθ=0 3)`

`Here b = damping coefficient, θ=angular displacement ,m = mass in kg`

To solve this eqn no. 3), we will break into two 1^st order ODES. We will make following assumptions.

`θ_1= θ 4)`

`(dθ_1)/dt= θ_2 5)`

When we substitute equation 4) & 5) in equation 3), we will have,

`(dθ_2)/dt+b/m θ_2+g/l sinθ_1=0 6)`

`(dθ_2)/dt = -b/m θ_2-g/l sinθ_1 7)`

`d/dt[〖(θ〗_1;θ_2)]=[θ_2;( -b/m θ_2-g/l sinθ_1)] 8)`

Here the purpose is to animate the motion of simple pendulum.

At 1^st We write the user defined function that will solve 2^nd order ODE by breaking it into two 1^st order differential equations. The function is as follow.

%User define function to solve 2nd order ode for pendulum%
function [dtheta_dt] = ode_fun3(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

In this code main function is defined with the initial condition and then it will be plotted for displacement and angular velocity varying wrt time.

The main program is as follow.

%MATLAB Program to simulate pendulum motion%
clear all
close all
clc 
b = 0.05 ; %Damping coefficient 
g = 9.81 ; %Gravitional constant 
m = 1.00 ; %Mass in kg 
l = 1.00 ; %length in meter 

%Initial displacement & velocity%
theta_0 = [0;5];

%Time span% 
t_span = linspace (0,20,500);

%Solving 2nd order ode using ode_pendulum external function% 
[t, results] = ode45(@(t,theta) ode_fun3(t,theta,b,g,l,m),t_span,theta_0);

%Procedure to plot displacement & velocity wrt  varying time %
figure(1)
plot(t,results(:,1))
hold on 
plot(t,results(:,2))
axis([0 20 -5 5]);
xlabel('Plot')
ylabel('Time')

%To plot simple pendulum%
%Values of theta w.r.t. time stored in column %
THETA  = results(:,1);

% Counter variable is intilised to 1,it will be incremented to capture
% frames further%
c = 1;
%fixed refference pointis considered to origin %
x0 = 0;
y0 = 0;
%For loop will be used to calculate the other end of the pendulum for all
%angles%
for i = 1:length(THETA)
    THETA1 = THETA(i);
    x1 = l*sin(THETA1);
    y1 = -l*cos(THETA1); %-VE Sign because of  swinging in 3rd & 4th quadrtant %
    
    figure(2);
    
    plot([x0,x1],[y0,y1],'color','r','linewidth',2)
    hold on 
    plot(x1,y1,'o','markerfacecolor','b')
    hold off 
    grid on 
    axis([-1.5 1.5 -1.5 0.5])
    pause(0.003)
    M(c) = getframe(gcf);
    c = c+1 ;
end 
%To get annimation of pendulum %
movie(M)
videofile = Videowriter('Simple_pendulum.avi','uncompressed avi');
open(videofile)
writeVideo(videofile,M)
close(videofile)

The plot showing variation of angular displacement & angular velocity w.r.t time.

The animation of simple pendulum is shown in below link.

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

From the video it can be seen that the motion is logarithmically decreasing.