Solving second order ODE for a simple pendulum using MATLAB

Aim: To observe the transient nature of the Simple Pendulum using MATLAB.

Theory: In this problem we consider a Simple Pendulum of mass ‘m’ attached to a string having length ‘l’ and the other end being fixed. ‘b’ is the Damping Coefficient and ‘g’ is acceleration due to gravity.

The mathematical equation for the motion of Simple Pendulum is given by an Ordinary Differential Equation;

`(d^2 Ɵ)/(dt^2 )+b/m*(dƟ)/dt+g/l*sin⁡θ`

In the given equation;

g = 9.81 m/s²; l =1m; m = 1kg; b = 0.05

Procedure: Steps to solve the problem:

Step1: We first consider Ɵ = Ɵ₁

Hence, `(dƟ)/dt=(dƟ1)/dt= Ɵ2`

Step2: We then find Ɵ2

`(d^2 Ɵ)/(dt^2 )=(d^2 Ɵ1)/(dt^2 )=d/dt* ((dƟ1)/dt)=(d /dt)**Ɵ2`

Step3: We then store it in a 2 columns in which the first column will consist of angular displacement values while the second columns of angular velocity values.

`d/dt*[(Ɵ1),(Ɵ2)]=[(Ɵ2),(-b/m*Ɵ2-g/l*sin⁡θ )]`

Now, we will write the code of the following problem in MATLAB in which we will use ode45 to solve the problem.

MATLAB code:

function dtheta_dt = odefunc(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

close all
clear all
clc
%Solving Second order ODE
%Inputs
b = 0.05;
g = 9.81;
l = 1;
m = 1;

%Initial Condition angular velocity of 3 rad/sec at time t = 0 sec
theta_0 = [0;3];

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


%Solve ODE
[t,results] = ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);

%plotting
hold on
plot(t,results(:,1))
plot(t,results(:,2))
xlabel('time')
ylabel('plot')
title('Variation of Angular displacement & Angular velocity w.r.t time');
theta2 = results(:,1);
ct =1;
for i = 1:length(theta2) 
    theta_x = theta2(i); %Theta_x keeps changing w.r.t time
     x0 = 0;
     y0 = 0;
     x1 = l*sin(theta_x);
     y1 = -l*cos(theta_x);
     
     %Plotting the pendulum
     figure(2)
     plot([-0.5 0.5],[0 0],'color','k','linewidth',3)
     hold on
     plot([x0 x1],[y0 y1],'color','b','linewidth',3)
     hold on
     plot(x1,y1,'o','markers',10,'markerfacecolor','r')
     hold off
     grid on
     axis([-2 2 -2 2]);
     pause(0.05)
     M(ct) = getframe(gcf);
     ct = ct +1;
end
movie(M)
videofile = VideoWriter('simple_Pendulum.AVI','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

FIGURE:

You Tube Link:

https://youtu.be/hK-DQw2E1dw

Errors or warnings in the program:

For ODE function we will get an error of:

Not enough input arguments.

Error in ode_func (line 4)
theta1 = theta(1)