Simulation of the Transient Behavior of a Simple Pendulum and its animation

Q) Write a program to simulate the transient behaviour of a simple pendulum and to create an animation of it\'s motion.

ANS:

Theory of Simple Harmonic Motion

In engineering mathematics, Ordinary Differential Equations (ODEs) are used to characterize the transient behavior of systems. In this case, we can do the same for a simple pendulum system. A pendulum is basically an object of small mass, m (called a bob) attached to a string of length, L which is fixed at a point. The pendulum undergoes simple harmonic motion when it is displaced from its initial position from rest i.e. when velocity is 0. When it is displaced from the initial position (referred to as the equilibrium position), it will swing sideways and continue to do so until due to energy losses, it will again come to rest back to its initial position.

The maximum extremities to which the pendulum travels is the amplitude, x (with respect to the equilibrium position). The pendulum travels from -x to x or vice versa. Any object in simple harmonic motion is characterized by two fundamental ideas:

• The acceleration of the object is proportional to its displacement (from the equilibrium position)
• the acceleration of the object acts in the opposite direction to the displacement

These 2 fundamentals are denoted mathematically below:

`a = - (omega^2*x)`where a denotes acceleration, omega denotes the anguar velocity, and x represents the displacement from the equilibrium position.

The pendulum is also going through an energy change between kinetic energy and potential energy. When the pendulum is at the equilibrium position, it has no acceleration and is at constant velocity and therefore has no resultant force. Because the acceleration is 0, the maximum velocity is always achieved at this position. Hence, at this point, the kinetic energy is maximum.

Similarly when the pendulum reaches the maximum amplitude, the acceleration is maximum and has no velocity. As the pendulum travels from the equilibrium position to maximum amplitude, it is decelerating as it loses its velocity. When it reaches the amplitude however, the velocity is not 0 but is changing rapidly as it is about to change direction and travel in the opposite direction. Acceleration is defined as the rate of change of velocity. So acceleration is maximum at the amplitudes and therefore, so is the resultant force as it changes direction. Because velocity is minimum at the amplitudes, this is where the potential energy is maximum.

Solving the ODE

The image above, shows the ODE that characterizes the motion of the pendulum and it is what we are trying to solve. The method of solving any ODE is to first deconstruct it in to a set of first order differential equations. As you can see, our ODE is a second order one and so it can be broken down in to two first order differential equations. The way this is done is shown below:

let `theta = theta_1`

1. `d theta / dt = d theta_1 / dt = theta_2`. This is the first equation and knowing this relation, we can make the assumption in the next line.
2. `d^2 theta/ dt^2 = d^2 theta_1/dt^2 = d /dt (d theta_1/dt) = d / dt (theta_2)`. This is the second equation.

Substituting `theta_1` in the place of `theta`in the second order ODE, we get another equation:

`d^2 theta_1/dt^2 + b/m*d theta_1/dt + g/L*sin(theta_1)`

Using the relation from the second step, we get:

`d theta_2/dt + b/m*theta_2 + g/L*sin(theta_1)`. Rearranging this, we get:

`d theta_2/dt = -b/m*theta_2 - g/Lsin(theta_1)`

Now, we can start with the programming.

Program

clear all
close all
clc

% Conditions of pendulum properties

b = 0.05;
g = 9.81;
l = 1;
m = 1;

As usual, we always start the script with the clear all, close all and clc codes just to reset all variables and conditions. Then, we move on to the conditions of the pendulum properties. As shown in the question, the pendulum needs to have a damping coefficient of 0.05 denoted by the variable b. The acceleration due to gravity is denoted by the variable, g and is set to 9.81m/s². The length of the string attached to the bob of the pendulum should be 1m as shown by the variable, l. Finally, the pendulum should have a mass of 1kg denoted by the variable, m.

% Initial Conditions for Displacement and Velocity 

theta_0 = [0;3];

When solving any ODE, we need to have a set of initial conditions or boundary conditions to help us to make a start in the solving process. In our case, we are asked to give the pendulum an initial displacement of 0 degrees and an initial velocity of 3 radians per second. These conditions are put in a column matrix and assigned to the variable, theta_0 which is appropriately chosen to show that these are the conditions when theta (i.e. the displacement) is 0. These conditions will be used in the solving process later on.

