SOLVING 2nd ORDER DIFFERENTIAL EQUATIONS USING MATLAB

PROJECT NAME – SOLVING 2^nd ORDER DIFFERENTIAL EQUATIONS USING MATLAB

2^nd order differential equation is-

Where,

b = damping coefficient.

m = mass of the body.

g = gravity

l = length

ODE’s are extremely important in engineering, they describe a lot of important phenomenon and solving ODE can actually help us in understanding these systems.

Couple of examples involve in 2^nd order ODE is simple pendulum

So if I solve this particular equation then I am specifically solving this for and also the rate at which particular value of changing with respect to time which is nothing but the velocity in addition to this same pendulum also describes the motion of a spring, mass & damper system so for e.g. we have a spring mass system & if I give this mass initial condition then that same ODE can be used to describe how the mass can move. This has a lot of application in vehicle dynamics where suspension system can be modelled as a set of mass, spring and dampers so we can actually solve the ODE and figure out how vehicle dynamics will be going to be?

The fundamental idea behind solving ODE’s both analytically and computational areas when you have higher order ODE’s then you can break down into set of first order ODE then you have only so what you are going to do is you will say -

Now, original ODE is, So in order to solve 2^nd order ODE one have to create 2 ODE’s And So, this is what it meant as breaking down 2^nc order ODE into a series of 1^st order ODE and that is why we did this transformation so, the idea is if you have third order ODE then break it down into 3 first order ODE’s.

This is called as solution vector because it contains two quantities as it is representes – on LHS there is so now we will solve it using programming language (it can be using any ODE solvers) either be MATLAB/python we need to provide equations in particular format and that is what we are going to looking at next write a programme in MATLAB-

Remember, while solving ODE conditions there are 2 first order ODE’s 1 is for velocity and 1 is for placement and to solve 1^st order ODE we need to provide initial conditions for displacement and we will provide initial conditions of pendulum i.e. 0 and initial velocity is 5 m/s and then we will integrate our ODE system to obtain the solution.

Ode 45 is a inbuilt solver in MATLAB so in general (@ t, theta) saying that it may input as t,theta and there are bunch of other constants but theta and t are both dependent and independent variables of our ODE system so we are saying that our ODE is afunction of (t, theta) and then we are defining our function ade_func(t, theta, b, g, l, m) these are input parameters that we are passing after that we will provide integration and t_span and provide time array to describe how integration needs to be done and then the initial conditions (theta_0) corresponds to start you time. So, ODE function (in a separate script of MATLAB we will write it as-

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

% ODE 45 gives results as 2 columns i- displacement  and ii- velocity by
% being converted into numerical solution.

I have defined theta 1 and theta 2 (theta 1 is our solution array) we are seeing the above derivation so in MATLAB programming combination of theta 1 and theta 2 as theta and now we are defining our first ODE and then second ODE and then we will combine them in an array separate it with 2 rows with the help of semicolon which means that it is going to be a column matrix. This is very important if you don’t use semicolon then you might get into an error and finally this function will give the initial time and the results since we are solving 2 ODE’s the results will be column vector and to be specific it has two columns and if you enter in workspace and open up your results you will see two spaces of the column the first refers to the displacement and second refers to the velocity and you can see first condition starts with zero with boundary conditions and initial conditions we have provided and 2^nd column (t) starts with 5 and now if you run the programme you will see plot –

clear all
close all
clc

b = 0.01;
g = 9.81;
l = 1;
m = 1;

theta_0 = [0; 5];
t_span = linspace(0,10,500);

[t, results] = ode45(@(t, theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);
plot(t, results(:,1))
hold on
plot(t,results(:,2))
xlabel('time')
ylabel('plot')

Fig – Displacement vs.Velocity

Now, to show the oscillation of pendulum

By initializing the inputs we will generate the time array using the linspace command starting with zero and ending with 20 with 150 points in between and ode45 command outputs resulted in theta and theta dot with time and on the right side equation of main variable is wrriten first and with another codes function is made with all the inputs and with the time span and it will integrate upon with the initial conditions first figure is generated by plotting time with all the rows and first column for theta and with the second column for theta dot and all the rows and the second figure is the animation of the pendulum by defining its x0 and y0 and initiate the counter and time array is constructed by using for loop to show the length of the time array and now x1 and y1 is defined by using geometric relation and uses output angular displacement in column 1 and uses counter i for row counter and now by using array brackets pendulum string & wall is plotted and by using symbol – o ball of the pendulum is plotted and defined the axis and in order to have only one plot on the screen after hold on is used to draw all the elements of the plot and the momentum of the pendulum is slowed down with the help of pause command to 0.07and now the frame is captured and counter is increased by 1 and now the movie is made in uncompressed format.

% Second ODE
clear all
close all
clc
b = 0.05      % Damping coefficient
g = 9.81      % Gravity (m/s^2)
l = 1            % Length of string attached to pensulum (meters)
m = 1          % Mass of the pendulum (kg)

% Initial conditions
% Angular displacement = 0 and Angular velocity = 3
theta = [0; 3] 
% Time points
t_span = linspace(0,20,150)
% ODE solver
[t, results] = ode45(@(t, theta) ode_func(t,theta,b,g,l,m),t_span,theta);
% Graph plotting
figure(1)
hold on 
plot(t,results(:,1))
plot(t,results(:,2))
xlabel('time(s.)')
ylabel (' angular displacement and velocity (rad.)')
legend('Angular displacement','Angular velocity')
grid on

figure(2)
hold on
x0 = 0
y0 = 0
ct = 1

% Pendulum oscillation with help of for loop
for i = 1:length(t)
  x1 = l*sin(results(i,1))
  y1 = -l*cos(results(i,1))
   
plot([x0,x1],[y0,y1],'color','b', 'linewidth', 6)
hold on
plot([-1.4,1.4],[0,0],'linewidth',10)
hold on
plot(x1,y1,'o', 'markersize',20, 'markerfacecolor','r')

axis([-1.4 1.4 -1.4 1.4])
grid on
hold off

pause (0.07)
M(ct)= getframe(gcf);
ct = ct + 1;
end

% Creating the animation

movie(M)
videofile = VideoWriter('vib.avi', 'uncompressed avi');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Fig. - Angular displacement vs. Angular velocity

Animation for oscillation of Pendulum -

LINK OF THE ANIMATION - https://youtu.be/oP4WS7LsNbk

Conclusion –

Using ODE45 Higher order ODE can be converted into first order ODE also trignometric as well as geometric calculation can be easily done by knowing the syntax of ODE45.

From the above graph and animation lowering the damping coefficient motion becomes more frequent as well as continuous so to increase the resistence of motion – increase the damping coefficient.