To simulate the motion of a simple pendulum by solving 2nd order ODE by Matlab/Octave programming.

AIM

To simulate the motion of a simple pendulum by solving its 2^nd order ODE (ordinary differential equation) using Matlab/Octave programming.

GOVERNING EQUATIONS

`(del^2theta)/(delt^2)+(b/m)*((deltheta)/(delt))+(g/l)*sintheta`

g = acceleration due to gravity in m/s^2

l = length of the pendulum in metres

m = mass of the ball/bob in kg

b = damping co-efficient

ODE WRT SIMPLE PENDULUM

ODE stands for ordinary differential equations. These equations are basically used to describe the transient behaviour of a system. These are differential equations containing one or multiple functions of one independent variable and their derivatives.

A simple pendulum in basic terms is a bob/mass suspended from a string/rod which can oscillate to and fro freely about a particular axis. When displacement is given externally, the pendulum swings to and fro until the force gets diminished.

When we consider the ODE of a simple pendulum, we consider a mass suspended by a spring of a particular length. The angle the string makes with the vertical with respect to time is also taken into consideration. We have seen the ODE of 2^nd order for a simple pendulum above. In Matlab/Octave programming, we cannot take up those equations directly into the program to

solve them. We need to break the differential equation up and then make the code accordingly. The final solution of the differential equation based on which the code is written is illustrated below:

Let `theta=theta1`from the main governing equation

The equation after partial differentiation becomes`(deltheta)/(delt)= (deltheta1)/(delt)=theta2`

Performing double differentiation, we get `(del^2theta)/(delt^2)=(del^2theta1)/(delt^2)=del/(delt)*((deltheta1)/(delt))`

From the original equation we get `(del^2theta1)/(delt^2)+(b/m)*(deltheta1)/(delt)+(g/l)*sintheta1=0`

Substituting the value of `theta1`, we get

`(deltheta2)/(delt)+(b/m)*theta2+(g/l)*sintheta1=0`

The two obtained Ordinary Differential Equation (ODE) are hence given as follows:

