Research Project Chapter 1-3 Sample –NUMERICAL METHODS FOR SOLVING DIFFERENTIAL EQUATIONS

CHAPTER ONE

  1. INTRODUCTION
    • Background of Study

Numerical analysis is a branch of mathematics that studies quantitative methods for solving complex equations using arithmetic operations, often so complex that it requires a computer to estimate the process of analysis (i.e. calculus) (Elzaki & Ezaki, 2011). according to numerical analysis is a valuable tool to generate and evaluate methods for calculating numerical information from a giving data set. Though therer are various types of differential equation which the numerical method can be used to solve, this project, focuses on “fractional fractional equations.” this follows the fact that In recent years, fractional calculus has attracted the attention of more well established scholars and many successful researchers in from various fields of science and engineering (Abell & Braselton 2016). according to Abell & Braselton (2016) one of the major advantages of fractional calculus is that fractional derivatives provide a superior approach to describing the memory and hereditary properties of a variety of materials and processes.

Several numerical methods using different types of derivative operators to solve different types of fractional measures have been proposed. The most commonly used is the Adomian decomposition method (ADM) (El-Kalla, 2011; Ahmed, 1998), General differential modification technique (GDTM) (Shaher, 2007; Zaid, 2008; Shaher, 2008), wave method (Li Zhu, 2012; Mingxu, 2012), difference method to (Deng, 2013), Laplace transformation method (Gupta, 2014), variation iteration method (Odibat, 2006), fractional variation transformation method (Arikoglu, 2009), operational strategy (Luchko, 1999; Li, 2013; Bengochea, 2014) and other methods (Zaid, 2009; Kai, 2004; Javidi, 2009). Moreover, orthogonal functions also play an important role in finding numerical solutions for FDE, such as Pulse Pulse function (Yi, 2013), Bernstein polynomial (Saadatmandi, 2014), replaced by Legendre polynomial (Khalil, 2014), Chebyshev wavelets (Zhu , 2012), Legendre Wave (Heydari, 2014), etc. in like manner, there are also several mathematical tools or integrated development environments (IDEs) (such as MATLAB, Python, Mathematica, and Maple) that provide efficiency, reliability, and simplicity – use the code to solve the differential equation numerically.  Yi  (2013) further noted that the concept of derivatives derived from variable sequences has also been introduced and several related practical application research papers have emerged.

Due to the breakdown of fractional orders in differential operators, FDE analysis solutions are usually difficult to obtain. As a result, various methods have been developed to provide numerical solutions for FDE, such as viscoelasticity and damping, wave propagation and propagation, and turbulence (Zhu, 2012). Regardless of the large number of definitions of unequal fractional derivatives that are widely used, and this study focuses on a specific form (called the Caputo derivative version) (Yao & Chen, 2013). Variant derivatives, which have traditionally been defined and studied by mathematicians, are derivatives of Riemann Liouville fractions but are not always the most appropriate definition for actual applications. When used in mathematical models, fractional Riemann-Liouville derivatives require initial conditions to be expressed in fractional integrals. Yuanlu, (2010) described it as a derivative without explicit physical interpretation, and therefore initial probabilistic values ​​are required to be immediately available. An alternative definition of fractional derivatives obtained after the detection and integration of interchangeable differences is the so -called Caputo derivative, which for sufficient functional variation has only the advantage of requiring initial conditions assigned to the overall derivative conditions of the order. The derivatives of this whole sequence represent the characteristics of the observable physical state, and therefore their value can be measured accurately. thus, in this project, the researchers will only use Caputo fractional displacement, thanks to its application to real models

Solving FDEs has been recognized by researchers as a daunty task due to the numerous technique available (sheun et al., 2019).  According to Yao & Chen (2013), solving FDEs is not an easy task but not saying it is impossible, as majority of the time that they are solved numerically by using some computational tools. The effect of such computational always result the problem of approximating the answer since they are very complex and they require an expertise in computational mathematics to write the code for execution (Yuanlu, 2010). For Elzaki & Ezaki, (2011) solving differential equation to some extent involves the use of different code which are mostly unavailable for a time dependent event and another for non-time dependent.

 

