Simulation of Simple Pendulum Kinematics using MATLAB

Logic -

The objective has been achieved by reducing the second order DE to two first order DE.

In this way, we can utilize the MATLAB solver ode45 to then solve both these equations.

The result is then a 2x1 column vector whose 1st column is the Angular displacement while the 2nd column is the angular velocity.

An animation of depicting the pendulum motion has then been created using the 'getframe' and 'VideoWriter' functions available in MATLAB. The time span of the simulation is 20 seconds.

Angular Displacement vs time and Angular velocity vs time curves have also been plotted.

From the graphs, it is clear that the velocity is maximum and minimum displacement which verifies that our solution is correct.

The link to the animation is -

https://youtu.be/Y8kAqelCiuQ

The MATLAB Codes are -

%Motion of pendulum with damping
clear all
close all
clc

%Constants
L = 1;      %Length in meters
g = 9.81;   %Acceleration due to gravity in m/s^2
m = 1;      %Mass of ball in kg
b = 0.05;   %Damping coefficient

t = [0 20];     %Time span. 20 seconds
theta0 = [0 3]; %Initial conditions.Angular displacement and Angular velocity

%Using External Function and parsing L,m,g and b in the function

[t,theta] = ode45(@(t,theta) ode_second_order(t,theta,L,g,m,b),t,theta0);

%Pendulum support
fixed_pt = plot(0,0,'Marker','.','MarkerSize',10,'MarkerEdgeColor','b','MarkerFaceColor','b');  
hold on

%Pendulum arm
arm = plot([0 0],[0 -L],'LineStyle','-','Color','k');
hold on

%Pendulum bob
bob = plot(0,-L,'Marker','o','MarkerSize',10,'MarkerEdgeColor','r','MarkerFaceColor','r');
hold on
ct = 1;
for m = 1:length(t)
    axis([-1 1 -1 0.5]);
    grid on
    x(m) = L*sin(theta(m,1));   %Bob x Location
    y(m) = -L*cos(theta(m,1));  %Bob y Location
    set(arm,'XData',[0 x(m)],'YData',[0 y(m)]);
    set(bob,'XData',x(m),'YData',y(m));
    drawnow;    %Updates output window at each iteration
    hold on
    M(ct) = getframe(gcf);   %Capture current frame using getframe command
    pause(0.01);
    ct = ct+1;  %Loop Counter
    title(sprintf('Simple Pendulum kinematics simulation for t = %0.2f seconds',t(end)));
end

%Crating an uncompressed file called Pendulum.avi and assigning it to
%variable video
video = VideoWriter('Pendulum.avi','Uncompressed AVI');

open(video);
writeVideo(video,M);
close(video);

close all
%Displacement and velocity plots
figure(1);
subplot(2,1,1)
plot(t,theta(:,1),'r-');
xlabel('Time(s)');
ylabel('Angular Displacement(rad)');
subplot(2,1,2)
plot(t,theta(:,2),'b-');
xlabel('Time(s)');
ylabel('Angular Velocity(rad/s)');

Supporting function -

function thetadot = ode_second_order(t,theta,L,g,m,b)

thetadot = zeros(2,1);
% Rate of change of angular displacement is angular velocity
thetadot(1) = theta(2);

thetadot(2) = -((g/L)*sin(theta(1))+((b/m)*theta(2)));

end

The result plot is -