SIMULATION OF MOTION OF A PENDULUM

• AIM: In Engineering, ODE is used to describe the transient behavior of a system. A simple example is a pendulum

The way the pendulum moves depends on Newton's second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "ball" w.r.t time.

Similarly, there are hundreds of ODEs that describe key phenomena and if we have the ability to solve these equations then we will be able to simulate any of these systems.

• GOVERNING EQUATION:

The `2^(nd)` order ODE for Pendulum is given by -

`(d^2theta)/dt^2 + b/m (d theta)/dt + g/l sintheta` ........................(I)

In the above equation,

g = gravity in m/`s^2` , l = length of the pendulum in m, m = mass of the ball in kg, b=damping coefficient.

• OBJECTIVE:

The objective is to write a program that solves the ODE. This ODE represents the equation of motion of a simple pendulum with damping.

To write a Program in Octave, following Input Values are considered to simulate the Pendulum Motion-

`l` = 1 m, `m` = 0.1 kg, `b` = 0.05, `g` = 9.81 m/`s^2`

Simulation of the motion is considered between 0-20 sec, for angular displacement=0, angular velocity=3 rad/sec at time t=0.

Program for this simulation is -

% PROGRAM TO SOLVE 2ND ORDER ODE
clear all
close all
clc

  % Inputs
   b = 0.05 ;  % damping factor
   g = 9.81 ; % Gravity 
   l = 1 ; % length of Pendulum
   m = 0.1 ; % mass
   
   % initial condition
   theta_0 = [0;3] ;  % [Initial displacement ; initial velocity]
   
   % time span
   t_span = linspace(0,20,500);
   
   % solving ODE
   [t, results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);
   

 subplot (5,4,[13,14,17,18])
   plot(t,results(:,1))
   hold on
   plot(t,results(:,2))
   legend('Displacement','Velocity')
   xlabel('Time')
   ylabel('Plot')
   
  
  subplot(5,4,[15,16,19,20])  
    plot(results(:,1),results(:,2))
    xlabel('Position-Radians')
    ylabel('Anguler 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));
   subplot(5,4,[1,2,3,4,5,6,7,8,9,10,11,12])  
     plot([-1 1],[0 0],'k','lineWidth',4)
     hold on;
     line([x0 x1],[y0 y1],'color','b','lineWidth',2) 
     hold on
     plot(x1,y1,'-o','markersize',20,'markerfacecolor','r')
     grid on ;
     axis([-1.5 1.5 -1.5 0.5]) ;   
     pause(0.0003)
     hold off;
     file_text = sprintf("Pendulum%d.png",ct)
     saveas (gca,file_text);
     ct = ct+1;
     
  
   
   
   end

To run this program, a function is written. The function is -

% Function to solve ODE

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] ;
  
  endfunction

This function is used in ode45. This command is used to solve the ODE equation. It solves a set of non-stiff Ordinary Differential Equations (non-stiff ODEs) with the well known explicit Dormand-Prince method of order 4. The function must accept two inputs where the first is time t and the second is a column vector of unknowns. In this project, the second column vectors are displacement and velocity of the bob of Pendulum.

Octave can display more than one plot in a single figure. The simplest way to do this is to use the subplot function to divide the plot area into a series of subplot windows that are indexed by an integer. In this Project, we divide the plot area into thee subplot windows, one is for - the motion of Pendulum, the second one is for - the plot of velocity and displacement with time & third one is for angular velocity with respect to angular position.

After successfully run the program, we get the plots and the position of the pendulum as a function of time.

During this project, some error has been faced but in the end, it is successfully overcome. The 1st ERROR raised during giving the range of 'i' which is a variable. The range of this value had given mistakenly from 0 to length of the vector, instead of 1 to the length of [results].

`2^(nd)` Error is related to the order of definition of inputs to the function. As ode45 is a order of time and theta, so function also writes according to this format.

After running the program, Command Prompt was opened to convert all the png files related to this project to make a gif file with the help of Magick Image software. After getting the gif file further it is converted to mp4 file to upload it to YouTube.

• CONCLUSION:

In this Project, we solve a `2^(nd)` order ODE by using ode45. From the output plots, we can observe the position and velocity of Pendulum, as we introduce damping, so after some time span, the pendulum stops. If we increase the value of damping then it will be stopped early. By varying the inputs, we can get different plots. In this project ode45, subplot functions are used. At the end Magick image software is used to stitch all the png files to a gif file.