SOLUTION OF SECOND ORDER ODE FOR SIMPLE PENDULUM

Second order differential equations are typically harder than first order. In most cases, students are only exposed to second order linear differential equations. A general form for a second-order linear differential equation is given by

`a(x)y\'\'(x) +b(x)y\'(x) +c(x)y(x) = f(x)`

One can rewrite this equation using operator terminology. Namely, one first defines the differential operator `L = a(x)D2 + b(x)D + c(x)`

The Function code in order to get recalled in the main program -->

Creating the function \'ode_function(~,theta,b,g,l,m)\'

• In a new script, a function is opened with return value being \'dtheta_dt\' as shown

function [dtheta_dt] = ode_function(~,theta,b,g,l,m)

end

• The name of the function is \'ode_function\', and the input values are ~,theta,b,g,l,m.
• In the body of the function,θ2\"> `theta2` and `(d theta)/dt`
  
  dθ2dt\">
  
  
  
  
  
  were calculated and are assigned as elements of the matrix.
• This matrix is returned to the main program.

Parameters considered are -

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

function [dtheta_dt] = ode_function(~,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];

end

The Main MATLAB code -->

• Initialize the parameters of the coefficients like damping factor, acceleration due to gravity, length of the pendulum and mass of the bob. And dependent variables like time interval, initial values of position and initial velocity.

• Call the function \'ode45\' and assign the output to two sets of values. In calling the function the main dependent variables are specified. And the function \'ode_function(~,theta,b,g,l,m)\', time interval and the initial position and velocity are given as inputs. The function \'ode_function(~,theta,b,g,l,m)\' returns a matrix.

• From the obtained values, a plot of displacement and velocity vs time is made which shows the rate of change of the quantities. The curve can be varied by changing the damping coefficient.

• Now, a counter is assigned to 1. A \'for\' loop is opened to plot every position of the pendulum for the specified time interval. A pointer \'i\' is assigned to the displacement column of the output matrix. So, whenever the \'for\' loop increments the value of the pointer, the successive displacement value of the pendulum is obtained.

• Using trigonometry, the exact position of the pendulum bob is calculated and stored in temporary variables as coordinates.

• A plot of [X, Y] is done and a marker is assigned as per requirements to represent a bob of the pendulum. This plot is stored in an array of the counter, and then the counter is incremented and the loop is closed.

• Outside the loop, a movie of the acquired plots is made by writing the following code.

clear all
close all
clc

b = 0.05; % Damping factor
g = 9.81; % Acceleration due to gravity
l = 1; % Length of the pendulum
m = 1; % Mass of the pendulum


% Initial position in radians and velocity
theta_0 = [0;3];

t_span = linspace(0,30,1000);

[t,y] = ode45(@(t,theta) ode_function(t,theta,b,g,l,m),t_span,theta_0);

figure(1),clf
plot(t,y(:,1));
hold on
plot(t,y(:,2));
xlabel(\'Time\');
ylabel(\'Parameter\');
legend(\'Angular Displacement\',\'Angular Velocity\');

ct = 1;
for i = 1:length(t_span);
    
    % position of the pendulum bob
    X(i) = cos(y(i,1)-pi/2); 
    Y(i) = sin(y(i,1)-pi/2);
    Xo = cos(y(1,1)-pi/2);
    Yo = l-sin(y(1,1)-pi/2);
    
  figure(2)
  clf   %Forgot to mention Initially hence error in the animation took place
     plot([Xo X(i)],[ 0 Y(i)],\'-\');
     hold on;
     plot(X(i),Y(i),\'marker\',\'o\',\'markersize\',15,\'markeredgecolor\',\'b\',\'markerfacecolor\',\'b\');
     axis([-2*l 2*l -2*l l]);
     set(gca,\'xtick\',[],\'ytick\',[])
     axis square
    
     M(ct) = getframe(gcf);
     ct = ct+1;
end

movie(M)
videofile = VideoWriter(\'Damped Harmonic Pendulum.avi\',\'Uncompressed AVI\');
open(videofile);
writeVideo(videofile,M);
close(videofile);

The output of the Program is given below -->

FIGURE 1 --

FIGURE 2 (ANIMATION) --

The Simulation of whole Matlab coding is shown below -->

There was a single error while writing the Code, which included the addition of \'\'cfl\'\' after the figure() code line.

Hence, as a result of this, the following result in the animation was obtained -->

Therefore, when the \'\'cfl\'\' was added to the line 33 in the main code, we obtained a refined and more sensible animation of the simple pendulum.

As shown,

Hence, the solution and the animation of the simple pendulum was obtained.