% Time Array for time period

t_span = linspace (0, 20, 500);

Because we are trying to characterize the transient behavior of the pendulum system, we know that the system is changing with respect to time. So, we need to know for how long we will characterize this system for. As required, we need to use a time duration of 20 seconds. So, an array is created using the \'linspace\' command from 0 to 20 taking 500 increments in between. This is just so that we get more data for showing the transient behavior and in the end, will give us a much more smooth and detailed animation to view. This array is assigned under the variable name \'t_span\'.

% Solving of ODE using solver

[t, results] = ode45(@(t, theta) ode_func(t, theta, b, g, l, m), t_span, theta_0);

Now that we have all the boundary conditions and parameters set up, we can move on to the solving of the ODE. Since we will be solving using MATLAB, we will make use of its in-built solvers to help us. MATLAB has the ode45 solver which is used for this purpose. The solver is used so that we obtain the results as an array of time in one column, and the result (the displacement) in another. The syntax of the code actually calls for a function, \'ode_func\' which contains the parameters of the time, displacement, damping coefficient, acceleration due to gravity, length and mass. The end of the line of code shows how the boundary conditions set up earlier come in to use.

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

The above code, shows how the function \'ode_func\' was created. Like all functions, this code was written in another script file. The first line shows the basic syntax for any function. The following lines, show how the first order differential equations were created from the second order ODE. The mathematics behind that is just translated to code for out use. Variables theta1 and theta2 were created for the assumptions which led to the simplification of the second order ODE. The last line shows the function name and also denotes the equation `d theta_2/dt = -b/m*theta_2 - g/Lsin(theta_1)`in matrix form.

If the code is executed, after all of this, we should get an array in the workspace section containing the results we want as shown below.

As you can see, the first row shows the time at 0 and the angular velocity at 3 radians per second which is what we wanted from the boundary conditions. The results go on for 500 iterations due to the time array we created from 0 to 20 seconds with 500 increments.

Plotting of Transient Behavior

% Plotting of pendulum behavior

