Simulating the Motion of a Simple Pendulum with Damping in MATLAB

PROGRAM IN MATLAB/OCTAVE TO SIMULATE THE PENDULUM MOTION

OBJECTIVE:-

1. To write a program in MATLAB/OCTAVE to solve the equation of motion of a simple pendulum (with damping).
2. To simulate the pendulum motion b/w 0-20sec, for angular displacement=0, angular velocity=`3 rad⁄sec`; at time t=0.

PROBLEM FORMULATION AND METHODOLOGY:-

The movement of a simple pendulum is governed by Newton’s second law of motion. When the second law of motion is written for the pendulum, we get a second-order ordinary differential equation (ODE); that describes the position of the ‘bob’ w.r.t time.

The ODE representing the equation of motion of a simple pendulum with damping is given by:-

`(d^2 θ)/(dt^2 )+(b/m)(dθ)/dt+(g/L)sinθ=0`

Where, acceleration due to gravity, `g=9.81m⁄s^2`

Length of the pendulum, `L=1m`

Mass of bob, `m=1kg`

Damping coefficient, `b=0.05`

`θ` = Angular displacement of the pendulum

`t` = time interval

To solve a second-order ODE in MATLAB/OCTAVE, we need to convert second-order ODE into two first-order ODE’s and then we can solve them simultaneously. To demonstrate this process, consider the dependent variable as theta (`θ`), and the independent variable is time (`t`).

`(d^2 θ)/(dt^2 )+(b/m)(dθ)/(dt)+(g/L)sinθ=0` …………………….. (I)

Let, angular displacement, `θ=x(1)` …………………….. (II)

And, angular velocity, `(dθ)/dt=(dx(1))/dt=x(2)` …………………….. (III)

Now, `(d^2 θ)/(dt^2 )=(d^2 x(1))/(dt^2 )=(dx(2))/dt` …………………….. (IV)

On substituting the values of (II), (III), and (IV) in equation (I), we get:-

`(dx(2))/dt+(b/m)x(2)+(g/L)sin⁡x(1)=0` OR `(dx(2))/dt=-(b/m)x(2)-(g/L)sin⁡x(1)` ………………… (V)

Since we have converted our second-order ODE equation (I) into two first-order ODE equations (III) and (V). Now, we’ll input the first-order ODE equations in MATLAB; in form of matrix:-

Procedure for developing the program:-

1. We’re going to make a function that takes the input variables `t,θ,b,m,L,g` and will generate the required matrix ‘theta_dot’.
2. We’ll call the created function in our program to get the solution for second-order ODE for the pendulum with damping. As an output, we’ll get the values for angular velocity and angular displacement. Then we’ll plot the graphs for these solutions.
3. Finally, we’re going to make an animation to show the pendulum motion; for 0-20sec timespan.

Creating function ‘pendulum_ode_func’ in MATLAB/OCTAVE:-

% creating function for pendulum ODE
 
function [theta_dot]=pendulum_ode_func(t,x)
 
% given, mass(m)=1kg, damping coefficient(b)=0.05, length of pendulum(L)=1m, acceleration due to gravity(g)=9.81m/(s^2)
m=1; b=0.05; L=1; g=9.81;
 
% converting into matrix form which contain 2 rows and 1 column
theta_dot=zeros(2,1);
 
% let, theta=angular displacement=x(1) 
% theta_dot(1)=angular velocity=x(2)
% theta_dot(2)=angular accleration
 
% Breaking second order ODE into two first order ODE's
%1st ODE:- rate of change of angular displacement=angular velocity
theta_dot(1)=x(2);
 
%2nd ODE:- rate of change of angular velocity=angular acceleration
theta_dot(2)= -(b/m)*x(2)-(g/L)*sin(x(1)); 
end

Creating a program to solve pendulum ODE and to show the pendulum motion:-

clear all
close all
clc
 
% initial and boundary condition
% time span 0-20 sec
% angular displacement=0, angular velocity=3
 
t=[0 20];
theta_O=[0;3];
 
% solving ODE using 'ode45'
 
[t,theta]=ode45('pendulum_ode_func',t,theta_O);
 
% plotting the solution
% plotting the graph of angular displacement and angular velocity
figure(1)
plot (t,theta(:,1),'linewidth',2)
grid on; hold on;
plot(t,theta(:,2),'linewidth',2,'color','r')
xlabel('time')
ylabel('amplitude')
legend('Angular displacement','Angular velocity')
 
% Making Animation
 
ct=1;    % loop counter
 
for i=1:length(theta(:,1))
 
% Creating Coordinates of the bob
 
x0=0;
y0=0;
x1=1*sin(theta(i,1));
y1=-1*cos(theta(i,1));
 
% Plotting components used for animation
 
animation = figure(2);
% Plotting horizontal support of pendulum
plot([-0.5 0.5],[0 0],'linewidth',10,'color','b')
hold on
grid on
% Plotting vertical support of bob (string of length 1m)
line([x0 x1],[y0 y1],'linewidth',5,'color','m')
hold on
% Plotting bob
plot (x1,y1,'marker','O','markersize',20,'markerfacecolor','c')
hold off
 
axis([-2 2 -2 1.5])
 
pause(0.003); % using Pause so that plots would get time to plot another entity
 
% Generating .png files and saving them
 
text=sprintf("pendulum%05d.png",ct);
saveas(animation,text);
ct=ct+1;
end

RESULTS:-

CONCLUSION

Since we can observe that, with damping, the amplitude of our motion slowly decreases over time, and oscillation is fading away with time. Hence, we can conclude that damping lowers the natural frequency of an object and also decreases the magnitude of the amplitude of waves.