MATLAB program to solve ODE and simulate motion of simple pendulum

Ordinary Differential Equation:

In mathematical sense, 'Ordinary Differential Equation' or 'ODE' is a differential equation containing one or more functions of one independent variable and derivatives of those functions.

Second order linear differential equations have a variety of applications in science and engineering. For eg. vibration of spring, electric circuit, etc.

In this project I am writing MATLAB program to solve ODE which represents the equation of motion of a simple pendulum with damping.

Mathematical Part:

Hence we have converted one 2nd order differential equation into two 1st order differential equations which can be solved using ode45 function of MATLAB.

MATLAB Code:

clear all
close all
clc
 
%Inputs
b = 0.05;        %Damping coefficient in 'N-s/m'
g = 9.81;        %Acceleration due to gravity in 'm/s^2'
l = 1;           % Length of pendulum in 'm'
m = 1;           % Mass attached to the string in 'kg'


%Initial conditions
theta_x = [0;3]; % Column matrix for initial conditions where initial_displacement = 0 m and initial_velocity = 3 m/s

%time span
time_span = linspace(0,20,500);   % Timespan from 0-20sec and dividing this span into 500 equal parts

% Solving ODE using MATLAB function ode45
[t,results] = ode45(@(t,theta) second_order_ode_func(t,theta,b,g,l,m), time_span, theta_x);

% plotting

figure(1)
plot(t,results(:,1),'color','g');     %plotting displacement vs time graph as we are only plotting first column of results matrix which is displacement
hold on
plot(t,results(:,2),'color','r');     %plotting velocity vs time graph as we are only plotting second column of results matrix which is velocity
xlabel('Time (sec)');
ylabel('Displacement(green) (m), Velocity(red) (m/s)');
theta_position = results(:,1);    % saving different position values in a separate array

ct = 1;
for i = 1:length(time_span)       % for loop to plot different positions of pendulum
    x0 = 0;                       % x co-ordinate of pivot point
    y0 = 1.5;                     % y co-ordinate of pivot point
    
    x1 = x0 - l*sin(theta_position(i));   % x co-ordinate of bob end , x = x0 - l*sin(theta)
    y1 = y0 - l*cos(theta_position(i));   % y co-ordinate of bob end , y = y0 - l*cos(theta)
    
    figure(2)
    plot([-1,1],[1.5,1.5],'color','black','linewidth',4); % Making top support
    hold on                           % to hold the previos plot so that plot don't get overwrite
    
    % plotting string by joining pivot point and instantaneous point of bob
    % bob is plotted using marker of size 15
    plot([x0,x1],[y0,y1],[x1,x1],[y1,y1],'o','MarkerSize',15,'MarkerFace','y','linewidth',3); 
    hold off                      % to plot one position of pendulum rod at a time otherwise all position gets plotted on same graph

    axis([-1,1,0,2]);     % setting x-axis limit as '-1 to +1' and y-axis limit as '0 to +2'
    pause(0.01);          % hold each plot for 0.01 sec

    N(ct) = getframe(gcf); % saving each frame in an array named as N
    ct = ct + 1; 

end

movie(N)                   % creating video file of animation
videofile = VideoWriter('Simple Pendulum','Uncompressed AVI');
open(videofile);           % opening video file
writeVideo(videofile,N);   % writing each frame in video file
close(videofile);          % close video file

Function to convert 2nd order ODE into two 1st order ODEs:

function [dtheta_upon_dt] = second_order_ode_func(t,theta,b,g,l,m)
    
    theta_1 = theta(1);
    theta_2 = theta(2);
    dtheta1_upon_dt = theta_2;
    dtheta2_upon_dt = -((b/m)*theta_2)-((g/l)*sin(theta_1));
    dtheta_upon_dt = [dtheta1_upon_dt ; dtheta2_upon_dt];
    
end

In function part I had used the same logic which I had shown in mathematical part above.

Logic to create animation :

To create animation I had used getframe command to save each figure showing a particular position of pendulum in an array. I had also used a counter 'ctr' to keep track of no. of frames. Then I had combined all these frames to make an animation video.

OUTPUT:

Following plot shows the variation of displacement of pendulum about mean position with time in green color.

Also variation of velocity with time can be seen in red color.

Animation video can be seen below:

By observing the plot and animation we can conclude that amplitude of vibration decreases with time. This is because of presence of damper in the system.

Errors faced while programming:

1)Only bob is coming on screen without rod:

I had plotted a particular position of pendulum but only bob is visible on the plot as shown below

I had written this code

I was using two separate plot commands for rod and bob.

To rectify this problem I had used a single plot command with two sets of co-ordinate data, and it worked perfectly.

2) All positions of pendulum gets plotted on the same plot:

After running the code I am getting simulation but all positions of pendulum are getting plotted on the same plot.

After some time I had realised that because of hold on command which is written above plot command all positions are getting printed.

To rectify this problem I had added hold off command after plot command.