plot (t, results (:,1), \'linewidth\', 2, \'color\', \'red\', \'DisplayName\', \'Displacement (rad)\')

hold on
plot (t, results (:,2), \'linewidth\', 2, \'color\', \'blue\', \'DisplayName\', \'Angular Velocity (rad/s)\')

legend 
title (\'Transient Behavior of Simple Pendulum\')
ylabel (\'Plot\')
xlabel (\'Time\')

Using the plot command, we basically want to plot the results from solving the ODE from the results array as seen in the image above. The code in this section will plot 2 curves. One will show the angular displacement and the angular velocity, both varying with time. Because the plot will show 2 curves on the same axes, I added a legend using the \'DisplayName\' command and gave each respective curve its required key, for ease of understanding. The \'hold on\' command is used so that the first curve remains while the second curve is also plotted. The legend command is used to activate the \'DisplayName\' code.

Each curve is given a line width of 2 and different line colours for differentiation. A title and labels are also added for presentation. The graph I obtained is shown below:

Animation of Transient Behavior

% Creating array for Displacement results

disp = results (:,1);

In order to create the animation for the transient behavior of a pendulum, we first need an array of the displacement values obtained in the results array. So, in this line of code, a variable called \'disp\' was created with contains the results array with just the displacement values as we need these to show how the pendulum is moving with time.

% Creating Animation for Pendulum Swing

figure (2)
ct = 1; 

Next, to start the animation, we need to use the \'counter\' command with syntax \'ct\' and assign it to a value of 1. This command will take frames of each displacement of the pendulum. The figure command is just so that the animation appears in a second figure window.

For Loop

A \'for loop\' will be used in order to iterate the different displacements of the pendulum for the given time duration. The code for this is shown below:

for i = 1 : length (disp)
    
    % Starting position of pendulum
    x0 = 0;
    y0 = 0;
    
    % Position of pendulum as the anguar diplacement changes
    x1 = l*sin(disp(i) - pi);
    y1 = l*cos(disp(i) - pi);

Using the \'disp\' array created before, we use it in the for loop statement so that MATLAB will iterate the different positions of the pendulum within the time period in the time array from 1 to the length of the array (which is 500). Next, I created variables x0 and y0 and equated them to 0 to denote the coordinates of the connection between the pendulum string and the fixed support. So this will give us (0,0).

But, we also need a set of coordinates for the each displacement that the pendulum takes which is denoted by (x1, y1). In other words, (x1, y1) tells us the instantaneous position of the pendulum. Using trigonometry, I assigned x1 and y1 as the sine and cosine components of the length of the pendulum string which was assigned as 1 at the start.

However, since the displacement is changing with time, we need to use the angles (displacements) from the \'disp\' array and link it to the \'for loop\'. This is so that each set of coordinates for (x1, y1) will be iterated as the for loop runs through each displacement value in the \'disp\' array. I had to include a \'- pi\' also. This would subtract 180 degrees (pi in radians) from the angles. If this was not included, the animation of the pendulum would be inverted.

% Creating of Fixed Pendulum Support
    plot( [-1 1], [0.06 0.06], \'linewidth\', 14, \'color\', \'black\');
    hold on

Next, the different parts of the pendulum need to appear in the animation. For this, I used the plot command to plot them. Starting with the pendulum\'s fixed support, I just had to give 2 sets of coordinates so that I could plot a horizontal line. Then I gave it a line width of 14 for a better look and a black colour. Once again, the \'hold on\' command is so that this line remains.

% Creation of Pendulum Rope 
    plot( [x0 x1], [y0 y1], \'linewidth\', 2, \'color\', \'red\' );

Then, for the string/rope, I used the plot command to plot a straight line between the fixed support (0, 0) to the instantaneous position of the pendulum (x1, y1). A line width of 2 and a red color was given to this line.

% Creation of Pendulum Bob
    plot(x1, y1, \'-o\', \'markersize\', 14, \'markerfacecolor\', \'black\')
    
    
    axis([-1.5 1.5 -2 0.5])

For the pendulum bob, the plot command was used to create one coordinate which is the end of the string. This is the (x1, y1) coordinate point, so that was plotted and given the marker type \'o\' for a round shape. A marker size of 14 was used for bigger size using the \'markersize\' command.

Also the \'markerfacecolor\' command was used to give the pendulum bob a black colour. The axis command was used to set up the axes for the whole animation and to ensure everything is visible and in frame.

% Compiling of frames for animation
    pause (0.04)
    M(ct) = getframe(gcf);
    ct = ct + 1;
    hold off
     
end

For the compiling of all the frames, we need to use the \'getframe\' command which was used previously in the robot arm animation program. The counter command \'ct\' is assigned to add 1 so that it cumulatively adds frames to each set of data in the \'disp\' array. Since, we already created the variable \'ct\' and assigned it to a value of 1, each frame will add do that successively. A pause value of 0.04 is used to control the speed of the animation i.e. frames per second. The \'hold off\' command is needed so that each frame is stored but not shown during the animation. We only need the frame at the current time value i.e. the instantaneous frame. Finally, the end command is to denote the end of the \'for loop\'.

Movie File Creation

% Creation of Movie File

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

Once, we have got all the frames we need for the animation, we can create the movie file for the animation and store it as a video file. The \'movie\' command is used for this and a video file. The \'VideoWriter\' command writes the frames stored to a video file called \'pendulum\' which is in the uncompressed AVI format.

The image below shows the pendulum animation at the end of the time cycle.

Errors

• The first error came in the function script. In the 7th line, I forgot to include brackets while writing the code for the variable, \'dtheta2_dt\'. I had also used the degree command which was not necessary, since all angles were taken in radians anyway. This was easily fixed.

• The second error came about with the compiling of the frames with the \'getframe\' command. I had entered \'getframes\' instead and this appeared as an error in the command window. This was also easily fixed.

• The third error came in the code for the creation of the movie file for the animation. I forgot the case sensitivity in the \'VideoWriter\' command and entered \'videoWriter\' by mistake. This was also easily rectified.

Conclusion

The transient behavior of the simple pendulum system was simulated and its animation was created using MATLAB. The video file of the animation is uploaded on YouTube and the link is provided below:

https://youtu.be/QFq4QH1-09k