Simulation of Transient Behaviour of Simple Pendulum and Creation of an Animation of Motion of Simple Pendulum using Matlab

Code for simulation of simple pendulum and creating animation of motion of simple pendulum:-

% Part 1
% Program for both simulation of transient behaviour of simple pendulum and
% creation of animation of motion
clear all
close all
clc

% Defining given inputs in the problem statement

% Damping Coefficient
b = 0.05;

% Gravitational Acceleration (m/s^2)
g = 9.81;

% Length of string of the pendulum (m)
L = 1;

% Mass of bob of the pendulum (kg)
m = 1;

% Inital values of angular displacement and angular velocity
theta = [0; 3];

% Defining time interval (second) 
t_interval = linspace(0, 20, 1000);

% Solving ODE by calling a function named as ode_function_simple_pendulum
[t, results] = ode45(@(t, theta) ode_function_simple_pendulum(t, theta, b, g, L ,m), t_interval, theta);

% Plotting the first plot which shows variation of angular displacement
% and angular velocity wrt time
figure (1)
hold on
plot(t, results(:, 1), \'color\', \'r\', \'linewidth\', 2)
plot(t, results(:, 2), \'color\', \'b\', \'linewidth\', 2)
ylabel(\'Parameter\')
xlabel(\'Time\')
legend(\'Angular Displacement\', \'Angular Velocity\')
grid on


% Plotting the second figure which shows animation
ct = 1; 

% Defining for loop for representing different positions of pendulum
for i=1:length(t_interval)
  
  %Defining the coordinates of fixed end of string of the pendulum
  X0 = 0;
  Y0 = 0;

  %Defining the coordinates of free end of string of the pendulum
  X(i) = sin(results(i, 1))*L;
  Y(i) = -cos(results(i, 1))*L;
  
  % Plot for animation
  figure (2)
  clf
  hold on
  
  % Plot for creating support 
  plot([-1 1], [0 0], \'linewidth\', 3, \'color\', \'r\')
  
  
  % Plot for creating string of the pendulum
  plot([X0 X(i)], [Y0 Y(i)],\'linewidth\', 2, \'color\', \'k\');
 
  
  % Plot for creating bob of the pendulum
  plot(X(i), Y(i), \'marker\', \'o\', \'MarkerSize\', 15, \'MarkerFaceColor\', \'m\', \'MarkerEdgeColor\', \'k\');
  grid on
  axis([-1.5 1.5 -1.5 1.0])
  
  % Displaying the angular displacement and angular velocity in animation
  Angular_Displacement = results(i, 1);
  Angular_Velocity = results(i, 2);
  txt1 = [\'Angular Displacement : \' num2str(Angular_Displacement) \' rad\'];
  txt2 = [\'Angular Velocity : \' num2str(Angular_Velocity) \' rad/s\'];
  text(-0.5, 0.5, 3, txt1, \'color\', \'r\')
  text(-0.5, 0.2, 3, txt2, \'color\', \'b\')
  
  % Defining initials for the creating animation
  M(ct) = getframe(gcf);
  ct = ct + 1;
end

% Syntax for creating animation of pendulum
movie(M)
videofile=VideoWriter(\'Motion of Simple Pendulum.avi\', \'Uncompressed AVI\');
open (videofile)
writeVideo (videofile, M)
close (videofile)


% Part 2
% Program to define function for Ordinary Differential Equation

function [dtheta_dt] = ode_function_simple_pendulum(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

Stepwise Explanation of Code:-

PART I - Main Program

1) First Step - Use of basic commands Description -

clear all % Used to clear the workspace

close all % Used to close all the plots

clc % Used to clear the command window

Output of the step - Workspace and command window gets closed. All the open plots get closed.

2) Second Step -Defining the inputs of the problem Description:-

b = 0.05 % Damping coefficient

g = 9.81 % Gravitational acceleration (meter/second^2)

L = 1 % Length of the string of the pendulum (meter)

m = 1 % Mass of bob of the pendulum (kg)

theta = [0; 3] % Column vector/matrix with first element of first row as angular displacement and first element of second row as angular velocity defined using concatenation operator

time_interval = linspace(0, 20, 1000) %Time interval of 20 seconds divided into 1000 parts and presented as a row vector

Output of the step:- All the inputs get defined and time interval is presented as row vector with 1000 columns usng linspace command.

3) Third Step - Solving of Ordinary Differential Equations using ode45 solver and by calling a function defined earlier Description:-

t % Time array as evaluation points used to perform integration from t_initial (i.e 0 sec) to t_final (i.e 20 sec)

results % Solution array having 1000 rows and 2 columns

ode45 % Solver in MATLAB used to solve ordinary differential equations (ode)

@(t, theta) % Used to defined that ode is a function of t and theta(array)

ode_function_simple_pendulum(t, theta, b, g, L, m) % Function defined earlier which is called upon to solve ode

t_interval % Time interval with inital value 0 second and final value as 20 seconds which are limits of integration

theta % Initial Conditions (i.e angular displacement and angular velocity)

Output of the step:- 't' gets generated as a column vector with 1000 rows (i.e 1000 intervals) and 'results' gets generated as a 1000x2 matrix (i.e matrix of 1000 rows and 2 columns).

4) Fourth Step - Plotting of variation of angular displacement and angular velocity wrt time Description:- Plot command is used to plot the variation of angular displacement and angular velocity wrt time using different colors. 'xlabel' and 'ylabel' commands are used to name the x and y- axis. 'legend' command is used to name the plot lines such as red color line represnts angular displacement and blue color line represents angular velocity. 'grid on' command is used to make the grid on.

Output of the step:- Plot of variation of variation of angular displacement and angular velocity wrt time gets created and x and y axis gets labelled and legends appeared in the plot and grid gets on.

5) Fifth Step:- Using for loop for changing the value of angle of inclination of string with the vertical axis and hence creating various positions of string and bob by defining the starting and end positions of the string and creating animation Description -

'for' loop % Used to vary the angle of inclination of string with vertical axis using variable 'i' that varies from 1 to length(t_interval) i.e 1000 which means 1000 iterations i.e different positions of string and bob of the pendulum are created

X0, Y0 % x and y coordinates of the fixed end of the string of the pendulum

X(i), Y(i) % x and y coordinates of the free end of string of the pendulum giving 1000 values as variable 'i' varies

Output of the step:- After using for loop 1000 different positions of the string gets defined for 1000 values of variable i.

6) Sixth Step:- Use of clf, hold on, plot and axis command inside the for loop Description -

clf % Use to clear the current figure window mainly used twhen animation is required to be done to print the single plot

plot([-1 1], [0 0] %Use of plot command to plot the setup containing fixed support

plot([X0 X(i)], [Y0 Y(i)]) % Use to plot the string and as value of 'i' changes in every iteration a different position is created

plot(X(i), Y(i)) % Use to plot the bob of the pendulum which has coordinates same as that of free end of the string

hold on % Used to have all the plots in the same figure

axis % Used to specify the limits of x and y axis

Output of the step:- Figure containing fixed support, string and bob of the pendulum is obtained for the particular iteration of the for loop.

7) Seventh Step:- Displaying the angular displacement and angular velocity in animation with values varying according to the position of the pendulum Desription -

Angular_Displacement = results(i, 1) %angular displacement is defined as first column of the solution array

Angular_Velocity = results(i, 2) %angular velocity is defined as second column of the solution array

txt1 and txt2 % name of the variables

num2str % command for converting numeric value to character array

text % command used to display values of angular displacement and angular velocity in the plot of animation of simple pendulum

Output of the step :- Angular displacement and angular velocity gets presented on the plot of animation showing variation of values according to the position of the pendulum.

8) Eighth Step - Capturing of frames Description:-

M(ct) = getframe(gcf) % To capture the current figure window/frame

ct % counter variable with initial value defined as 1

Output of the step - Frames are captured for each iteration and hence 1000 frames are captured for 1000 different values of variable 'i'.

9) Ninth Step - Termination of for loop Description:-

end % Used to terminate the particular iteration after the value of ct i.e counter variable gets incremented

Output of the step:- The 'for' loop gets terminated.

10) Tenth Step - Defining syntax for the creation of animation Description:-

movie (M) - to play the movie defined by the frames captured using 'getframe(gcf)' function

videofile - variable to define the video name and format using 'VideoWriter' function

Output of the step :- After this, using 'VideoWriter', animattion video file is defined and named and then using 'open', 'writeVideo' and 'close' comands the video file gets generated.

PART II - Program for defining function named as ode_function_simple_pendulum

1)First Step - Definition of the function
Description:-
Function named as ode_function_simple_pendulum and its input variables i.e time array(t), initial conditions(theta), damping coefficient(b), gravitational acceleration(g), length of the string(L) and mass of bob of the pendulum(m) and its ouput variable i.e array(dtheta_dt) are defined.

Output of the step:- Function definition gets created.

2)Second Step - Defining input parameters
Description:-

theta1 = theta(1) % First element of first row of the array named as theta and that element is angular displacement

theta2 = theta(2) % First element of second row of the array named as theta and t that element is angular velocity

dtheta1_dt = theta2 % Derivative of variable theta1

dtheta2_dt % Derivative of variable theta2

Output of the step:- Other parameters are defined.

3)Third Step - Calculation of dtheta_dt
Description:-

dtheta_dt = [dtheta1_dt; dtheta2_dt] %array of derivatives calculated using all the input variables

Output of the step :- Array containing values of angular displacement and angular velocity gets calculated using input parameters and at the different values of timei.e t.

4)Fourth Step - Terminating the function
Description:-

end % Keyword is used to end the function

Output of the step:- Function gets terminated.

Output of the program:- Simulation of simple pendulum i.e plot showing variation of angular displacement and angular velocity wrt time gets obtained and also animation of motion of simple pendulum gets created.

Plot showing simulation -

Animation -