Simulation of the transient behaviour of a simple pendulum using MATLAB
Program to simulate transient behaviour of a simple pendulum using MATLAB
Aim: To program to simulate the transient behaviour of a simple pendulum and to create an animation of it's motion using MATLAB.
Given:
L = length of the pendulum = 1 m
m = mass of the ball = 1 Kg
b = damping coefficient = 0.5
g = acceleration due to gravity = 9.81 m/s^2
At time t = 0 sec, angular displacement = 0 rad, angular velocity = 3 rad/sec.
This above ODE represents the equation of motion of a simple pendulum with damping.
Program:
close all
clear all
clc
%input values(given)
% damping coefficient = b
b = 0.05;
%mass of the bob = m (kg)
m = 1;
%acceleration due to gravity = g (m/s^2)
g = 9.81;
%lenght of the pendumlum = L (m)
L = 1;
% initial displacement = 0 and initial velocity = 3 m/s at time = 0 sec.
%initial condition
Theta = [ 0 ; 3 ];
%time_constraints(given)
%time considered
time_span = linspace(0,20,100);
%solving 2nd order ODE and calling function : ode_func
[t , results] = ode45(@(t,theta) ode_func(t, theta, b, g, L, m),time_span,Theta);
figure(1);
angular_displacement = results(:,1);
angular_velocity = results (:,2);
figure(1)
% plotting graph for displacement and velocity of the bob w.r.t. time
plot(t,[angular_displacement angular_velocity])
xlabel('Time (sec)');
ylabel('Theta (rad and rad/sec)');
legend('ang. displacement (rad)','ang. velocity (rad/sec)')
title('Theta vs Time')Note: 'ode45' is an in-build function in MATLAB. The MATLAB function 'ode45' is based on a combination of fourth- and fifth-order Runge-Kutta methods.
This above program gives graph of theta1(ang. displacement) vs time and thata2(ang. velocity) vs time.
Function:
function [dtheta_dt] = ode_func(t,theta, b, g, L, m)
theta1 = theta(1);
theta2 = theta(2);
dtheta1_dt = theta2;
dtheta2_dt = -(b/m)*theta2 - (g/L)*sin(theta1);
dtheta_dt = [dtheta1_dt; dtheta2_dt];
end
Explanation:
1) This function takes the all the given values b, g, L, m as input.
2) theta(1) ie angular displacement gets stored in theta1 and theta(2) ie angular velocity gets stored in theta2.
3) Now the given 2nd Order ODE gets converted into two 1st order ODE. Two 1st Order ODE are
4) These two 1st Order ODE are then stored in an array 'dtheta_dt'.
This above function takes the given values as input to the function and convert the 2nd order ODE into two 1st order ODE (ie if there is nth order ODE this function converts it into 'n' number of 1st order ODE.)
Program to create the animation of motion of a simple pendulum:
close all
clear all
clc
%input values(given)
% damping coefficient = b
b = 0.05;
%mass of the bob = m (kg)
m = 1;
%acceleration due to gravity = g (m/s^2)
g = 9.81;
%lenght of the pendumlum = L (m)
L = 1;
% initial ang. displacement = 0 and initial ang. velocity = 3 m/s at time = 0 sec.
%initial condition
Theta = [ 0 ; 3 ];
%time_constraints(given)
%time considered
time_span = linspace(0,20,100);
%solving 2nd order ODE and calling function : ode_func
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,L,m),time_span,Theta);
angular_displacement = results(:,1);
angular_velocity = results (:,2);
ct = 1;
figure(1);
for i = 1:length(time_span)
dis(i) = results(i,1);
% x and y co-ordinates of bob
x(i) = L * sin(dis(i));
y(i) = -L * cos(dis(i));
xc = x(i);
yc = y(i);
figure(ct);
hold on ;
%taking the hinged support at (0,0)
%plotting a line between the hinged support and bob of the pendulum
plot([ xc 0],[ yc 0],'linewidth',4,'color','b');
plot([-0.6 0.6],[0 0],'linewidth',3,'color','y');
% plotting bob
plot(xc,yc,'O','markersize',25,'MarkerFaceColor','r');
%setting the axis
axis ([-1 1 -1.2 0.2]);
%adding one to counter for increasing the frames
pause(1.0);
grid on;
M(ct)=getframe(gcf);
ct = ct + 1;
end
%creating a movie
movie(M);
videofile = VideoWriter('simple pendulum.avi','Uncompressed AVI');
open(videofile);
writeVideo(videofile,M);
close(videofile)YouTube link for animation: https://www.youtube.com/watch?v=BBkE1A15q_Y
Explanation of above program:
1) Firstly, all the given parameters ie b, m, g, L are stated. Then the initial angular displacement and initial angular velocity are stored in column vector called Theta.
2) The motion of pendulum is studied from 0 sec to 20 sec.
3) The command time_span = linspace(0,20,100) gives equally spaced 100 values between 0 and 20 and get stored in time_span row array.
4) Now, the function ode45 is used. First the variables are defined t and theta then ode_func is called which passes the array 'dtheta_dt' then it is integrated from 0 to 20 sec for all those 100 values and initial values of ang. displacement and ang. velocity (ie. Theta array) is passed.
5) The function 'ode45' gives the 't' and 'results'. 't' is transpose of 'time_span' and 'results'
is 100*2 array in which first column represents the theta1 (ang. displacement) and second column represents the theta2 (ang. velocity) each corresponding to each value in 'time_span'.
6) For creating the animation, counter (ct) is used and is set to ct = 1.
7) Every thime the value in 'time_span' changes the position of pendulum changes, thus the position of pendulum should be captured. So for every change of position, new image should be stored.
8) So 'for' loop is used. The loop is executed for 100 times. Then displacement of bob is stored from results. X and Y co-ordinate of bob is plotted and assumming the pendulum is hinged at (0,0). From bob to the hinged support, the line is drawn by using 'linewidth' command. And spherical bob is drawn by 'O' and 'MakerFaceColor' command.
9) The X-axis is bounded from -1 to 1 and Y-axis is bounded from -1.2 to 0.2.
10) Every time the snap is stored there is pause of 0.8 sec.
11) Snap gets stored in 'M' and ct is increamented by 1.
12) After the loop is executed for 100 times ie 100 images are stored in 'M', movie(M) is executed.
13) VideoWriter is command in MATLAB which is used for creating the video of motion of pendulum. It is stored in same folder where the program is saved.
Refrences:
1) https://www.tutorialspoint.com/matlab/
2) https://in.mathworks.com/help/matlab/