Motion of a Pendulum Simulated in Matlab- Second Order ODE

Motion of a Pendulum Simulated in Matlab- Second Order ODE

Obective

The objective is to write a program in matlab that solves a second order ordinary differential equation, which represents the equation of motion of a simple pendulum with damping. And to simulate in an animation form, the damped motion of a pendulum for 20 seconds.

The equation for a damped pendulum motion is:

where,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

In matlab, direct second order ordinary differential equations cannot be computed, therefore it needs to be reduced into two separate first order ordinary differential equations. Which may then be coded using an inbuilt matlab function called \'ode45\'.

Two separate equations were found to be

Equation 1:

Equation 2:

The above equations were used in the following code using matlab.

The initial conditions were taken to be 0 rad/s for angular velocity and 3 rad/s for angular displacement.

Other variables were taken to be

L=1 metre,

m=1 kg,

b=0.05,

g=9.81 m/s^2.

clear all
close all
clc

%time span
t_span=linspace(0,20,500);
%initial conditions
theta0=[0 3];

%solving first order derivative using ode45
[t,theta]=ode45(\'odefun_second_order\',t_span,theta0);

% plot of time vs angular displacement and angular velocity 
figure(1)
plot(t,theta)
axis ([0 20 -5 5])
hold on
xlabel(\'Time(s)\')
ylabel(\'Angular Velocity and Displacement(rad/s)\')
legend(\'Angular Velocity\',\'Angular Displacement\')
title(\'Angular displacement and velocity vs Time\')

% Pendulum Animation 
figure(2)
O=[0 0];
axis(gca,\'equal\');
axis([-1.5 1.5 -1.5 1]);
title(\'Pendulum Motion Animation- A R Shaikh\')

for i=1:length(t_span)
    L=1;
    % Motion equation
    P=L*[sin(theta(i,1)) -L*cos(theta(i,1))];
    %Pendulum line
    pendulum_line=line([O(1) P(1)],[O(2) P(2)]);
    %Point at which line is hinged
    hinge_point=viscircles(O,0.01);
    %Ball that is suspended
    ball=viscircles(P,0.05);
    
    pause(0.0001);
      if i<length(t_span)
          delete (pendulum_line);
          delete (hinge_point);
          delete (ball);
      end

end

File to represent the ode45 function in above code is given below.

Please note the name of the below code file should be \'odefun_second_order\' and it is required to correctly run the above code file.

function [thetadot]=odefun_second_order(t,theta);
L=1;
m=1;
b=0.05;
g=9.81;
thetadot=zeros(2,1);
thetadot(1)=theta(2);
thetadot(2)=(-b/m)*theta(2)-(g/L)*sin(theta(1));
end

Results

The plot (figure 1-Angular displacement and velocity vs time) shows damping of both angular displacement and angular velocity with time. Although, 20s does not seem to be long enough for the pendulum to come to become stationary. For the pendulum to stop completely, either of the two things should happen. The simulation time should either be increased or damping coefficient should have been higher. If the damping coefficient was 0.5 instead of 0.05, the result would nearly display complete halt of the pendulum within the same time period as seen in the figure below.

The animation generated by matlab has been uploaded on youtube and can be found on the link given below.

https://youtu.be/pvavvkbxHmk

It can be concluded that the pendulum has not becoeme stationary within the twenty second time that was computed. By looking at the above plot generated for damping coefficient b=0.5, it can be inferred that the pendulum in the animation would become stationary (vertically positioned in our case as shown in the figure below) had the damping coefficient was 0.5.