The purpose of this project is to provide illustrative and basic principles of some FDE solutions, so as to be able to do a brief exercise on numerical FDE solutions using some MATLAB procedures and to discuss some non-trivial problems related to effective implementation. methods such as permanent memory treatment, solving equations involving implicit methods, etc; at the same time explaining the MATLAB routines used to solve various FDEs, and concluding by solving some numerical examples.

The significance of this project is to provide and easy-to-use step by step method on how these FDEs can be solved using some routines coded in a MATLAB software. It will move away from solving problems that will not be practically possible to something that is visible and possible. This will cut away the amount of time and resources spent on solving an FDEs that will result in an approximated solution to an event driven solution that will be just good enough for that real-world problem.

This research will explore the standard way of solving FDEs and relate to it by using one of the available frequently used mathematical computational tool to solve some examples that depends on different events. One of the limitations of solving fractional differential equations numerically is that it requires an experienced or expertise in computational mathematics with reasons being that, one faces some non-trivial difficulties, mainly related to the presence of a persistent memory (which makes the computational procedure to be extremely slow and expensive), to the low-order accuracy of the majority of the methods, to the not always straightforward computation of the coefficients of several schemes, and so on. With no particular expertise in computational mathematics, it is will be difficult to code routines that will solve the problems, instead, This research will rely on efficient and already tested routines.

CHAPTER TWO

LITERATURE REVIEW

2.1 Introduction

This section contains literature pertaining to calculus, application of differentiation and the numerical methods of solving Fractional Differential Equations (FDEs). It gives a detailed view of how FDEs are solved numerically.

2.2.1 Functions of Single Variable and their Graphs  

Functions are essential ingredients in the study of calculus. Their absence means an incomplete study of this branch of mathematics which is more or less the care of many aspects of mathematics.

Calculus by its basis is the notion of correspondence for example, the surface area (A) of a sphere relates to its radius by the formula

A= 4πr2

The sphere volume (V) of a given mass of gas is related to the pressure (P) of the gas.

V=

These example give us an idea what a function is. Conclusion from them could imply that if the value of one quantity, say y depends on the value of another quantity, say x, then for every value of x  there corresponds one and only the value one of y.

On this basis we say that y is a function of x. thus, A is a function of r and V is a function of p as above.

Though, the same element of y may correspond to different element of x. for example, two different books may have the same number of pages.

Definition 2.0: A function f from a set x to a set y is a correspondence that assigns to each element x of X a  unique element y of y.

The set of values of X is called the domain of the function. The element y is called the image X under f and is denoted by f (x). the range of the function consists of all images of elements of X.

We can represent these explanations by means of diagram.

Values of f at each of the points indicated ( or the images under f ) cam be obtained by more evaluation

Example2.1   f (x)= x3 +4x-3

F (2)= 23 + 4 (2) -3

=8+8-3

=16-3

=13

F(-3)   =(-3)3+4(-3)-3

=-27-12-3

=-42

F (c) = c3+4 (c) -3

=c3+ 4c-3

Definition 2.1 .A function f is said to be defined at x or f(x) exists if x is in the domain. Put  differently, a function is defined for a certain value b of  x, if a definite value of f (b) corresponds to the value of x.

On the other hand, f is undefined at x if x is not in the domain of f or at a given value of x the function is meaningless, since division by zero is not a valid operation.

Example2.2 : for what value of x is the function y = 5x-5 defined? What is the domain and co- domain?

Solution: y=5x-5 is defined for every value of x since for every value of x, we obtain one value for y.

The domain are the real numbers and the co-domain is the set of all the real numbers.

Y or f depends on the value of x hence the domain is called the DEPENDENT VARIABLE while

X is called the INDEPENDENT VARIABLE.

DEFINITION 2.2  a function f form x to y is one to one function. (injective)

Example2.3 : determine if the function f is one to one if : (i) f (x) =2x+9,  (ii) f(x) =3x2 at x = 3 and 4 for (i) and at x=-2 and 2 for (ii)

Solution: (i) if a≠ b then with a =3 and b = 4 we have 2.3+9 ≠ 2.4+9 or 15≠17 hence f is a one to one function

