Motion of pendulum using lsode

AIM:

To write a program in MATLAB/Octave to solve second order ordinary differential equation that represents the motion of the simple pendulum.

GOVERNING EQUATION:

Motion of a simple pendulum with damping is represented by following second order ordinary differential equation :

`(d^2 θ)/(dt^2 ) +b/m dθ/dt +g/l sin⁡θ`

b= Damping coefficient (kg/s)

m= mass of bob (kg)

g= acceleration due to gravity (m/`s^2`)

l= length of the string (m)

OBJECTIVES:

1. To write a program in MATLAB/Octave to solve second order ordinary differential equation that represents the motion of the simple pendulum.
2. To simulate the motion of the simple pendulum and to create its animation.

THEORY:

A simple pendulum consists of relatively massive object hung by a string from a fixed support. It typically hangs vertically in its equilibrium position. The massive object is affectionately referred to as the pendulum bob. When the bob is displaced from equilibrium and then released, it begins its back and forth vibration about its fixed equilibrium position. The motion is regular and repeating, an example of periodic motion.

Motion of a simple pendulum with damping is represented by following second order ordinary differential equation :

`(d^2 θ)/(dt^2 ) +b/m dθ/dt +g/l sin⁡θ ...............eq.1`

In order to solve this second order ordinary differential equation in MATLAB/Octave, It is converted into two first order differential equation, the procedure for which is described below:

Let `θ= θ_1`

Therefore,

`dθ/dt = (dθ_1)/dt = θ_2` (Assume) ........................eq.2

Now,

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

Replacing the values from eq.2 & eq.3 in primary differential eq. (i.e. eq.1), we get:

`(dθ_2)/dt + b/m θ_2 +g/l sin⁡θ_1`

Rearranging it, we get:

`(dθ_2)/dt = - b/m θ_2 - g/l sin⁡θ_1`...................eq.4

Therefore the primary second order differential equation gets converted into two first order differential equations:

1. `(dθ_1)/dt = θ_2`
2. `(dθ_2)/dt = - b/m θ_2 - g/l sin⁡θ_1`

The above two equations can be represented as:

`d/dt [[θ_1], [θ_2 ]] = d[θ] =[[θ_2] , [-(bθ_2)/m- (g sin⁡θ_1 )/l]]`

Displacement values `(θ_1)` and angular velocities `(θ_2)` of the pendulum can be obtained by solving above equations.

PROGRAM: To write a program in MATLAB/Octave to solve second order ordinary differential equation that represents the motion of the simple pendulum.

PROGRAMMING LANGUAGE: Octave

%solving ODE for pendulum

clear all
close all
clc

%input parameters ( damping coefficient, acceleration due to gravity,length of 
%pendulum, mass of bob respectively)
b=0.05;
g=9.81;
l=1;
m=0.1;

%initial condition
% Providing the initial values for displacement and angular velocity
%displacement = 0 ; angular velocity = 5
theta_0=[0;5];

%time_span
time_span=linspace(0,10,500);

%Solving ode
%calling ode_func
[results]=lsode('ode_func',theta_0,time_span);

%plotting the variatin of angular displacement and angular velocity wrt time
figure(1)

%plotting the variation of displacement of the pendulum with respect to time
plot(time_span, results(:,1),'color','b');
hold on
grid on;

%plotting the variation of angular velocity with respect  to time
plot(time_span, results(:,2),'color','r');
hold off

%labelling of Axes
ylabel('angular velocity and angular displacement');
xlabel('Time')

%placing a legend on graph
legend('angular displacement' , 'angular velocity');

%defining the variable for frame counting and initialized the value 1
ct=1;

%Dispalcement= θ_1= an array containing first coloumn and all rows of 'results'.
displacement=results(:,1);

%using 'for' loop in order to iterate different positions of pendulum
for i= 1:length(displacement)
THETA=displacement(i)

%Plotting the Pendulum
%length of pendulum
l=1;

%Defining the Fixed and Free end of the pendulum
x0=0 ; y0= 0 ;
x1= l*sin(THETA);
y1= -l*cos(THETA);

%plotting
figure(3)

%plotting of horizontal line containing the pivot point of pendulum
plot([-0.5 0.5],[0 0],'linewidth',5,'color','b');
hold on;

%plotting the string of pendulum
plot([x0 x1],[y0 y1],'linewidth', 3,'color','g');
hold on;

%plotting the bob of pendulum
plot(x1,y1,'marker','.','markersize',50,'markerfacecolor','r');
hold off

%axis generation
axis([-1.5 1.5 -1.5 0]);
%Turning on the grid
grid on;

%Freezing the frames for every iteration 
pause(1)

%change the file name with every iteration
 file_text=sprintf(“ode_for_pendulum%d.png”,ct);
 
%saving the plot for every iteration
 saveas(gca,file_text);
  ct=ct+1;
 
 end

CALLING FUNCTION:

function [dtheta_dt]=ode_func(theta,t)

