Pendulum solver using Octave

Introduction:

Pendulum is the most basic system used to describe vibration systems. Vibrations can be classified into Free and Forced vibrations. Mass, Spring constant and Damping form the three key elements to define the free vibration system. This could be conveniently represented by means of an ordinary differential equation as shown below.

Where b referes to damping coefficient, m refers to mass and (mg/l) forms the spring constant.

Similar to this, engineering problems can be easily converted into system of ordinary differential equations and its solution present wide array of application.

ODE45 is the solver that is widely used across the scientific computing tools to solve these ode equations. The ode45 expects the function which forms a differential matrix for the given problem.

Syntax for ode45: [time, output] = ode45(fcn, time array,init condition)

A second order ODE can be reduced to system of two first order ODE. In general, nth order ODE can be reduced to system of n first order ODE, as shown below. This information is coded in the 'fcn' which is fed to ode45 solver.

In the code, after intial cleanup and after passing input values, intial conditions for the solver is defined and time array is constructed. ODE matrix information is fed into the funtion "ode_func". The results are then obtained and plotted.

Animation

In Octave, video file cannot be obtained directly from coding, however an array of image file can be obtained.

'ct' is used as the counter in the code, each frame is extracted by utilising the For loop. counter is advanced before ending the for loop. 'saveas' and 'sprintf' commands are used in writing the image file in the directory. External means can then be employed to stitch the image files into video or a movie. Plotting command is effectively utilised to create the plot, marker size is changed to resemble that of a ball mass at the end of the string. Similar change in plot properties results in formation of a vertical wall.

Octave code:

%This is coded by Sweekar, contact [email protected] for queries
%This is an intellectual property of Sweekar
%Ordianry differential equations define many aspects of engineering and solving ODEs lead to answers to some interesting problems
%This program is coded to solve ODE that defines motion of an damped oscillating pendulum
%However, with the inputs used, output is extremely underdamped system.
%Zeta vaue or ratio of damping to critical damping coefficient is 7E-03

%Cleanup
clear all
close all
clc

%inputs
c = 0.05; %Damping coefficient
g = 9.81;
l = 1; %length of string
m = 1; %mass of ball

%initial condition
theta_0 = [0;3];

%Time array for integration
t_span = linspace(0,20,500);

%ODE45 SOlVER
[t,results] = ode45(@(t,theta)ode_func(t,theta,c,g,l,m),t_span,theta_0);

%Plotting results
subplot(2,2,3)
plot(t,results(:,1),'color','b')
xlabel('Time(s)')
hold on
plot(t,results(:,2),'color','r')
hold on
title('Angular displacement and Angular velocity vs Time')
legend('Angular displacement','Angular velocity')
grid on
subplot(2,2,4)
plot(results(:,1),results(:,2))
xlabel('Angular displacement(rad)')
ylabel('Angular velocity(rads-1)')
title('Angular displacement vs Angular velocity')
grid on
hold on


%Coding for the animation of simple pendulum motion
theta = results(:,1);
ct = 1;
for i=1:length(theta)
  thetaj = results(i,1);
  x1 = l*sin(thetaj);
  y1 = -l*cos(thetaj);
  subplot(2,2,[1,2])
  plot([0 x1],[0 y1],'linewidth',3,x1,y1,'.','markersize',40,[-1.2 1.2],[0 0],'linewidth',5,'color','black')
  axis([-1.5 1.5 -1.5 0.5])
  title('Simple pendulum')
  xlabel('x axis')
  ylabel('y axis')
  pause(0.03)
  Anime1 = sprintf("pendulum%d.png",ct)
  saveas(1,Anime1)
  ct = ct+1;
end

function [dtheta_dt] = ode_func(t,theta,c,g,l,m)
  
  theta1 = theta(1);
  theta2 = theta(2);
  dtheta1_dt = theta2;
  dtheta2_dt = -(c/m)*theta2-(g/l)*sin(theta1);
  dtheta_dt = [dtheta1_dt;dtheta2_dt];
 end

Results:

Animated video can be found in this link: https://youtu.be/6BGcrnqQf94

Discussion and Report:

1) System represents extremely underdamped system, zeta = 7e-03

2) System of nth order ODE can be reduced to n first order ode and could be solved using ode45

3) Simple harmonic motion of a pendulum can be represented in form of a circle

4) When the damping is added. free vibration of a pendulum can be represented by means of a system of spirals as shown in Angular displacement vs angular velocity

5) Angular velocity vs anguar displacement spirals inwards owing to damping which is equivalent to reduction of radius of circle which is due to dissipation of energy from the system.

6) Animated video can be found in this link https://youtu.be/6BGcrnqQf94