Simulation of a simple pendulum motion with damping using MATLAB
Aim:
To simulate the motion of a simple pendulum and create an animation using MATLAB
Introduction:
(a) Motion of a simple pendulum:
In Engineering, ODE is used to describe the transient behavior of a system. A simple example is a pendulum.
A simple pendulum is a system consisting of a mass at one end of string while the other end is fixed to a support. It is initially in vertical position, which is the equilibrium. The mass is called 'bob' and when the pendulum is displaced by an angle `theta` and released, the force of gravity pulls it back towards its resting position. Its momentum causes it to overshoot and come to an angle `-theta` (if there are no frictional forces), and so on. The restoring force along the motion of the pendulum due to gravity is `-mg sin theta`. Thus, according to Newton's second law, the mass times the acceleration must equal `-mg sintheta`. This gives rise to periodic motion.
Source : www.physicsclassroom.com
The above image is a simple description of the pendulum motion, where
l is the length of the pendulum,
m is the mass of the pendulum.
(b) Governing equation:
The motion of the simple pendulum can be studied using a second order differential equation as follows:
`(d^2theta)/(dt^2) + (b/m)*(d theta)/(dt) + (g/l)* sin theta = 0`
Here,
`theta` = Angular displacement in radians
`b` = Damping co-efficient
`m` = Mass of the pendulum in kg
`g` = Acceleration due to gravity in `m/(s^2)`
`l` = Length of the pendulum in metres
(c) Solving a differential equation in MATLAB:
In order to solve an ordinary differential equation, MATLAB uses a function called ode45. But to solve a second order differential equation, we need to convert it into two First order differential equations.
(1) Let `theta = theta_1`
On differentiating the above relation twice,
`(d theta)/(dt) = (d theta_1)/(dt) = theta_2`
`(d^2theta)/(d t^2) = (d^2 theta_1)/(dt^2) = d/dt((d theta_1)/dt)`
`d/dt((d theta_1)/dt) = d/dt(theta_2)`
Hence, `(d theta_1)/(dt) = theta_2`
The above equation is the first differential equation.
(2) `(d^2theta)/(dt^2) + (b/m)*(d theta)/(dt) + (g/l)* sin theta = 0`
Based on the considerations made, the above relation will become
`(d^2theta_1)/(dt^2) + (b/m)*(d theta_1)/(dt) + (g/l)* sin theta_1 = 0`
`(d theta_2)/(dt^2) + (b/m)*(theta_2)+ (g/l)* sin theta_1 = 0`
`(d theta_2)/(dt^2) = - (b/m)*(theta_2) - (g/l)* sin theta_1`
The above equation is the second differential equation.
Hence, the motion of the simple pendulum can now be visualized as:
`d/dt[ (theta_1),(theta_2)] = [(theta_2), (-b/m theta_2 - g/l sin theta_1)]`
(d) Initial conditions:
The pendulum is displaced with angular displacement of zero and angular velocity of 3 rad/sec.
Other initial parameters are :
- Length l = 1.5 m,
- Mass m = 1 kg,
- Damping co-efficient b = 0.5,
- Acceleration due to gravity = 9.81 m/s^2.
MATLAB Program:
%Solving second order ODE pendulum motion
close all
clear all
clc
%Inputs
%Damping co-efficient
b = 0.5;
%Acceleration due to gravity in m/s^2
g = 9.81;
%Length of the pendulum in metres
l = 1.5;
%Mass of the bob in kg
m = 1;
%Initial conditions
%Angular displacement is zero
%Angular Velocity is 3 rad/sec
theta_1= [0;3];
%time - 0 to 20 seconds
time = linspace(0,25,500);
%Solving ODE using ode45 function
[t, result] = ode45(@(t,theta) ode_function(t,theta,b,g,l,m), time,theta_1);
%Plotting
figure(1)
plot(t,result(:,1))
hold on
plot(t,result(:,2))
xlabel('Time in seconds')
ylabel('Plot')
legend('Displacement','Velocity')
theta_2 = result(:,1);
%Counter variable
k = 1;
for i = 1:length(theta_2)
theta_3 = theta_2(i);
x0 = 0;
y0 = 0;
%x co-ordinate of the bob
x1 = l*sin(theta_3);
%y co-ordinate of the bob
y1 = -l*cos(theta_3);
figure(2)
plot([-1 1],[0 0],'linewidth',2);
hold on
plot([x0 x1],[y0 y1],'color','b','linewidth',2)
hold on
plot(x1,y1,'o','markers',20,'markerfacecolor','r')
hold off
grid on
axis([-2 2 -2 1])
pause(0.3)
%To get and store the frames
M(k) = getframe(gcf);
k = k+1;
end
%To create a video
movie(M)
videofile = VideoWriter('Simulation of pendulum motion.avi','uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)
Program for ode_function:
%ode_function
function [dtheta_dt] = ode_function(t,theta,b,g,l,m)
theta_1 = theta(1);
theta_2 = theta(2);
%First ODE
dtheta1_dt = theta_2;
%Second ODE
dtheta2_dt = -(b/m)*theta_2 - (g/l)*sin(theta_1);
dtheta_dt = [dtheta1_dt; dtheta2_dt];
endExplanation:
The initial conditions of the pendulum are given as per the problem statement.
An array is used to define the Angular displacement and Angular velocity and is assigned to a variable 'theta_1'.
The time is defined using 'linspace' command from 0 to 25 seconds.
Since the motion of a pendulum is a 'second order differential equation' , to solve it, we use a function called 'ode45', which integrates a system of differential equations from '0 to 25 seconds' with the initial conditions specified under 'theta_1'
A function called 'ode_function' is defined to calculate the displacement and velocity during the movement of the pendulum for the specified time and is assigned to a variable called 'result'.
From the calculated values, a plot of 'Displacement and Velocity against Time' is drawn as shown in the 'Results'.
To animate the simulation, 'For loop' is used. To get frames 'Counter variable k' is defined and for iterations, a variable 'i' is used.
At any instant of time, the position of the bob is given by '(x1,y1) = (l *sin(theta_3), -l*cos(theta_3))' .
Every movement of the pendulum is captured as a separate frame using 'getframe' command and using the 'movie' command, a video is created.
Results:
The above plot is the result of first plot command used in the code.
Explanation:
The initial displacement is 0 and angular velocity is 3 rad/sec and hence the starting point of the curves are from the respective values.
Due to the damping co-efficient, the oscillation attains equilibium within 25 seconds.
Legends are used to indicate the respective curves.
The oscillation has taken a sinusoidal waveform with respect to time.
The animation of my simulation was YouTube and the link is as follows:
Conclusion:
Motion of a simple pendulum based on the given conditions was programmed using MATLAB and the simulation was animated. A detailed explanation of the code and results was given with the help of output plots. The animation was uploaded to YouTube and the link is also shared in the report.