%input parameters
b=0.05;
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

The calling function ‘ode_func’ is defined by using following equations:

1. `(dθ_1)/dt = θ_2`
2. `(dθ_2)/dt = - b/m θ_2 - g/l sin⁡θ_1`

The above two equations can be represented as:

`d/dt [[θ_1], [θ_2 ]] = d[θ] =[[θ_2] , [-(bθ_2)/m- (g sin⁡θ_1 )/l]]`

Displacement values `(θ_1)`and angular velocities `(θ_2)`of the pendulum can be obtained by solving above equations.

These equations are already discussed in detail in the theory part of this project.

PROCEDURE:

Input parameters which include damping coefficient, acceleration due to gravity, length of pendulum, and mass of bob of pendulum are provided.

%solving ODE for pendulum

clear all
close all
clc

%input parameters ( damping coefficient, acceleration due to gravity,length of 
%pendulum, mass of bob respectively)
b=0.05;
g=9.81;
l=1;
m=0.1;

Initial condition for displacement and angular velocity is set to be 0 and 5 respectively and time span of the motion is defined as array which contains 500 values starting from 0 and ends on 10.

%initial condition
% Providing the initial values for displacement and angular velocity
%displacement = 0 ; angular velocity = 5
theta_0=[0;5];

%time_span
time_span=linspace(0,10,500);

Equation is solved by using ‘lsode’ solver and results (i,e angular displacements and angular velocities) are obtained and variation of angular displacement and angular velocity with respect to time is plotted on a single graph.

%Solving ode
%calling ode_func
[results]=lsode('ode_func',theta_0,time_span);

%plotting the variatin of angular displacement and angular velocity wrt time
figure(1)

%plotting the variation of displacement of the pendulum with respect to time
plot(time_span, results(:,1),'color','b');
hold on
grid on;

%plotting the variation of angular velocity with respect  to time
plot(time_span, results(:,2),'color','r');
hold off

%labelling of Axes
ylabel('angular velocity and angular displacement');
xlabel('Time')

%placing a legend on graph
legend('angular displacement' , 'angular velocity');

‘for’ loop is defined in order to iterate different position of pendulum and also a variable is defined for frame counting and its value has been initialized as 1.

%defining the variable for frame counting and initialized the value 1
ct=1;

%Displacement= θ_1= an array containing first column and all rows of 'results'.
displacement=results(:,1);

%using 'for' loop in order to iterate different positions of pendulum
for i= 1:length(displacement)
THETA=displacement(i)

Plotting of pendulum by consecutively plotting its parts which includes horizontal line containing the pivot point for pendulum, length of the pendulum, and bob of the pendulum.

%Plotting the Pendulum
%length of pendulum
l=1;

%Defining the Fixed and Free end of the pendulum
x0=0 ; y0= 0 ;
x1= l*sin(THETA);
y1= -l*cos(THETA);

%plotting
figure(3)

%plotting of horizontal line containing the pivot point of pendulum
plot([-0.5 0.5],[0 0],'linewidth',5,'color','b');
hold on;

%plotting the string of pendulum
plot([x0 x1],[y0 y1],'linewidth', 3,'color','g');
hold on;

%plotting the bob of pendulum
plot(x1,y1,'marker','.','markersize',50,'markerfacecolor','r');
hold off

%axis generation
axis([-1.5 1.5 -1.5 0]);
%Turning on the grid
grid on;

Pause between frames is set to be 1 sec and then ‘sprintf’ command is used for changing the file name with every iteration and ‘saveas’ command for saving the plot for every iteration

%Freezing the frames for every iteration 
pause(1)

%change the file name with every iteration
 file_text=sprintf(“ode_for_pendulum%d.png”,ct);
 
%saving the plot for every iteration
 saveas(gca,file_text);
  ct=ct+1;
 
 end

All the output images (png.) are then converted into gif type and then clubbed together into an animation file by using ‘image magick’ software.

ERRORS ENCOUNTERED:

ERROR 1

This error occurred because of providing wrong syntax for ‘lsode’. According to the syntax, initial conditions (theta_0) should be provided before the time duration (time_span) and can be resolved by doing the mentioned adjustment.

ERROR 2

Parse error encountered while using ‘axis’ command which is resolved by providing suitable brackets.

OUTPUT GRAPH :

The graph shows the variation of angular displacement and angular velocity with respect to time. It can be clearly seen that at time equals to zero, angular displacement is zero and angular velocity is five radian per second which also justifies our initial condition. Also because of damping, the magnitude of angular displacement and angular velocity decreases with time.

ANIMATION:

The link of the animation showing the motion of pendulum is attached below:

https://www.youtube.com/watch?v=QNv61Go5G40

CONCLUSION:

• Second order ordinary differential equation that represents a motion of simple pendulum is solved using the program written in Octave.
• Graph showing the variation of angular displacement and angular velocity with respect to time is plotted.
• Simulation of motion of simple pendulum is done through the animation file.