(ii) if a≠ b in x then with a = -2 and b =2 we have 3 (-2)2 = 3(2)2 or     f (-2) = f(2)

Therefore 12 =12 or f(-2) =f (2)

By this f is not a one to one function

Definition 2.3 :

(a) If f is a function from x to x or if f (x) =x for every x, that is every element x maps into itself then f is an identity function on x. For example

F (x) =x

f(i)=1 etc

(b) A function f is a constant function it there is a fixed element b such that f(x)=b for every x in the domain. for example f (x) = 1  x

2.2.2 Graph of a Function

This aspect of calculus deals with the sketching of graph of a function.

The tread in this topic is to make a table of values for a specified range of x to obtain values f (x). From this we can plot the pair of points by choosing a suitable scale

Example2.4  sketch the graph of f if

Solution: the table below list some points x, f(x) in the graph.

X½1345
F (x)o13

Plotting the points yield

Definition 2.4

Where the coefficient ao, a1, a2,………an are real numbers and the exponents are non negative integer polynomial of degree 3 is given below.

F (x) = ax3 + bx2 + cx+d. a0

F(x) = anxn nR

Where a, b, c and d are constant and R is a set of real numbers.

If n = 0

F (x) = a0

If n = 5

F (x)=a5x5+ a4x4 + a3 x3+ a2 x2+a1x +a0

f (x) =

Where 2x2 + 3x + 1and 5x3 + 6x+7 and polynomials

An algebraic function is given below

                f(x)= 1+4x2 +

However there is one exception where the variables x and y are not given separately. Thus, a function this form is said to be implicitly defined if neither x nor y is expressed as an explicit function of the other.

2.2.3 The Rate of Change of A Function (Gradient Of A Function )

In this section we shall endeavor to investigate and observe what happens, when a function changes along its curve and observe the average rate of change. This average rate of change could equally be considered as the gradient or the slope of that function.

Thus, given a function y=x, we observe that as x is changing y is equally changing.

In this light, the rate of change of a given function could be measured in between two points by taking the average rate of change.

Consequently, it x, and x2 are two points in x and y1 and y2 in y, then the rate of change is the ratio of the increase in y to that of x hence

Rate of change =

Example 2.5 find an expression for the average rate of change of the function

 Solution. By the formula obtained above for the average rate of change of a function, the average rate of y=3x+4 is

for y =2x2 we have

By this expression we can obtain different values at different intervals.

2.3 Differentiation as a limit of rate of change of elementary function As we rightly said in the study of rate of change of a function, the gradient of a line is measured by taking the ratio of the increase in Y and the increase in x in moving from one point to another on the line. Such that if (x1, y1,) and (x2, y2,) are two point on a line, the gradient is then taken as

 

 

 

 

         (x, y)

 

 

Infact the gradient of a line can be seen as the slope of that line. That is, the rising or falling of the line. It is quite easy to find the gradient or slope of a straight line since this will always remain constant all along the line. This is always equal to the tangent of the angle made with the positive direction of the x-axis the gradient of a straight line is always equal to tan

 

 

 

 

 

,

Curve

However, this task is not quite easy in a curve, since the curve changes direction at each point. This difficulty brings us to the idea of a tangent line, which has been a problem of  the great Greek Scientist Archimedes 287-212) B.C. The problem has been that of finding a unique tangent line (if one exist) to a given curve a given point on the curve.

Examples of typical tangent line at different points are shown below.

Diagram

 

 

 

 

 

 

 

 

One might ask, what has tangent line got to do with gradient? Infact, the gradient at a point on curve is the gradient of the tangent line at that point.

As a result of this difficulty, we develop a means of finding a gradient function which will be for a particular curve.

The gradient function is a function is a function from which we can find the gradient at any given point.

The means is in line with the study of limits.

Suppose we want to find the gradient of a function f(x) = x2 at the point 2. when x = 2, y = x2 = 4

We can approach this by finding it act a point N(2, 4) where x=2 i.e the gradiant of the tangent at that point, then we take another pointy R(3, 9) the gradient point of the curve drawing

