Simulate Simple Harmonic Motion of A Pendulum By Solving Second Order Differential Equation (Newtons 2nd Law) Using Matlab/Octave Program.

Abstract:- The work is focused on transforming Newtons 2nd Law of motion, which is in the form of a second-order differential equation into two ordinary differential equations (ODE) & writing a code using "ODE45" or "LSODE" function in MATLAB or Octave (depending on the programing language used) respectively to simulate the harmonic motion of a pendulum. This ODE represents the equation of motion of a simple pendulum with damping.

In the above equation,

                                                                    g = gravity in m/s2,  L = length of the pendulum in m,  m = mass of the ball in kg,  b= damping coefficient.

                                                                    g= 9.81 m/s2,  L= 1 metre,  m= 1 kg, b=0.05.  

                                                                    Simulate the motion between 0-20 sec, for angular displacement=0, angular velocity=3 rad/sec2 at time t=0.

Method:- The above-given equation is in the form of a second-order differential equation & in order to solve the equation in MATLAB or Octave we need to convert this equation into an ordinary differential equation (ODE). MATLAB allows the solution of first-order differential equations by using a special function called "ODE45" or "LSODE" in the case of older versions of Octave GUI. These special functions use the Runge-Kutta method with a variable time step for efficient computation. But we are given a second-order differential equation that represents the motion of the simple pendulum. To use these special functions we need to find a way to convert the second-order differential equation into two first-order differential equations.  

We have, `(d^2theta)/(dt^2)+b/m(d theta)/(dt)+g/Lsin theta=0`

Now consider, `theta = theta1`  &  `(d theta)/(dt) = (d theta1)/dt = theta2`

We can write `(d^2 theta)/dt` as `d/dt (d theta)/dt => d/dt theta2`

Substituting this in the equation we get,  `(d theta2)/dt + b/m theta2 + g/L sintheta1=0`

Further, we can rearrange,   `(d theta2)/dt =- b/m theta2 -g/L sintheta1`

where `theta1 =` angular displacement     &     `theta2 =` angular velocity

Now that we have sucesfully transformed the second order differential equation into first order differential equation we will need the geneal command for the function ODE45, which is given by 

                                                                                                                      [t,results] = ode45 (@fname, time span, initial condition)

Function to solve the first-order ODE

function [dtheta_dt] = ode_function_pendulum(t, theta, b, m, L, g)
                 
                 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

Main Code

clear all
close all
clc

%input parameters

  m = 1; %mass of the bob
  b = 0.05; %damping coefficient
  g = 9.81; %acceleration due to gravity
  L = 1; %length of the wire
  ct = 1; %counter variable

%initial condition

  theta_0 = [0;3]; 

  
%duration for oscillation

  t_span = linspace(0,20,500);
  
%Solve the differential eqation by calling the function

  [t, results] = ode45(@(t, theta) ode_function_pendulum(t, theta, b, m, L, g), t_span, theta_0);
  
%plot the Variation of angular displacement & angular velocity w.r.t time

  figure (1)
  hold on
  plot (t,results(:,1))
  plot (t,results(:,2))
  xlabel ('x-axix')
  ylabel ('y-axix')
  grid on
  title ('Variation of angular displacement & angular velocity w.r.t time', 'color', 'b')
 
 theta2 = results(:,1);
 
 for i = 1:length(theta2);
   theta_x = theta2(i);
   
   x0 = 0;
   y0 = 0;
   
   x1 = L*sin(theta_x);
   y1 = -L*cos(theta_x);
 
%plot the pendulum motion

 figure (2)
 plot ([-0.5 0.5], [0 0], 'color', 'b', 'linewidth', 5);
 hold on
 plot ([x0 x1], [y0 y1], 'color', 'r', 'linewidth', 2)
 hold on
 plot (x1,y1, 'o', 'markers', 15, 'markerfacecolor', 'y')
 hold off
 grid on
 axis ([-1.5 1.5 -1.5 0.5]);
 pause (0.002)
 M(ct) = getframe (gcf);
 ct = ct+1;
 
 
 end

 movie (M)
 videofile = VideoWriter ('pendulum.avi', 'Uncompressed avi');
 open (videofile)
 writeVideo (videofile, M)
 close (videofile)
 

Observation:- 

Reference 

https://in.mathworks.com/help/symbolic/simulate-physics-pendulum-swing.html