Simulation of Pendulum Motion by using Octave

AIM: TO WRITE A PROGRAM TO SOLVE SECOND ORDER ODE AND TO SIMULATE THE PENDULUM MOTION BY USING OCTAVE.

OBJECTIVE:

1. To write a program to solve second order (ODE) and to simulate the pendulum motion
2. To understand transient behavior of pendulum
3. To make an animation

THEORETICAL BACKGROUND:

an ordinary differential equation (ODE) is a differntial equation containing one or more functions of one independent variable and the derivatives of those functions, In Engineering ODE are used to describe transient behavior of system example – Newton second law, simple pendulum

To solve second order ordinary differential equations ( ODE) we need to convert second order ODE into two first order ODE ,suppose if we have to solve third order ODE then we need to convert that into three first order ODE.

Let’s say given equation below is second order differential equation

`(d^2 θ)/(dt^2 )+b/m (dθ)/dt+g/L sin⁡θ=0`

To convert this equation into two first order ODE

Let’s say

`(dθ)/dt=(dθ_1)/dt=θ_2`

`(d^2 θ)/(dt^2 )=(d^2 θ_1)/(dt^2 )=d/dt ((dθ_1)/dt)=d/dt (θ_2 )`

So by substituting

`(dθ_2)/dt+b/m θ_2+g/L sin⁡θ=0`

`(dθ_2)/dt=-b/m θ_2+g/L sin⁡θ`

`(dθ_1)/dt=θ_2`

So we can represent those two first order ODE in matrix format or we can say solution vector

``d/dx [θ_1/θ_2 ]= [θ_2/(-b/m θ_2+g/L sin⁡θ )]`

In the above equation,

g = gravity in [ m/s²] ,

L = length of the pendulum in [m],

m = mass of the ball in [kg],

b = damping coefficient.

INTIAL INPUT:

Length of the pendulum = 1 [meter]

Mass of the ball = 1 [kg]

Damping coefficient = 0.01

Acceleration due to gravity = 9.81 [m/s²]

MAIN PROGRAM :

Programming language used: High level Programming Language using Octave software

clear all
close all
clc

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

% initial condition
theta = [0;3];
% time points
t = linspace (0,20,500);
% to solve ode
[results] = lsode('ode_fun',theta,t)

figure(1);
plot(t,results(:,1),'color','b');
hold on;
grid on;
plot(t,results(:,2),'color','r');
xlabel('Time');
ylabel('Plot');
legend('Displacement','velocity');
 ct=1;
 for i=1:length(results(:,1))
 x0=0;
 y0=0;
 x1= l*sin(results(i,1));
 y1= -l*cos(results(i,1));

 figure(2)
 plot([-1,1],[0,0],'linewidth',6,'color','g');
 axis([-1.5 1.5 -1.5 0.5]);
 line([x0,x1],[y0,y1],'linewidth',2,'color','b');
 hold on
 plot(x1,y1,'o','markersize',10,'markerfacecolor','r','markeredgecolor','r');
 hold off;
 pause(0.05)
 
 text=sprintf("animation_pendulum%05d.png",ct);
 saveas(2,text)
 ct=ct+1;
 end

EXPLANATION:

1. First, we have closed all files, cleared all workspace, cleared screen.
2. Define constant that we going to use in program for example length of pendulum, mass of ball, damping coefficient and acceleration due to gravity.
3. Then we have defined the initial condition for displacement and velocity.
4. We have created array of time span by using “linspace” command for which we are going to see the variation in displacement and velocity.
5. In octave we have used “lsode” function to solve two first order ordinary differential equation
6. Then we plot variation of angular displacement with time span and variation of velocity with time span.
7. While plotting we have used “hold on” and “grid on” command to see variation of displacement and velocity in figure window.
8. We have labeled the plot by using command “ xlabel ” ,” ylabel ” .
9. To create the animation, we have defined counter as 1.
10. We have used for loop which will run for 500 times because length of “length(results(:,1) ”
11. We know that in animation one point of pendulum will be fixed so that we have defined coordinate of that point as (0,0).
12. We have found out coordinates of second point which will varry and define as x1,y1.
13. Then we have made one line by using line command to which fixed point will be pivoted.
14. After that we have defined how will our axis is going to varry.
15. Then we plot the line which joins fixed point and (x1, y1).
16. We have used hold on command to see the line and sphere that we are going to make in animation.
17. We have made sphere by using marker command and center of the sphere (x1, y1) as reference and given diameter and size and color as input.
18. We have pause code to see variation in displacement and velocity.
19. We saved all files to get animation using “saveas” command.
20. We have done increment in counter for next png file and ended the for loop.

SUPPLEMENTARY PROGRAM :

Programming language used: High level Programming Language using Octave software

function [dtheta_dt] = ode_fun(theta,t)
b=0.1;
g=9.81;
l=1;
m= 0.1;

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

EXPLANATION:

1. We have made function to convert second order ordinary differential equation to two first order differential equation.
2. We are going to give theta and time as input argument to function “ode_fun”.
3. We have defined constant that we going to use in function for example length of pendulum, mass of ball, damping coefficient and acceleration due to gravity.
4. And used the process of converting secondary order differential equation into two first order ordinary differential equation as explained above in theoretical background.
5. And ended the loop and this function is used in “lsode” function which solves this two first order ordinary differential equation.

COMMAND FOR ANIMATION:

EXPLANATION FOR ANIMATION IN OCTAVE:

1. run main program we will get 500 png files .
2. open command prompt, you will see working directory.
3. to change directory use cd followed by directory of all png files.
4. once your directory changed, try checking files in that location by using dir and hit enter.
5. now use magick “*.png output.gif” command to get animation.

OUTPUT:

1. plot of time verses angular displacement and velocity.

2. Animation link:

https://drive.google.com/open?id=1QpfCJ5kZ7RPn3xSyCY8aRsYCqgE4lWoX

CONCLUSION:

We have understood the variation of angular displacement and velocity with respect to time span, we have written code for solving second order differential equation .

RESULT: Simulation of pendulum motion has been done successfully.