Motion of a simple pendulum with Damping

Aim: To write a program in MATLAB to simulate the simple pendulum motion between 0-20 seconds considering an initial angular displacement =0,initial angular velocity=3 rad/s at time t=0. The other parameter values are given and considered according to that : length of string l=1m,mass of the ball m=1kg,damping coefficient b=0.05 and gravity g=9.81 m/s^2.

Theory:

A simple pendulum basically consists of a mass m which in this case is the ball of 1kg,hanging from a string of length l which is 1 m and fixed at a pivot point on the other side. The way the pendulum moves depends on Newtons second law. When this law is written down,we get a second order ODE that describes the position of ball w.r.t time. In the code, we have broken down the second order ODE into two first order ODE\'s and then proceeded,as the solver ode45 can solve first order ODE\'s only. ODE\'s are used to describe the transient behaviour of a system which in this case is a simple pendulum system.

Procedure:

1)Let us first do the function definition portion: Here ,we have created a function ode_func that takes t,theta,b,g,l,m as inputs and outputs dtheta_dt which is an array storing every single value of displacement and velocity. As dtheta_dt=[dtheta1_dt;dtheta2_dt] where dtheta1_dt represents the first ODE and dtheta2_dt represents the second ODE after being converted to two first order ODE\'s. theta1=theta(1) is where initial angular displacement condition is specified and similarly,theta2=theta(2) is where initial angular velocity condition is specified. The \'function\' keyword is used to start a function and \'end\' keyword ends the function.

function [dtheta_dt]= ode_func(t,theta,b,g,l,m)
   theta1=theta(1) %initial condition for dispalcement
  
   theta2=theta(2) %initial condition for velocity
   %Apply first equation
   dtheta1_dt= theta2;
   %Apply second equation
   dtheta2_dt= -(b/m)*theta2 - (g/l)*sin(theta1);
   %Output of the function
   dtheta_dt=[dtheta1_dt;dtheta2_dt];
end

2)In this part,we basically do the function calling portion and get the solution for both ODE\'s,create plots and then create an animation of the pendulum motion.

i)Using theta_0=[0;3]; the initial angular displacement and angular velocity values are specified which is 0 and 3 rad/s respectively. It consists of two rows and one column.

ii)t_span is an array containing 200 different values between 0-20 seconds. It basically indicates the time,for which time integration would be performed in order to get the solution for the ODE\'s provided inside function definition portion ,by calling them this way:

[t,results]=ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0); where ode45 is a solver in matlab that can solve first order ODE\'s only and hence the conversion to two first order ODE\'s is necessary. The solution for both the ODE\'s is computed where , the first column would represent the solution for first ODE and bring out all angular displacement values,and similarly the second column is the solution for second ODE and outputs all angular velocity values.

iii)We then proceed to plot the curve that indicates displacement vs time and velocity vs time. When we say : plot(t,results(:,1)) we mean to plot the x-axis as time array and results(:,1) represents y-axis which is all the rows of the first column and is nothing,but the angular displacements. Similarly,plot(t,results(:,2)) represents all velocity values vs time array. Both displacement and velocity vs time are represented in a single figure window.

iv)We then define a second figure window by using figure(2) where we use \"for loop\" to extract each value of time with their respective displacements which is the solution and then animate it. We start with : for i=1:length(t) means the loop will run for each value of time upto 200 different values which are present here. t contains 200 rows and 1 column.

We then go on to define the coordinates of the single link which is being used here.

x0=0,y0=0 is the first point of the link(origin) and the second point of the link is defined by:

x1,y1 where x1=sin(results(i,1)) and y1=-cos(results(i,1)). One important thing to note here is that since the length of the pendulum makes an angle with y-axis and not x-axis while in motion,hence there is an inversion in the formula for x1,y1.

v)We then use plot([x0 x1],[y0 y1],\'linewidth\',3) to plot the link whose coordinates we have previously defined. We provide a linewidth of 3 so that the line is clearly visible. We use the hold on command to give it a hold and then we plot the ball using marker command at the end of the link i.e at point x1,y1 and use a green color to represent it by using the command : plot(x1,y1,\'o\',\'MarkerFaceColor\',\'g\'). We use the hold off command as well to maintain smooth functionality.

vi)We have seen an out of frame result when the axis provided is manual. Hence,checking the values of time and displacements from workspace,we can provide a suitable axis to keep things within our frame.

vii)We then proceed to label the plot using xlabel and ylabel and provide suitable understanding of the plot using legend command.

viii)M(ct) is nothing but an array which stores the frames we have captured using getframe(gcf) and ct is the counter variable taking care of all the frames to be captured in the animation. We then create a movie,name the video file pendulum_motion,open the video file ,write into the videofile and then close it finally.

The animation has been created and the youtube link for it is:

https://youtu.be/SQMA25t2POU

The above picture is one of the frames which have been captured during the output is running.

Errors:

1)When we use the formula : x1= cos(results(i,1)) and y1=-sin(results(i,1)) the output goes out of frame and makes an angle in the opposite direction as shown below:

In order to fix this: x1=sin(results(i,1)) and y1=-cos(results(i,1)) . This keeps the pendulum motion like the one shown in the figure and within our frame.

2) When we use just 50 values of time,the displacement vs time graph is not smooth . Using more points like here we have used 200 values,the curves are much more smoother.

%Solving 2nd order ODE\'s
clear all
close all
clc
%Assigning values to constants in simple pendulum equation
%We have checked that increasing the damping coefficient \'b\' by a factor of
%10 i.e b=0.2,dampens out the vibration really well
%damping coefficient,gravity,length of pendulum,mass of bob
b=0.05;
g=9.81;
l=1;
m=1;

%Initial conditions
%First value is initial displacement 0,second value is initial velocity 3
theta_0 = [0;3];
%Time points
t_span=linspace(0,20,200);
%For loop for extracting each value of t_span upto 20
%t_span1 = single values,whereas t_span= array

 %Solving the ODE
[t,results]= ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);
 plot(t,results(:,1))
 hold on
 plot(t,results(:,2))
 %--------------------------------------------%
 
 figure(2)
 ct=1;
 
for i=1:length(t)


%Co-ordinates of the origin(first point) 
x0=0;
 y0=0;
 %Co-ordinates of the second point
 x1= sin(results(i,1));
 y1= -cos(results(i,1));
 
 plot([x0 x1],[y0 y1],\'linewidth\',3)
 
 
 hold on
 p=plot(x1,y1,\'o\',\'MarkerFaceColor\',\'g\')
 hold off
 axis ([-1 1 -1 0])

 
  xlabel(\'Time\')
 ylabel(\'Disp or Velocity\')
 legend(\'Displacement\',\'ball\')
 pause(0.3)
 M(ct)=getframe(gcf);
 ct=ct+1;
 
 end

 
 
 movie(M)
 videofile=VideoWriter(\'pendulum_motion1.avi\',\'Uncompressed AVI\');
 open(videofile)
 writeVideo(videofile,M)
 close(videofile)

 

Inference:

The animation for a simple pendulum motion under damping is created and the youtube link has been provided above as well.