Simulating the transient behaviour of a simple pendulum and creating an animation of its motion

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 \"bob\" with respect to time.

The equation of motion of a simple pendulum with damping is as follows:

where,

g = acceleration due to gravity in m/s^2

l = length of the pendulum in m

m = mass of the ball in kg

b = damping coefficient

Matlab Code to plot the graph of Angular displacement and Velocity against time and to create an animation of the pendulum motion:

clear all
close all
clc

% Inputs

b = 0.05;                         % damping coefficient
g = 9.81;                         % gravity constant
l = 1;                            % length of the pendulum from the fixed point
m = 1;                            % mass of the ball
theta_0 = [0;3];                  % initial conditions (Angular Displacement , Velocity)             
t_span = linspace(0,20,500);      % time plots from 0 to 20 seconds
[t, results] = ode45(@(t,theta) ode_function(t,theta,b,g,l,m), t_span, theta_0); % Calling ODE function

% Plotting the graph

figure(1)                         
plot(t,results(:,1),\'color\',\'r\')
hold on 
plot(t,results(:,2),\'color\',\'b\')  
xlabel(\'Time\')
ylabel(\'Angular Displacement and Velocity\')
legend(\'Angular Displacement\', \'Angular Velocity\')
axis([0 10 -5 5])

% Plotting the position of pendulum at any given time

x0 = 0;                             % fixed points 
y0 = 0;
frames_ct = 1;                      % variable defined to count number of frames


for i = 1:length(t_span)            % for loop used to plot all values of x1 and y1
    x1 =  l*sin(results(i,1));      % coordinates of bob on x-axis  
    y1 = -(l*cos(results(i,1)));    % coordinates of bob on y-axis
    
   % Plotting and creating an Animation

   figure(2)    
   plot([-1,1],[0,0],\'linewidth\',2,\'color\',\'g\')
   hold on 
   plot([x0 x1],[y0 y1],\'linewidth\',3,\'color\',\'b\')
   plot(x1,y1,\'o\',\'markersize\',30,\'markerfacecolor\',\'r\',\'markeredgecolor\',\'r\')
   hold off
   axis([-1 1 -1.5 0.5])
   pause(0.001)
   M(frames_ct) = getframe(gcf);
   frames_ct = frames_ct + 1;
end
    
% Creating a video file to save and run Animation

movie(M)
v = VideoWriter(\'Pendulum Simulator.avi\',\'Uncompressed AVI\');
open(v)
writeVideo(v,M)
close(v)

Converting the 2nd order ODE to 1st order ODEs:

Let,

`theta = theta1`

`(d theta) / dt = (d theta_1)/dt = theta_2`

`(d^2 theta)/dt^2 = (d^2 theta_1)/dt^2 = (d theta_2)/dt`

We get,

`(d theta_1)/dt = theta_2` (Angular Displacement values)

`(d theta_2)/dt = -(b/m)theta_2 - (g/l)sintheta_1` (Angular Velocity values)

Matlab Code to define function to perform ODE calculation which is called in the main program:

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

Output:

1. Graph of Angular Displacement and Angular Velocity against Time:

2. Animation of motion of simple pendulum with damping: