clear all
close all
clc
%% PART 1: INPUTS:
% given data:
b=0.05; % damping Coefficient.
g=9.81; % Acceleration due to gravity(N/m^2).
l=1; % Pendulum length(m).
m=1; % pendulum Mass(Kg).
% initial conditions
theta_0 = [0;3];
% time points
t_span = linspace(0,20,500);
%% PART 2: Solving the differential Equation of the pendulum:
% solve ODE
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);
ct=1;
for i=1:length(results)
% Pendulum fixed end(origin).
x0=0;
y0=0;
% Pendulum free end.
a=l*sin(results(i)); % X-coordinate of Pendulum free end.
b=-l*cos(results(i)); % Y-coordinate of Pendulum free end.
plot([-1 1],[0 0],'linewidth',4,'color','b');
hold on;
line([x0 a],[y0 b],'linewidth',2,'color','r'); % Pendulum Rod Line(or Length)
hold on;
plot(a,b,'o','markers',20,'markerfacecolor','g') %Pendulum Bob.
axis([-2 2 -2 2]);
xlabel('X-Axis')
ylabel('Y-Axis')
M(ct)=getframe(gcf);
hold off
ct=ct+1;
end
%% PART 3: Angular Position, Velocity Plots, and Animation of Pendulum Motion.
% Animation of Pendulum Motion.
movie(M)
videofile = VideoWriter('pendulum animation.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)
% Angular Position and Velocity plot against time(t).
figure(2)
plot(t,results(:,1));
hold on
plot(t,results(:,2));
xlabel('time')
ylabel('plot')
legend('Angular position','Angular Velocity');
Explanation of the code:
The code is divided into 3 parts, namely:
PART 1: INPUTS:
This part contains all the given data, as follows:
a) Line 7, damping coefficient(b) = 0.05.
b) Line 8, Acceleration due to gravity(g) = 9.81.
c) Line 9, pendulum length(l) = 1.
d) Line 10, pendulum mass(m) = 1 Kg.
Line 13, initial values for the position and velocity of the pendulum bob are named and stored as ‘theta_0’.
Line 18, ‘t_span’ is a row vector for time, its value increases from 0 to 20 sec linearly divided into 500 timesteps.
These timesteps will be used later in the code to calculate pendulum angular position and velocity.
PART 2: Solving the differential Equation of the pendulum:
Line 24, the MATLAB inbuilt function ‘ode45’ is called which will solve the differential equation of the pendulum and return as output the values of the pendulum position and velocity ‘[t, results]’ at different time values taken up from vector ‘t_span’.
The function is called by providing 3 inputs into it.
First is ‘ode_func’, it is a user-defined function which in the present case will solve the differential equation of the pendulum.
The code for ‘ode_func’ is provided below:
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
Line 1, function output is ‘dtheta_dt’, followed by function name and after that the input variables necessary for its execution.
Line 3-6, are the formulation of the differential equation and solution of these equations as ‘dtheta_dt’, which is a column vector storing values of the velocity and acceleration.
Second, is the time-span for which the velocity and acceleration equation will be integrated by ode45.
Third is the initial or the start values of displacement and velocity which are needed to start the integration.
Line 26, the program counter (ct) is defined and set to1.
Line 27, for-loop is defined, for-loop counter running from 1 to length of the ‘results’ matrix obtained by ode45.
Line 31, 32 represent the coordinates of the fixed end of the pendulum.
Line36, ‘a’ is the X-coordinate of the free end of the pendulum (Angular Position). This value is the first column of the results matrix from ode45. For the corresponding value of ‘i’, the value of ‘a’ (Angular position) is drawn the ith row of the ‘results’ matrix.
Line 37, Y-coordinate of the free end (b), calculated in the same way as the X-coordinate. However, Y-coordinate actually is in negative Y-direction so (-1) is multiplied in the beginning.
Line 39, plotting a line from (-1,0) to (1,0), which will be acting as a fixed datum.
Line 40, ‘hold on’ will prevent the above-generated line from being overwritten.
Line 42, pendulum rod is plotted with the fixed end as (x0,y0) on the mid-point of the blue line. Free end as (a, -b) distance between fixed and free end being fixed equal to ‘l’ (equal to 1 m).
Line 43, ‘hold on’ to prevent overwriting.
Line 45, creates the pendulum bob at the free end using ‘markers’ command. In this, ‘a’ and ‘b’ are respective coordinates for the bob ‘o’ is for picking the shape 20 is the bob size, and color is filled using ‘markerfacecolor’.
Line 46, the extent of both axes is fixed from -2 to 2 units respectively.
Line 47, 48 label the axes X-axis and Y-Axis respectively.
Line 49, ‘M(ct)’ is a variable that will store the plots generated as its elements, corresponding to different values of ct using the ‘getframe(gcf)’ command.
Line 50, hold command is turned off, it is necessary to turn it off because if not done all the frames of pendulum motion will get superimposed.
Line 51, ‘ct’ is incremented by 1 for a new iteration of the loop.
PART: 3: Angular Position, Velocity Plots, and Animation of Pendulum Motion:
Line 58, the function ‘movie(M)’ is called which requires ‘M’ as input.
Line 59, videofile is created using the function ‘VideoWriter’, the name and the profile type of the video are mentioned in the brackets.
Line 60, ‘open’ command opens the videofile.
Line 61, commands to write in all the frames in the row vector, M, into the videofile.
Line 62, closes the video file after completion of the video writing process.
Line 66, a new figure is created with the name ‘figure2’.
Line 67, plots the variation of the first column of ‘results’ matrix (Angular Position) w.r.t time (t).
Line 69, plots the variation of the second column of ‘results’ matrix (Angular Velocity) w.r.t time t).
Line 68, superimposes both the plots in a single figure for better understanding.
Line 70, 71 provide labels for respective axes.
Line 72, distinguishes both the plots by providing legend to respective plots, red color for Angular Position and blue color for Angular Velocity.
There was one error encountered while writing this code, which was not giving the ‘hold off’ command at line 50. This led to the superimposition of all the plots of the position of the pendulum in the animation.
The error screenshot is also attached for better understanding.
The error is rectified by adding the 'hold off' command before the program counter is (ct) incremented.
Results of this code are as follows:
a) Animation of the motion of the pendulum during the time-span provided in the code.
b) Plot of the variation of the angular position of the pendulum bob, the blue curve, and the variation of the angular velocity of the bob the red curve w.r.t the time. The image of the plot is attached below.
The link to the video of the animation of the pendulum motion is attached below:
https://youtu.be/Stx7XuatuV0