Next we let R assume a new position closer to N and calculate the gradient then. We will observe that the approximation will get better until R is. Very close to N and it will be clear that the gradient approaches the value 4.

This process of allowing R get close N is known as the limiting process and we say that as R tends to N the gradient of NR approaches the limiting value of 4 hence we say the gradient of the tangent at N is this value 4.

However, a general method of finding the gradient function as follows:

Suppose we want to find the gradient function of Y=x3 at the point N(x, x3). Next take another point R with x-coordinate an increase in x. thus the ordinate of R will have a corresponding increase y and so its y – coordinate y+

Now at N, Y=x3 and at R.

The gradient of NR=

As R→ N, . The limiting value will be 3x2 for meaningfulness, we rewrite the limiting value of  (read dee y/dee x) and thus obtain

This limiting value of  was due to the work of Leibnitz.

We shall at times use f1(x) to denote  . Thus the gradient function of y = f(x) is  is the limiting value of  as x →0.

Note, any letter or symbol can be used for the increment in X and Y as such one should not be surprise when letters or symbols such as  are used.

Also this method of finding the gradient function from which the gradient at a point is obtained is known as the 1st principle of finding the gradient.

Example 3.1:

Find the gradient function of 2x2.

Solution: let y = 2x2

Take a point N coordinate (x, y) i.e. (x, 2x2) on the curve, next take another point R whose X coordinate is then y coordinate of R 2(x+x)2. corresponding increase in y2 is y+y

Gradient of the curve NR=

2.4 differentiation as a limit of rate of change of a function

The procedure for finding  or f1(x) as explained above is called differentiation. Often at times, we call it the derivation of y with respect to x or the differential coefficient. The  measures the way or how often a given function changes. Hence we say it measures the rate of the change of the function y when viewed with x.

 Thus when we say the rate of change of 3x2 is 6x. we mean the gradient of the tangent or the derivative of 3x2 is 6x.

Example 3.2: find the rate of change of y = 4x2 and the gradient of the tangent to the curve y = 4x2 at the point 2.

Solution:

In terms of

Hence the gradient of the tangent at the point is f1(2) = 8.2= 16

Example 3.3: find the derivative of y = 3x2. from the first principle.

Solution: y = 3x2 …………………….. (1)

With  as increments in x and y respectively.

Subtract (i) from (ii)

Example 3.4

Find the derivative of f(x) = sin x

Solution: let y = sinx.

Take as an increment in x and y respectively.

Then:

We now apply the trigonometrical formula:

2.5 Rules of Differentiation

In this section we shall endeavour to look at the various rules each with examples.

Rule 1: The derivative of a constant is zero.

For clearness, we use the idea of limits. Let y = c, we wish to show that f1(x) = 0 since any value of f is C it follows that C + x = C for all x

Hence

Rule 2:

The power Rule.  If n is a +ve integer then derivative of xn is nxn-1

Example 3.5: Find the derivative of 5x2

Solution: let

Example 3.6 find   if y = 7x4

Thus if y = axn

=naxn-1

By this, we have the theorem.

Theorem: the derivation of this sum or difference of two or more functions is equal to the derivative of the individual function and their sums or differences taken afterwards.

Example 3.7

Find f(x)1 if f(x) = 10x2 + 9x-4

Solution

Rule 3: The Product Rule

The product rule for differentiations of x, then the derivative of y with respect to x is done by taking the derivative of s while keeping t constant plus the derivatives of t while s is kept constant.

The t is;

Y = st ………………….(i)

Take  as increment is s and t respectively thus producing a change  in y

Then:

Example 3.8:

Differentiate

g(x) = (x3-7) (2x2+3)

let g(x) = (x3-7) (2x2+3)

When compared with the rule above, hence product rule.

Theorem: If n is positive integer then the derivative of x-n is –nx-n-1

Example 3.9: differentiate

Solution y = (3s)-4

Example 3.1.0  Differentiate

2.6 Review of Related Works

