Simulation of transient behaviour of simple pendulum

Simulation of Transient behaviour of simple pendulum

Aim: To create a simulation of transient behvaiour of simple pendulum using matlab program.

Theory involved

Simple pendulam is having a string which is attached to a fixed end and a bob tied at the other end of string and it moves in Simple harmonic motion.Newton\'s second law is used to solve the forces acting on pendulum bob and that leads to formation of a differential equation which is of second order.The The ordinary differential equation is having two variables one is independant variable which is time and the other is dependant variable called angular displacement.The first differential of the angular displacement will give us the angular velocity of the pendulam.

b = damping constant; m = mass of the bob; g = acceleration due to gravity; L = length of the pendulam

For solving the differential equation,we need to break this differential equation into a system of some first order differential equations.

Here the second order differential equation can be break into two diffrential equations of first order.

(dθ2/dt)=-(b/m)θ2-(g/l)sin(θ1)....................................... (1)

(dθ1/dt) = θ2.......................................... (2)

This can be combined into a single matrix form.

Matlab Program

clc
clear all
close all
%damping constant
b = .05;
%accleration due to gravity
g = 9.81;
%length of the pendulam
l = 1;
%mass of the bob
m = 1;
%Intial conditions
theta_0 = [0;3];
%Integration time span
t_span =  linspace(0,20,500);

[t,results] = ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);
%plotting
figure(1)
%angular displacement vs time
plot(t,results(:,1))
hold on
%angular velocity vs time
plot(t,results(:,2))
xlabel('time')
ylabel('results')
legend('angular displacement','velocity')

%theta values
theta = results(:,1);
ct =1;
for i = 1:length(theta)
    THETA = theta(i);
    
    x0 = 0;
    y0 = 0;
    x1 = l*sin(THETA);
    y1 = -l*cos(THETA);
   figure(2)
   %plotting
    plot([x0 x1],[y0,y1],[-1 1],[0 0],'linewidth',3,'color','r')
    hold on
    %plotting the bob 
    plot(x1,y1,'o','MarkerSize',20,'color','k','Markerfacecolor','k')
    hold off
    axis([-2 2 -2 1]) 
    pause(.03)
    
    M(ct) = getframe(gcf);
    ct = ct+1;
end
movie(M)
videofile = VideoWriter('pendulam.aVI','uncompressed aVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Function Template

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];

Code Explanation of function template

1. A function template in name ode_func is made which can be called when needed.Function have the output dtheta_dt OR (dtheta/dt) and have the inputs theta,t,b,g,l and m.
2. Then we define the Theta1 and theta2 which is to be used in next step for defining the differential equations.
3. Then we define the two first order differential equations we made from the actual second order differential equation.

Code explanation of actual program

1. As usual the matlab program is started with clear all,close all,clc,for clearing the workspace,closing all plots and to empty the command window.
2. Then we define the inputs such as l,b,m and g.We also typed in the objective above the code using % sign.
3. Then we are defining the intial conditions for the system.Here the intial conditions include the intial angular displacement and intial angular veocity,that means the condition of the system at time t =0.And we defined it as theta_0.
4. Next to define our integration time span,How far we need to intergrate our differential equation.That is provided next and here we are integrating over a time of 20 seconds.And we need to create an animation of the pendulum.So we are using linspace command to define the t_span so that we can have as many values of time as required in between these 20 seconds with same interval.
5. Then we need to solve our differential equations.So we are using an inbuilt solver called ODE45 here to solve our equations.Outputs we required are t and the displacement and velocity.The displacement and velocity is obtained in colum matrix,since we write the code in such a way that the dtheta_dt is a colum matrix.
6. While using ode45 we need to first give the variables which is both independant and dependant variables,Here it is t and theta.Then we will be calling the function odc_func and given its inputs.Then we provide the intergration time span and then we enter the intial conditions.
7. Running the program now will give as the value of displacement and velocity in a colum matrix called results.
8. Then we define the figure(1),which means the first figure.
9. Then we use plot command to plot the curves between diplacement and velocity versus time.For diplacement vs time plot we need to use results(:,1) which means all the values in first colum of all rows.Similarly results(:,2) for velocity values.
10. x label and ylabel commands are used to label the axis.
11. Then we used the legend command to clearly differentiates the graphs as displacement and velocity.

After all this done we need to move on with creating the animation of pendulum.For that we need to set the fixed point of the pendulum,which is the orgin.And the end point of string should have a bob and it should move according to the angular displacement.So we will be creating images of the string with bob at differant instances and we will then stitch these plots using movie command in matlab.We use counter to store the images in an array.

1. We first define the angular displacement seperately as theta.
2. Then we use a for loop in which i is varying from 1 to length of theta, which is basically the number of rows in results matrix.So the program will take the values of theta one by one.
3. Then we define THETA as theta(i).So we can use the THETA in cosine and sine functions while entering the end point of the pendulum.
4. Then we inputs the fixed point as x0,y0 as 0,0. And we define the x1 and y1 which are the end points of the pendulum.The y cordinate must be on negative side since the pendulum is upside down.
5. Then we mention the figure 2 and we use plot command to plot the string of the pendulam and the base in which the pendulum is hung.The base is basically the x axis in which the y values are zero.And we also provide color to the plot.
6. We use hold on command since we need plot command to show both the bob and the string of pendulum.
7. Then we plot the bob using the marker and we specify the end point of pendulum to plot the marker and we given size of marker as 10 and we also gives the marker face color as black and marker color also as black
8. We don\'t need the plot command to show all the plot of bob at differant values of i,so we use hold off command to prevent this.
9. Then we specify the axis using axis command.
10. For taking the pause in between, the pause command is used.
11. M(ct) is used to store the figure in counter into a movie array.And we gives the code to increment the counter.
12. We provide the end for the for loop.
13. Movie command make the movie with the figures stored in the M array.
14. Then we save the video file in the name pendulam.avi.And then write the file.Then we close the file.

ERRORS FACED

This happened when i mistakely provide x1 as y1 and y1 as x1.And also i forget the negative sign while defining y1.When everything correctly entered the problem gets solved.

This happens when hold on command is not given ,as plot command only shows the latest command as outputs.When hold on command is given ,the problem is soved.

Here hold off command after plotting the bob is not given and this problem occured ,then i give the hold off command and the problem is solved.

Results

The plot of diplacement and velocity vs time is plotted.As damping constant inctreases the oscillations gets damped easily and we gets the diplacement as zero in graph, if the oscillations gets damped.

This is the graph we obtained as the result of the above problem.Here displacement is not completely damped but the amplitude is reducing as time passes due to damping factor.

Animation is obtained and the youtube link of the animation is provided below.

<iframe src="https://www.youtube.com/embed/NI2Cx0uPw7k" frameborder="0" allow="accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>