``(del/(delt))*theta1=theta2` and `((del)/(delt))*theta2=(-b/m)*theta2-(g/l)*sintheta`

OBJECTIVES

The main objective of this project is to obtain the simulated motion or animation of a simple pendulum through Matlab/Octave programming. The code can also be used to obtain the graph of angular velocity and displacement as a function of time. The code has two parts, one part where the code for the plots and animation is written and the other part where the function with respect to the 2^nd order differential equation is called from.

EXPLANATION OF THE CODE

The following program is written in OCTAVE.

The code consists of two portions, the main code and the code consisting of the function.

1) Code used for calling the function

function[theta_final]=diffeqnfunc(theta,t) – The function command is used to call the values in another code. ‘diffeqnfunc’ is the name of the saved file which must be the same as used in the function. Theta_final is the final simplified form of the ODE based on which the main code is written.

g=9.81;

b=0.05;

l=1;

m=1; - These are the standard values of the parameters as mentioned in the question.

theta1=theta(1);

theta2=theta(2);

d_dt_Theta1 = theta2;

d_dt_Theta2 = -((b/m)*theta2)-((g/l)*sin(theta1));

theta_final=[d_dt_Theta1;d_dt_Theta2];

end - These statements are used for the simplification of the second order ODE as explained in the ODE portion of the report. Theta1 and theta2 are the variables that are used to store the two values of the input, theta. D_dt_Theta1 and d_dt_Theta2 are the variables that store the final simplified form of the equation in the form of a matrix so that the function can be called, each parameter can be separated and the simulation can be done easily. Theta_final is used to store the values in a matrix form. The end command is used to close the function.

2) Code for simulation

g=9.81;

b=0.05;

m=1;

l=1; - These are the input standard values as given in the question.

theta=[0;3];

time=linspace(0,20,100); - These statements are used to define the limits of angular velocity(theta) and time. The values of time are taken from 0 to 20 seconds with 100 equal spaced values in between using the linspace command. The angular displacement is taken from 0 to 3.

[Final_Result]=lsode('diffeqnfunc',theta,time); - This statement is the most important statement in this code. This statement is used to call the ODE function from the ‘diffeqnfunc’ code and implement it here in order to get the solution of the ODE. The ain use of the ‘lsode’ statement is to return the value of the corresponding option when one argument is given. We can also set option parameters for the ‘lsode’ function when called with two arguments.

figure(1)

plot(time,Final_Result(:,1),'r') – This statement is used to plot the result of angular displacement as a function of time. The lsode statement enables us to consider only the first part of the row matrix, ie, the angular displacement part and plot a graph with respect to time.

hold on – Used to plot both graphs in a single sheet.

plot(time,Final_Result(:,2),'b') - – This statement is used to plot the result of angular velocity as a function of time. The lsode statement enables us to consider only the first part of the row matrix, ie, the angular velocity part and plot a graph with respect to time.

xlabel('Time') – Specifying values contained in X-axis

ylabel('angular displacement/velocity') – Specifying values contained in Y-axis

legend('Angular Displacement','Angular Velocity') – Denoting the colors on the graph with their respective representation.

c=1; - Dummy variable used for repetitive iterations to save all plots for animation.

A=Final_Result(:,1); A is a variable which is used to store the displacements with respect to time. Here, there are 100 values of time as mentioned in the linspace command. Hence, there will be 100 positions of the bob.

for i=1:length(A) – The for loop runs from 1 to the total number of values held by A which amount to 100.

B=A(i); - Stores each value of A in B

x0=0;

y0=0;

x1=l*(sin(B));

y1=-l*(cos(B)); - Specifying the position of the pendulum with respect to change in their angular positions. The calculation of the length of the string is also done.

figure(2)

plot([-0.5 0.5],[0 0],'color','b') – Plotting the base of the pendulum on which the string is attached. The X-axis varies from -0.5 to 0.5 and the Y-axis is 0.

hold on – Used to get all plots in a single frame/graph.

plot([0,0],[0,0],'o','color','r')

hold on

plot([x0 x1],[y0 y1],'color','g') – Plotting the string with its length as calculated above.

hold on

plot(x1,y1,'o','markers',50,'color','r') – Used to make a bob that is attached to the string. The ’o’ command is used to create a circular bob with a size of 50.

hold off

axis([-1.5 1.5 -1.5 1]) – Specifying the axis within which the simulation must occur.

pause(0.02) -Denotes the speed of motion of the frames of motion.

textfile=sprintf('animation%05d.png',c); - Used to save each plot of motion as a png file with the name ‘animation’ followed by 5 digit number in order which is denoted by %05d . This is the function of the sprintf command.

saveas(gca,textfile); - saves each plot after getting its current axes.

c=c+1; - Used to plot the next position graph for the next iteration.

End – Used to close the function.

CODE

The code in written in octave and consists of two parts, one function code and another main code.

Code used as a function to be called in the main code

% Function for obtaining the simplified version of the 2nd order ODE

function[theta_final]=diffeqnfunc(theta,t)
g=9.81;
b=0.05;
l=1;
m=1;

% Simplification of the formula

theta1=theta(1);
theta2=theta(2);
d_dt_Theta1 = theta2;
d_dt_Theta2 = -((b/m)*theta2)-((g/l)*sin(theta1));

% Arranging the values in form of a matrix

theta_final=[d_dt_Theta1;d_dt_Theta2];
end

Main code used for simulating the motion of simple pendulum

% Standard inputs as given in the question

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


% Specifying limits of angular velocity and time
theta=[0;3];
time=linspace(0,20,100);

% Calling the function 

[Final_Result]=lsode('diffeqnfunc',theta,time);

% Plotting the figures

figure(1)
plot(time,Final_Result(:,1),'r')
hold on
plot(time,Final_Result(:,2),'b')
xlabel('Time')
ylabel('angular displacement/velocity')
legend('Angular Displacement','Angular Velocity')

% Animation Plot

c=1;
A=Final_Result(:,1);
for i=1:length(A)
  B=A(i);
  
  x0=0;
  y0=0;
  x1=l*(sin(B));
  y1=-l*(cos(B));
  
  figure(2)
  plot([-0.5 0.5],[0 0],'color','b')
  hold on
  plot([0,0],[0,0],'o','color','r')
  hold on
  plot([x0 x1],[y0 y1],'color','g')
  hold on
  plot(x1,y1,'o','markers',50,'color','r')
  hold off
  axis([-1.5 1.5 -1.5 1])
  pause(0.02)
  

  % Saving each plot to obtain final animation
textfile=sprintf('animation%05d.png',c);
saveas(gca,textfile);
c=c+1;
end

ERRORS

The following error occurs when the marker size for the bob is not given in proportion with the size of the string used. This error can be rectified by either increasing the size of the string to a proportionate value or by reducing the size of the bob to the prescribed one.

RESULTS

The above graph shows the variation of both angular velocity and angular displacement as a function of time. The legend command used in the code indicates the parameter represented for each color.

The YouTube link for the pendulum animation is given below:

https://youtu.be/P9HDyeyNANs