The study by Garrappa (2018) on the Numerical Solution of Fractional Differential Equations was carried out to examine how to solve FDEs analytically and with the aid of a computational tool (MATLAB) to solve equations numerically. The research work also outlines the challenges that is faced while trying to solve FDEs which could be time dependent, lack of built it routines to enable equations be solved faster. One of the major problems it pointed out was that most of the mathematical computational tools used today fail to provide a robust and reliable functions to solve FDEs and in most times, the researcher has to write his/her own code depending on which mathematical computational of his/her choice. Some of the tools that are popularly used in solving FDEs includes but not limited to these alone: Analytica, FreeMat, MATLAB, LabVIEW and TK Solver. A quick google search on any these tools will show the features and what problems they can solve best.

For the sake clarity, some preliminary materials on fractional calculus were also read and properly referenced in Zayernour & Karniadakis. (2014). Some of the most useful definitions in fractional calculus were recalled and their properties were also used in the subsequent section of the work. As the starting point for introducing Riemann-Liouvile (RL) integral; for a function y(t) ∈ L1 ([t0, T]) (as usual, L1 is the set of Lebesgue integral functions), the RL fractional integral of order α > 0 and origin at t0 is defined as (Li, Qian, & Chen, 2011):

y(t) = α-1 y(τ)dτ.          (1)

It provides a generalization of the standard integral, which, indeed, can be considered a particular case of the RL integral (1) when  = 1. The left inverse of is the RL fractional derivative:

y(t) : = y(t) = m-α-1 y(τ)dτ (2)

where m = is the smallest integer greater or equal to  and Dm, y(m) or dm/dtm denotes the standard integer-order derivative.

An alternative definition of the fractional derivative, obtained after interchanging differentiation and integration in Equation (2), is the so-called Caputo derivative, which, for a sufficiently differentiable function, namely for y ∈ Am ([t0, T]) (i.e., y(m-1) absolutely continuous), is given by:

y(t) : = y(t) = m-α-1 y(m)(τ)dτ (3)

They observe that  y(t) is a left inverse of the RL integral, namely

y = y

but not its right inverse, since

(4)

where Tm-1 [y; ](t) is the Taylor polynomial of degree m -1 for the function y(t) centred at , that is

y(k)(

The two definitions (2) and (3) are interrelated, and indeed, by deriving both sides of Equation (4) in the RL sense, it is possible to observe that

 = y(t) – Tm-1

and, consequently:

y(t) =

Observe that in the special case 0 < a < 1, the above relationship becomes:

y(t) =

clearly showing how the Caputo derivative is a sort of regularization of the RL derivative at t0. Another feature that justifies the introduction of the Caputo derivative is related to the differentiation of constant function; indeed, since:

in some applications, it is preferable to deal with operators where the inverting derivative is zero, as in the case of the Caputo derivative.

However, one of the most important applications of Caputo derivatives is in FDE. Unlike FDE derived from RL, which is initialized by a non-integer derivative, the initial value problem for FDE (or FDE system) derived from Kaputo can be formulated as:

)

y() = y0,y’(t0 ) = y0(1),…,y(m-1)( = y0(m-1) (6)

where f (t, y) is assumed to be continuous and

y0(1),…,y(m-1)( = y0(m-1)

 is the derivative value specified at t0. Clearly, the beginning of an FDE with an integer derivative value is more useful because it has a clearer physical meaning with respect to fractional derivatives. The application on both sides of equation (6) of the integral RL ogether with equation (4) leads to a rearrangement of the FDE from the point of view of the weak Volterra integral equation (VIE):

y(t) + α-1f(τ ,y(τ))dτ (7)

The integral formula (7) is undoubtedly useful because it allows the use of theoretical and numerical results already available for this VIE class to study and solve FDE.

They emphasize the non -local nature of FDE: the presence of real forces in the kernel does not make it possible to separate the solution of equation (7) at any point as the solution at some previous point tn – h plus the increment term related to the interval [tn – h, tn], as is common with ODEs (Li, Qian, & Chen, 2011)

CHAPTER THREE

METHODOLOGY

3.0 Introduction

 One of the most important ways of solving or dealing with FDEs is looking at as a delay fractional differential equation because it consists of one or some of its arguments solved at a time which differs by any fixed number or variable called time or lag. Theoretically has it that, using delay fractional differential equation can be considered as a general approach to the theory of FDEs and applying the forementioned Linear Multistep Methods (LMM’s) to evaluate delay fractional differential equations of the retarded and can be used to derive some of the well-known numerical methods to solve delay fractional differential equations. The same technique will be applied here in this section but before that, I will outline the important derivatives that we have in FDEs. For one to be able to solve FDEs, one need the understanding of fractional calculus, how to take fractional derivatives of some basic functions and some integrals. All these concepts add up to bringing out a correct solution to any FDEs, there are many fractional integrals in FDEs, out of which, one is the most outstanding which is the   Riemann-Liouvile integral. There are also varies of fractional derivatives but these are the most famous of them all namely: Riemann-Liouvile derivative and the Caputo derivative. With these three fractional operatives and some other elementary derivatives, one can solve an FDE successfully. There are many methods of solving these equations, one can use the Laplace transform method, Power series method, Fourier Transform method etc. In this case I will stick to the Laplace transform method.

 

3.1 Laplace transform of the function y(t)

For y(t) of any exponential order  then  exist for all

Res > . The Laplace transform of y(t) will be denoted by L{ y(t)} and this is defined as

L{ y(t)} =  (s) =

From the above, I can say that y(t) = { (s)} is the inverse Laplace transform of  (s).

Laplace transform of is given by

 L{} =

More useful pairs of Laplace transform of the function) are

L { )} = (Re(s) >| a |)

 {} = )}

