Solution of Second Order Ordinary Differential Equation using MATLAB

AIM: To solve Second Order Ordinary Differential Equation using MATLAB.

OBJECTIVES:

1. To understand and describe Ordinary Differential Equation.
2. To develop a programme to solve Second Order Ordinary Differential Equation using MATLAB.
3. To simulate motion of Pendulum by solving Second Order Ordinary Differential Equation using MATLAB.

GOVERNING EQUATIONS FOR THE PROJECT UNDERTAKEN:

ODE are very important in engineering they describe a lot of phenomenon and solving those ODEs can actually help us understand these systems. Examples that involves second order ODE includes simple pendulum.

Second order differential equation

`("d" ^"2" "θ" )/("d" "t" ^"2" ) "+" "b" /"m" "" "dθ" /"dt" "+" "g" /"L" "" "sin" ⁡"θ" "= " 0`

Where,

b = damping coefficient

m = mass of the ball [Kg]

g = acceleration due to gravity [m/sec^2]

L = length of pendulum [m]

If we solve this particular equation for theta and also rate at which theta values changing with respect to time which is nothing but the velocity. In addition to this the same pendulum can also describe the motion of spring mass damper system.

The idea of solving ODE both analytically and computationally is when you have higher order ODE than break it out into set of first order ODEs. So the first thing you will notice here is this is actually a second order ODE because it has exponent 2.

Let,

`"θ= " "θ" _"1"`

`"dθ" /"dt" "=" ("d" "θ" _"1" )/"dt" "= " "θ" _"2" ……………………………………………..(1)`

`("d" ^"2" "θ" )/("d" "t" ^"2" ) "=" ("d" ^"2" "θ" _"1" )/("d" "t" ^"2" ) "=" "d" /"dt" "" (("d" "θ" _"1" )/"dt" )"=" ("d" "θ" _"2" )/"dt" " "`

`("d" ^"2" "θ" )/("d" "t" ^"2" ) "+" "b" /"m" "" "dθ" /"dt" "+ " "g" /"L" "" "sin" ⁡"θ" "=0 "……………………………….(2)`

Put

`"θ=" "θ" _"1" in equation (2)`

`("d" ^"2" "θ" _"1" )/("d" "t" ^"2" ) "+" "b" /"m" "" ("d" "θ" _"1" )/"dt" "+" "g" /"L" "" "sin"⁡"θ" _"1" "=0"`

`("d" "θ" _"2" )/"dt" "+" "b" /"m" """θ" _"2" "+" "g" /"L" "" "sin" ⁡"θ" _"1" "=0"`

`("d" "θ" _"2" )/"dt" " = -" "b" /"m" "" "θ" _"2" "-" "g" /"L" "" "sin" ⁡"θ" _"1" ……………………………….(3)`

So we have two single order DE (1 & 3) that is nothing but breaking down a second order ODE into two single order ODEs. Similarly, if we have a third order ODE you can break it out into three single order ODEs.

`("d" "θ" _"1" )/"dt" "=" "θ" _"2"`

`("d" "θ" _"2" )/"dt" " = -" "b" /"m" "" "θ" _"2" "-" "g" /"L" "" "sin" ⁡"θ" _"1"`

Writing above two equations in Matrix form

There are two first order ODE one for displacement and other for velocity

• When you solve first ODE you need to give initial condition for displacement so here we are saying that the initial displacement of the pendulum is ‘0’ and initial velocity is 3 rad/s.
• Then we are going to be integrating the ODE system to obtain solution and for which we need to define the value for time so I’m saying that t_span is my time array goes from 0 to 20 seconds with 500 values in between.

MAIN PROGRAMME

Software used: MATLAB R2020a trial version

Following is the program used to simulate the motion of pendulum between 0-20 sec, for angular displacement = 0, angular velocity=3 rad/sec at time t=0.

Program includes comments which explains commands used in program and various other parameters.

% To Solve Ordinary Differential Equation of Second Order.

clear all
close all
clc

% inputs
b = 0.05;                           % damping coefficient
g = 9.81;                           % mass of the ball              (kg)
l = 1;                              % acceleration due to gravity   (m/sec^2)
m = 0.1;                            % length of pendulum            (m)


% initial conditions                % initial condition corresponding to starting time
theta_0 = [0;3];                    % angular velocity=3 rad/sec at time t=0.                  

% time points
t_span = linspace(0,20,300);        % time array which says how far integration needs to be done


% solve ODE
[t, results] = ode45(@(t,theta) ode_func(t, theta, b,g,l,m), t_span, theta_0);
% ode45      : It is an inbuilt solver in MATLAB
% @(t,theta) : Main functions of ODE
% ode_func   : separate MATLAB program 

figure(1)

%subplot(5,4,[13,14,17,18]);
plot(t,results(:,1),'color','b')                    % result is a two 500 row & two column array (:,1) represents all rwo only first column
hold on                                             % to avoid overwrite of plot
plot(t,results(:,2),'color','r')                    % result is a two 500 row & two column array (:,1) represents all rwo only second column
ylabel('Amplitude')                                 % Y-axis title
xlabel('Time (s)')                                  % X-axis title
legend('Displacement','Velocity')
grid on
hold off

figure(2)
% subplot(5,4,[15,16,19,20])
plot(results(:,1),results(:,2))
xlabel('Position')                                 % X-axis title
ylabel('Angular Velocity')                         % Y-axis title


ct=1;
for i=1:length(results(:,1))
    x0 = 0;                                         %to locate pivot point
    y0 = 0;
    x1 = l*sin(results(i,1));                       % coordinates of bob trajectory
    y1 = - l*cos(results(i,1));
    figure(3)
    %subplot(5,4,[1,2,3,4,5,6,7,8,9,10,11,12]);
    plot([-1,1],[0,0],'linewidth',4,'color','k');   % line on which pivot is loacted
    hold on
    line([x0,x1],[y0,y1],'linewidth',2,'color','b');
    hold on
    plot(x1,y1,'-o','markers',20,'markerfacecolor','r'); % ploting thr bob by using marker & facecolor
    grid on
    axis([-1.5 1.5 -1.5 0.5])
    pause(0.03);                                    % To make animation smoother using pause of (0.3) which time required
                                                    % to represent next change
    hold off

    q(ct) = getframe(gcf);                          % gcf get current figure window
                                                    % after getting individual frame we are going to save them in an array
                                                    %first saving the frames than icnreasing the counter
    
    ct = ct+1;                                      % ct reresents each and every frame
    
end
% at the end of for loop to generate the movie using movie(M) with all those frmae
movie(q)

% we need to create an avi file where we are going to store all the frame
% we need to provide file name and format
videofile = VideoWriter('pendulum3.avi','Uncompressed AVI');

%once we created a file we need to open it by using following command
open(videofile)

% to write the fames which we have captured
writeVideo(videofile,q)
close(videofile)

Programme for function used in main programme

function [ dtheta_dt ] = ode_func(t,theta,b,g,l,m)

theta1 = theta(1);
theta2 = theta(2);
dtheta1_dt = theta2;                                        % first ODE
dtheta2_dt = (-((b/m)*theta2) - ((g/l) * sin(theta1)));     % second ODE
dtheta_dt = [dtheta1_dt; dtheta2_dt];

end

Output: Amplitude vs Time

Output 2:

Output 3: Motion of Pendulum using SUBPLOT COMMAND in MATLAB

Video of Pendulum Motion