Simulation Of Damped Simple Pendulum

Aim - To write a matlab code for simulating and determining the angular displacement and angular velocity of a damped simple pendulum.

Background - The governing differential equation can be derived by applying Newton's second law i.e. `∑F_i=ma_i` to the pendulum bob. The resulting equation will be a second order non linear ordinary differntial equation.

`(d^2 theta)/(dt^2) + (b/m)*((d theta )/(d t)) + (g/l)*sin(theta) = 0`

where `theta`= angular displacement

m = mass of the bob.

b = damping coefficient.

l = length of string.

g = acceleration due to gravity.

Procedure - The above second order differential equation has to converted into a set of 2 first order differential equations because matlab inbuilt function ode45 uses Ranga-Kutta method to solve the equations, which is valid only for first order ordinary differential equation.

Converting the single second order differential equation into 2 first order equations -

For our convenience

`theta= y(1)`

`\dot theta = y(2)`

now we can write the equation as

`dot y(1) = y(2)`

`dot y(2) = -(g/l)*sin(y(1)) - (b/m)*y(2)`

writing a matlab function for the above governing differential equation.

% Governing equation for animating the simple pendulum

function yout = simple_pendulum_func(t,y,m,l,g,b)
yout = [y(2);
       (-g/l)*sin(y(1))-(b/m)*y(2)];
end

For the given problem we are assuming

mass = 1 kg

length = 1 m

acceleration due to gravity = 9.81 m/s^2

damping coefficient = 0.05

Code for simulating the pendulum:

% Animating the simple pendulum
clear all
close all
clc
%% Initializing the constants
m = 1; % mass 
l = 1; % length 
g = 9.81; % acceleration due to gravity 
b = 0.05; % damping coefficient 

%% Solving the differential equation
t = [0 20];     % Time span
y0 = [0; 3];    % Initial conditions
[t,y] = ode45(@(t,y)simple_pendulum_func(t,y,m,l,g,b),t,y0);

%% Plotting angular_displacement and angular_velocity vs time
figure(1)
plot(t,y(:,1),'b','linewidth',1.2)  % angular_displacement
hold on
plot(t,y(:,2),'r','linewidth',1.2)  % angular_velocity
title('simple pendulum')
xlabel('time')
legend('Angular Displacement \theta','Angular Velocity \omega ')


%% Animating 
figure(2);
for i = 1:length(t)
    clf
    plot([-0.5 0.5],[0 0],'linewidth',10,'color',[0.5 0.5 0.3])
    hold on
    plot(l*sin(y(:,1)),l*-cos(y(:,1)),':y')
    axis([-1 1 -1.2 0.2])
    title(['pendulum at time ',num2str(t(i)),' sec'])
    plot([0 0],[-1.2 0.2],'--k')
    xl = l*sin(y(i,1));  % horizontal location of bob
    yl = -l*cos(y(i,1)); % vertical location of bob
    plot([0 xl],[0 yl],'linewidth',3,'color',[0.75 0.75 0.75]) % pendulum string
    plot(xl,yl,'mo','markersize',20,'markerfacecolor','m')     % bob
    quiver(xl,yl,0.2*y(i,2)*cos(y(i,1)),0.2*y(i,2)*sin(y(i,1)),'linewidth',2)
    shot(i) = getframe(gcf);
end


%% Saving the animation
videofile = VideoWriter('forward_kinematics.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,shot)
close(videofile)

Results - figure(1) shows the variation of angular displacement and angular velocity with time. We can observe that as we progress in time, due to the damping effect the amplitude goes on reducing.

Link for the animation

https://youtu.be/l0rb-n-TlKA