Riemann-Liouville Fractional Derivative and Riemann-Liouville Fractional Integral

Let α ∈ R+ and n = [α], where [.] is the greatest integer function. The Riemann-Liouville Fractional Derivative of order α is denoted by Dα a is defined as

y(t) = m-α-1 y(m)(τ)dτ, for t0 ≤ t ≤ b

Also, the Riemann-Liouville Fractional Integral of order α is denoted by   and is defined as thus:

y(t) = α-1 y(τ)dτ.

Let’s consider these two examples
1.     C =

 C = ) =  ,

Where C is any constant

  1.   C =

 C = ) =  ,

Where C is any constant

3.2 Power Function of Riemann-Liouville Fractional Derivative and

 The Riemann-Liouville Fractional Derivative of order α of a Power function is given as

tv = ,  V > -1,

and the Riemann-Liouville Fractional Integral of order α of a Power function is given by

tv = ,  V > -1,

Let’s consider another example

  1. 0t2= ,  =    =
  2. 0t2= ,  =    =

Properties of Riemann-Liouville Fractional Derivatives

For any function y(t) that is continuous at t ≥ a, then integration of the arbitrary real order has the following properties

  1. [ y(t)] =  y(t)
  2. [ y(t)] = y(t),  for
  3. [(t) + ] = (t) + ]

3.3 Laplace transform of the Riemann-Liouville fractional derivative

The Laplace transform of Riemann-Liouville fractional derivative is of the order α > 0 is given as

L[ y(t); s] = sαY(s) –

Substituting for α = 1,  = 1, k = 1, a = 3 in property 1 of the Riemann-Liouville fractional derivative, it will yield

L{)} =

{} =

Below are the parallels of some Laplace transform

3.4 MATLAB Analysis

Moving over to MATLAB, there are a couple of functions that are prewritten by MATHWORKS: the sole owners of MATLAB software, they have a great community that enables individuals, groups and organizations to collaborate and work together as a team. They also have a host of libraries in order to simplify complex situations or calculations to a simple function. You can access those functions by importing any of the desired or concern library(ies) into your script and call those functions to do what to use them for. Some can have return type while some without. You are free to ask questions in the community and someone with an idea or solution will definitely reach out to you via the community and see how to solve that problem, by so doing they have because one of the best computational used around the world, though you to subscribe in other to get your license yet there is a trial version that you can use to do all that you want to do without any limitation. Students also have an extended period of license. This project was on version R2015a on a computer with the following specifications to enable the software to run its calculations smoothly without lagging. 8GB DDR3 RAM, 3.0GHz processor speed, Running Windows 10 Operating system and 500HDD ROM storage. Below is the overview of the MATLAB work environment.

Screenshot of the MATLAB script for a single euqtion

Below is the set of equations that will be solved with  the MATLAB script

