Week 3 - Solving second order ODEs

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 'I' and the other end being fixed. 'b' is the Damping Coefficient and 'g'

is the acceleration due to gravity.

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

g = 9.81 m/s2; I =1m; m 1kg; b = 0.05

Procedure: Steps to solve the problem:

Step1: We first consider = 0

Step 2: We then find θ2

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

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:

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) theta 1 = theta(1)

Animation Link: https://www.youtube.com/watch?v=dY8IMJLegDc&feature=youtu.be