To solve the ODE of a second order to study a motion of a simple pendulum

1.AIM:

To solve the ODE of the simple pendulum and to plot the curve of various parameter by using MATLAB coding

1. Time vs angular displacement and Time vs angular velocity
2. angular displacement vs angular velocity

To simulate the motion of a pendulum for time t= 0 to 20 second

To make the animation of the motion of pendulum

2.GOVERNING EQUATIONS:

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

In the above equation,

g = gravity in `m/s^2`

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

Given data,

g = 9.81 `m/s^2`

L = 1 m,

m = 1 kg,

b=0.05,

angular displacement = 0 m and angular velocity = 3 `(rad)/sec` at time t= 0 second.

3.BODY OF THE CONTENT

A simple pendulum consists of a mass attached at the end of the string and the other end of the string is fixed as shown in fig

The very first step to solve the above ODE is to convert the second order differential equation into first order differential equation

`theta`=`theta1`
`theta1`= `(d(theta))/dt`

`(d(theta1))/dt` = `theta2`

`(d(theta2))/dt` = `-(b/m)*theta2 - (g/l)*sin(theta1)`

The program is completed into two parts. First part which is function to solve the ordinary differential equation and second part is for plotting and making the animation of the simple pendulum.

MATLAB CODE

First part of the program is function which is discussed below

% ode_func is the function which want t, theta, b,g, l , m as an input
%where these are input parameter
function [dtheta_dt] = ode_func(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

Second part of the program is discussed below

clear all
close all
clc
% defining values of damping coefficient(b), acceleration due to
% gravity(g),length (l), mass(0.1)
b= 0.05;
g= 9.81;
l= 1;
m= 1;
%initial condition theta_0 
theta_0 = [0,3];       %[initial displacement; initial velocity]

% time span
t_span = linspace(0,20,500);    %spliting t=0 to 20 into 500 parts by using linspace command

%solve ode
%ode45 is a versatile ODE solver 
%@(t,theta) variable used in the ODE
% ode_func is the function, which is call by this ode45
%  t_span is the time span
% theta_0 is the initial condition
[t, result]= ode45(@(t, theta) ode_func(t,theta,b,g,l,m), t_span,  theta_0);
%defining figure 1
%plot for angular displacement and time for time 0 to 20 seconds
figure(1)
%plot formation under figure 1
plot(t, result(:,1)) %all the values of displacement for time 0 to 20s

%hold on command is used to hold the plot
hold on
%plot formation under figure 1
%plot for angular velocity and time for time 0 to 20 seconds
plot(t, result(:,2)) %all the values of angular velocity for time 0 to 20s 
%defining the x axis by xlabel as time and y axis by ylabel as plot
xlabel('time')
ylabel('plot')
%legend used for detailing plot
legend('angular displacement','angular velocity')

%defining figure 2
figure(2)
%plot for angular displacement and angular velocity for time 0 to 20
%seconds
plot(result(:,1),result(:,2),'linewidth',1,'color','r')
%hold off command is used to off the plot just like end
hold off

%starting of animation of pendulum 
%defining counter fot ct=1
ct=1;
%begining of for loop from i to length of t_span
for i=1:length(t_span)
    %%defining figure 3
    figure(3)
    %initial condition of the string 
    x0=1;
    y0=1;
    
    %figure(3);
    %defining the coordinates of end point as x1 and y1
    x1= 1- l*sin(result(i));
    y1= 1- l*cos(result(i));
    %plot of line on which one end of string is fixed
    plot([2.5 -0.5],[1 1],'linewidth', 3,'color','b')
    hold on
    %figure(3);
  %plot for the string which is fixed at x0, y0 and end at x1, y1 
    plot([x0 x1],[y0 y1])
   % figure(3);
   %defining the staring and end point of the axis i.e x and y axis
    axis([-1 3 -1 3])
   % making a ball like structure at the end of string
   %'o' is used to define the shape at the end of string
   %'MarkerSize',20 is used to give size of the shape used
   %'MarkerFaceColor','r' is used to give color on the face of the shape
   %used
   %'MarkerEdgeColor','b' is used to give color at the edge of the shape
   plot(x1,y1,'o','MarkerSize',20,'MarkerFaceColor','r','MarkerEdgeColor','b')
   hold off
   %pause time 0.3 sec
   pause(0.3);
    %getting a frame for different position for the pendulum and stored in
    %M(ct)
    
    M(ct)= getframe(gcf);
    %increament of counter
    ct=ct+1;
    
   
end


%command for generating animation file
movie(M)
%defining the name and format of the animation
videofile= VideoWriter('practice.avi','uncompressed AVI');
% open the video file to store the data
open(videofile)
%writing the data in the video
writeVideo(videofile,M)
%closing of video file
close(videofile)

Error Obtain While Programming

1. Because of calling function ode_fun instead of ode_func.

2. Because of not defining the shape of bob which later I use ‘o’ to give the shape of bob.

5.CONCLUSION

The plot obtain after solving the ODE is

1. Time vs angular displacement and Time vs angular velocity

2. angular displacement vs angular velocity

3. To make the animation of the motion of pendulum

Youtube link