(t) = (t-1)

(t) = (t) – (t-1) + (t – 0.5)

(t) = (t) – (t)|

  t ≤ 0; (t) = 1, (t) = 0, (t)|= -1

There are two lags in this set of equations, 1 and 0.5

(t-1) represent the history of the equation

Below is the MATLAB code to solve the above set of equations

clc;

clear;

close all;

lags = [1 0.5];

tf = 5;

sol = dde23(@ddefunc, lags, @yhist, [0 tf]);

t = linspace(0, tf, 100);

y = deval(sol, t);

figure;

plot(t,y);

function yp = ddefunc(t, y, YL)

    yl1 = YL(:,1); %lag: 1

    yl2 = YL(:,2); %lag: 2

    yp = [yl1(1)

          y(1) – yl1(1) + yl2(2)

          y(2) – y(3)];

end

function y = yhist(t)

    y = [1 0 -1]’;

end

3.5 Explanation of the code in the script above
Each complete statement in MATLAB ends with a semicolon
clc: This is used to clear all the text from the command window and it results in a clear screen
clear: This removes all variables from the workspace and frees up the system memory
close all: It closes all open MATLAB figure windows and clear all data stored to a variable.
lags: this variable is defined to hold the lags or the delay of the equation to be solved
yp: is a the user defined function that is created to handle the set of the equations
figure: this used to bring out a visible presentation of the solution graphically
plot(t,y): this is used to plot the value or the solution with the delay. This will be the solution for the set of equations entered
t,y: this variables are used to parse the delay and solution of the equation
yl1: this is used to handle the first delay in the equation
yl1: it is used to handle the second delay in the equation
yhist: this is used to hold the history of the equation
dde23: this the built-in function in MATLAB, that is designed to solve delay differential equation. It is a very powerful function that will accept the delay of the equation(s), history of the equation and
sol: this is the variable that holds the solution of the solved equation
ddefunc: this represent the model of the equation
tf: this represent the time span of the equation, this can be adjusted to see how the model will be likle for such a time.
The solution to the solved equation with the graphs will be presented in the next chapter for more explanation, adjustment will be made to some variable to see how they perform with time is being adjusted.

REFERENCE

Abell, M. L., & Braselton, J. P. (2016). Differential equations with Mathematica. Academic Press.

Ahmed. M. A. (1998), Nonlinear functional differential equations of arbitrary orders , Nonliear Analysis: Theory, Methods & Applications, Elsevier, vol. 33, no.2, pp. 181-186.

Arikoglu, I. (2009), Solution of fractional integro-differential equations by using fractional differential transform method’’, Chaos Solitons Fractals. 40 521-529.

Bengochea G. (2014), Operational solution of fractional differential equations, Applied Mathematics Letters. 32 48-52.

