Simulation of a pendulum by solving a second order ordinary differential equation using MATLAB

OBJECTIVES - To write a MATLAB program to solve a second-order ODE and simulate the pendulum motion

THEORY - The Differential equation is a mathematical equation that relates some functions with its derivatives. In application the function represents a physical quantity and the derivatives represent its rate of change and the differential equation represents the relation between the two. The highest order derivative sets the order of the equation.

ODE refers to an ordinary differential equation in which the highest power of equation is one. ie, The order of the equation is one.

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

The way the pendulum moves depends on the Newtons second law. 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.

Similarly, there are hundreds of ODEs that describe key phenomena and if you can solve these equations then you will be able to simulate any of these systems.

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.

INPUTS

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

angular velocity= 20 rad/sec

angular displacement = 0

time = 20seconds

FUNCTION PROGRAM

function [dtheta_dt] = func_ode(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

 

Initially, the 2nd order ODE representing the pendulum motion is brought down to ODE. The function program is used to create a matrix dtheta_dt which consists of of elements dtheta1_dt = theta2 and dtheta2_dt. Later this function program will be called in the main program to obtain the displacement and angular velocity values.

MAIN PROGRAM

clear all
close all
clc

%inputs
b = 0.05;
g = 9.81;
l = 1;
m = 1;
%intitial condition
theta_0 = [0;3];
%time points
t_span = linspace(0,20,500);
%solve ODE
[t,results] = ode45(@(t,theta) func_ode(t,theta,b,g,l,m),t_span,theta_0);
plot(t,results(:,1))
hold on
plot(t,results(:,2))
xlabel('time')
ylabel('Angular displacement vs Angular velocity')
hold off
grid on

%plotting pendulum
figure

cp = 1;
for i = 1:length(t_span)
    x_init = 0;
    y_init = 0;
    x_fin = l*sin(results(i,1));
    y_fin = -l*cos(results(i,1));
plot([x_init x_fin],[y_init y_fin],[-0.5 -0.5],[0 0],'linewidth',3)
hold on
plot(x_fin, y_fin, 'Marker','O', 'MarkerSize',30,'MarkerFaceColor','b')
axis([-1.5 1.5 -1.5 0.5])
grid on 
pause(0.01)
hold off
M(cp) = getframe(gcf);
cp = cp+1;

end

%video file
movie(M);
Videofile = VideoWriter('pendulum_animation_ode.avi','Uncompressed AVI');
open(Videofile);
writeVideo(Videofile,M);
close(Videofile);

STEPS

1. Input all the parameters including damping coefficient =0.05 , g=9.81 m/s2, length of string = 1m, mass of pendulum= 1kg.

2. Set the speed of the pendulum = 3rad/sec and the time span to be 20 seconds.

3. Call the function program and using the ode45 solver solve the function which includes the differential equation governing the pendulum motion.

4. Plot the results with time in the x-axis and angular velocity and displacement in the y-axis.

5. Using a FOR loop plot the marker for the given pendulum motion.

6. Create a videofile in AVI format using Videofile command.

OUTPUT

Angular velocity and displacement vs Time

Animation link

https://youtu.be/JLYCrj1lR1s

ERRORS

I made a mistake of using shift+enter command to change the 'm' to 'M' in movie command but it also changed the 'm' in mass to 'M' which caused this error.

This error arose because I forgot to use the hold off command.

CONCLUSION

The pendulum motion using MATLAB solver ode45 is achieved.