Simulation of the Transient Behaviour of a Simple Pendulum by solving ODEs

Simulation of the Transient Behaviour of a Simple Pendulum by solving ODEs

Introduction:

In mathematics, an ordinary differential equation (ODE) is a differential equation containing one or more functions of one independant variable and the derivatives of those functions. The term ordinary is used in contrast with the term partial differential equation which may be with respect to more than one independent variable.

In Engineering, ODE is used to describe the transient behavior of a system. A simple example is a pendulum

The way the pendulum moves depends on the Newtons second law. It states that the acceleration of the object is directly proportional to the net force acting on the object and inversly proportional to the mass of the object

When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the \"ball\" w.r.t time.

In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

Similarly, there are hundreds of ODEs that describe key phenomena and this ODE represents the equation of motion of a simple pendulum with damping.

Inputs given:

L=1 m

m=1 kg

b=0.05

g=9.81 m/s2

Time period = 0-20 sec

Angular displacement=0

Angular velocity=3 rad/sec at time t=0.

Solving the above second order ODE:

• d^2theta/dt^2 + b/m * dtheta/dt + g/l sind(theta) = 0
• theta = theta1
• dtheta/dt = dtheta1/dt = theta2
• d^2theta/dt^2 = d^2theta1/dt^2 = d/dt * dtheta/dt = d/dt * (theta2)

In equation sub the values,

by that we are converting the second order ODE into two first order ODE

• d^2theta1/dt^2 + b/m * dtheta1/dt + g/l * sind(theta1) = 0
• dtheta2/dt + b/m * theta2 + g/l * sind(theta1) = 0
• d/dt [theta1; theta2] = [theta2; -b/m * theta2 - g/l * sind(theta1)]

The above matrix is the solution vector from this we can solve the program.

Program - 1:

%Solving the second order ODE
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/1)*sin(theta1);
dtheta_dt = [dtheta1_dt; dtheta2_dt];
end

Code Explanation :

• In the above program we are using the function command for giving the inputs and solving the second order ODEs equation.
• And we are solving the solution vector dtheta1_dt and dtheta2_dt and determining it in the main program.

Program - 2:

clear all
close all
clc

%Input values
b = 0.05;
g = 9.81;
l = 1;
m = 1;

%Initial conditions
theta0 = [0;3];

time = linspace(0,20,500);

%Calculation of ODE
[a,b] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),time,theta0);

%Plotting
figure(1)
plot(a,[b(:,1) b(:,2)])
xlabel(\'Time in sec\')
ylabel(\'Angular Displacement & Angular velocity(rad/sec)\')
legend(\'Angular Displacement\' , \'Angular Velocity(rad/sec)\')
axis([0 22 -3 3])
hold on

%Simulation of Pendulum
theta = b(:,1);
c = 1;

for i = 1:length(theta)
    theta1 = theta(i);
    
    x0 = 0;
    y0 = 0;
    
    x1 = l*sin(theta1);
    y1 = -l*cos(theta1);
    
    figure(2)
    
    plot([-1.4 1.4],[0 0],\'linewidth\',5,\'color\',\'g\')
    hold on
    plot([x0 x1],[y0 y1],\'color\',\'b\',\'linewidth\',3);
    hold on
    
    plot(x1,y1,\'O\',\'markers\',30,\'color\',\'r\',\'markerfacecolor\',\'r\');
    hold off
    
    axis([-1.5 1.5 -1.5 0.5])
    pause(0.003)
    M(c) = getframe(gcf);
    c = c+1;
end
movie(M)
videofile = VideoWriter(\'pendulum.avi\',\'uncompressed AVI\');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Code Explanation:

• Initially we are describing the Input values, Initial conditions of angular displacement and angular velocity in a row matrix and then we are giving the time period using the linspace command
• And then we are solving the ODE equation using ode45 command and by determining the input values from function program and also adding the time span and initial conditions respectively.
• Finally we are plotting the time vs angular displacement and time vs angular velocity graph and we are using (:) operator for getting the result of 1st & 2nd column respectively.
• And then we are labelling the graph and adding legend command to define the angular displacement and angular velocity and also adding axis command to determine the x&y axis size.
• Next step is to simulate the pendulum therefore we are getting the angular displacement value and assuming it as theta and we are using the for loop to repeat the set of action again and again.
• Then in the for loop we are assigning the theta1 as length of theta and also assigning the x0,y0,x1,y1 values.
• And then we are plotting the bar were the string is tied and then we are plotting the string and finally we are plotting the ball.
• The ball is plotted by using the marker command and inserting the circle market style to it then giving the markerfacecolor command to fill the circle.
• Then we are assigning the counter value as 1 initially and each time the for loop is being executed the value gets incremented.
• Then to create a movie, we have to save the plotting window into a array variable. And then movie is generated by plotting window with all the frames.
• And the video file is created with the ‘avi’ format and we are opening it and writing the frames we captured from the plotting window to it.

Errors:

The error occured due to usage of (:) instead of (;) as it should be row matrix and its corrected and the program is executed.

Output1:

The plot represents change in displacement & change in velocity with respect to the time.

Output2:

Simple Pendulum motion Animation

Result:

Hence the given second order ordinary differential equation is solved and the graph has been plotted for change in displacement & change in velocity with respect to the time and then the Simple Pendulum motion has been simulated.