Deng, W.H. (2013), High order finite difference WENO schemes for fractional differential equations, Applied Mathematics Letters. 26(362-366.

EI-Kalla I. L. (2011),  Error estimate of the series solution to a class of nonlinear fractional differential equations,  Commun. Nonlinear Sci. Numer. Simulate, Elsevier, vol. 16, no.3, pp. 1408-1413.

Elzaki, T. M., & Ezaki, S. M. (2011). On the Elzaki transform and ordinary differential equation with variable coefficients. Advances in Theoretical and Applied mathematics6(1), 41-46.

Garrappa, R. (2018). Numerical solution of fractional differential equations: A survey and a software tutorial. Mathematics6(2), 16.

Gohari, H., Barari, A., & Kishawy, H. (2018). An efficient methodology for slicing NURBS surfaces using multi-step methods. The International Journal of Advanced Manufacturing Technology95(9), 3111-3125.

 Gupta, S. D. (2014), Numerical study for systems of fractional differential equations via Laplace transforms, Journal of the Egyptian Mathematical Society, doi:10.1016/j.joems.2014.04.003.

 Heydari, M.R. (2014), Legendre wavelets method for solving fractional partial differential equations with Dirichlet boundary conditions, Applied Mathematics and Computation. 234 267-276.

Javidi, A. (2009), Modified homotopy perturbation method for solving system of linear Fredholm integral equations, Mathematical and Computer Modelling, Elsevier, vol. 50, no. 1-2, pp. 159-165.

Kai J. (2004), Multi-order fractional differential equations and their numerical solution, Applied Mathematics and Computation, Elsevier, vol. 154, no.3, pp. 621–640.

 Khalil, H. R. A. (2014), A new method based on Legendre polynomials for solutions of the fractional two-dimensional heat conduction equation, Computes & Mathematics with Applications. 671938-1953.

Li Zhu, Qibin F. (2012), Solving fractional nonlinear Fredholm integro-differential equations by the second kind Chebyshev wavelet,  Communications in Nonlinear Science and Numerical Simulatation, Elsevier, vol. 17, no.6, pp. 2333-2341.

Li, C., Qian, D., & Chen, Y. (2011). On Riemann-Liouville and caputo derivatives. Discrete Dynamics in Nature and Society2011.

 Li, W. (2013),  Solving Abel’s type integral equation with Mikusinski’s operator of fractional order, Advances in Mathematics Physics. Vol. 2013, Article ID 806984, 4 pages.

Liu, Y., Roberts, J., & Yan, Y. (2018). Detailed error analysis for a fractional Adams method with graded meshes. Numerical Algorithms78(4), 1195-1216.

Lubich, C. (1986). Discretized fractional calculus. SIAM Journal on Mathematical Analysis17(3), 704-719.

Luchko, R. (1999), An operational method for solving fractional differential equations with the Caputo derivatives, Acta Mathematica Vietnamica. 24(2) 207-233.

Mingxu Y. (2012), Haar wavelet operational matrix method for solving fractional partial differential equations,  Computer Modeling in Engineering &Sciences, Tech Science, vol. 88, no.3, pp. 229-244.

Odibat, S. (2006), Application of variational iteration method to nonlinear differential equations of fractional order, Int. J. Nonlinear Sci. Numer. Simul. 7 27-34.

Rackauckas, C., & Nie, Q. (2017). Differentialequations. jl–a performant and feature-rich ecosystem for solving differential equations in julia. Journal of Open Research Software5(1).

Rezaiee-Pajand, M., & Sarafrazi, S. R. (2010). A mixed and multi-step higher-order implicit time integration family. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science224(10), 2097-2108.

 Saadatmandi A. (2014), Bernstein operational matrix of fractional derivatives and its applications, Applied Mathematical Modeling. 38 1365-1372.

Shaher M. Z (2008), Generalized differential transform method: Application to differential equations of fractional order”, Applied Mathematics and Computation, Elsevier, vol. 197, no.2, pp. 467-477.

Shaher M. Z. (2007),  Generalized differential transform method for solving a space and time-fractional diffusion-wave equation , Physics Letters A, Elsevier, vol.370, no. 5-6, pp. 379-387.

Yao, K., & Chen, X. (2013). A numerical method for solving uncertain differential equations. Journal of Intelligent & Fuzzy Systems25(3), 825-832.

Yi, M.X. (2013), Block pulse operational matrix method for solving fractional partial differential equation, Applied Mathematics and Computation. 221  121-131.

Young, A. (1954). Approximate product-integration. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences224(1159), 552-561.

Yuanlu, L. I. (2010). Solving a nonlinear fractional differential equation using Chebyshev wavelets. Communications in Nonlinear Science and Numerical Simulation15(9), 2284-2292.

Zaid. O. (2009), An algorithm for the numerical solution of differential equations of fractional order, J. Appl. Math. Inform, Elsevier, vol. 26, no. 1-2, pp. 15–27.

Zaid. O. M (2008),  A generalized differential transform method for linear partial differential equations of fractional order , Applied Mathematics Letters, Elsevier, vol. 21, no. 2, pp. 194-199.

Zayernouri, M., & Karniadakis, G. E. (2014). Exponentially accurate spectral and spectral element methods for fractional ODEs. Journal of Computational Physics257, 460-480.

 Zhu, Q.B. (2012), Solving fractional nonlinear Fredholm integro-differential equations by the second kind Chebyshev wavelet, Commun Nonlinear Sci Numer Simulat. 17 2333-2341.