Solving 2nd order ODE of a simple pendulum using MATLAB

INTRODUCTION

A Simple Pendulum is a dynamic system. A mass-spring-damper system (MSD) can be considered as a discretised model of any dynamic system. So a simple pendulum can describe the motion of a spring, mass and damper system. If the mass is given an initial condition then the ODE of the motion of pendulum can be used to describe how this mass is going to move.

The mass of the dynamic system is lumped into a single point mass in the MSD system. The inertial effect of the dynamic system is related through this lumped mass.

The stored energy/internal energy of the dynamic system is modelled as a one-dimensional spring in the MSD system. The spring is able to store energy inside when it is stretched or compressed from its original length.

The energy dissipated out of the dynamic system is modelled through a one-dimensional damper in the MSD system. The viscous damper, for instance, is able to dissipate energy as heat outside the dynamic system.

These three components, mass, spring and the damper can model any dynamic response situation in a general sense.

In this project the 2^nd order ODE of a simple pendulum is solved.

Mathematically the linear dynamic equation of equilibrium for a pendulum with dampling is represented as:

(d^2 θ)/dt+  b/m  dθ/dt+  g/l  sinθ=0

To solve this ODE computationally, the higher order ODE must be broken down into a set of first order ODEs.

Let θ = θ₁

dθ₁/dt = θ₂

d²θ/dt² = dθ₂/dt

With these, the series of first order ODEs are

(dθ_1)/dt= θ_2       &       (dθ_2)/dt=  (-b)/m θ_2-g/l sinθ_1

So the second order ODE is broken down into a set of two first order ODEs and now a solution can be obtained.

PROGRAMMING IN MATLAB

clear all
close all
clc
 
b = 0.05; %The damping coefficient
g = 9.81; %Acceleration due to gravity
l = 1; %Length of the string
m = 1; %Mass of the ball
 
% Initial Conditions(Angular displacement & Angular velocity)
theta_0 = [0 3];
 
% Time points(0 to 20 seconds)
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);
 
figure(1)
plot(t, results(:,1))
hold on
plot(t, results(:,2))
xlabel(\'time\')
ylabel(\'Plot\')
 
theta = results(:,1);
 
ct = 1;
 
for i=1:length(theta)
    Theta = theta(i);
    
    x0 = 0;
    y0 = 0;
    
    x1 = l*sin(Theta);
    y1 = -l*cos(Theta);    
    
    figure(2)
    %Plotting the string of the pendulum
    plot([x0 x1],[y0 y1],\'linewidth\',3)
    hold on
    
    %Plotting the ball of the pendulum
    plot(x1,y1,\'o\',\'MarkerSize\',15,\'MarkerFaceColor\',\'r\') 
    hold off 
    
    axis([-1 1 -1.5 0])
    grid on
    pause(0.001)
    M(ct) = getframe(gcf);
    ct = ct + 1;
end  
%Creating a movie
movie(M);
videofile = VideoWriter(\'pendulum.avi\', \'Uncompressed avi\');
open(videofile)
writeVideo(videofile,M)
close(videofile)

STEP WISE EXPLANATION

STEP 1: The input are defined values with the following commands.

b = 0.05; %The damping coefficient
g = 9.81; %Acceleration due to gravity
l = 1; %Length of the string
m = 1; %Mass of the ball

STEP 2: This command defines the initial condition angular displacement (0) and angular velocity (3) of the pendulum

theta_0 = [0 3];

STEP 3: The time array is next defined as,

t_span = linspace(0,20,500);

The motion is stimulated between 0 to 20 seconds with 500 values in between.

STEP 4: Solving the ODE

The command is as follows:

\'[t, results] = ode45(@(t, theta) ode_func(t,theta, b, g, l, m), t_span, theta_0);‘

The output matrix gives two results, the time‘t’ and the ‘displacement and velocity’ in ‘results’ which is column vector which holds the results of the two first order ODEs.

Explanation:

• ode45 is an inbuilt solver in MATLAB
• (@(t, theta)says our ODE is a function of ‘t’ and ‘theta’
• ode_func(t,theta, b, g, l, m)is a function which takes the input values t, theta, b, g, l and m.
• t_span, theta_0 are the initial inputs which are to be integrated in the ODE

The function ‘ode_func(t,theta, b, g, l, m)’

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

Explanation: At first the function is assigned a name: ‘ode_func’ which take the inputs, ‘t’, ’theta’, ‘b’, ‘g’, ‘l’, ‘m’ and gives the output ‘dtheta_dt’. ‘theta1’ and ‘theta2’ take values of ‘theta(1)’ and ‘theta(2)’.

The first order ODEs and the defined,

dtheta1_dt = theta2;

dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1);

These two are later incorporated in an array called,

dtheta_dt = [dtheta1_dt; dtheta2_dt];

This function is then called in the main program to solve the two first order ODEs.

STEP 5: Once the ODE is solved it can be plotted. This plot includes two curves of ‘velocity and displacement’ vs ‘time’.

figure(1)
plot(t, results(:,1))
hold on
plot(t, results(:,2))
xlabel(\'time\')
ylabel(\'Plot\')

STEP 6: Plotting the pendulum motion and creating an animation.

• A new variable ‘theta’ is defined which includes the results array.

theta = results(:,1);

• ‘for’ loop is used where a new variable ‘i’ takes values from 1 to length of ‘theta’. Each value that ‘theta(i)’ takes is assigned to ‘Theta’

for i=1:length(theta)
    Theta = theta(i);  

• The following coordinates are defined, ‘x0’ and ‘y0’ are assigned a value zero and ‘x1’ and ‘y1’ are calculated.

x0 = 0;
y0 = 0;
   
x1 = l*sin(Theta);
y1 = -l*cos(Theta);​

· The following commands create a new figure to plot lines between the points ‘x0’,’x1’,’y0’,’y1’and a bob using the ‘marker’ command. The axes of the plot are defined and grid is set. The ‘pause’ command determines the pause time between displaying two plots.

figure(2)
%Plotting the string of the pendulum
plot([x0 x1],[y0 y1],\'linewidth\',3)
hold on
    
%Plotting the ball of the pendulum
plot(x1,y1,\'o\',\'MarkerSize\',15,\'MarkerFaceColor\',\'r\') 
hold off 
    
axis([-1 1 -1.5 0])
grid on
pause(0.001)

· The ‘getframe’ command is used to create an array of frames in a variable ‘M’ and these frames are incremented in steps on ‘1’ using a constant ‘ct’.

M(ct) = getframe(gcf);
ct = ct + 1;

STEP 7: The frames are then used by a command ‘movie’ to create a movie. Further Inbuilt commands are used to save, open, write and close the movie.

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

Here is the video link to the plot obtained: https://youtu.be/sb-VI7EQRcY

RESULTS

The second order ODE for the motion of a simple pendulum is solved giving plots of ‘angular displacement and angular velocity’ vs ‘time’ and an animation is created.