MATLAB program to simulate the transient behaviour of a simple pendulum

AIM : To write a MATLAB program to simulate the transient behavior of a simple pendulum and to create an animation of its motion.

BACKGROUND : 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 Newton’s 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 we have the ability to solve these equations then we will be able to simulate any of these systems.

Our objective is to write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping.

`frac{d^2theta}{dt^2}+frac{b}{m}frac{d theta}{dt}+frac{g}{L}sin theta=0`

g = gravity in m/s²,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

INPUT VARIABLES :

L=1 metre, m=1 kg, b=0.05, g=9.81 m/s².The motion is simulated between 0-20 seconds.

INITIAL CONDITIONS :

At time t=0, angular displacement=0 ; angular velocity= 3 rad/sec.

CALCULATIONS : In order to solve the above ODE, the following calculation is done to get the solution through MATLAB. The basic idea involved here is to break down a single `n^{th}`order ODE into 'n' lower order ODEs and representing them in matrix form to obtain the solution vector.

let, `theta=theta1`

`frac{d theta}{dt}=frac{d theta1}{dt}=theta2` let

hence, `frac{d^2 theta}{dt^2}=frac{d^2 theta1}{dt^2}=frac{d}{dt}(frac{d theta1}{dt})=frac{d}{dt}theta2`

Using these expressions in the governing ODE,

`frac{d^2 theta1}{dt^2}+frac{b}{m}frac{d theta1}{dt}+frac{g}{L}sin theta1=0`

`frac{d theta2}{dt}+frac{b}{m}theta2+frac{g}{L}sin theta1=0`

hence, `frac{d theta2}{dt}=-frac{b}{m}theta2-frac{g}{L}sintheta1`

Representing in matrix form,

`frac{d}{dt}((theta1),(theta2))= [[theta2,0],[-frac{b}{m}theta2,-frac{g}{L}sin theta1]]`

SOLUTION ALGORITHM

1. The input variables will be entered as mentioned in the problem.
2. The time period will be split into 500 parts between limits of 0-20 seconds as mentioned in the problem.
3. The displacement and velocity as a function of time will be calculated in a separate function file and it will be called in the main program file.
4. First, the plotting will be done to represent displacement and velocity as a function of time with displacement and velocity taken along y-axis and time taken along x-axis.
5. Then an animation plot will be created to simulate the movement of pendulum.
6. The horizontal bar, vertical bar of the pendulum and the pendulum bob will be created.
7. The position of the pendulum at particular time will be captured and finally all the frames will be combined.
8. The output will be obtained in a video file to represent the motion in the form of an animated video.
9. The entire script will be written in MATLAB.

MATLAB CODE(for main file) :

%program to solve a second order ODE
%ODE represents motion of a simple pendulum with damping

%input variables

%gravity in m/s^2
g= 9.81;

%length of the pendulum in metres
L=1;
 
%damping coefficient
b= 0.05;

%mass of the ball in kg
m= 1;

%initial condtions for angular displacement and angular velocity
theta_init= [0;3];

%time period over which the motion is to be simulated is between 0 to 20
%seconds
t_period= linspace(0,20,500);

%the results are calculated in a separate function file 
%calling the function in the main file
[t,result]= ode45(@(t,theta)ode_solve(t,theta,b,g,L,m),t_period,theta_init);

%plotting the results
figure(1)
hold on
plot(t,result(:,1))
plot(t,result(:,2))
grid on
xlabel('time')
ylabel('displacement&velocity')
title('Displacement and Veolcity v/s Time')
legend('displacement','velocity')


%creating an animation out of the plot
for i= 1:length(t_period)
    x0=0;
    y0=0;
    THETA1= result(i,1)
     x1= L * sin(THETA1)
     y1= -L * cos(THETA1)
      figure(2)
%creating the horizontal bar
     plot([-1 1],[0 0],'Linewidth',2,'color',[0.3452 0.4567 0.5632])  
     hold on
%creating the vertical bar of the pendulum     
     plot([x0 x1],[y0 y1],'linewidth',3,'color','b')
%creating the pendulum bob     
     plot(x1,y1,'marker','o','MarkerSize',10,'MarkerFaceColor','r')
     axis([-1.5 3 -1.5 0.4])
     grid on
     hold off
%capturing the frame     
     M(i)= getframe(gcf);
    end
%creation of animated video clip
movie(M)
videofile= VideoWriter('pendulum_simu.avi', 'Uncompressed AVI');
%creating an empty video file of AVI format
open(videofile)
%writing the content in the video file
writeVideo(videofile,M)
%closing the file
close(videofile)

MATLAB CODE(for function file) :

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

OUTPUT : The following output is obtained after running the code :

* The animation of the output can be viewed from the following youtube link:

https://youtu.be/9IoFPPS5-5k

ERRORS :

1. The error occurred because of typing mistake and hence was rectified.

2. The error occurred because I was not using hold off command and hence it was retaining the current plot and axes. With the use of hold off command, the axes were reset to default resulting in a continuous motion.

CONCLUSION :

• The same model can be used to solve other order ODE and simulate any other system.
• Solving the ODEs is extremely helpful especially while doing vehicle body dynamics.

NOTE :

• ode45 is an inbuilt MATLAB solver that is used to solve the ODE here.