1
CHAPTER I
INTRODUCTION
Differential equation
A Differential Equation is an equation with a function and one or more of its derivatives.
A differential equation in mathematics is an equation that connects one or more unknown
functions and their derivatives[1]. In Applications, a Differential equation typically defines a
relationship between functions to describe physical quantities and derivatives to indicate the rates
at which those quantities change. Due to the prevalence of these relationships, differential
equations are widely used in many fields, including engineering, physics, economics, and
biology. The solutions that satisfy the differential equations and the solutions' characteristics are
the primary goals of studying differential equations [2]. In applied mathematics, engineering, and
physics, differential equations are very significant and helpful, and a lot of mathematical and
numerical equipment has been developed to solve differential equations[3]. We are in a time of
incredible advancement. Engineers can build robots, physicists can explain how waves,
pendulums, and chaotic systems move, and humans can communicate wirelessly across a
massive global network[4]. However, differential equations are the deep and enigmatically
potent forces underlie these modern marvels [4].
They are remarkably adept at predicting our environment. In addition to solving problems
with radioactive decay, continuous compound interest, flow, cooling, and heating, orthogonal
trajectories, fluid mechanics, circuit design, heat transfer, population or conservation biology,
2
seismic waves, and in the field of medicine where differential equations are used. Fractional
differential equation is also used to describe the exponential growth and decay, species
population growth, changes in investment return over time, and bank interest[4]. They are
employed in chemistry to simulate chemical reactions and compute radioactive half-lives. They
describe in economics to identify the best investment plans. This also describes physics’ motion
of waves, pendulums, and chaotic systems. As they relate to the temperature of objects and their
environment, they are also utilized in physics in conjunction with Newton's Second Law of
Motion and the Law of Cooling. These equations are employed in engineering to describe the
motions of electricity. Differential equations are also used in software development to
comprehend how computer hardware relates to electrical engineering or applied physics. They
are also utilized in gaming features to simulate character velocity. They are a crucial component
of computer vision and graphics models because they are fundamental tools for defining the
physical universe's nature [4].
In calculus, which was developed by Newton and Leibniz, is when differential equations
were first introduced. Chapter 2 of Methodus fluxional et Serierum Infinitarumas was published
in 1671[5]. In this book, three distinct differential equation types were listed by Isaac Newton,
𝑑𝑦
𝑑𝑥=𝑓(𝑥)
𝑑𝑦
𝑑𝑥=𝑓(𝑥,𝑦)
𝑥1𝜕𝑦
𝜕𝑥1+𝑥2𝜕𝑦
𝜕𝑥2=𝑦
(1)
The variables are x, 𝑥1, and 𝑥2. The function is 𝑓. It addresses the non-uniqueness of answers
while using infinite series to solve these and other examples.
3
There are different types of differential equations. Distinct differential equations are often
employed, including ordinary, partial, linear, nonlinear, homogeneous, and heterogeneous
equations[6].
Ordinary differential equations
An equation with ordinary derivatives is an ordinary differential equation. A differential
equation with variables also includes a derivative of the dependent variable with respect to the
independent variable[7]. A differential equation having one or more functions of one
independent variable and their derivatives is known as an ordinary differential equation (ODE) in
mathematics[8]–[11]. Ordinary differential equations are utilized as opposed to partial
differential equations, which may be with respect to multiple independent variables[10], [11].
Some of the ordinary differential equations are listed below,
𝑑𝑦
𝑑𝑥=sin(𝑥)
𝑑2𝑦
𝑑𝑥2+𝑘2𝑦=0
𝑑2𝑦
𝑑𝑡2+𝑑2𝑥
𝑑𝑡2=𝑥
(2)
Partial differential equations
An equation with two or more independent variables, an unknown function that depends
on those variables, and partial derivatives of the unknown function with respect to the
independent variables is referred to as a partial differential equation (PDE). The largest
derivative included determines the partial differential equation's order. To formally formulate
and assist in solving physical and other issues involving functions of several variables, such as
the propagation of heat or sound, fluid flow, elasticity, electrostatics, electrodynamics, etc.,
4
partial differential equations are used[12], [13]. Listed below are a few partial differential
equations[14]:
𝑢𝑥=𝜕𝑢
𝜕𝑥
𝑢𝑥𝑥 =𝜕2𝑢
𝜕𝑥2
𝑢𝑥𝑦 =𝜕2𝑢
𝜕𝑦𝜕𝑥=𝜕
𝜕𝑦(𝜕𝑢
𝜕𝑥)
(3)
Where 𝑢(𝑥,𝑦) is a function of two variables. 𝑥 and 𝑦 are independent variables[14].
Linear differential equation
A differential equation described by a linear polynomial of the unknown function and its
derivatives is called a linear differential equation[15]. The following equation is the ideal form of
a linear differential equation[15]:
𝑎0(𝑥)𝑦+𝑎1(𝑥)𝑦′+𝑎2(𝑥)𝑦′′………+𝑎𝑛𝑦(𝑛)=𝑏(𝑥)
(4)
Where, 𝑎0(𝑥)……𝑎𝑛(𝑥) and 𝑏(𝑥) are arbitrary linear differentiable functions, and the
consecutive derivatives of an unidentified function of the variable x are 𝑦′………𝑦(𝑛)[15].
When the linear differential equation is plotted on the graph, it shows a straight line.
Differentiable functions are used as coefficients in a linear combination of fundamental
differential operators to create a linear differential operator. Thus, a linear operator has the
following form in the univariate situation[16]:
𝐿=𝑎0(𝑥)+𝑎1(𝑥)𝑑
𝑑𝑥+𝑎2(𝑥)𝑑2
𝑑𝑥2+⋯……+𝑎𝑛𝑑𝑛
𝑑𝑥𝑛
(5)
5
Where L is a linear differential operator, and 𝑎0(𝑥)……𝑎𝑛(𝑥) are differentiable functions. The
nonnegative integer n also serves as the operator's order[15], [16].
Nonlinear differential equation
A nonlinear differential equation is a differential equation that is not linear in the
unknown function and its derivatives. Nonlinear equations have terms with dependent variables
with indexes higher than one and contain multiples of their derivatives, or these terms can have a
maximum degree of two or more. When the nonlinear differential equation is plotted on the
graph, it forms a curved line on the graph[17]. There are a minimal number of techniques to
solve the exact solutions for nonlinear differential equations, and those that exist often rely on
the equation's symmetry. Nonlinear differential equations can display chaotic behaviors over
very long-time scales. Even the most basic questions regarding the extendibility, uniqueness, and
existence of solutions to nonlinear differential equations, as well as the well-suitable of initial
and boundary value problems for nonlinear partial differential equations. These are challenging
issues, and their solution in certain circumstances is regarded as a significant advancement in
mathematics. However, one would anticipate that the differential equation would have a solution
if it were an adequately constructed representation of an essential physical process.[6], [18].
Linear differential equations typically approximate nonlinear equations. These estimates are only
reliable in certain circumstances. For modest amplitude oscillations, the harmonic oscillator
equation, for instance, approximates the nonlinear pendulum equation[6], [18]. Some of the
Nonlinear differential equations are listed below[17]:
𝑦𝑑𝑦
𝑑𝑥+𝑦2=0,
𝑑sin(𝑦)
𝑑𝑥 +𝑥𝑦𝑑𝑦
𝑑𝑥=𝑠𝑖𝑛𝑥,
(6)
6
(𝑑2𝑦
𝑑𝑥2)3+𝑙𝑛𝑥
𝑥𝑦+1
𝑦2𝑑𝑦
𝑑𝑥=𝑐𝑜𝑠𝑦
Basis set
A basis is a group of vectors that creates every component of the vector space and is
composed of linearly independent vectors[19]. If every member of vector V can be expressed
distinctly as a finite linear combination of elements of B, then the set B of vectors in vector space
V is referred to as a basis in mathematics. The components or coordinates of the vector with
respect to the set of B components are known as this linear combination's coefficients[20], [21].
The term "basis vectors" is used to describe basis components. In the same way that every
member of the vector V is a linear combination of elements of B, a set of B components is a
basis set if and only if each of its elements is linearly independent. A basis is, in other words, a
linearly independent spanning set. Multiple bases can exist in a vector space, but each basis has a
fixed number of elements, known as the vector space's dimension[20], [21].
Polynomial and Polynomial equation
One of the fundamental ideas of algebra is a polynomial equation[22]. The sum of a finite
number of terms, each term is the result of a constant coefficient and one or more variables
raised to a positive integer exponent, can be used to represent a polynomial or an expression
made up of variables and coefficients and involving the operations of addition, subtraction,
multiplication, and non-negative integer exponentiation of variables solely is referred to as a
polynomial[23]–[25]. An equation with a polynomial set to zero is called a polynomial
equation[26]. The following is a list of several polynomial equation equations[24]:
𝑎𝑛𝑥𝑛+𝑎𝑛−1𝑥𝑛−1+⋯…………+𝑎2𝑥2+𝑎1𝑥+𝑎0=∑𝑎𝑘𝑥𝑘=0,
𝑛
𝑘=0
𝑥3+2𝑥𝑦𝑧2−𝑦𝑧+1=0.
(7)
7
CHAPTER II
THEORETICAL BACKGROUND
There are several methods to solve fractional differential equations. But we choose to use
a modified Bernstein polynomial with the help of the Galerkin method[27]. B-poly can be
efficiently coded into any computational and symbolic language like Mathematica, MATLAB,
etc. It is also easy to code and easy to implement the initial condition. The B-poly method does
not use a grid between interval points; it’s a grid-less calculation process. Our calculations show
that the modified B-poly method leads to more correct results for all nonlinear, linear, and
multidimensional fractional differential equations or (Integer order or factional
differential equation (IOFDE)) [28], [29].
Bernstein polynomial
A Bernstein polynomial is a polynomial created by linearly combining Bernstein basis
polynomials, a branch of mathematics known as numerical analysis. Sergei
Natanovich Bernstein is honored by the idea's name[30], [31]. The universal formula for
the nth-degree Bernstein-polynomials is [32],
𝐵𝑖,𝑛(𝑥)=(𝑛𝑖)𝑥𝑖(𝑅−𝑥)𝑛−𝑖
𝑅𝑛,0≤𝑖≤𝑛
(8)
For 𝑖=0,1,………,𝑛, where (𝑛𝑖) are the binomial coefficients, which are given by [32],
8
(𝑛𝑖)= 𝑛!
𝑖!(𝑛−𝑖)!,
(9)
Polynomials are defined to form a full basis across the interval [0, R], and R is the broadest range
over which this is possible. For instance, 𝐵2,5(𝑥)=(5
2)𝑥2(1−𝑥)3=10𝑥2(1−𝑥)3. There
exist nth-degree polynomials of order (n+1). We conveniently set, 𝐵𝑖,𝑛(𝑥)=0, for 𝑖<0 or 𝑖>
𝑛. A straightforward Mathematica or Maple program may be utilized to produce any non-zero
polynomial of any degree of nth support throughout the interval. Generally, the boundary
conditions of the topic under inquiry are connected to the first and final polynomials[32].
The B-polynomials during this period may also be produced using a recursive definition,
allowing for the writing of the i-th nth-degree B-polynomial[32]:
𝐵𝑖,𝑛(𝑥)=(𝑅−𝑥)
𝑅𝐵𝑖,𝑛−1(𝑥)+𝑥
𝑅𝐵𝑖−1,𝑛−1(𝑥)
(10)
The nth-degree B-polynomials' derivatives are polynomials of degree n-1 and are
provided by[32],
𝑑𝐵𝑖,𝑛(𝑥)
𝑑𝑥 =𝑛
𝑅(𝐵𝑖−1,𝑛−1(𝑥)−𝐵𝑖,𝑛−1(𝑥))
(11)
For any real x falling inside the range [0, R], it is easily demonstrated that each B-
polynomial is positive and that the total of all B-polynomials equals unity, or ∑𝐵𝑖,𝑛(𝑥)
𝑛
𝑖=0 . It is
simple to demonstrate that every given degree nth polynomial may be enlarged in terms of a
linear combination of the basis functions[30]–[32]:
9
𝑃(𝑥)=∑𝐶𝑖𝐵𝑖,𝑛(𝑥),𝑛≥1
𝑛
𝑖=0
(12)
This is referred to as a degree n-Bernstein polynomial or a polynomial in Bernstein form. 𝐶𝑖 are
also known as Bézier coefficients or Bernstein coefficients[31].
Caputo’s Fractional Differential-Order Operator
The following is an explanation of the fractional-order derivative of Caputo:
𝐷𝛾𝑓(𝑥)=𝐽𝑚−𝛾𝐷𝑚𝑓(𝑥)=1
𝛤(𝑚−𝛽)∫(𝑥−𝑡)𝑚−𝛾−1𝑓(𝑚)(𝑡)𝑑𝑡,
𝑥
0
𝑓𝑜𝑟𝑚−1<𝛾≤𝑚,𝑚∈𝑁,𝑥>0,𝑓∈𝐶−1
𝑚,
(13)
where in Equation (13), 𝑓(𝑥) is the fractional variable and 𝐷𝛾 is Caputo's fractional operator or
fractional derivative. A constant's derivative by Caputo has a value of zero. In other words,
𝐷𝛾𝐶=0 and C is a constant, the fractional derivative of the Bernstein-poly 𝐷𝛾𝑥𝛼 is given by,
𝐷𝛾𝑥𝛼={0𝑓𝑜𝑟𝛼∈𝑁0𝑎𝑛𝑑𝛼<[𝛾]
𝛤(𝛼+1)𝑥𝛼−𝛾
𝛤(𝛼+1−𝛾)𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒.
(14)
The fractional function is denoted by the letter α, while the derivative is characterized by the
letter 𝛾.
Bhatti polynomial (Fractional-Order Bernstein-Poly Basis): The fractional, 𝐵𝑖,𝑛(𝛼,𝑥) in
its generalized form. B-polys are defined in terms of variable x across the range [0, R] [33], [34],
𝐵𝑖,𝑛(𝛼,𝑥)=∑𝛽𝑖,𝑘
𝑛
𝑖=0 (𝑥
𝑅)𝛼𝑘.
(17)
10
Fractional-order parameter 𝛼 denotes the degree of the Bhatti-fractional poly. In Eq. (17),
every n-value has (n + 1) fractional-order B-polynomials associated with it. The factor 𝛽𝑖,𝑘 is
defined as follows in Eq.
𝛽𝑖,𝑘 =(−1)𝑖−𝑘(𝑛
𝑘)(𝑘𝑖).
(18)
The definition of this binomial coefficient is (𝑛𝑖)= 𝑛!
𝑖!(𝑛−𝑖)!, The non-zero fractional
polynomials might be generated by Mathematica or Maple software using a straightforward
prewritten algorithm that supports any value of n across an interval. The first and last
polynomials in the basis set are often connected to the problem's boundary conditions.
Generation of approximate solution: The process takes benefit of the continuous and the unitary
property of the B-Polys [35]. The new development is successfully applied to partial differential
equations in two variables on a closed interval [0, R]. The B-Polys definite integration matrix
elements are converted into a big operational matrix providing flexibility to include initial and
boundary conditions for the problems at hand. Note that the set of B-Polys formed a complete
basis set of polynomials. The same set of B-Polys has been used to solve the partial differential
equations. For example, the desired result of the partial differential equation may be expressed
in terms of the B-Poly basis set, i.e.
𝑈(𝑥,𝑡)=∑𝑎𝑖(𝑡)𝐵𝑖,𝑛(𝑥)
𝑛
𝑖=0 ,
(19)
where 𝑎𝑖(𝑡) is the ith variable coefficient of the linear mixture in Eq. (19). In Eq. (19), we impose
the initial conditions on variable x. The coefficients 𝑎𝑖(𝑡) are a function of variable (t) and
11
𝐵𝑖,𝑛(𝑥) is nth degree B-Poly in variable x. Furthermore, the coefficients 𝑎𝑖(𝑡) can be developed
in terms of the constant coefficients 𝑏𝑘
𝑖 and the B-Polys as polynomials in t, such as,
𝑎𝑖(𝑡)=∑𝑏𝑘
𝑖𝐵𝑖,𝑛(𝑡).
𝑛
𝑖=0
(20)
Over the interval [0, T]. In Eq. (20), initial conditions on variable t can be imposed. We
present estimate solutions to several partial differential equations using a complete set of B-Polys
of degree n, explaining the procedure and comparing the graphs of the results obtained in the
following sections. To avoid repetition, we will reference our earlier work [28], [29], [32], [36]–
[38], where we have laid the groundwork for generalizing a complete set of B-Poly basis sets
that are employed to estimate solutions to a variety of differential equations. Graphs of B-Polys
were also supplied to show unique characteristics of the B-Poly basis set for calculating the
solutions [35], [39]–[42]. The paper addresses the IOFDE equation in two variables (x, t). As
mentioned above, for the application process, we first integrate internal products with respect to
x and second, integrate with respect to t to convert the IOFDE equation into the operational
matrix with a nonzero determinant which is then inverted to determine the desired solution from
Eq. (19) and Eq. (20). The initial conditions can be directly imposed on both equations to start
the process for determining the desired solution of the IOFDE. For one variable, PDE, we have
presented error analysis in many published papers [35], [40], [41], [43], [44]. Therefore, in the
following section, we shall provide converged solutions of the IOFDEs with the number of
polynomials used in both variables x and t.
12
It is possible to think of the generalized fractional-order B-polynomials 𝐵𝑖,𝑛(𝛼,𝑥) used to
extend the unknown one-variable dependent function 𝑦(𝑥). The following equation represents an
approximation to the nonlinear fractional differential equation,
𝑦(𝑥)=∑𝑏𝑗𝑖𝐵𝑖,𝑛(𝛼,𝑥)+𝑓(𝑥).
𝑛
𝑖,𝑗=0
(15)
A j-th and n-degree fractional-order B-poly as a fractional-order parameter across an
interval and 𝑓(𝑥) as the initial condition of a B-poly in the variable x. The Bernstein expansion
coefficients 𝑏𝑗𝑖 in Equation (15) reflect the variables that are determined using the Galerkin
method of minimization. Fractional differentiation may be carried out by employing Caputo's
derivative property as a linear operator.
𝐷𝑥𝛾(∑𝑏𝑗𝑖𝐵𝑖,𝑛(𝛼,𝑥)+𝑓(𝑥)
𝑛
𝑖,𝑗=0 )=∑𝑏𝑗𝑖(𝐷𝑥𝛾(𝐵𝑖,𝑛(𝛼,𝑥)))
𝑛
𝑖,𝑗=0 .
(16)
13
CHAPTER III
SOLVING LINEAR PARTIAL DIFFERENTIAL EQUATIONS USING BERNSTEIN
POLYNOMIAL BASES
This paper aims to broaden processes for solving hyperbolic partial differential equations
in two variables (x, t) using a mixture of B-Poly basis, operational matrix, and Galerkin method.
The process takes benefit of the continuous and the unitary property of the B-Polys [35].
The new development is successfully applied to partial differential equations in two variables
on a closed interval [0, R]. The B-Polys definite integration matrix elements are converted into a
big operational matrix providing flexibility to include initial and boundary conditions for
the problems at hand. Note that the set of B-Polys formed a complete basis set of polynomials.
The same set of B-Polys has been used to solve the partial differential equations. For
example, the desired result of the partial differential equation may be expressed in terms of the
B-Poly basis set, i.e.
𝑈(𝑥,𝑡)=∑𝑎𝑖(𝑡)𝐵𝑖,𝑛(𝑥),
𝑛
𝑖=0
(19)
where 𝑎𝑖(𝑡) is the ith variable coefficient of the linear mixture in Eq. (19). In Eq. (19), we impose
the initial conditions on variable x. The coefficients 𝑎𝑖(𝑡) are a function of variable (t) and
𝐵𝑖,𝑛(𝑥) is nth degree B-Poly in variable x. Furthermore, the coefficients 𝑎𝑖(𝑡)can be developed
in terms of the constant coefficients 𝑏𝑘
𝑖 and the B-Polys as polynomials in t, such as,
14
𝑎𝑖(𝑡)=∑𝑏𝑘
𝑖𝐵𝑖,𝑛(𝑡),
𝑛
𝑖=0
(20)
over the interval [0, T]. In Eq. (20), initial conditions on variable t can be imposed. We present
estimate solutions to several partial differential equations using a complete set of B-Polys of
degree n, explaining the procedure and comparing the graphs of the results obtained in the
following sections. To avoid repetition, we will provide references to our earlier work[35], [40],
[45], where we have laid down the groundwork for the generalization of a complete set of B-
Poly basis sets that are employed to estimate solutions to a variety of differential equations.
Graphs of B-Polys were also supplied to show unique characteristics of the B-Poly basis set to
be utilized in estimating the solutions [35], [39]–[42]. This paper addresses the Hyperbolic
partial differential (HPD) equation in two variables (x, t). As mentioned above, for the
application process, we first integrate internal products with respect to x and second, integrate
with respect to converting the HPD equation into an operational matrix with a nonzero
determinant which is then inverted to determine the desired solution from Eq. (19) and Eq. (20).
The initial conditions can be directly imposed on both equations to start the process for
determining the desired solution of the HPD. For one variable, PDE, we have presented error
analysis in many published papers [35], [40], [41], [43], [44]. Therefore, in the following
section, we shall provide converged solutions of the HPDEs with the number of polynomials
used in both variables x and t.
15
Computation of results on a B-Poly basis
A process to estimate solutions, 𝑈 (𝑥, 𝑡), of the 2-Dimensional (2-D) partial differential
equations, considered a linear mixture of B-Polys in the variables x and t with initial condition
(U (x, 0) = f(x)) built-in is given as,
𝑈(𝑥,𝑡)=∑𝑎𝑖(𝑡)𝐵𝑖(𝑥)+𝑓(𝑥)
𝑛
𝑖=0 .
(21)
With the application of this mixture, Eq. (21), to a PDE, we convert it into a matrix
formation, and in terms of variables, 𝑎𝑖(𝑡). For the sake of simplicity, we shall drop subscript n
in the following sections. A second estimation of variable coefficients 𝑎𝑖(𝑡) is considered by
expanding and imposing an initial condition on the equation,
𝑎𝑖(𝑡)=∑𝑏𝑘
𝑖𝐵𝑘(𝑡)
𝑛
𝑖=0 .
(22)
We employ Galerkin method[46] and inversion of the operational matrix to solve the 2-D
partial differential equations. We also provide several examples to show successful applications
of the proposed process for solving HPD equations. Consider a more general 2-dimensional
HPDE of the form,
𝛼𝑑𝑈(𝑥,𝑡)
𝑑𝑡 +𝛽𝑑𝑈(𝑥,𝑡)
𝑑𝑥 +𝛾𝑈(𝑥,𝑡)=0
(23)
with constants 𝛼, 𝛽, 𝛾, and non-homogenous initial condition 𝑈(𝑥,0)=𝑓(𝑥) at 𝑡=0. Let’s
substitute the desired solution of Eq. (21) into Eq. (5), which returns,
∑𝛼𝑎𝑖(𝑡)𝐵𝑖(𝑥)
𝑛
𝑖=0 +∑𝛽𝑎𝑖(𝑡)𝐵𝑖′(𝑥)
𝑛
𝑖=0 +∑𝛾𝑎𝑖(𝑡)𝐵𝑖(𝑥)
𝑛
𝑖=0 =−𝛽𝑓′.
(24)
16
Here prime (′) and dot (·) denote the differentials with respect to x and t, respectively.
We may multiply the above equation with another B-Poly from the set and integrate it over the
interval [0, R] to generate,
∑[𝛼𝑎𝑖(𝑡)⟨𝐵𝑖(𝑥)|𝐵𝑗(𝑥)⟩+𝛽𝑎𝑖(𝑡)⟨𝐵𝑖′(𝑥)|𝐵𝑗(𝑥)⟩+𝛾𝑎𝑖(𝑡)⟨𝐵𝑖(𝑥)|𝐵𝑗(𝑥)⟩]
𝑛
𝑖=0 =−𝛽⟨𝑓′|𝐵𝑗(𝑥)⟩
(25)
Where matrix elements of the integrals of B-poly products are defined,
𝑀𝑖𝑗 =⟨𝐵𝑖(𝑥)|𝐵𝑗(𝑥)⟩=∫ 𝐵𝑖(𝑥)𝐵𝑗(𝑥)𝑑𝑥
𝑅
0
𝑁𝑖𝑗 =⟨𝐵𝑖′(𝑥)|𝐵𝑗(𝑥)⟩=∫ 𝐵𝑖′(𝑥)𝐵𝑗(𝑥)𝑑𝑥
𝑅
0𝑎𝑛𝑑𝐹𝑗=−⟨𝛽𝑓′|𝐵𝑗(𝑥)⟩
(26)
we obtain the matrix equation,
∑[𝛼𝑎𝑖(𝑡)𝑀𝑖𝑗+𝛽𝑎𝑖(𝑡)𝑁𝑖𝑗+𝛾𝑎𝑖(𝑡)𝑀𝑖𝑗]
𝑛
𝑖=0 =−𝐹𝑗
(27)
In the next step, we build 𝑎𝑖(𝑡) in terms of B-polys as presented in Eq. (22); we attain,
∑∑[𝛼𝑏𝑘
𝑖𝐵𝑘(𝑡)𝑀𝑖𝑗+𝛽𝑏𝑘
𝑖𝐵𝑘(𝑡)𝑁𝑖𝑗+𝛾𝑏𝑘
𝑖𝐵𝑘(𝑡)𝑀𝑖𝑗]
𝑛
𝑖=0
𝑛
𝑘=1 =−𝐹𝑗
(28)
Again, multiplying both sides of the above equation with another 𝐵𝑙(𝑡) from the set and
integrating it in the interval 𝑡 ∈ [0, 𝑇], we obtain a simple representation in terms of the matrix
elements,
∑∑𝑏𝑘
𝑖[𝛼𝑈𝑘𝑙𝑀𝑖𝑗+𝛽𝑁𝑖𝑗𝑉𝑘𝑙+𝛾𝑀𝑖𝑗𝑉𝑘𝑙]
𝑛
𝑖=0
𝑛
𝑘=1 =𝑊𝑗𝑙, 𝑙=𝑗=0,…..𝑛,
(29)
where the elements of matrices in the above Eq. (29) are defined
by,
17
𝑈𝑘𝑙 =⟨𝐵𝑘(𝑡)|𝐵𝑙(𝑡)⟩=∫ 𝐵𝑘(𝑡)𝐵𝑙(𝑡)𝑑𝑡
𝑅
0
𝑉𝑘𝑙 =⟨𝐵𝑘(𝑡)|𝐵𝑙(𝑡)⟩,
𝑊𝑗𝑙 =⟨𝐹𝑗|𝐵𝑙(𝑡)⟩
(30)
The precise solution of Eq. (23) is found, reference [47],
𝑈(𝑥,𝑡)=𝑔(𝛼𝑥−𝛽𝑡)𝑒−𝛾
𝛽𝑥
(31)
The above exact result may be further streamlined subject to the initial conditions. We
shall use the Galerkin Technique[46] to estimate the result of the partial differential equations
(PDE) in both variables (x, t). In the following section, we will utilize the proposed process in
four examples to illustrate how it works for assessing solutions of the HPD equations.
Example 1. We consider an equation which is attained presuming 𝛼 = 2, 𝛽 = 1, and 𝛾 = 0 in the
Eq. (23),
2𝑑𝑈(𝑥,𝑡)
𝑑𝑡 +𝑑𝑈(𝑥,𝑡)
𝑑𝑥 =0.
(32)
Whose solution we are looking for in the intervals 0 ≤ x ≤ 2 and 0 ≤ t ≤ 2, with initial condition.
𝑈 (𝑥, 0) = 𝑓(𝑥) = 𝑥3, at 𝑡 = 0. An estimate solution to Eq. (9) may be written via Eq. (21),
𝑈(𝑥,𝑡)=∑𝑎𝑖(𝑡)𝐵𝑖(𝑥)+𝑥3, 𝑛
𝑛
𝑖=0 ≥1
(33)
The 𝑎𝑖(𝑡) is the 𝑖𝑡ℎ variable coefficient in the development above Eq. (33). After substituting Eq.
(33) into Eq. (32), we shall come at Eq. (24), which would now appear like,
18
∑[2𝑎𝑖(𝑡)𝑀𝑖𝑗+𝑎𝑖(𝑡)𝑁𝑖𝑗]
𝑛
𝑖=0 =𝐹𝑗
(34)
Again, we make another approximation to the 𝑎𝑖(𝑡) coefficients using the expansion
provided in Eq. (22). By replacing Eq. (22) with Eq. (33) with a little bit of oversimplification,
and from the variational property with respect to the coefficients, we may attain an expression
given in Eq. (25) which streamlines to the following expression,
∑∑𝑏𝑘
𝑖[2𝑈𝑘𝑙𝑀𝑖𝑗+𝑁𝑖𝑗𝑉𝑘𝑙]
𝑛
𝑖=0
𝑛
𝑘=1 =𝑊𝑗𝑙, 𝑙=𝑗=0,…..𝑛,
(35)
The matrix elements are clearly provided,
𝑀𝑖𝑗 =∫ 𝐵𝑖(𝑥)𝐵𝑗(𝑥)𝑑𝑥,𝑁𝑖𝑗 =∫𝐵𝑖′(𝑥)𝐵𝑗(𝑥)𝑑𝑥
2
0
2
0
𝑈𝑘𝑙 =∫𝐵𝑘(𝑡)𝐵𝑙(𝑡)𝑑𝑡
2
0, 𝑉𝑘𝑙 =∫𝐵𝑘(𝑡)𝐵𝑙(𝑡)𝑑𝑡
2
0
𝑊𝑗𝑙 =𝐹𝑗∫𝐵𝑙(𝑡)𝑑𝑡
2
0𝑎𝑛𝑑𝐹𝑗=−∫𝑓′(𝑥)𝐵𝑙(𝑥)𝑑𝑥
2
0
(36)
This algorithm leads to an (𝑛 + 1)2 by (𝑛 + 1)2 system of equations 𝐴𝐵 = 𝑊, in the unknown
variables 𝐵0, 𝐵1, 𝐵2, … 𝐵𝑛 where the matrix A is given by,
𝐴 = 2 𝑀 𝑈 + 𝑁 𝑉
(37)
Clearly, the HPD equation is converted into a large matrix whose inverse provides
specific values of the unknown coefficients 𝑏𝑘
𝑖 of the linear mixture in Eq. (22) by solving
the B = A-1W. Then the estimated solution is constructed from the product of these
coefficients 𝑎𝑖(𝑡) and B-Poly basis set, 𝐵𝑖(𝑥) in the Eq. (32). Before solving the matrix
equation 𝐴 𝐵 = 𝑊, we also employ initial conditions on the matrix by deleting rows and
19
columns of matrices A and W defined in the Eq. (35). The precise solution for the Eq.
(27) is obtained analytically applying the initial condition via the Eq. (26) which is given,
𝑈𝑒𝑥𝑎𝑐𝑡(𝑥,𝑡)=(𝑥−𝑡2)3,𝑤ℎ𝑒𝑛𝑡=𝑥,𝑤𝑒𝑔𝑒𝑡𝑈(𝑥)=(𝑥
2)3.
(38)
In this example, both the exact solution in Eq. (32) and the estimated solution after
ignoring small terms in Eq. (33) are equivalent. To find the numerical solution of the differential
Eq. (32), we only used 𝑛=3-degree B-Polynomials. In Fig. 1, we submit a plot of the absolute
difference between the estimate and exact solutions supplied in Eq. (38) and Eq. (39) at 𝑡 = 𝑥 and
at 𝑡 ≠ 𝑥. Both solutions overlap, showing no noticeable differences. In example 1, the absolute
difference is so tiny that it reaches the order of 10-15, which suggests that the estimated solution is
in superb agreement with the exact one. This kind of precision was achieved only with 𝑛=3
polynomials; see Fig. 1(a). Also, In Fig. 1 on the right, we have provided the 3-D graphs of the
exact and approximate solutions for comparisons. The analytic solution with 𝑡 = 𝑥 is provided
below:
𝑈(𝑥)=.125𝑥3𝑎𝑡𝑡=𝑥,𝑎𝑛𝑑
𝑈 (𝑥, 𝑡) = 𝑥3 + 𝑡2(1.8652 × 10−14 + 0.7500 𝑥 + 9.9920 × 10−15 𝑥2 − 1.3323 × 10−15 𝑥3) +
𝑡3(−0.1250 + 3.0309 × 10−14𝑥 − 7.6605 × 10−15 𝑥2 − 3.6082 × 10−16 𝑥3) + 𝑡 (−1.0658 ×
10−14 + 6.3949 ×10−14𝑥 − 1.5000𝑥2 + 1.9651 × 10−14𝑥3)
(39)
20
Figure 1. A narrative of the absolute difference between exact and estimated solutions is illustrated on the
left for t = x. A 3-D graph U (x, t) of both exact and estimated solutions is shown on the right. To resolve
the differential Eq. (32), only n = 3-degree B-Polynomials were used. A narrative of the absolute
difference between exact and estimated solutions is illustrated on the left for t = x. A 3-D graph 𝑈 (𝑥, 𝑡) of
both exact and estimated solutions is shown on the right in intervals 𝑥 ∈ [0, 2] and 𝑡 ∈ [0, 4] and. To
resolve the differential Eq. (32), only n = 3-degree B-Polynomials were used. Results of different closed
intervals, for example, 𝑥 ∈ [0, 2] and 𝑡 ∈ [0, 4], provided similar accuracy of the order of 10-14. The
approximate solution to Eq. (32) for different intervals is given below in both variables. The 3-D graph is
also included in Fig. 1(b).
(a)
(b)
21
U (x, t) = 𝑥3 + 𝑡2𝑥 (0.7500 − 1.5987 × 10−14𝑥 + 1.4988 × 10−15𝑥2) + 𝑡3(−0.1250 + 2.9976 ×
10−15𝑥 + 8.3267 × 10−17𝑥2 + 6.9389 × 10−18𝑥3) + 𝑡 (−2.1316 × 10−14 + 4.7962 × 10−14𝑥 −
1.5000 𝑥2 +1.6986 × 10−14𝑥3)
(40)
For simplicity, in the following examples, we will consider the integration intervals the
same for both variables x and t.
Example 2. Let us consider another HPD Equation with slightly different initial conditions. We
substitute 𝛼 = 1, 𝛽 = 1, 𝛾 = 0, and the initial condition in the general hyperbolic Eq. (23). So that
the HPDE in two dimensions is given,
2𝑑𝑈(𝑥,𝑡)
𝑑𝑡 +𝑑𝑈(𝑥,𝑡)
𝑑𝑥 =0.
(41)
with the initial condition U (x, 0) = f(x) = sin (x) at t = 0. The result of Eq. (41) is sought
in the intervals 0 ≤ x ≤ 2 and 0 ≤ t ≤ 2. We shall follow the same process as in Example 1 to
determine the resolution of Eq. (41) in the preferred regions, with the initial condition at t = 0.
According to Eq. (21), the estimated solution is,
𝑈(𝑥,𝑡)=∑𝑎𝑖(𝑡)𝐵𝑖(𝑥)+𝑆𝑖𝑛(𝑥), 𝑛
𝑛
𝑖=0 ≥1
(42)
Substituting Eq. (42) in terms of the basis set of B-Polys of degree n=9, we will reach Eq. (35),
as in Example 1, only the elements of column matrix 𝑊𝑗ℓ will change too,
𝑊𝑗𝑙 =𝐹𝑗∫𝐵𝑙(𝑡)𝑑𝑡
2
0𝑎𝑛𝑑𝐹𝑗=−∫cos(𝑥)𝐵𝑗(𝑥)𝑑𝑥
2
0
(43)
The exact solution of Eq. (41) is,
𝑈(𝑥,𝑡)=𝑆𝑖𝑛(𝑥−𝑡2),𝑤ℎ𝑒𝑛𝑡=𝑥,𝑤𝑒𝑔𝑒𝑡𝑈(𝑥)=sin(𝑥
2).
(44)
This can be directly attained from Eq. (31) using 𝛼 = 2, 𝛽 = 1, 𝛾 = 0, and involving the
initial condition 𝑓(𝑥) = sin(𝑥). Also, the estimated solution to Eq. (41) is provided in Eq. (45),
22
which is contrasted with the exact solution at t = x. The 3-D graph of both solutions (exact and
estimated) is also supplied, and the contrast is shown in Fig. 2 on the right. A narrative of the
absolute difference between the estimated and exact solutions at t = x is exhibited in Fig. 2. The
precision is of the order of 10-11 with only n = 9-degree B-Polys.
𝑈(𝑥) = 0.5 𝑥 − 0.0208 𝑥3 + 0.0003 𝑥5 − 1.5501 × 10−6 𝑥7 + 5.3823 × 10−9𝑥9− 1.2232 ×10−11
𝑥11 + 1.9603 × 10−14 𝑥13 − 2.3337 × 10−17𝑥15+ 2.1450 × 10−20 𝑥17 − 1.5680 × 10−23 𝑥19 (−0.5
+ 0.25𝑥2 − 0.0208𝑥4) + 𝑡2(−0.125𝑥) + 𝑡3(0.0208) + sin (𝑥)
(45)
Figure 2. A narrative of the absolute difference between exact and estimated solutions is depicted
on the left for t = x. A 3-D graph
U
(𝑥,𝑡) of both exact and estimated solutions is given on the
right. To estimate the result of the differential Eq. (41), only n = 9-degree B-Polynomials were
utilized.
The complete solution in terms of both variables x and t is provided below:
𝑈 (𝑥, 𝑡) = 𝑡9(−9.5824 × 10−9 + 4.4224 × 10−9𝑥) + 𝑡2(−4.5121 × 10−9 − 0.125𝑥 − 1.4593 ×
10−6𝑥2 + 0.0208𝑥3 − 1.4117 × 10−5𝑥4 − 0.0010 𝑥5 − 1.6643 × 10−5𝑥6 +3.3430 × 10−5𝑥7
−2.5473×10−6𝑥8) + 𝑡8(1.1863 × 10−9 + 1.3059 × 10−7𝑥 − 3.9801 × 10−8𝑥2) + 𝑡3(0.0208 +
(46)
23
4.8643 × 10−7𝑥 − 0.0104 𝑥2 + 9.4115 × 10−6𝑥3 +8.5177 × 10−4𝑥4 + 1.6644 × 10−5𝑥5 −
3.9001 × 10−5𝑥6 + 3.3964 × 10−6𝑥7) +𝑡7(1.4947 × 10−6 + 1.4860 × 10−7𝑥 − 1.0447 × 10−6𝑥2
+ 2.1227 × 10−7𝑥3)+𝑡4(−3.7801 × 10−8 + 0.0026 𝑥 − 3.5293 × 10−6𝑥2 − 4.2588 × 10−4𝑥3 −
1.0402 ×10−5𝑥4 + 2.9251 × 10−5𝑥5 − 2.9718 × 10−6𝑥6) + 𝑡6(−1.8569 × 10−8 − 2.1294 ×
10−5𝑥 − 1.0402 × 10−6𝑥2 + 4.8752 × 10−6𝑥3 − 7.4296 × 10−7𝑥4) + 𝑡5(−2.6053 × 10−4
+7.0586 × 10−7𝑥 + 1.2774 × 10−4𝑥2 + 4.1609 × 10−6𝑥3 − 1.4626 × 10−5𝑥4 +1.7831 ×
10−6𝑥5) + 𝑡 (−0.5000 + 2.4115 × 10−8𝑥 + 0.2500𝑥2 + 1.9457 × 10−6𝑥3 − 0.0208 𝑥4 + 1.1294
× 10−5𝑥5 + 6.8142 × 10−4𝑥6 + 9.5106 × 10−6𝑥7 − 1.6715 × 10−5𝑥8+1.1321 × 10−6𝑥9) +
sin(𝑥)
Example 3. Consider another more complicated HPD equation to show that the proposed
process is extremely beneficial to estimating a solution to an anticipated accuracy. The
hyperbolic differential equation is,
3𝑑𝑈(𝑥,𝑡)
𝑑𝑡 +2𝑑𝑈(𝑥,𝑡)
𝑑𝑥 +𝑈(𝑥,𝑡)=0.
(47)
With the initial condition 𝑈(𝑥,0)=𝑠𝑖𝑛(𝑥) at t = 0. Eq. (47) has an exact answer under
the above initial condition, which we achieved after substituting 𝛼 = 3, 𝛽 = 2, and 𝛾 = 1 in the
Eq. (31), 𝑈(𝑥,𝑡)=sin(𝑥−2
3𝑡)𝑒−1
3𝑡. We utilized a basis set of 9 B-Polys to estimate the
solution of Eq. (47) and contrasted it with the exact result of this equation. The same process
was used to seek the estimated solution as given in Examples 1 and 2. The procedure provides
matrix elements of the equation 𝐴 𝐵 = 𝑊 as follows,
∑∑𝑏𝑘
𝑖[3𝑈𝑘𝑙𝑀𝑖𝑗+2𝑁𝑖𝑗𝑉𝑘𝑙+𝑀𝑖𝑗𝑉𝑘𝑙]
𝑛
𝑖=0
𝑛
𝑘=1 =𝑊𝑗𝑙,
(48)
24
The matrix elements of Eq. (48) are defined in Eq. (29). The only change in the matrix 𝐹𝑗
is,
𝐹𝑗=−∫cos(𝑥)𝐵𝑗(𝑥)𝑑𝑥
2
0
(49)
We note that in the above examples considered, the inverse of matrix A in the equation A
B = W is acquired after imposing the initial condition on the matrix A elements at t = 0 and 𝑥 =
0. Finally, the equation A B = W is resolved for the unknown coefficients to construct the
anticipated solution of the differential equation (47). The final estimated solution to Eq. (47) is
provided in Eq. (49) over the intervals 0 ≤ 𝑥 ≤ 2 and 0 ≤ 𝑡 ≤ 2. Both the exact and the estimated
results of the partial differential Eq. (47) were equated. The absolute difference between the
exact and estimated analytic solution of Eq. (47) over the ranges is exhibited in Fig. 3. The
absolute discrepancy between the solutions is of the order of 10-11 with the usage of only n = 9-
degree B-Poly basis. However, the anticipated accuracy of the numeric solution of the
differential equation depends on the size of the basis set chosen and the degree of the
polynomials. The larger the basis set, the better the estimated solution's accuracy. However, the
drawback is that the larger the size of the matrix, the larger the CPU time required to invert the
matrix. The estimated solution of Eq. (47) for t = x is shown employing the set of n = 9-degree
polynomials in Eq. (50). We also present a 3-D graph of the solution in Fig. 3 on the right,
which is overlapping.
𝑈(𝑥) = −0.6667 𝑥 − 0.1111 𝑥2 + 0.1790 𝑥3 + 2.2030 × 10−6 𝑥4 − 8.4753 × 10−3 𝑥5 +2.1878
× 10−5 𝑥6 + 1.9175 × 10−4 𝑥7 + 3.4687 × 10−6 𝑥8 − 4.0656 × 10−6 𝑥9 + 3.1480 ×10−7 𝑥10
−2.3201 × 10−8 𝑥11 + 4.8104 × 10−9 𝑥12 − 1.8274 × 10−10 𝑥13 − 8.2700 ×10−11 𝑥14 − 1.9815
(50)
25
× 10−13 𝑥15 + 5.0809 ×10−13 𝑥16 + 2.3476 × 10−14 𝑥17 + 3.0402 ×10−16 𝑥18 + sin(𝑥)
Figure 3: illustrates a narrative of the fundamental difference between exact and estimated
solutions on the left for t = x. A 3-D graph U (x, t) of both exact and estimated results is given on
the right. To solve the differential Eq. (20), only n = 9-degree B-Polynomials were exploited.
The full approximate solution in terms of both variables x and t is provided here:
𝑈 (𝑥, 𝑡) = 𝑡(−0.6667 − 0.3333 𝑥 + 0.3333 𝑥2 + 0.0556𝑥3 − 0.0278𝑥4 − 0.0028 𝑥5 +8.9890 ×
10−4𝑥6 + 8.5163 × 10−5𝑥7 − 2.4701 × 10−5𝑥8 + 1.0249 × 10−6𝑥9 + 𝑡 (0.2222 − 0.1666𝑥 −
0.1111𝑥2 +0.0278𝑥3 + 0.0092 𝑥4 − 0.0013𝑥5 − 0.00035𝑥6 + 5.1675 ×10−5𝑥7 + 1.0422 ×
10−6𝑥8 − 1.7081 × 10−7𝑥9) + 𝑡3(−0.0123 − 0.0036𝑥 + 0.0062𝑥2 +6.3164 × 10−4𝑥3 − 5.4982
× 10−4𝑥4 − 9.4638 × 10−6𝑥5 + 1.2332 × 10−5𝑥6 +3.3176 × 10−7𝑥7 − 4.7288 × 10−8𝑥8 −
1.5816 × 10−9𝑥9) + 𝑡5(8.49109 × 10−5 +2.2629 × 10−4𝑥 − 4.7331 × 10−5𝑥2 − 3.3259 ×
10−5𝑥3 + 2.7273 × 10−6𝑥4 +1.269 × 10−6𝑥5 + 4.1096 × 10−8𝑥6 − 6.0505 × 10−9𝑥7 − 3.8576 ×
10−10𝑥8 − 5.8558 × 10−12𝑥9) + 𝑡7(1.64103 × 10−6 − 1.6285 × 10−6𝑥 − 7.6756 × 10−7𝑥2
(51)
26
+1.6234 × 10−7𝑥3 + 6.1308 × 10−8𝑥4 + 2.8908 × 10−9𝑥5 − 3.5597 × 10−10𝑥6 − 3.6431 ×
10−11𝑥7 − 1.1130 × 10−12𝑥8 − 1.1089 × 10−14𝑥9) + 𝑡8(−1.4845 × 10−8 +1.1387 × 10−7𝑥 +
1.4398 × 10−8𝑥2 − 1.28718 × 10−8𝑥3 − 2.41918 × 10−9𝑥4 − 2.0921 × 10−11𝑥5 + 1.8858 ×
10−11𝑥6 + 1.3434 × 10−12𝑥7 + 3.4565 × 10−14𝑥8 +3.0402 × 10−16𝑥9) + 𝑡6(−2.6438 × 10−5 −
2.4396 × 10−6𝑥 + 1.2006 × 10−5𝑥2 +6.7251 × 10−7𝑥3 − 7.2374 × 10−7𝑥4 − 8.4726 × 10−8𝑥5 +
1.9274 × 10−9𝑥6 +5.7452 × 10−10𝑥7 + 2.3230 × 10−11𝑥8 + 2.7765 × 10−13𝑥9) + 𝑡4(0.0013 −
0.0014𝑥 −6.6280 × 10−4𝑥2 + 2.5132 × 10−4𝑥3 + 4.2979 × 10−5𝑥4 − 9.2357 × 10−6𝑥5 −
1.1320 × 10−6𝑥6 + 2.8312 × 10−8𝑥7 + 5.0502 × 10−9𝑥8 + 1.0544 × 10−10𝑥9) +𝑡2(0.01235 +
0.0679𝑥 − 0.0062𝑥2 − 0.0113𝑥3 + 4.6301 × 10−4𝑥4 + 6.1347 × 10−4𝑥5 −4.16579 × 10−5𝑥6 −
7.5944 × 10−6𝑥7 + 2.258 × 10−7𝑥8 + 1.8979 × 10−8𝑥9)) + sin(𝑥)
Example 4. The final example we consider is that of a hyperbolic differential equation,
3𝑑𝑈(𝑥,𝑡)
𝑑𝑡 +2𝑑𝑈(𝑥,𝑡)
𝑑𝑥 +8𝑈(𝑥,𝑡)=0.
(52)
With initial condition 𝑈(𝑥,0)=𝑐𝑜𝑠(𝑥). The Eq. (52) is obtained by exchange of 𝛼 = 3, 𝛽 = 2,
and 𝛾 = 8 in the general Eq. (23). The precise solution subjected to the initial condition 𝑈 (𝑥, 0) =
cos (𝑥) is,
𝑈(𝑥,𝑡)=cos(3𝑥−2𝑡)𝑒−8
3𝑡
(53)
Expanding the solution of Eq. (52) in the estimated form via Eq. (21),
𝑈(𝑥,𝑡)=∑𝑎𝑖(𝑡)𝐵𝑖(𝑥)+𝑐𝑜𝑠(𝑥),
𝑛
𝑖=0
(54)
Substituting Eq. (54) into the Eq. (52) generates an equation in the form of Eq. (24).
27
Furthermore, the coefficients 𝑎𝑖(𝑡) in Eq. (24) are developed in terms of B-Polys as a function of
(t) in order to determine the coefficients of Eq. (22). This process leads to an articulation
provided in Eq. (48). All the matrix elements remain the same as in Example 3, except the
elements of 𝑊𝑗ℓ which are furnished in the matrix,
𝑊𝑗𝑙 =𝐹𝑗∫𝐵𝑙(𝑡)𝑑𝑡
2
0𝑎𝑛𝑑𝐹𝑗=−∫sin(𝑥)𝐵𝑗(𝑥)𝑑𝑥
2
0
(55)
Eq. (54) coefficients were attained by employing the Galerkin method[46] and inverting
the matrix A. The values of these coefficients were manipulated to determine the coefficients
𝑎𝑖(𝑡) in terms of the B-Polys as a function of (t) and the result to Eq. (52) was analyzed from Eq.
(52). Only n = 9, B-Polys of degree 9 were used to estimate the solution of the Eq. (54). The
errors were of the order of 10-11 between the exact and the estimated solutions of the Eq. (52).
Fig. 4 on the left shows a narrative of the absolute discrepancy between the estimated and the
exact results when t = x is substituted in the solutions. The absolute error between the two
solutions improved as the number of B-Polys was systematically increased from 2 through 9-
degree polynomials. A 3-D graph is also supplied in Fig. 4 on the right to show the comparison
of both solutions. Again, this example demonstrates the validity of the process employed in the
two-dimensional HPD equations. The estimated solution of Eq. (52) for t = x is submitted using
the set of n = 9-degree polynomials in Eq. (56).
𝑈(𝑥) = −0.3333 𝑥 + 0.5000 𝑥2 + 0.0123 𝑥3 − 0.0437 𝑥4 + 0.00014 𝑥5 + 0.001 𝑥6+ 3.4590 ×
10−7 𝑥7 − 2.5591 × 10−5 𝑥8 + 4.2254 × 10−7 𝑥9+ 1.5033 × 10−7 𝑥10 + 1.3902 × 10−8 𝑥11 +
1.3342 × 10−9𝑥12− 1.1149 × 10−9 𝑥13 − 8.4884 × 10−13 𝑥14 + 2.0058 × 10−11 𝑥15 + 1.6700 ×
10−12 𝑥16 + 4.9942 × 10−14 𝑥17 + 5.114 × 10−16 𝑥18 + cos (x)
(56)
28
Figure 4: A narrative of the absolute difference between exact and estimated solutions is depicted
on the left for t = x. A 3-D graph 𝑈(𝑥,𝑡)of both exact and estimated solutions is given on the
right.
To solve the differential Eq. (52), only n = 9-degree B-Polynomials were employed. The
full estimated solution in terms of both variables x and t is provided here.
𝑈 (𝑥, 𝑡) = 𝑡(−0.3333 + 0.6666 𝑥 + 0.1667 𝑥2 − 0.1111𝑥3 − 0.0139𝑥4 + 0.0055𝑥5 + 4.7358 ×
10−4𝑥6 − 1.3862 × 10−4𝑥7 − 6.3152 × 10−6𝑥8+ 1.7240 × 10−6𝑥9 + 𝑡 (−0.1667 − 0.2222𝑥 +
0.0833𝑥2 + 0.0370𝑥3− 0.00693 𝑥4 − 0.00188 𝑥5 + 2.4453 × 10−4𝑥6 + 3.9945 × 10−5𝑥7−
4.1195 × 10−6𝑥8 − 2.8733 × 10−7𝑥9) + 𝑡3(−0.0036 + 0.0123𝑥+ 0.0018𝑥2 − 0.0021 𝑥3 −
1.3895 × 10−4𝑥4 + 9.9155 × 10−5𝑥5+ 4.0766 × 10−6𝑥6 − 1.6169 × 10−6𝑥7 − 1.3392 ×
10−7𝑥8− 2.6605 × 10−9𝑥9) + 𝑡5(2.2406 × 10−5 − 8.6208 × 10−5𝑥− 1.1067 × 10−4𝑥2 + 1.3628
× 10−5𝑥3 + 8.7271 × 10−6𝑥4− 2.6396 × 10−7𝑥5 − 2.6283 × 10−7𝑥6 − 2.4669 × 10−8𝑥7−
8.5023 × 10−10𝑥8 − 9.8505 × 10−12𝑥9) + 𝑡7(−1.7547 × 10−6− 1.7167 × 10−6𝑥 + 7.6125 ×
10−7𝑥2 + 3.1042 × 10−7𝑥3− 3.1053 × 10−8𝑥4 − 1.7725 × 10−8𝑥5 − 2.0325 × 10−9𝑥6− 1.0041
× 10−10𝑥7 − 2.2535 × 10−12𝑥8 − 1.8652 × 10−14𝑥9)+ 𝑡8(1.2966 × 10−7 +2.4218 × 10−8𝑥 −
(57)
29
5.8362 × 10−8𝑥2− 8.7093 × 10−9𝑥3 + 3.3734 × 10−9𝑥4 + 9.0003 × 10−10𝑥5+ 8.1697 ×
10−11𝑥6 + 3.4564 × 10−12𝑥7 + 6.8594 × 10−14𝑥8+ 5.1139 × 10−16𝑥9) + 𝑡6(−3.1649 × 10−6 +
2.6014 × 10−5𝑥+ 1.8144 × 10−6𝑥2 − 4.0460 × 10−6𝑥3 − 3.3892 × 10−7𝑥4 + 1.5531 × 10−7𝑥5 +
2.8446 × 10−8𝑥6 + 1.8045 × 10−9𝑥7+ 4.8622 × 10−11𝑥8 + 4.6704 × 10−13𝑥9) + 𝑡4(−0.00141 −
0.00130 𝑥+ 7.0807 × 10−4𝑥2 + 0.000212𝑥3 − 5.6912 × 10−5𝑥4 − 9.8825 × 10−6𝑥5+ 1.2366 ×
10−6𝑥6 + 2.5062 × 10−7𝑥7 + 1.2120 × 10−8𝑥8+ 1.7736 × 10−10𝑥9) + 𝑡2(0.0679 − 0.0124𝑥 −
0.0339𝑥2 + 0.0020𝑥3+ 0.00285𝑥4 − 0.0001.1811 × 10−4𝑥5 − 8.9306 × 10−5𝑥6+ 2.8852 ×
10−6𝑥7 + 1.0324 × 10−6𝑥8 + 3.1926 × 10−8𝑥9) + cos (𝑥)
Results and Discussions
In this article, we have given a broad description of the 2-D algorithm to show how the
B-Polynomial basis may be employed to provide highly accurate results of the 2-D Hyperbolic
Partial Differential Equations (HPDEs). To describe the application of this process, we laid
down a 2-D groundwork. We utilized it to solve four examples of the HPD equations with
various initial conditions enforced on the results. In each of the four examples worked out, we
were capable of contrasting the exact and the estimated solutions achieved using the Galerkin
method [46] in two variables (x, t), and an agreement was discovered to desired accuracy of 10-11
as exhibited in Figs. 1–4. We have noticed that increasing the number of B-Polys in the
estimated solutions increases the precision of the results [35]. All computations and analytic
integrations over the intervals were conducted utilizing Wolfram Mathematica's symbolic
program version 11 [48]. Comparisons between the exact and the estimated solutions were
shown in the 2-D and 3-D graphs in Figs. 1-4. In each case, the precision of the results was
increased by boosting the number of B-polys in the basis set. Furthermore, errors were also
30
contrasted with the exact results of the PDEs when a parameter (t) was set equal to x and
precision was shown to be better than the order of 10-11 for all examples considered.
The current process, which utilizes the continuous B-Polynomials, may offer great
potential for solving linear and nonlinear examples of 2-D problems in other disciplines,
particularly in physics. More recently, several authors[43], [44], [49], [50] have used B-polys
techniques to construct an operational matrix to solve a variety of differential equations in one-
dimensional variables. For the first time, we successfully extended the method to solve HPDEs
in two variables. We had already published KdV and nonlinear Burger equations results on a B-
polynomial basis. Our method worked very well for solving both equations using an operational
matrix scheme [43], [44]. The only difference is that we had discretized the time variable using
the fourth-order Runge-Kuttta [44] method while the spatial variable was expanded in terms of
the B-polynomial basis. First, we would like to have the present valuable procedure published.
By implementing the current method, we plan to publish some examples of nonlinear PDEs
separately, such as the Burger equation. All calculations were performed using the symbolic
Mathematica Code [48]. Typically, the CPU time used to perform all the calculations for each
example is about 32 seconds, except for example 1, the CPU time used was about 2.5 seconds.
This article shows that based on the B-Poly technique, the process will give relevant
results for 2-D partial differential equations after the initial conditions are enforced on the
operational matrix. It is a powerful tool that we may utilize to surmount the difficulties
associated with complex systems of differential equations where there are no exact solutions
available, particularly in two-variable differential equations. It has also been established to be
effective in producing precise results and could be quickly executed in various disciplines.
31
CHAPTER IV
APPROXIMATE SOLUTIONS OF NONLINEAR PARTIAL DIFFERENTIAL EQUATIONS
USING B-POLYNOMIAL BASES
Muhammad I. Bhatti*, Md. Habibur Rahman, and N. Dimakis
University of Texas Rio Grande Valley, Edinburg Texas, 78539
*Corresponding author: [email protected]
Abstract
A multivariable technique has been incorporated for guesstimating solutions of Nonlinear
Partial Differential Equations (NPDE) using bases set of B-Polynomials (B-polys). To
approximate the anticipated solution of the NPD equation, a linear product of variable
coefficients 𝑎𝑖(𝑡) and B-polys 𝐵𝑖(𝑥) has been employed. Additionally, the variable quantities in
the anticipated solution are determined using the Galerkin method to minimize errors. Before
the minimization process takes place, the NPDE is converted into an operational matrix equation
that yields values of the undefined coefficients in the expected solution when inverted. The
nonlinear terms of the NPDE are combined in the operational matrix equation using the initial
guess and iterated until converged values of coefficients are obtained. A valid converged
solution of NPDE is established when an appropriate degree of B-poly basis is employed, and the
initial conditions are imposed on the operational matrix before the inverse is invoked. However,
the accuracy of the solution depends on the number of B-polys of a certain degree expressed in
multidimensional variables. Four examples of NPDE have been worked out to show the efficacy
and accuracy of the 2-dimensional B-poly technique. The estimated solutions of the examples are
32
compared with the known exact solutions, and an excellent agreement is found between them. In
calculating the solutions of the NPD equations, the currently employed technique provides a
higher-order precision compared to the finite difference method. The present technique could be
readily extended to solving complex partial differential equations in multivariable problems.
Keywords: Nonlinear partial differential equations, Partial differential equations, B-poly basis
set, Hyperbolic partial differential equations, Differential equations
Introduction
Using a B-polynomial basis set, one can solve very complicated 2D partial differential
equations, which could appear in the fields of physics, engineering, chemistry, and computer
science[51], [52]. The flawless integration and differentiation are the nature of B-polys that aid
in using symbolic programming languages such as Mathematica or Maple. Over any closed
interval, B-ploys are smooth functions that provide the basis to represent an arbitrary function to
desirable correctness[53]–[55]. The specific details and properties of the B-polys are provided in
our previous work in Refs. [35], [40], [41]. In the earlier years, various methods were employed
to solve linear and nonlinear differential equations, including fractional-order differential
equations[56]–[64]. In the articles ref. [53], [55], the authors have solved various differential
equations employing a concoction of the operational matrix and B-poly basis set. In the earlier
progression, the authors used the B-poly technique to calculate solutions of the single-variable
differential equations Ref. [35]. In our recently extended task, two-variable dependent
Hyperbolic Partial Differential (HPD) equations have been solved using the B-ploy bases[65]. In
the year 2011, the B-ploys were expressed in terms of the Legendre basis that has been used to
solve the linear differential equations [66].
33
In the previous work, the authors successfully applied a similar progression to the linear
partial differential equations and reported highly accurate solutions [67]. In the present work, we
aim to use an extended version of the technique for solving two variables (x, t) nonlinear partial
differential (NPD) equations employing the B-poly basis, operational matrix, and Galerkin
method. It is well known that the continuous and unitary properties of B-ploy help to determine
semi-analytic and, in some cases, exact solutions much quicker way in terms of CPU time [53].
The newly designed progression has effectively solved two variables' NPD equations in closed
intervals, such as [0, R] and [0, T]. A complete basis set of polynomials in two variables (x, t) in
terms of the product of B-poly sets have been utilized to figure out the solution of the NPDE.
The technique is explained step by step and applied to a more general 2-dimensional
NPDE,
𝛼𝑑2𝑦(𝑥,𝑡)
𝑑𝑥2+𝛽𝑦(𝑥,𝑡)𝑑𝑦(𝑥,𝑡)
𝑑𝑥 +𝛾𝑑𝑦(𝑥,𝑡)
𝑑𝑡 =𝑔(𝑥,𝑡).
(58)
Where 𝛼, 𝛽, 𝛾, could be constants or variables. A desired solution of the NPDE is
expressed as a linear combination of B-poly basis set as follows.
𝑦(x,t)=∑𝑎𝑖(𝑡)𝐵𝑖,𝑛(𝑥),
𝑛
𝑖=0
(59)
where 𝑎𝑖(𝑡) is the i-th expansion unknown coefficient in equation (59) that is a function of
variable t. In equation (59), we impose the initial conditions on variables (x, t). The 𝐵𝑖,𝑛(𝑥) is the
n-th degree B-poly in variable 𝑥 from the basis set. Furthermore, the coefficients 𝑎𝑖(𝑡) could be
expressed in terms of the constant coefficients 𝑏𝑗𝑖 and the B-polys 𝐵𝑗,𝑚(𝑡) in variable 𝑡 that could
have same or a different set of polynomials over the interval [0, T], such as,
34
𝑎𝑖(𝑡)=∑𝑏𝑗
𝑖
𝑚
𝑗=0 𝐵𝑗,𝑚(𝑡).
(60)
In equation (60), the coefficients 𝑎𝑖(𝑡) could be subjected to initial and boundary
conditions if required. We plan to present results by using a complete set of B-polys of the
different degree to a few nonlinear partial differential equations. In the following sections, the
procedure is applied to solve equations, and the results of the NPD equations are compared with
available 2D exact solutions. An excellent agreement has been found between exact and
estimated solutions. In the previous work, following the notation in References[53]–[55], [67],
the work has been protracted to include complete B-poly bases sets that were involved to
approximate results of a variety of differential equations, Refs. [67], [68].
The current technique is applied by substituting the approximate solution, Eq. (59), into
the NPDE (58) to separate out inner products in variables 𝑥 and𝑡. Both sides of the equations
are multiplied by a product of B-polys, 𝐵𝑚(𝑥)𝐵𝑛(𝑡), and integration over the closed intervals is
carried out. The inner products of B-polys are multiplied to form an operational matrix with a
nonzero determinant. Finally, the inverse of the operational matrix is carried out to determine the
unknown coefficients of the linear combination. The coveted solution of the NPDE is assembled
with the initial condition imposed on the operational matrix equation. In the following sections,
we shall explain the process of how to find an appropriate solution, present plots of the solutions,
and calculate semi-analytic solutions for each of the four examples considered in this paper.
Comparisons between exact and approximate solutions will be made in section 2, and finally, the
error analysis of the final example will be presented in section 3.
35
Computations of solutions of NPD equations
As mentioned in the previous section, details of the B-ploys will be left out to avoid
duplication of the formulas. Furthermore, information on how to generate sets of B-poly basis
and construct a solution from the sets are provided in refs.[32], [45], [69]. To make things
simpler, once the degree of the B-poly basis set is chosen, we may neglect subscript 𝑛
representing the degree in B-polys 𝐵𝑖,𝑛(𝑥) . The current technique to estimate solution, 𝑦(𝑥, 𝑡), of
the two-dimensional differential equation (58) is outlined in this section. The solution is
considered as a combination of the variables, 𝑎𝑖(𝑡)𝑎𝑛𝑑𝐵𝑖(𝑥). The NPDE is solved by adding an
initial condition in Eq. (60), employing the Galerkin method [70], and taking the inverse of the
operational matrix for calculating the unknown variables, 𝑎𝑖(𝑡), see refs. [32], [45], [69]. The
approximate solution with the initial condition (𝑦(𝑥, 0) = 𝑓(𝑥)) using the Galerkin method [70] is
given by,
𝑦(𝑥,𝑡)=∑𝑎𝑖(𝑡)𝐵𝑖(𝑥)
𝑛
𝑖=0 +𝑓(𝑥).
(61)
Implementing this solution, Eq. (61), into a nonlinear partial differential equation (58),
we may transform NPDE into an operational matrix by calculating the inner products of B-poly
in both variables x and t. Below, we have provided four examples of the NPD equations, which
are solved using the proposed technique and the initial condition 𝑦(𝑥,0)=𝑓(𝑥) at 𝑡=0.
Putting Eq. (61) into the second order NPD equation (58), we have an equation,
36
𝛼∑𝑎𝑖(𝑡)𝐵𝑖′′(𝑥)+𝛼𝑓′′(𝑥)
𝑛
𝑖=0
+𝛽(∑𝑎𝑗(𝑡)𝐵𝑗(𝑥)
𝑛
𝑗+𝑓(𝑥))(∑𝑎𝑖(𝑡)𝐵𝑖′(𝑥)
𝑛
𝑖+𝑓′(𝑥))
+𝛾∑𝑎𝑖(𝑡)𝐵𝑖(𝑥)
𝑛
𝑖=𝑔(𝑥,𝑡).
(62)
Where do (∙) and prime (′) denote derivatives with respect to t and x, respectively. The
equation (62) can be further simplified by moving some of the terms that do not depend on
unknown coefficients 𝑎𝑖(𝑡) to the right-hand side of the equation,
∑𝑎𝑖(𝑡)[𝛼𝐵𝑖′′(𝑥)
𝑛
𝑖=0 +𝛽𝐵𝑖′(𝑥)∑𝑎𝑗(𝑡)𝐵𝑗(𝑥)+𝛽𝑓(𝑥)𝐵𝑖′(𝑥)+𝛽𝑓′(𝑥)𝐵𝑖(𝑥)]
𝑛
𝑗=0
+𝛾∑𝑎𝑖(𝑡)𝐵𝑖(𝑥)=𝑔(𝑥,𝑡)−𝛽𝑓′(𝑥)𝑓(𝑥)
𝑛
𝑖=0 −𝛼𝑓′′(𝑥).
(63)
The expansion of 𝑎𝑖(𝑡)=∑𝑏𝑗𝑖𝑛
𝑗=0 𝐵𝑗(𝑡)𝑎𝑛𝑑𝑎𝑘(𝑡)=∑𝑏𝑙𝑘𝑛
𝑙=0 𝐵𝑙(𝑡) coefficients can be
used to convert Eq. (63) in terms of constants, 𝑏𝑗𝑖. After multiplying both sides of the Eq. (63)
with the product of B-polys 𝐵𝑚(𝑥)𝐵𝑛(𝑡) on both sides and integrating with respect to t and x
over the intervals 𝑡∈[0,𝑇] and 𝑥∈[0,𝑅], the Eq. (63) is transformed into Eq. (64),
37
∑∑𝑏𝑗𝑖[
𝑛
𝑖=0
𝑛
𝑗=0 〈𝛼𝐵𝑖′′(𝑥)+𝛽𝑓(𝑥)𝐵𝑖′(𝑥)+𝛽𝑓′(𝑥)𝐵𝑖(𝑥)|𝐵𝑚(𝑥)〉⟨𝐵𝑗(𝑡)|𝐵𝑛(𝑡)⟩
+𝛽∑𝑏𝑙𝑘𝑔𝑖𝑘𝑚
′
𝑛
𝑘,𝑙 𝑔𝑗𝑙𝑛+𝛾⟨𝐵𝑖(𝑥)|𝐵𝑚(𝑥)⟩⟨𝐵𝑗(𝑡)|𝐵𝑛(𝑡)⟩]
=⟨⟨𝑔(𝑥,𝑡)−𝛽𝑓′(𝑥)𝑓(𝑥)−𝛼𝑓′′(𝑥)|𝐵𝑚(𝑥)⟩|𝐵𝑛(𝑡)⟩.
(64)
The integrals and inner products of B-polys including the nonlinear terms are given
below:
⟨𝐵𝑚(𝑥)|𝐵𝑛(𝑥)⟩=∫ 𝐵𝑚(𝑥)𝐵𝑛(𝑥)𝑑𝑥
𝑅
0
⟨𝐵𝑚(𝑡)|𝐵𝑛(𝑡)⟩=∫ 𝐵𝑚(𝑡)𝐵𝑛(𝑡)𝑑𝑡
𝑇
0
⟨⟨𝐵𝑖(𝑥)|𝐵𝑚(𝑥)⟩|𝐵𝑛(𝑡)⟩=∬𝐵𝑖(𝑥)𝐵𝑚(𝑥)𝐵𝑛(𝑡)𝑑𝑥𝑑𝑡.
𝑅,𝑇
0
𝑔𝑖𝑘𝑚
′=∫𝐵𝑖′(𝑥)𝐵𝑘(𝑥)𝐵𝑚(𝑥)𝑑𝑥
𝑅
0
𝑔𝑗𝑙n=∫𝐵𝑗(𝑡)𝐵𝑙(𝑡)𝐵𝑛(𝑡)𝑑𝑡
𝑇
0
(65)
The terms 𝑔𝑖𝑘𝑚&𝑔𝑗𝑙𝑛are factors that exhibit nonlinear terms due to the appearance of
the coefficients 𝑏𝑙𝑘in front of them. Using the initial condition and including the nonlinear terms,
we can carry out calculations of the coefficients 𝑏𝑙𝑘. The converged coefficients are determined
after a few iterations of the nonlinear terms in the above Eq. (64). Sometimes, the initial guess to
the nonlinear terms may be given zero in the first iteration. We also employ the Galerkin method
[71] to minimize the error in the solution of the nonlinear partial differential equations. In this
method, the coefficients in Eq. (61) are minimized by increasing or decreasing to a certain
degree the number of polynomials in the process.
38
We may consider another type of general nonlinear PDE of the form in which a nonlinear
term appears at a different place in the equation,
𝛼𝑑2𝑦(𝑥,𝑡)
𝑑𝑥2+𝛽𝑑𝑦(𝑥,𝑡)
𝑑𝑥 +𝛾𝑦(𝑥,𝑡)𝑑𝑦(𝑥,𝑡)
𝑑𝑡 =𝑔(𝑥,𝑡).
(66)
Notice that now the nonlinear term appears in the third term as opposed to the second
term in Eq. (58). Using similar steps as described above, we can change the above equation (66)
into an operational matrix,
∑∑𝑏𝑗𝑖[
𝑛
𝑖=0
𝑛
𝑗=0 ⟨𝛼𝐵𝑖′′(𝑥)+𝛽𝐵𝑖′(𝑥)|𝐵𝑚(𝑥)⟩⟨𝐵𝑗(𝑡)|𝐵𝑛(𝑡)⟩
+𝛾⟨𝑓(𝑥)𝐵𝑖(𝑥)|𝐵𝑚(𝑥)⟩⟨𝐵𝑗(𝑡)|𝐵𝑛(𝑡)⟩+𝛾∑𝑏𝑙𝑘
𝑛
𝑘,𝑙=0 𝑔𝑖𝑘𝑚𝑔𝑗𝑙𝑛]
=⟨⟨𝑔(𝑥,𝑡)−𝛽𝑓′(𝑥)−𝛼𝑓′′(𝑥)|𝐵𝑚(𝑥)⟩|𝐵𝑛(𝑡)⟩.
(67)
Where, 𝑔𝑗𝑙𝑛&𝑔𝑖𝑘𝑚are the nonlinear terms which exhibit nonlinearity via the
coefficients 𝑏𝑙𝑘. The results of the Eq. (67) are further updated at each iteration when it is
subjected to the initial condition and initial guess for the nonlinear term. The error in the solution
of NPDE is minimized using the Galerkin method [71]. The current technique is applied to four
nonlinear differential examples. It is demonstrated that the method is suitable for finding an
accurate solution to NPDE.
First Example. Consider an NPDE which is obtained presuming parameters 𝛼=0,𝛽=𝑡, 𝛾=
−1𝑎𝑛𝑑𝑔(𝑥,𝑡)=𝑥 in the general Eq. (66). The NPDE in two dimensions is given,
𝑡𝑑𝑦
𝑑𝑥−𝑦𝑑𝑦
𝑑𝑡=𝑥.
(68)
39
The exact solution of Eq. (68) is 𝑦𝑒𝑥𝑎𝑐𝑡(𝑥,𝑡)=(𝑥−𝑡). For evaluating numerical
solutions using initial condition y(x,0)=f(x)=𝑥 in intervals0≤𝑡≤1and 0≤𝑥≤1,
an approximate solution to Eq. (68) may be assumed as 𝑦(𝑥,𝑡)=∑𝑎𝑖(𝑡)𝐵𝑖(𝑥)
𝑛
𝑖=0 +𝑥. After
substituting this expansion into Eq. (68), we get a similar equation as Eq. (67), which
restructures to the following equation,
∑𝑏𝑗𝑖[
𝑛
𝑖,𝑗=0 ⟨𝐵𝑖′(𝑥)|𝐵𝑚(𝑥)⟩⟨𝑡𝐵𝑗(𝑡)|𝐵𝑛(𝑡)⟩−⟨𝑥𝐵𝑖(𝑥)|𝐵𝑚(𝑥)⟩⟨𝐵𝑗(𝑡)|𝐵𝑛(𝑡)⟩
−∑𝑏𝑙𝑘𝑔𝑖𝑘𝑚𝑔𝑗𝑙𝑛
𝑛
𝑘,𝑙=0 ]=⟨⟨𝑥−𝑡|𝐵𝑚(𝑥)⟩|𝐵𝑛(𝑡)⟩.
(69)
Where, 𝑔𝑖𝑘𝑚&𝑔𝑗𝑙𝑛 are the nonlinear terms which exhibit nonlinearity via the unknown
coefficients 𝑏𝑙𝑘. This algorithm leads to an (𝑛+1) by (𝑛+1) system of equations𝑋𝐵=𝑊, in
unknown variables {𝑏1
1,𝑏2
1,𝑏3
1,…,𝑏1
2,𝑏2
2,𝑏3
2,…,}, which are elements of matrix B. The matrices
𝑋&𝑊 in terms of inner products of B-polys are provided below,
𝑋=⟨𝐵𝑖′(𝑥)|𝐵𝑚(𝑥)⟩⟨𝑡𝐵𝑗(𝑡)|𝐵𝑛(𝑡)⟩−⟨𝑥𝐵𝑖(𝑥)|𝐵𝑚(𝑥)⟩⟨𝐵𝑗(𝑡)|𝐵𝑛(𝑡)⟩
−∑𝑏𝑙𝑘𝑔𝑖𝑘𝑚𝑔𝑗𝑙𝑛
𝑛
𝑘,𝑙=0 ,
𝑊=⟨⟨𝑥−𝑡|𝐵𝑚(𝑥)⟩|𝐵𝑛(𝑡)⟩=∬(𝑥−𝑡)𝐵𝑚(𝑥)𝐵𝑛(𝑡)𝑑𝑥𝑑𝑡.
𝑅,𝑇
0
(70)
The nonlinear partial differential Eq. (70) is converted into an operational matrix form X,
whose inverse is multiplied by the column matrix W to yield values of the unknown coefficients
𝑏𝑗𝑖 by solving the equation 𝐵=𝑋−1𝑊. Before constructing the operational matrix equation,
𝑋𝐵=𝑊, initial conditions are imposed by deleting rows and the corresponding columns of Eq.
(70) so that the solution vanishes at t=0 and x=0. The resulting approximate solution is gathered
40
from the B-poly basis set via Eq. (60) and Eq. (61). After two iterations, the converged values of
the coefficients 𝑏𝑗𝑖 were attained, {0, -1, 0, -1}. With the application of this technique, the
solution of Eq. (68) is obtained, which is given for intervals 𝑡∈[0,1] and x ∈[0,1],
𝑦(𝑥,𝑡)=𝑡(−1.0+0.×10−30𝑥)+𝑥≈𝑥−𝑡.
(71)
As you may note, the above result is very accurate. To solve the NPD Equation (68), the
B-polys of degree n=1 have been utilized in both variables x and t. The B-poly basis sets used
are {1-t, t} and {1-x, x}, which gave a 4x4 operational matrix by multiplying both sets. We have
presented 3D plots of the exact and estimated results of Eq. (71) for comparison; see Fig. 5,
which shows the exact agreement between both solutions at the level of machine precision. Note
that when 𝑡=𝑥 is replaced in Eq. (71), the error can be seen in the order 10−17 representing the
high quality of the resolution in one-dimension x. The solution is essentially zero when t = x is
replaced. In this example, the absolute error between the solutions is negligible, showing that
both solutions are in perfect agreement.
Figure 5: The illustrations of the contrast between exact (sol) and approximate (fx) solutions
are presented on the left-hand for t = x replaced in the solutions of Eq. (68). This shows
complete overlap of the solutions. On the right-hand side, a 3D plot of the absolute error
between exact and estimated solutions is also shown in the intervals 𝑥∈[0,1] and 𝑡∈[0,1].
The graph shows the accuracy of the numerical results is of the order of 10−17. This kind of
accuracy occurred with an only n=1-degree polynomial basis set.
fx
sol
41
Second Example. Another case of the NPDE with a distinct nonlinear term is considered. We
replace 𝛼=0,𝛽=1,𝛾=1, 𝑔(𝑥,𝑡)=𝑥−𝑡−1in the general NPD equation (58). This
equation has a nonlinear term associated with the first term as opposed to the nonlinearity term
associated with the second term of the first example. We want to show that the current technique
can handle nonlinear terms associated to any term in the NPDE. The second example adds
another level of difficulty to be considered with nonhomogeneous terms on the right-hand side of
Eq. (72). The NPDE in two dimensions is given by,
𝑦𝑑𝑦
𝑑𝑥+𝑑𝑦
𝑑𝑡=𝑥−𝑡−1.
(72)
We are searching for a solution of Eq. (72) in the intervals 0≤𝑡≤1and, 0≤𝑥≤
1with initial condition y (x, 0) = 𝑓(𝑥)=𝑥 at t = 0. The exact solution to Eq. (72) is known,
𝑦𝑒𝑥𝑎𝑐𝑡(𝑥,𝑡)=(𝑥−𝑡). An approximate solution to Eq. (72) may be assumed to be 𝑦(𝑥,𝑡)=
∑𝑎𝑖(𝑡)𝐵𝑖(𝑥)
𝑛
𝑖=0 +𝑓(𝑥). In the assumed solution 𝑦(𝑥,𝑡), the coefficient 𝑎𝑖(𝑡) is the i-th
coefficient in the expansion which depends on t. By replacing the assumed solution into Eq. (14)
and multiplying both sides of Eq. (72) with the product of B-polys 𝐵𝑚(𝑥)𝐵𝑛(𝑡), we can
separately perform integration over both variables t and x in the intervals 𝑡∈[0,𝑇] and 𝑥∈
[0,𝑅], respectively. Applying an additional approximation to the coefficients 𝑎𝑖(𝑡)=
∑𝑏𝑗𝑖𝑛
𝑖=0 𝐵𝑗(𝑡), the Eq. (60) and with substitution of the term 𝑔(𝑥,𝑡)=𝑥−𝑡−1, we attain Eq.
(73) given below:
∑𝑏𝑗𝑖[
𝑛
𝑖,𝑗=0 ⟨𝑓(𝑥)𝐵𝑖′(𝑥)+𝑓′(𝑥)𝐵𝑖(𝑥)|𝐵𝑚(𝑥)⟩⟨𝐵𝑗(𝑡)|𝐵𝑛(𝑡)⟩+ ∑𝑏𝑙𝑘𝑔𝑖𝑘𝑚
′𝑔𝑗𝑙𝑛
𝑛
𝑘,𝑙=0
+⟨𝐵𝑖(𝑥)|𝐵𝑚(𝑥)⟩⟨𝐵𝑗(𝑡)|𝐵𝑛(𝑡)⟩]
=⟨⟨𝑥−𝑡−1−𝑓(𝑥)𝑓′(𝑥)|𝐵𝑚(𝑥)⟩|𝐵𝑛(𝑡)⟩.
(73)
42
Where, 𝑔𝑖𝑘𝑚
′&𝑔𝑗𝑙𝑛 are the nonlinear terms linked via the unknown coefficients 𝑏𝑙𝑘. The
Eq. (72) leads to an (𝑛+1) by (𝑛+1) system of equations𝑋𝐵=𝑊 , in the unknown variables
{𝑏1
1,𝑏2
1,𝑏3
1,…,𝑏1
2,𝑏2
2,𝑏3
2,…,},elements of matrix B, where the matrices 𝑋&𝑊 are given,
𝑋=⟨𝑓(𝑥)𝐵𝑖′(𝑥)+𝑓′(𝑥)𝐵𝑖(𝑥)|𝐵𝑚(𝑥)⟩⟨𝐵𝑗(𝑡)|𝐵𝑛(𝑡)⟩+ ∑𝑏𝑙𝑘𝑔𝑖𝑘𝑚
′𝑔𝑗𝑙𝑛
𝑛
𝑘,𝑙=0
+⟨𝐵𝑖(𝑥)|𝐵𝑚(𝑥)⟩⟨𝐵𝑗(𝑡)|𝐵𝑛(𝑡)⟩
𝑊=⟨⟨𝑥−𝑡−1−𝑓(𝑥)𝑓′(𝑥)|𝐵𝑚(𝑥)⟩|𝐵𝑛(𝑡)⟩.
(74)
The operational matrix X is provided in Eq. (74), whose inversion yields values of
unknown coefficients 𝑏𝑗𝑖 through solving the equation 𝐵=𝑋−1𝑊. Before the inverse of matrix
X is called for, initial conditions are imposed by deleting the appropriate rows and columns of
the matrices X and W to make the solution vanish at t=0 and x=0. Using a variational property
with respect to the coefficients and carrying out several iterations of the nonlinear terms, the
converged solution of the Eq. (72) is given,
𝑦(𝑥,𝑡)=𝑡(−1.0+0.×10−27)+𝑥≈𝑥−𝑡.
(75)
43
The values of the coefficients are listed as 𝑏𝑗𝑖={0,−1,0,−1} which are needed to
construct the numerical solution via Eq. (60). The expected solution, in the closed intervals, e.g.,
𝑡∈[0,1]𝑎𝑛𝑑𝑥∈[0,1], has the accuracy of the order of 10−27. Again, only n = 1-degree B-
polys were required for solving the differential Eq. (74) in both variables (x, t).
Figure 6: The illustrations of the contrast between precise (sol) and approximate (fx) results are
presented when t = x is replaced in the solutions of the Eq. (72), on the left-hand of Fig. This
graph shows accuracy in the single variable x. On the right-hand side of the Fig, 3D plot of the
absolute error between both solutions is also given over the intervals 𝑥∈[0,1] and 𝑡∈[0,1].
The graphs show the high precision of the numerical results. This kind of accuracy occurred
when n=1-degree polynomials were used in variables (x, t) for approximating the solution of Eq.
(72).
Third Example. Next, we consider a more complex second-order NPD equation with a distinct
nonlinear term. By replacing parameters 𝛽=−1,𝛼=𝜇=1, 𝛾=𝑥𝑡3, 𝑎𝑛𝑑𝑔(𝑥,𝑡)=2𝜇𝑡2 in
the general NPD equation (58), the NPD equation becomes,
𝜇𝑑2𝑦
𝑑𝑥2−𝑦𝑑𝑦
𝑑𝑥+𝑥𝑡3𝑑𝑦
𝑑𝑡=2𝜇𝑡2.
(76)
44
Again, we are searching for the solution in the intervals 0≤𝑥≤1and 0≤𝑡≤1.
Here we are going to use a slightly different initial condition, y (1, t) = f (t) = 𝑡2. The desired
solution for Eq. (76) can be assumed, 𝑦(𝑥,𝑡)=∑𝑎𝑖(𝑡)𝐵𝑖(𝑥)
𝑛
𝑖=0 +𝑡2. The unknown coefficient
𝑎𝑖(𝑡) is the i-th coefficient as a function of variable t in the evolution of expression, 𝑦(𝑥,𝑡). By
inserting this approximate solution into Eq. (76) and multiplying both sides of Eq. (76) with the
product of B-polys 𝐵𝑚(𝑥)𝐵𝑛(𝑡), we perform integration separately over both variables x and t in
the intervals 𝑡∈[0,𝑅] and 𝑥∈[0,𝑇], respectively. Additional approximation to the coefficients
𝑎𝑖(𝑡)=∑𝑏𝑗𝑖𝑛
𝑖=0 𝐵𝑗(𝑡)was also used to obtain the following equation,
∑𝑏𝑗𝑖[
𝑛
𝑖,𝑗=0 ⟨𝐵𝑖′′(𝑥)|𝐵𝑚(𝑥)⟩⟨𝐵𝑗(𝑡)|𝐵𝑛(𝑡)⟩−⟨𝐵𝑖′(𝑥)|𝐵𝑚(𝑥)⟩⟨𝑡2𝐵𝑗(𝑡)|𝐵𝑛(𝑡)⟩
−∑𝑏𝑙𝑘𝑔𝑖𝑘𝑚
′𝑔𝑗𝑙𝑛+⟨𝑥𝐵𝑖(𝑥)|𝐵𝑚(𝑥)⟩⟨𝑡3𝐵𝑗(𝑡)|𝐵𝑛(𝑡)⟩
𝑛
𝑘,𝑙=0 ]
=⟨⟨(2𝑡2−2𝑥𝑡4)|𝐵𝑚(𝑥)⟩|𝐵𝑛(𝑡)⟩.
(77)
Where, 𝑔𝑖𝑘𝑚
′&𝑔𝑗𝑙𝑛 are the nonlinear terms which exhibit nonlinearity via the unknown
coefficients 𝑏𝑙𝑘. This algorithm leads to an (𝑛+1)by (𝑛+1) order of matrix equation𝑋𝐵=
𝑊, where unknown coefficients are {𝑏1
1,𝑏2
1,𝑏3
1,…,𝑏1
2,𝑏2
2,𝑏3
2,…,}, that represent elements of
matrix B, and the matrices 𝑋𝑎𝑛𝑑𝑊 are given,
𝑋=⟨𝐵𝑖′′(𝑥)|𝐵𝑚(𝑥)⟩⟨𝐵𝑗(𝑡)|𝐵𝑛(𝑡)⟩−⟨𝐵𝑖′(𝑥)|𝐵𝑚(𝑥)⟩⟨𝑡2𝐵𝑗(𝑡)|𝐵𝑛(𝑡)⟩
−∑𝑏𝑙𝑘𝑔𝑖𝑘𝑚
′𝑔𝑗𝑙𝑛+⟨𝑥𝐵𝑖(𝑥)|𝐵𝑚(𝑥)⟩⟨𝑡3𝐵𝑗(𝑡)|𝐵𝑛(𝑡)⟩,
𝑛
𝑘,𝑙=0
𝑊=⟨⟨(2𝑡2−2𝑥𝑡4)|𝐵𝑚(𝑥)⟩|𝐵𝑛(𝑡)⟩=∬(2𝑡2−2𝑥𝑡4)𝐵𝑚(𝑥)𝐵𝑛(𝑡)𝑑𝑥𝑑𝑡.
𝑅,𝑇
0
(78)
45
The NPDE is transformed into an operational matrix equation𝑋𝐵=𝑊, where the
inversion of the matrix X yields values of the unknown coefficients 𝑏𝑗𝑖 in Eq. (77) by solving
equation 𝐵=𝑋−1𝑊. The approximate result is constructed from the product of the B-poly basis
set 𝐵𝑖(𝑡) and the coefficients 𝑎𝑖(𝑡). The converged values of the coefficients for this example
are,𝑏𝑗𝑖={0,0.×10−28,0.×10−27,0,0.×10−28,0.×10−27,0,0.×10−27,1.0}. The initial
condition is imposed on the matrices X and W defined in Eq. (78) by deleting the first row and
the corresponding first column because the solution must be zero at t=0 and x=0. This approach
delivered a valid solution of the Eq. (76), and the estimated result is provided,
𝑦(𝑥,𝑡)=𝑡(0.×10−27+0.×10−27𝑥+0.×10−27𝑥2)
+𝑡2(0.×10−27+0.×10−26𝑥+1.0𝑥2)≈𝑥2𝑡2.
(79)
The numerical result in Eq. (79) is compared with the exact solution𝑦𝑒𝑥𝑎𝑐𝑡(𝑥,𝑡)=𝑥2𝑡2.
Obviously, the approximate solution has high precision in the closed intervals 𝑡∈[0,1] and 𝑥∈
[0,1] that matched the exact solution when contributions of the small terms of the order
~10−27were discounted. This kind of accuracy shown in Eq. (79) was achieved with n = 2-
degree B-polys in both variables (x, t) which gave a 9x9 dimension operational matrix. To
observe the accuracy level of the solution in only variable x, we replaced t=x in the approximate
solution (fx) and the exact solution (sol) of Eq. (76), we essentially found overlays of graphs
shown on the left-hand side of Fig. 7. We have also shown 3D graphs of both the exact and the
estimated solutions on the right of Fig. 7. It is observed that there is no appreciable difference
between the graphs of both solutions as the error is very small and hence indicating the technique
works for calculating the solution of the NPDE.
46
Figure 7: Graphs of the approximate (f x) and the exact (sol) results are presented on the figure's
left for t = x, and both solutions overlap in one dimension. On the right, a 3D plot of exact and
approximate solutions is provided over the intervals t ∈ [0,1] and x ∈ [0,1], showing the
numerical results' accuracy because both graphs overlapped, showing no appreciable difference
between them.
Fourth Example. We shall now consider 4th example of the NPDE by substituting parameters
𝛼=1,𝛽=−1,𝛾=−1,𝑎𝑛𝑑𝑔(𝑥,𝑡)=− 3𝑥
(2𝑡+1)2 Into the general NPD equation (58). So, the
NPDE that we want to solve looks,
𝑑2𝑦(𝑥,𝑡)
𝑑𝑥2−𝑦(𝑥,𝑡)𝑑𝑦(𝑥,𝑡)
𝑑𝑥 −𝑑𝑦(𝑥,𝑡)
𝑑𝑡 =− 3𝑥
(2𝑡+1)2.
(80)
The exact solution of Eq. (80) is well-known, 𝑦𝑒𝑥𝑎𝑐𝑡 =3𝑥
(2𝑡+1). We are seeking a
solution in the closed intervals 0≤𝑡≤1and 0≤𝑥≤1 by applying initial condition at t = 0,
𝑓(𝑥)=𝑦(𝑥,0)=3𝑥, 𝑓′(𝑥)=3,𝑎𝑛𝑑𝑓′′(𝑥)=0. As we have done in the previous examples,
we would substitute the presumed solution Eq. (61) into the NPD Eq. (80) in order to convert it
into matrix form. The projected solution to Eq. (80) can be noted as 𝑦(𝑥,𝑡)=∑𝑎𝑖(𝑡)𝐵𝑖(𝑥)
𝑛
𝑖=0 +
3𝑥,where in this expansion𝑎𝑖(𝑡) is the i-th coefficient, which depends on variable t. Putting this
approximate solution into Eq. (80) and multiplying both sides of Eq. (80) with the product of B-
fx
sol
0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
47
polys 𝐵𝑚(𝑥)𝐵𝑛(𝑡), we can perform integration separately over both variables t and x in the
intervals 𝑡∈[0,𝑇] and 𝑥∈[0,𝑅], respectively. Additional approximation to the coefficients
𝑎𝑖(𝑡)=∑𝑏𝑗𝑖𝑛
𝑖=0 𝐵𝑗(𝑡) can also be applied to obtain the following equation,
∑𝑏𝑗𝑖[
𝑛
𝑖,𝑗=0 ⟨𝐵𝑖′′(𝑥)−3𝑥𝐵𝑖′(𝑥)−3𝐵𝑖(𝑥)|𝐵𝑚(𝑥)⟩⟨𝐵𝑗(𝑡)|𝐵𝑛(𝑡)⟩
−∑𝑏𝑙𝑘𝑔𝑖𝑘𝑚
′𝑔𝑗𝑙𝑛−⟨𝐵𝑖(𝑥)|𝐵𝑚(𝑥)⟩⟨𝐵𝑗(𝑡)|𝐵𝑛(𝑡)⟩
𝑛
𝑘,𝑙=0 ]
=⟨⟨ −3𝑥
(2𝑡+1)2+9𝑥|𝐵𝑚(𝑥)⟩|𝐵𝑛(𝑡)⟩.
(81)
Where, 𝑔𝑖𝑘𝑚
′&𝑔𝑗𝑙𝑛 are the nonlinear terms that exhibit nonlinearity via the unknown
coefficients 𝑏𝑙𝑘. This process leads to an (𝑛+1)by (𝑛+1) order of equation𝑋𝐵=𝑊, in
terms of the unknown variables 𝐵={𝑏1
1,𝑏2
1,𝑏3
1,…,𝑏1
2,𝑏2
2,𝑏3
2,…,}, where the operational matrix
𝑋 and the column matrix 𝑊 are,
𝑋=⟨𝐵𝑖′′(𝑥)−3𝑥𝐵𝑖′(𝑥)−3𝐵𝑖(𝑥)|𝐵𝑚(𝑥)⟩⟨𝐵𝑗(𝑡)|𝐵𝑛(𝑡)⟩
−∑𝑏𝑙𝑘𝑔𝑖𝑘𝑚
′𝑔𝑗𝑙𝑛−⟨𝐵𝑖(𝑥)|𝐵𝑚(𝑥)⟩⟨𝐵𝑗(𝑡)|𝐵𝑛(𝑡)⟩
𝑛
𝑘,𝑙=0 ,
𝑊=⟨⟨ −3𝑥
(2𝑡+1)2+9𝑥|𝐵𝑚(𝑥)⟩|𝐵𝑛(𝑡)⟩
=∬(−3𝑥
(2𝑡+1)2+9𝑥)𝐵𝑚(𝑥)𝐵𝑛(𝑡)𝑑𝑥𝑑𝑡.
𝑇,𝑅
0
(82)
It is tricky to choose the number of B-polys in both variable x and t for this example
because the same degree set of B-polys cannot be chosen in both variable x and t. We decided to
choose n=1 degree of B-polys in x-variable and 𝑛=13-degree of B-polys in t-variable. This
48
choice gave us a total of 28 B-polys basis set because n starts from 0. So, a large set of B-polys
were used to calculate the solution of Eq. (80). The composition of the exact solution shows that
there is 1-degree of polynomials present in x, {1-x, x} and a series of the polynomials present in
variable t. This is why we chose the approximate solution in this manner to minimize the error in
the solution. The result converged after 15 iterations to the desired accuracy. Converged values
of the 28 coefficients of the expansion for this example are also provided, 𝑏𝑗𝑖=
{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,−0.46153,−0.76925,−1.00687,−1.19203,
−1.34601,−1.47242,−1.58160,−1.67462,−1.75631,−1.82801,−1.89174,−1.94872,
−1.99999}, which are the elements of matrix B needed to build the desired solution via Eq. (82).
Again, the accuracy and the quality of the results depend on the number of B-polys and the
degree of polynomials used in the expansion of Eq. (82). On the left side of Figure 4, when t=x
is substituted in both exact and the approximate solutions, we plotted a graph for comparison in
linear dimension x. The graph shows the absolute error is of the order of ~10−8 between
solutions. As the number of B-poly increases from 2nd-degree through 13-degree polynomials in
variable t, the error between results (exact and approximate) was further decreased. This trend
will be presented in the error analysis section. On the right-hand side of Figure 8, we also present
a 3D plot of both solutions, which exhibits the accuracy of results in both variables x and t. The
converged approximate solution of example 4 is,
𝑦(𝑥,𝑡)=(3.0−5.99995𝑡+11.99753𝑡2−23.942006𝑡3+47.24151𝑡4−
89.71297𝑡5+156.15497𝑡6−234.51855𝑡7+286.60518𝑡8−270.34520𝑡9+
186.54281𝑡10−87.84899𝑡11+25.09446𝑡12−3.26879𝑡13)𝑥.
(83)
49
Graphs and excellent results of all examples indicate that the current technique works for
solving a variety of linear [72] and nonlinear partial differential equations.
Figure 8: On the left-hand side, a plot of the absolute error between exact and estimated solutions
are depicted for t = x. The graph shows that the error is of the order of 10-8 in one dimension. A
3D graph of the absolute error for exact and estimated solutions is also given on the right, which
shows the accuracy of the numerical results of the order of 10-8 over the intervals t ∈ [0,1] and x
∈ [0,1]. To achieve this level of accuracy in the solution, a total of 28 B-polys basis set was
involved, which is a product of 2 B-polys in variable x and 14 B-polys in variable t.
Error Analysis
It should be noted that the calculations based on B-polys are performed without a grid
representation in the intervals. Numerical results depend on the number of B-polys and the
degree of the polynomials chosen. In this section, we would like to show that accuracy starts to
improve as the number of B-poly increases. Both solutions, exact and approximate, begin to
overlap, and the precision is achieved after a certain number of iterations of the nonlinear terms.
We are going to present an error analysis for the approximate and exact solutions of the fourth
50
example. The same error analysis can be performed on other NPD equations. As you have seen
in example 4, only a set of two B-polys of degree n=1 was used in variable x because by
increasing the size of the B-poly basis set in variable x did not help improve the error. However,
by increasing the set of B-polys in variable t greatly helped improve the accuracy of the results
of the NPDE. We can observe from the graphs that the absolute error decreases between the
solutions as we steadily increase the set of B-polys basis in variable t. We are presenting two
graphs for n=4 and n=9 degrees set of B-poly in variable t, while n=1-degree of B-poly set was
unchanged. These graphs are depicted in Figs. 9 and 10 showing how the error steadily decreases
as we increase the number of B-polys from n=4 through n=13 degrees. The desired agreement
between exact solutions started to converge as n is increased for the inclusion of additional B-
polys in variable t, see Figs. 9-10. All the calculations are carried out without a grid on the
intervals of integration. When n= 4 degrees of B-polys are used and calculations are iterated 15
times, the absolute error is found to be of the order of 10-3 which is shown in Fig. 9. We have
also presented the 3D graph of the absolute error in terms of two variables (x, t) in Fig. 10. In
Fig. 10, n = 9 degrees of polynomials were used to show that the error further decreased to the
order 10-6. We have already shown a graph of the absolute error between the solutions, when
n=13 degrees of B-polys were chosen, see Fig. 8. Clearly, the error is systematically decreased as
the number of B-polys are increased in the calculations and each time results were iterated 15
times to achieve convergence of the solution. It indicates that the technique works to provide a
converged solution that is comparable with the exact solution. The CPU time for the calculation
significantly increased as we included a larger set of B-polys in the calculations.
51
Figure 9: A 3D plot of the absolute error between approximate and exact solutions of example 4
with n = 4-degree of B-polys is shown in two variables (x, t). The approximate solution has not
converged yet.
Figure 10: A 3D plot of the absolute error between approximate and exact solutions of example 4
with n = 9-degree of B-polys is shown in both variables (x, t). The approximate solution has not
converged as we need to increase the number of B-poly in variable t for better accuracy in the
solution.
Results and discussions
In the present work, we have applied 2-dimensional B-poly basis set technique to solve 4
examples of the NPD equations subjected to various initial conditions. Furthermore, A broader
52
explanation of the 2D algorithm has been given to calculate precise results of the NPDE. We
were approximating solutions employing the Galerkin method [71] in two variables (x, t) and the
converged results have been found which are presented in Figures 5–8. Normally, these solutions
converged after 10 to 15 iterations. For the first three examples, the approximate solutions were
so precise that the solutions matched the exact solutions after ignoring the tiny contributions and
a few iterations of nonlinear terms. With the increasing number of B-polys in the estimated
solution, we have also observed that the accuracy of the solutions[53] increases. In the first two
examples, we have used n = 1-degree set of polynomials in both variables, and in the third
example, n = 2-degree B-polys were used to obtain converged results. For the fourth example,
only 1-degree polynomials in variable x and 13-degree polynomials in variable t were used to
calculate the results. For example, 4, the B-poly basis set contained a total of 28 polynomials that
were used to calculate the solution. We have also displayed graphs of the absolute error between
the precise and the approximated results in Figs. 5-8. In every case, the accuracy of the solutions
has been different because different B-poly basis sets have been used. When variable t is set
equal to x, the absolute errors between exact and approximate solutions have been compared.
The accuracy was shown to be more significant for the converged solutions, as shown in Figs. 5-
8. Our method worked very well for solving nonlinear and linear differential equations using an
operational matrix scheme [71] as shown by the data and graphs presented in this paper.
Wolfram Mathematica symbolic program version-12[48] was used to perform all analytic
integrations and computations over the closed intervals for both variables (x, t).
Using B-poly basis sets, our method may display enormous possibilities for solving
nonlinear and linear 2D problems in physics and other disciplines. Recently, many authors [71],
[73]–[76] have formed an operational matrix utilizing B-polys techniques to solve one-
53
dimensional partial differential equations. We have successfully expanded this method to solve
the two-dimensional NPDE. We also presented a detailed error analysis for the last example,
showing how the error can be minimized as the number of B-polys increases in the desired
solution. To get a converged solution of the NPD Eq. (80), only 15 iterations were used. The
CPU time for executing examples 1-3 took only about 3 minutes, while for example 4, it took 10
minutes of CPU time as it involved higher degrees of B-polys and 15 iterations to converge.
In this article, we have shown an extended version of the technique[32] to solve NPD equations,
based on the B-poly basis set that provided in some cases exact and suitable solutions of 2D
nonlinear partial differential equations. This technique is superb for solving the problems related
to a complex system of nonlinear partial differential equations where there are no known
solutions that exist. As demonstrated in this paper, we may employ this particular method in
solving 2D linear and nonlinear partial differential equations.
54
CHAPTER V
TECHNIQUE TO SOLVE LINEAR FRACTIONAL DIFFERENTIAL EQUATIONS USING
B-POLYNOMIALS BASES
MUHAMMAD I. BHATTI*, MD. HABIBUR RAHMAN
University of Texas Rio Grande Valley, Edinburg Texas, 78539
*Corresponding author: [email protected]
Abstract
A multidimensional modified fractional-order B-polys technique has been implemented
for finding solutions to linear fractional-order partial differential equations. To calculate
the results of the linear Fractional Partial Differential Equations (FPDE), the sum of the
product of fractional B-polys and the coefficients is employed. Moreover, the minimization of
error in the coefficients is found by employing the Galerkin method. Before applying the
Galerkin method, the linear FPDE is transformed into an operational matrix equation that is
inverted to provide the values of the unknown coefficients in the approximate solution. A
valid multidimensional solution is determined when an appropriate number of basis sets and
fractional order of B-polys are chosen. In addition, initial conditions have to be applied to the
operational matrix to seek proper solutions in multidimensions.
The technique has been applied to 4 examples of linear FPDEs, and the agreements
between exact and approximate solutions are found to be excellent. The current technique can be
55
expanded into finding multidimensional fractional partial differential equations in other areas,
such as physics and engineering fields.
Keywords: Fractional B-Polynomials (B-Ploy), partial fractional differential equations,
Multidimensional formalism
Introduction
In real-world scientific phenomena, most problems follow either linearity or nonlinearity
in the systems. In different fields, for example, engineering, computer science, and
chemistry[77]–[79], fractional-order differential equations emerge more often. Most physical
phenomena are described by differential systems, which are integral-order systems. Many
systems could be expressed with the help of the fractional differential equation [80]–[85]. Due to
the real-world problems' materials, chemical properties, memory, and genetic characteristics,
many physical problems follow fractional dynamical behavior [85]–[87]. The partial fractional-
order differential equations are becoming a great tool for modeling physical phenomena [85]. For
this reason, there has been an urgent need to find a solution to the fractional problems. However,
there has been difficulty in obtaining accurate analytical or numerical results for most fractional-
order differential model equations. There is a need for a suitable technique for finding the
solutions to fractional-order differential problems, linear and nonlinear. Our current paper aims
to apply a technique to resolve multivariable linear fractional-order differential equations.
Nonlinear fractional-order partial differential equations will be considered in another future
work. In recent years many authors have used various numerical and analytical procedures to
unravel fractional-order differential equations, such as the modified simple equation approach
[88], [89], the variational iteration procedure [90], the Lagrange characteristic approach [91],
Adomian decomposition method [92], the finite difference procedure [93], the differential
56
transformation method [94], the finite element technique [95], the fractional sub equation
procedure [96], the (𝐺′/𝐺)-expansion method [97], first integral approach [98], and the fractional
complex transform technique[99]. Every method has its own pros and cons.
In this study, we are going to implement the modified fractional-order Bhatti polynomial
(B-poly) technique [53]–[55], [67], [100], [101] that is significantly capable of solving varieties
of multivariable linear fractional-order differential equations. We choose fractional-order B-poly
due to its well-defined basis set and precision [100]. With these basis sets, it can be demonstrated
that an arbitrary function can be represented to the desired accuracy and directly differentiable
over a closed interval. In several papers [53]–[55], [67], [100], [101] using the B-poly basis of
fractional-order and a generalized Galerkin method, the authors were able to find the solutions of
the fractional-order partial differential equations.
In the earlier work [53]–[55], [67], [100], [101], the authors have used a similar technique
to find the solution of ordinary nonlinear and linear multidimensional differential equations. The
current study focuses on linear fractional-order differential equations using the generalized
Galerkin method [53] and the B-poly basis of fractional order. This technique has the unique
advantage of the unitary partition property and the continuity of the generalized fractional-order
B-polys over an interval [0, R] which are differentiated seamlessly. With the help of fractional-
order B-polys, a fractional-order differential equation is transformed into an operational matrix
using matrix formalism that provides greater flexibility to apply boundary as well as initial
conditions on the operational matrix. The current study seeks solutions to four examples of linear
fractional-order partial differential equations using the fractional-order B-poly technique.
Employing Caputo’s fractional-order derivative definition, the derivatives of the fractional-order
B-polys are taken. The following sections present an analytical formalism to employ Caputo’s
57
fractional-order derivative on the polynomials. We present the process to create fractional-order
basis sets and develop an algorithm to resolve various linear fractional-order partial differential
equations. We apply this technique in four examples. Finally, we shall present an error analysis
on one of the fourth examples.
Caputo’s Fractional differential-order operator
The explanation of the fractional-order derivative of Caputo is provided as [79],
𝐷𝛾𝑓(𝑥)=𝐽𝑚−𝛾𝐷𝑚𝑓(𝑥)=1
𝛤(𝑚−𝛽)∫(𝑥−𝑡)𝑚−𝛾−1𝑓(𝑚)(𝑡)𝑑𝑡,
𝑥
0
𝑓𝑜𝑟𝑚−1<𝛾≤𝑚, 𝑚∈𝑁, 𝑥>0, 𝑓∈𝐶−1
𝑚.
(84)
Where 𝐷𝛾 are Caputo’s fractional operator, and fractional derivative in Caputo’s sense is
𝐷𝛾𝑓(𝑥), Eq. (84). Caputo’s derivative of a constant is zero, i.e., 𝐷𝛾𝐶=0 and a fractional
derivative of the polynomial 𝐷𝛾𝑥𝛼 is given by,
𝐷𝛾𝑥𝛼={0𝑓𝑜𝑟𝛼∈𝑁0𝑎𝑛𝑑𝛼<[𝛾]
𝛤(𝛼+1)𝑥𝛼−𝛾
𝛤(𝛼+1−𝛾)𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒.
(85)
Here α denotes the order of the fractional function. The unknown two-variable dependent
function 𝑈(𝑥,𝑡) is expanded as a product of two generalized fractional-order B-polynomials,
𝐵𝑗,𝑚(𝛼,𝑡)𝐵𝑖,𝑛(𝛼,𝑥), which may be considered as an approximate outcome of the FPD equation
represented by,
𝑈(𝑥,𝑡)= ∑𝑏𝑗𝑖𝐵𝑗,𝑚(𝛼,𝑡)𝐵𝑖,𝑛(𝛼,𝑥)
𝑛
𝑖,𝑗=0
(86)
58
Where, 𝐵𝑗,𝑚(𝛼,𝑡) is a j-th and m-degree fractional-order B-poly in variable 𝑡𝑜𝑟𝑥 with α
as a fractional-order parameter over an interval. The expansion coefficients 𝑏𝑗𝑖 in Eq. (86) are the
set of variables that are determined in the Galerkin scheme of minimization. Using Caputo’s
derivative property as a linear operator, we can perform fractional differentiation,
𝐷𝑥𝛾(∑𝑏𝑗𝑖𝐵𝑗,𝑚(𝛼,𝑡)𝐵𝑖,𝑛(𝛼,𝑥)
𝑛
𝑖,𝑗=0 )=∑𝑏𝑗𝑖(𝐵𝑗,𝑚(𝛼,𝑡)𝐷𝑥𝛾(𝐵𝑖,𝑛(𝛼,𝑥)))
𝑛
𝑖,𝑗=0
(87)
In the following section, we shall briefly mention the generalized fractional-order B-
Polys basis and some of their properties that could be useful to determine a solution to the linear
fractional-order partial differential equation.
Fractional-order B-Poly basis
The generalized form of fractional B-polys 𝐵𝑖,𝑛(𝛼,𝑥) in terms of variable x or t over an
interval [0, R] or [0, T] are defined in Refs. [53], [102],
𝐵𝑖,𝑛(𝛼,𝑥)=∑𝛽𝑖,𝑘
𝑛
𝑖=0 (𝑥
𝑅)𝛼𝑘.
The fractional-order parameter α represents the fractional degree of the B-poly. There are
(n + 1) fractional-order B-polynomials associated with any n value noted in the below equation.
The factor 𝛽𝑖,𝑘 in below equation is defined as,
𝛽𝑖,𝑘 =(−1)𝑖−𝑘(𝑛
𝑘)(𝑘𝑖).
where this binomial coefficient is defined as, (𝑛
𝑘)= 𝑛!
𝑘!(𝑛−𝑘)!. For convenience, if i < 0 or
i > n, we can set 𝐵𝑖,𝑘(𝛼,𝑥)=0. Mathematica or Maple software could be employed to develop
all the non-zero fractional polynomials using a simple code prewritten with any value of n
supported over an interval. The boundary conditions of the problem are generally associated with
59
the first and last polynomial of the basis set. As an example, when n=10 and fractional order
𝛼=1
2,5
3,9
4 The corresponding basis sets of B-polys are chosen in the above equation, plotted in
Figs. 11. Graphs of these fractional-order B-polys show how these B-polys add up to 1 at any
given point, x. Such type of sets of B-polys may be used to represent an arbitrary function at
higher accuracy.
(a) (b) (c)
Figure 11: (a) For n= 10, there is a total of 11 fractional B-polys of order 𝛼=1/2. (b) For n=
10, there is a total of 11 fractional B-polys of order 𝛼=5/3. (c) For n= 10, there is a total of 11
fractional B-polys of order 𝛼=9/4. The graphs of the set of 11 B-polynomials are presented in
the region x = [0, 10]. The 𝐵𝑖(1
2,𝑥), 𝐵𝑖(5
3,𝑥), 𝐵𝑖(9
4,𝑥), and 𝑥 are dimensionless quantities.
60
Table 1: For different values of 𝛼 (order of fractional-polynomials) and 𝛾 (order of the fractional
differential equation), the table below shows fractional polynomial basis sets with n = 1, gives
two B-polys and the corresponding derivatives. The symbol 𝛤 represents the Gamma function.
𝛼
𝛾
n
Basis set
Caputo’s Derivative
of Basis set (Eq. (2))
1/2
1/2
1
{1−√𝑥,√𝑥}
{−√𝜋
2,√𝜋
2}
3/4
3/4
1
{1−𝑥3 4
⁄,𝑥3 4
⁄}
{−Γ(7
4),Γ(7
4)}
5/3
5/3
1
{1−𝑥5 3
⁄,𝑥5 3
⁄}
{−Γ(8
3),Γ(8
3)}
5/4
5/4
1
{1 −𝑥5 4
⁄,𝑥5 4
⁄}
{−Γ(9
4),Γ(9
4)}
9/4
9/4
1
{1 −𝑥9 4
⁄,𝑥9 4
⁄}
{−Γ(13
4),Γ(13
4)}
9/5
9/5
1
{1−𝑥9 5
⁄,𝑥9 5
⁄}
{−Γ(14
5),Γ(14
5)}
Technique for approximating solutions
We exploit a technique to seek practical solutions to fractional-order partial differential
equations using the Galerkin method [47] and the generalized fractional-order B-poly basis set.
We intend to apply Caputo’s fractional derivative to the fractional-order B-ploys. Examples of
Caputo’s derivatives of the B-polys are provided in the last column of Table 1. Using the recent
technique, we transform the fractional-order linear partial differential equation into an
operational matrix, and the initial conditions and boundary conditions are applied to the
operational matrix. The presumed approximate solution, Eq. (86), is substituted into the
fractional-order differential equation, and by-products are separated in terms of integral products
61
in both variables x and t. Finally, both sides of the fractional equation are multiplied with
fractional B-polys basis elements, 𝐵𝑚(𝛼,𝑥)𝐵𝑛(𝛼,𝑡). The integrations are carried out using the
symbolic program Mathematica over the closed intervals [0, R] and [0, T]. For example, the
integration of the two fractional B-polys is given in the closed symbolic formula,
𝑚𝑖,𝑗 =(𝐵𝑖,𝑛(𝛼,𝑥),𝐵𝑗,𝑛(𝛼,𝑥))=∑𝛽𝑖,𝑘
𝑛
𝑘=𝑖 (𝑥
𝑅)𝛼𝑘∑𝛼𝑖,𝑘
𝑛
𝑙=𝑗 (𝑥
𝑅)𝛼𝑙 𝑅
(𝑘+𝑙)𝛼,
(88)
Caputo’s derivative defined in Eq. (85) is applied to the fractional B-ploy basis set, leading to the
following closed results,
𝐷𝑥𝛾(𝐵𝑖,𝑛(𝛼,𝑥))=∑𝛼𝑖,𝑘
𝑛
𝑘=𝑖 𝐷𝑥𝛾(𝑥
𝑅)𝛼𝑘 =∑𝛽𝑖,𝑘
𝑅𝛼𝑘
𝑛
𝑘=𝑖 𝛤(𝛼𝑘+1)
𝛤(𝛼𝑘+1−𝛾)𝑥𝛼𝑘−𝛾
𝑑𝑖,𝑗
(𝛾)(𝑥)=(𝐷𝑥𝛾𝐵𝑖,𝑛(𝛼,𝑥),𝐵𝑗,𝑛(𝛼,𝑥))=⟨𝐷𝑥𝛾𝐵𝑖,𝑛(𝛼,𝑥)|𝐵𝑗,𝑛(𝛼,𝑥)⟩
= ∑ 𝛽𝑖,𝑘
𝑛
𝑘=𝑖,𝑙=𝑗 𝛽𝑗,𝑘 𝛤(𝛼𝑘+1)
𝛤(𝛼𝑘+1−𝛾)𝑅1−𝛾
((𝑘+𝑙)𝛼+1−𝛾).
(89)
And the integrals of some arbitrary functions are given,
𝐹(𝑥,𝑡)=(𝑓(𝑥,𝑡),𝐵𝑖,𝑛(𝛼,𝑥))=∑𝛽𝑖,𝑘
𝑅𝛼𝑘
𝑛
𝑘=𝑖 ∫ 𝑓(𝑥,𝑡)𝑥𝛼𝑘𝑑𝑥,
𝑅
0
𝑊𝑚,𝑛 =∬𝑓(𝑥,𝑡)𝐵𝑚(𝛼,𝑥)𝐵𝑛(𝛼,𝑡)𝑑𝑥𝑑𝑡.
𝑅,𝑇
0
(90)
With the help of these analytic formulas, Eqs. (88-90), the operational matrix is
constructed. The inverse of the operational matrix is required to find out the unknown
coefficients 𝑏𝑗𝑖 of the linear combination in Eq. (86). In the next section, we will describe our
62
technique and how to obtain a desirable result of the linear fractional-order partial differential
equation. The technique will be employed in 4 examples to demonstrate that it works
appropriately for approximating the accurate solutions. We shall show how the inverse of the
operational matrix is calculated using the symbolic program Mathematica 12.0 [48]. Plots of the
approximate as well as exact solutions will be presented for the purpose of making comparisons.
Also, the absolute error analysis of the fourth example will be introduced to show that enlarging
the basis set of the fractional B-polys increases the accuracy of the solution. Similarly, the error
investigation can be carried out for other examples considered in this study. In the following
section, for simplicity, we would like to drop off subscripts 𝑛 and 𝑚 from the fractional B-polys,
𝐵𝑖,𝑛(𝛼,𝑥)=𝐵𝑖(𝛼,𝑥)and𝐵𝑗,𝑚(𝛼,𝑡)=𝐵𝑗(𝛼,𝑡).
Example 1: Let us introduce a linear partial fractional-order differential equation of the form,
2𝑑𝛾𝑈(𝑥,𝑡)
𝑑𝑡𝛾+𝑑𝛾𝑈(𝑥,𝑡)
𝑑𝑥𝛾=0.
(91)
The ideal solution of Eq. (91) is 𝑈𝑒𝑥𝑎𝑐𝑡(𝑥,𝑡)=(𝑥𝛾−𝑡𝛾/2). A numerical solution is
sought out in the intervals 0≤𝑥≤1&0≤𝑡≤1usinginitialconditionU(x,0)=f(x)=
𝑥𝛾. The assumed solution, 𝑈(𝑥,𝑡)=∑𝑏𝑗𝑖𝐵𝑗(𝛼,𝑡)𝐵𝑖(𝛼,𝑥)
𝑛
𝑖,𝑗=0 +𝑥𝛾, is substituted into the Eq.
(91) and the result is presented below,
2𝑑𝛾
𝑑𝑡𝛾(∑𝑏𝑗𝑖𝐵𝑗(𝛼,𝑡)𝐵𝑖(𝛼,𝑥)
𝑛
𝑖,𝑗=0 +𝑥𝛾)+𝑑𝛾
𝑑𝑥𝛾(∑𝑏𝑗𝑖𝐵𝑖(𝛼,𝑡)𝐵𝑖(𝛼,𝑥)
𝑛
𝑖,𝑗=0 +𝑥𝛾)=0
(92)
63
Caputo’s derivative operator is applied to Eq. (92). The product of fractional B-polys
𝐵𝑚(𝛼,𝑥) 𝐵𝑛(𝛼,𝑡) from the basis set is multiplied on both sides of the Eq. (92) and the
integration on both variables is calculated over the intervals using a symbolic program. This
operation provides the following equation.
∑𝑏𝑗𝑖[2⟨𝐵𝑖(𝛼,𝑥)|𝐵𝑚(𝛼,𝑥)⟩
𝑛
𝑖,𝑗=0 ⟨𝐷𝑡𝛾𝐵𝑗(𝛼,𝑡)|𝐵𝑛(𝛼,𝑡)⟩
−⟨𝐷𝑥𝛾𝐵𝑖(𝛼,𝑥)|𝐵𝑚(𝛼,𝑥)⟩⟨𝐵𝑗(𝛼,𝑡)|𝐵𝑛(𝛼,𝑡)⟩]
=⟨⟨−𝑓𝛾(𝛼,𝑥)|𝐵𝑚(𝛼,𝑥)⟩|𝐵𝑛(𝛼,𝑡)⟩.
(93)
Where 𝑓𝛾(𝑥)=𝐷𝑥𝛾(𝑥𝛾)=𝛤(𝛾+1)with 𝛼=𝛾. The current technique leads to a
system of (𝑛+1)×(𝑛+1) equations. The elements of matrix 𝐵=
{𝑏1
1,𝑏2
1,𝑏3
1,…,𝑏1
2,𝑏2
2,𝑏3
2,…,} are the unknown constants that are involved in those equations.
After further simplification, the right-hand side column matrix 𝑊 and the matrix elements of
operational matrix 𝑋 in terms of inner products of B-polys are given,
𝑋𝑚,𝑛 =∑[2⟨𝐵𝑖(𝛼,𝑥)|𝐵𝑚(𝛼,𝑥)⟩⟨𝐷𝑡𝛾𝐵𝑗(𝛼,𝑡)|𝐵𝑛(𝛼,𝑡)⟩
𝑛
𝑖,𝑗=0 −⟨𝐷𝑥𝛾𝐵𝑖(𝛼,𝑥)|𝐵𝑚(𝛼,𝑥)⟩⟨𝐵𝑗(𝛼,𝑡)|𝐵𝑛(𝛼,𝑡)⟩]
𝑊𝑚,𝑛 =⟨⟨−𝑓𝛾(𝑥)|𝐵𝑚(𝛼,𝑥)⟩|𝐵𝑛(𝛼,𝑡)⟩=−∬𝛤(𝛾+1)𝐵𝑚(𝛼,𝑥)𝐵𝑛(𝛼,𝑡)𝑑𝑥𝑑𝑡.
𝑅,𝑇
0
(94)
The partial fractional-order differential equation (91) is now transformed into a matrix
equation 𝑋𝐵=𝑊. By deleting rows and the corresponding columns of the equation (94), the
initial conditions are imposed on the operational matrix equation 𝑋 and the corresponding matrix
W so that the solution vanishes at t = 0 and x = 0. The operational matrix X has been coded in
64
the symbolic language Mathematica to determine its inverse. The inverse matrix was multiplied
by the column matrix W to yield values of the unknown coefficients 𝑏𝑗𝑖. The emerging estimated
result is composed of the linear combination of the B-poly basis set via Eq. (86). The process
provides a valid approximate solution 𝑈(𝑥,𝑡) of the Eq. (91) using B-polys of fractional-order
𝛼=1
2 and fractional differential order of 𝛾=1
2 in Eq. (91),
𝑈(𝑥,𝑡)=𝑥γ+𝑡γ(−0.5+0.×10−30𝑥𝛾)≈𝑥𝛾−𝑡𝛾/2.
(95)
From the above result, it is noted that the approximate solution is very accurate. We have
experimented with different values of fractional order 𝛾 of the differential equation while
keeping the same order 𝛾=𝛼 of the fractional polynomials basis set; the results remain the same
with various values of 𝛾𝑎𝑛𝑑𝛼. To solve the fractional-order partial differential equation (91),
we choose n = 1 and 𝛼=1
2 order B-poly basis set {1−√𝑡,√𝑡} and {1−√𝑥,√𝑥} in variables t
and x, respectively. The corresponding coefficient values of the constant we obtained are {0,
−20
81Γ(9
4), 0, −16
81Γ(9
4)}. The Caputo’s derivative of the fractional B-poly basis set is {−√𝜋
2,√𝜋
2}.
A 3D plot of the estimated and exact results of equation (91) is presented in Fig. 12 for
comparison, which presents an excellent agreement between both results at the level of machine
accuracy. Note that when t = x is substituted into Eq. (95), the absolute error can be observed in
the order of 10−17 exhibiting the great aspect of constancy in one-dimension x. In the example,
the absolute error between the exact and approximate results shows that both results are of
excellent reliability. The absolute error in the 3D graph is also presented on the right-hand side in
65
Fig. 12. The 3D graph displays that the absolute error in the converged solution is of the order of
10−17.
Figure 12: A 1D plot of the absolute error between approximate (fx) and exact (sol) solutions is
depicted on the left-hand for t = x changed in the solution, Eq. (95). The one-dimensional graph
shows the overlap of both results. On the right-hand side, a 3D plot of the absolute error between
approximate and exact results is also presented in the intervals 𝑡∈[0,1] and 𝑥∈[0,1]. The
figure represents the consistency of the numerical solution is of the order of 10−17. This kind of
accuracy occurred with only two fractional B-polynomials in the basis set.
Example 2: Consider another example of a fractional-order linear partial differential equation
with different initial conditions U(x,0)=f(x)=𝐸𝛼,1(𝑥𝛼),
2𝑑𝛾𝑈(𝑥,𝑡)
𝑑𝑡𝛾+𝑑𝛾𝑈(𝑥,𝑡)
𝑑𝑥𝛾=0.
(96)
The ideal solution to equation (96) is 𝑈𝑒𝑥𝑎𝑐𝑡(𝑥,𝑡)=𝐸𝛼,1(𝑥𝛼−𝑡𝛼/2). The function 𝐸𝛼,𝛽(𝑧), is
called the Mittag-Leffler function [103] that and is described as 𝐸𝛼,𝛽(𝑧)=∑𝑍𝑘
𝛤(𝑘𝛼+𝛽)
∞
𝑘=0 . In the
66
summation of the Mittag-Leffler function, we have only kept k =15 in the summation of terms.
So, the numerical solution's accuracy will most likely depend on the number of terms we would
keep in the summation of Mittag-Leffler function. According to Eq. (86), an estimated solution
of Eq. (96) using the initial condition may be assumed as 𝑈(𝑥,𝑡)=∑𝑎𝑖(𝛼,𝑡)𝐵𝑖(𝛼,𝑥)
𝑛
𝑖=0 +
𝐸𝛼,1(𝑥𝛼). After substituting this expression into the partial fractional differential Eq. (96), a
numerical solution is pursued in the intervals 0≤𝑥≤1 and 0≤𝑡≤1. The Galerkin method,
[53], [71], is also applied to the presumed solution to obtain,
2𝑑𝛾
𝑑𝑡𝛾(∑𝑏𝑗𝑖𝐵𝑗(𝛼,𝑡)𝐵𝑖(𝛼,𝑥)
𝑛
𝑖,𝑗=0 +𝐸𝛼,𝑙(𝑥𝛼))
+𝑑𝛾
𝑑𝑥𝛾(∑𝑏𝑗𝑖𝐵𝑗(𝛼,𝑡)𝐵𝑖(𝛼,𝑥)
𝑛
𝑖,𝑗=0 +𝐸𝛼,𝑙(𝑥𝛼))=0.
(97)
Caputo’s fractional derivative is applied to Eq. (97) as well as the product of fractional B-
polys 𝐵𝑚(𝛼,𝑥) 𝐵𝑛(𝛼,𝑡) from the basis set is multiplied on both sides of Eq. (97). The resulting
integration of both variables (t and x) is calculated over the intervals 0≤𝑥≤1 and 0≤𝑡≤1,
respectively. After further simplification of Eq. (97), we obtain,
∑𝑏𝑗𝑖[2⟨𝐵𝑖(𝛼,𝑥)|𝐵𝑚(𝛼,𝑥)⟩
𝑛
𝑖,𝑗=0 ⟨𝐷𝑡𝛾𝐵𝑗(𝛼,𝑡)|𝐵𝑛(𝛼,𝑡)⟩
−⟨𝐷𝑥𝛾𝐵𝑖(𝛼,𝑥)|𝐵𝑚(𝛼,𝑥)⟩⟨𝐵𝑗(𝛼,𝑡)|𝐵𝑛(𝛼,𝑡)⟩]
=⟨⟨−𝑓𝛾(𝛼,𝑥)|𝐵𝑚(𝛼,𝑥)⟩|𝐵𝑛(𝛼,𝑡)⟩.
(98)
Where the fractional-order derivative of the Mittag-Leffler function 𝑓𝛾(𝛼,𝑥)=
𝑑𝛾
𝑑𝑥𝛾(𝐸𝛼,1(𝑥𝛼))=𝐸𝛼,1(𝑥𝛼), with 𝛼=𝛾 is used. The current technique leads to a system of
equations of (𝑛+1) × (𝑛+1) equations. This system of equations may be summarized in the
67
matrix equation 𝑋𝐵=𝑊, where the elements of matrix 𝐵={𝑏1
1,𝑏2
1,𝑏3
1,…,𝑏1
2,𝑏2
2,𝑏3
2,…,} are
the unknown constants. The right-hand side column matrix elements of 𝑊 and the matrix
elements of operational matrix 𝑋 in terms of inner products of B-polys are given as,
𝑋𝑚,𝑛 =∑[2⟨𝐵𝑖(𝛼,𝑥)|𝐵𝑚(𝛼,𝑥)⟩⟨𝐷𝑡𝛾𝐵𝑗(𝛼,𝑡)|𝐵𝑛(𝛼,𝑡)⟩
𝑛
𝑖,𝑗=0 −⟨𝐷𝑥𝛾𝐵𝑖(𝛼,𝑥)|𝐵𝑚(𝛼,𝑥)⟩⟨𝐵𝑗(𝛼,𝑡)|𝐵𝑛(𝛼,𝑡)⟩],
𝑊𝑚,𝑛 =⟨⟨−𝑓𝛾(𝛾,𝑥)|𝐵𝑚(𝛼,𝑥)⟩|𝐵𝑛(𝛼,𝑡)⟩=∬𝑓𝛾(𝛾,𝑥)𝐵𝑚(𝛼,𝑥)𝐵𝑛(𝛼,𝑡)𝑑𝑥𝑑𝑡.
𝑅,𝑇
0
(99)
By deleting rows and the corresponding columns of the equation (99), the initial
condition has been imposed on the operational matrix 𝑋, to ensure the solution vanishes at x = 0
and t = 0. The operational matrix X has been inverted using Mathematica symbolic program and
multiplied with the column matrix W for solving the equation 𝐵=𝑋−1𝑊 to yield values of the
unknown coefficients 𝑏𝑗𝑖 .The emerging estimated result is composed of the B-poly basis set and
the coefficients via Eq. (86). The process provides approximate solution 𝑈(𝑥,𝑡) of the Eq. (96)
using B-polys of fractional-order 𝛼=1
2and fractional differential order of 𝛾=1
2 . The final
approximate solution is provided below:
68
𝑈(𝑥,𝑡)=1.0+𝑡11 2
⁄(−0.0000016542−0.0000019138√𝑥)+𝑡6(2.9981×10−7
+3.8261×10−7√𝑥)+𝑡5(0.0000081056+0.0000091828√𝑥
+0.000008138𝑥)+1.1283√𝑥+1.0𝑥+0.75225𝑥3 2
⁄+0.5𝑥2
+0.3009𝑥5 2
⁄+0.16666𝑥3+0.085972𝑥7 2
⁄+0.041666𝑥4
+0.019104𝑥9 2
⁄+0.0083333𝑥5+0.0034736𝑥11 2
⁄+0.0013888𝑥6
+0.0005344𝑥13 2
⁄+0.00019841𝑥7+0.000071253𝑥15 2
⁄
+𝑡9 2
⁄(−0.000037295−0.000042104√𝑥−0.000037314𝑥
−0.000028069𝑥3 2
⁄)+𝑡4(0.00016275+0.00018365√𝑥
+0.00016276𝑥+0.00012243𝑥3 2
⁄)+𝑡7 2
⁄(−0.00067165
−0.00075788√𝑥−0.00067165𝑥−0.00050525𝑥3 2
⁄
−0.00033582𝑥2)+𝑡3(0.0026041+0.0029385√𝑥+0.0026041𝑥
+0.001959𝑥3 2
⁄+0.001302𝑥2+0.0007836𝑥5 2
⁄+0.00043403𝑥3)
+𝑡5 2
⁄(−0.0094032−0.01061√𝑥−0.0094032𝑥−0.0070736𝑥3 2
⁄
−0.0047016𝑥2−0.0028294𝑥5 2
⁄−0.0015672𝑥3
−0.00080841𝑥7 2
⁄)+𝑡2(0.03125+0.035262√𝑥+0.03125𝑥
+0.023508𝑥3 2
⁄+0.015625𝑥2+0.0094032𝑥5 2
⁄+0.0052083𝑥3
+0.0026866𝑥7 2
⁄+0.001302𝑥4+0.00059703𝑥9 2
⁄
+0.00026041𝑥5)+𝑡3 2
⁄(−0.094032−0.1061√𝑥−0.094032𝑥
−0.070736𝑥3 2
⁄−0.047016𝑥2−0.028294𝑥5 2
⁄−0.015672𝑥3
−0.0080841𝑥7 2
⁄−0.003918𝑥4−0.0017964𝑥9 2
⁄−0.0007836𝑥5
−0.00032663𝑥11 2
⁄−0.0001306𝑥6−0.00005025𝑥13 2
⁄)+𝑡(0.25
+0.28209√𝑥+0.25𝑥+0.18806𝑥3 2
⁄+0.125𝑥2+0.075225𝑥5 2
⁄
+0.041666𝑥3+0.021493𝑥7 2
⁄+0.010416𝑥4+0.0047762𝑥9 2
⁄
+0.0020833𝑥5+0.0008684𝑥11 2
⁄+0.00034722𝑥6
+0.0001336𝑥13 2
⁄+0.000049603𝑥7)+√𝑡(−0.56419
−0.63662√𝑥−0.56419𝑥−0.42441𝑥3 2
⁄−0.28209𝑥2
−0.16976𝑥5 2
⁄−0.094032𝑥3−0.048504𝑥7 2
⁄−0.023508𝑥4
−0.010778𝑥9 2
⁄−0.0047016𝑥5−0.0019597𝑥11 2
⁄−0.0007836𝑥6
−0.0003015𝑥13 2
⁄−0.00011194𝑥7−0.0000402𝑥15 2
⁄)
(100)
69
To solve the fractional order partial differential equation (96), we have chosen n = 15 and
𝛼=1
2 order B-poly basis set in both variables t and x. The corresponding fractional-order B-poly
basis set is given in terms of variable x, {1−15√𝑥+105𝑥−455𝑥3 2
⁄+1365𝑥2−
3003𝑥5 2
⁄+5005𝑥3−6435𝑥7 2
⁄+6435𝑥4−5005𝑥9 2
⁄+3003𝑥5−1365𝑥11 2
⁄+455𝑥6−
105𝑥13 2
⁄+15𝑥7−𝑥15 2
⁄,15√𝑥−210𝑥+1365𝑥3 2
⁄−5460𝑥2+15015𝑥5 2
⁄−30030𝑥3+
45045𝑥7 2
⁄−51480𝑥4+45045𝑥9 2
⁄−30030𝑥5+15015𝑥11 2
⁄−5460𝑥6+1365𝑥13 2
⁄−
210𝑥7+15𝑥15 2
⁄,105𝑥−1365𝑥3 2
⁄+8190𝑥2−30030𝑥5 2
⁄+75075𝑥3−135135𝑥7 2
⁄+
180180𝑥4−180180𝑥9 2
⁄+135135𝑥5−75075𝑥11 2
⁄+30030𝑥6−8190𝑥13 2
⁄+1365𝑥7−
105𝑥15 2
⁄,455𝑥3 2
⁄−5460𝑥2+30030𝑥5 2
⁄−100100𝑥3+225225𝑥7 2
⁄−360360𝑥4+
420420𝑥9 2
⁄−360360𝑥5+225225𝑥11 2
⁄−100100𝑥6+30030𝑥13 2
⁄−5460𝑥7+
455𝑥15 2
⁄,1365𝑥2−15015𝑥5 2
⁄+75075𝑥3−225225𝑥7 2
⁄+450450𝑥4−630630𝑥9 2
⁄+
630630𝑥5−450450𝑥11 2
⁄+225225𝑥6−75075𝑥13 2
⁄+15015𝑥7−1365𝑥15 2
⁄,3003𝑥5 2
⁄−
30030𝑥3+135135𝑥7 2
⁄−360360𝑥4+630630𝑥9 2
⁄−756756𝑥5+630630𝑥11 2
⁄−
360360𝑥6+135135𝑥13 2
⁄−30030𝑥7+3003𝑥15 2
⁄,5005𝑥3−45045𝑥7 2
⁄+180180𝑥4−
420420𝑥9 2
⁄+630630𝑥5−630630𝑥11 2
⁄+420420𝑥6−180180𝑥13 2
⁄+45045𝑥7−
5005𝑥15 2
⁄,6435𝑥7 2
⁄−51480𝑥4+180180𝑥9 2
⁄−360360𝑥5+450450𝑥11 2
⁄−
360360𝑥6+180180𝑥13 2
⁄−51480𝑥7+6435𝑥15 2
⁄,6435𝑥4−45045𝑥9 2
⁄+135135𝑥5−
225225𝑥11 2
⁄+225225𝑥6−135135𝑥13 2
⁄+45045𝑥7−6435𝑥15 2
⁄,5005𝑥9 2
⁄−30030𝑥5+
75075𝑥11 2
⁄−100100𝑥6+75075𝑥13 2
⁄−30030𝑥7+5005𝑥15 2
⁄,3003𝑥5−15015𝑥11 2
⁄+
30030𝑥6−30030𝑥13 2
⁄+15015𝑥7−3003𝑥15 2
⁄,1365𝑥11 2
⁄−5460𝑥6+8190𝑥13 2
⁄−
5460𝑥7+1365𝑥15 2
⁄,455𝑥6−1365𝑥13 2
⁄+1365𝑥7−455𝑥15 2
⁄,105𝑥13 2
⁄−210𝑥7+
105𝑥15 2
⁄,15𝑥7−15𝑥15 2
⁄,𝑥15 2
⁄}
70
To get the B-polys basis set in variable t, we replace x = t. We have verified that as we
enlarge, the number of fractional-order B-polys sets, the accuracy of the numerical solution
increases. A 3D plot of the absolute error between the estimated and exact solution of equation
(96) is presented for comparison in Fig. 13, which illustrates the reliability between both results
at the level of 10−6. Note that when 𝑡=𝑥 is substituted in the solution Eq. (100), the absolute
error can be observed in the order of 10−6 in one-dimension. The absolute error between the
approximate and exact solutions is found to be in good agreement.
Figure 13: The absolute error plot between approximate (fx) and exact (sol) results is depicted on
the left-hand side for t = x (1D plot). This graph is obtained when x = t is substituted in Eq.
(100). On the right-hand side, a 3D plot of the absolute error between exact and estimated
solutions is also presented in the intervals 𝑥∈[0,1] and 𝑡∈[0,1]. Both graphs show the
efficiency of the numerical solutions is of the order of 10−6. This kind of accuracy observed
when n=16 and 𝛼=1
2 order polynomial basis set was used.
Example 3: Consider another partial fractional differential equation with an initial condition
depending on the generalized fractional-order sine function,
71
2𝑑𝛾𝑈(𝑥,𝑡)
𝑑𝑡𝛾+𝑑𝛾𝑈(𝑥,𝑡)
𝑑𝑥𝛾=0.
(101)
A numerical solution is pursued using initial condition 𝑈(𝑥,0)=𝑓(𝑥)=𝑠𝑖𝑛𝛾(𝑥𝛾) in
intervals 0≤𝑥≤1&0≤𝑡≤1. The generalized definitions of sine and cosine function [79]
are given,
𝑠𝑖𝑛𝛾(𝑥𝛾)=∑(−1)𝑘𝑥2𝑘
Γ(2𝑘𝛾+1)
∞
𝑘=0 𝑎𝑛𝑑𝑐𝑜𝑠𝛾(𝑥𝛾)=∑ (−1)𝑘𝑥(2𝑘+1)
Γ(2𝑘𝛾+𝛾+1)
∞
𝑘=0 .
In this example, the assumed approximate solution contains an initial condition with a
generalized sine function defined above. The efficiency of the numerical result is based on the
number of terms that are kept in the summation of the generalized function. For this example, we
kept k = 15 terms in the summation of the generalized sine function.
The exact solution of Eq. (101) is 𝑈𝑒𝑥𝑎𝑐𝑡(𝑥,𝑡)=𝑠𝑖𝑛𝛾(𝑥𝛾)𝑐𝑜𝑠𝛾(−𝑡𝛾2
⁄)+
𝑐𝑜𝑠𝛾(𝑥𝛾)𝑠𝑖𝑛𝛾(−𝑡𝛾2
⁄). The approximated solution 𝑈(𝑥,𝑡)=∑𝑏𝑗𝑖𝐵𝑗(𝛼,𝑡)𝐵𝑖(𝛼,𝑥)
𝑛
𝑖,𝑗=0 +
𝑠𝑖𝑛𝛾(𝑥𝛾) is substituted into the fractional-order differential Eq. (101) and the result is given
below:
2𝑑𝛾
𝑑𝑡𝛾(∑𝑏𝑗𝑖𝐵𝑗(𝛼,𝑡)𝐵𝑖(𝛼,𝑥)
𝑛
𝑖,𝑗=0 +𝑠𝑖𝑛𝛾(𝑥𝛾))
+𝑑𝛾
𝑑𝑥𝛾(∑𝑏𝑗𝑖𝐵𝑗(𝛼,𝑡)𝐵𝑖(𝛼,𝑥)
𝑛
𝑖,𝑗=0 +𝑠𝑖𝑛𝛾(𝑥𝛾))=0
(102)
The Caputo’s derivative is applied to the above expression and as well as the product of
fractional B-polys 𝐵𝑚(𝛼,𝑥) 𝐵𝑛(𝛼,𝑡) from the basis, sets are multiplied on both sides of Eq.
72
(102). The integration of both variables (x and t) is carried out over the intervals, respectively,
and after further simplification, Eq. (102) may be written in the form,
∑𝑏𝑗𝑖[2⟨𝐵𝑖(𝛼,𝑥)|𝐵𝑚(𝛼,𝑥)⟩
𝑛
𝑖,𝑗=0 ⟨𝐷𝑡𝛾𝐵𝑗(𝛼,𝑡)|𝐵𝑛(𝛼,𝑡)⟩
−⟨𝐷𝑥𝛾𝐵𝑖(𝛼,𝑥)|𝐵𝑚(𝛼,𝑥)⟩⟨𝐵𝑗(𝛼,𝑡)|𝐵𝑛(𝛼,𝑡)⟩]
=⟨⟨−𝑓𝛾(𝛼,𝑥)|𝐵𝑚(𝛼,𝑥)⟩|𝐵𝑛(𝛼,𝑡)⟩.
(103)
The fractional-order derivative of the sine function is 𝑓𝛾(𝛾,𝑥)=𝑑𝛾
𝑑𝑥𝛾(𝑠𝑖𝑛𝛾(𝑥𝛾))=
𝑐𝑜𝑠𝛾(𝑥𝛾)with 𝛼=𝛾. The technique leads to a system of (𝑛+1) × (𝑛+1) operational matrix.
This system of equations may be summarized in the matrix equation of the form 𝑋𝐵=𝑊,
where the elements of matrix 𝐵={𝑏1
1,𝑏2
1,𝑏3
1,…,𝑏1
2,𝑏2
2,𝑏3
2,…,} are the unknown constants.
After further simplification, the right-hand side column matrix 𝑊 and the matrix elements of
operational matrix 𝑋 in terms of inner products of B-polys are given as,
𝑋𝑚,𝑛 =∑[2⟨𝐵𝑖(𝛼,𝑥)|𝐵𝑚(𝛼,𝑥)⟩⟨𝐷𝑡𝛾𝐵𝑗(𝛼,𝑡)|𝐵𝑛(𝛼,𝑡)⟩
𝑛
𝑖,𝑗=0 −⟨𝐷𝑥𝛾𝐵𝑖(𝛼,𝑥)|𝐵𝑚(𝛼,𝑥)⟩⟨𝐵𝑗(𝛼,𝑡)|𝐵𝑛(𝛼,𝑡)⟩]
𝑊𝑚,𝑛 =⟨⟨−𝑓𝛾(𝛾,𝑥)|𝐵𝑚(𝛼,𝑥)⟩|𝐵𝑛(𝛼,𝑡)⟩=∬−cos𝛾(𝑥𝛾)𝐵𝑚(𝛼,𝑥)𝐵𝑛(𝛼,𝑡)𝑑𝑥𝑑𝑡.
𝑅,𝑇
0
(104)
To construct an appropriate solution to the Eq. (101), the partial fractional order differential
equation (101) is transformed into a matrix equation 𝑋𝐵=𝑊. By deleting rows and the
corresponding columns of the equation (104), the initial conditions are imposed on the
operational matrix equation, 𝑋 so that the solution vanishes at t = 0 and x = 0. The operational
matrix X is programmed in the symbolic language Mathematica to determine its inverse. The
73
inverse of the matrix X is multiplied by the column matrix W and by solving the matrix equation
𝐵=𝑋−1𝑊 values of the unknown coefficients 𝑏𝑗𝑖 are determined. The resulting approximate
solution is composed of the product of the B-poly basis set and the expansion coefficients using
Eq. (86). The technique provides approximate solution 𝑈(𝑥,𝑡) of the Eq. (101) using 16 number
of B-polys with n = 15, fractional-order 𝛼=1
2 and fractional differential-order of 𝛾=1
2 . The
approximate solution is provided here,
74
𝑈(𝑥,𝑡)=−2.06×10−8𝑡7+𝑡13 2
⁄(−1.2767×10−7+1.0904×10−7√𝑥)
+𝑡6(−1.7825×10−7+7.8386×10−7√𝑥)+𝑡11 2
⁄(0.000001411
+9.5536×10−7√𝑥−0.0000034747𝑥)+𝑡5(−3.8062×10−7
−0.0000075878√𝑥−0.0000040623𝑥+0.000012541𝑥3 2
⁄)
+𝑡9 2
⁄(−0.000037658+0.0000019308√𝑥+0.000030833𝑥
+0.000014011𝑥3 2
⁄)+𝑡4(−2.4231×10−7+0.00018537√𝑥
−0.0000074639𝑥−0.00010117𝑥3 2
⁄−0.000040623𝑥2)
+𝑡7 2
⁄(0.00067153+0.0000011199√𝑥−0.00067793𝑥
+0.00002317𝑥3 2
⁄+0.00027749𝑥2+0.00010088𝑥5 2
⁄)
+𝑡3(−4.3795×10−8−0.0029379√𝑥−0.0000038482𝑥
+0.0019773𝑥3 2
⁄−0.000059711𝑥2−0.00064749𝑥5 2
⁄
−0.00021666𝑥3+0.00045867𝑥7 2
⁄)+𝑡5 2
⁄(−0.0094032
+1.7793×10−7√𝑥+0.0094015𝑥+0.000010452𝑥3 2
⁄
−0.0047455𝑥2+0.00012975𝑥5 2
⁄+0.0012949𝑥3
+0.00040354𝑥7 2
⁄−0.00080267𝑥4)+𝑡2(0.035262√𝑥
−5.2405×10−7𝑥−0.023503𝑥3 2
⁄−0.000023089𝑥2
+0.0094911𝑥5 2
⁄−0.00023884𝑥3−0.0022199𝑥7 2
⁄
−0.00064998𝑥4+0.0012231𝑥9 2
⁄−0.00038577𝑥5)
+𝑡3 2
⁄(0.094032−0.094032𝑥+0.0000011862𝑥3 2
⁄+0.047007𝑥2
+0.000041811𝑥5 2
⁄−0.015818𝑥3+0.00037072𝑥7 2
⁄
+0.0032374𝑥4+0.00089676𝑥9 2
⁄−0.0016053𝑥5
+0.00048385𝑥11 2
⁄−0.000011479𝑥6−0.000011505𝑥13 2
⁄)
+𝑡(−0.28209√𝑥−1.4172×10−8𝑥+0.18806𝑥3 2
⁄
−0.0000020962𝑥2−0.075212𝑥5 2
⁄−0.000061572𝑥3
+0.021693𝑥7 2
⁄−0.00047769𝑥4−0.0039466𝑥9 2
⁄
−0.0010399𝑥5+0.001779𝑥11 2
⁄−0.00051436𝑥6
+0.000011743𝑥13 2
⁄+0.000011357𝑥7)+√𝑡(−0.56419
+0.56419𝑥+2.406×10−8𝑥3 2
⁄−0.28209𝑥2
+0.0000028469𝑥5 2
⁄+0.094015𝑥3+0.000071676𝑥7 2
⁄
−0.023727𝑥4+0.00049429𝑥9 2
⁄+0.0038849𝑥5
+0.00097829𝑥11 2
⁄−0.0016053𝑥6+0.00044663𝑥13 2
⁄
−0.0000098398𝑥7−0.0000092045𝑥15 2
⁄)+√𝑥(1.1283
−0.75225𝑥+0.3009𝑥2−0.085972𝑥3+0.019104𝑥4
−0.0034736𝑥5+0.0005344𝑥6−0.000071253𝑥7
+0.0000083828𝑥8).
(105)
75
From the above result, it is noted that the desired approximate solution converged and
reached the desired accuracy. To find the solution of fractional-order partial differential Eq.
(101), we used the same fractional-order B-poly basis set as example 2. With the higher number
(n) of fractional B-polys, the higher order of accuracy is achievable at the expense of computer
CPU time. A 3D plot of the estimated and exact results of Eq. (101) is presented in Fig. 14 for
the purpose of comparison. The plot shows an excellent agreement between both solutions at the
level of 10−7. Note that when t = x is substituted in Eq. (105), the absolute error can be
observed at the same level as 10−7exhibiting the significant aspect of constancy in one-
dimension x.
Figure 14: A 1D plot of the absolute error between approximate (fx) and exact (sol) solutions is
presented on the left-hand side when t = x is replaced in Eq. (105). The plot shows the desired
error is smaller. On the right-hand side, a 3D plot of the absolute error between approximate and
exact results is also presented in the intervals 𝑡∈[0,1] and 𝑥∈[0,1]. The figure represents the
efficacy of the numerical solutions is of the order of 10−7. This kind of accuracy occurred with
n=15 number of fractional order B-poly basis set in x variable, and the same set of basis set was
used in t variable.
76
It is further noted that from the traditional trigonometric identity, we know that.
𝑠𝑖𝑛(𝑥+𝑡)=sin(𝑥)cos(𝑡)+cos(𝑥)sin(𝑡),
(106)
However, in Ref. [103] the authors state this kind of trigonometry identity does not hold
in fractional calculus. In this example, we have computationally proven that the above identity is
no longer valid in fractional calculus, i.e.
𝑠𝑖𝑛𝛾(𝑥𝛾+(−𝑡𝛾
2))≠𝑠𝑖𝑛𝛾(𝑥𝛾)𝑐𝑜𝑠𝛾(−𝑡𝛾
2)+𝑐𝑜𝑠𝛾(𝑥𝛾)𝑠𝑖𝑛𝛾(−𝑡𝛾
2)
(107)
For further verification, we have plotted both sides of the identity Eq. (107), and they
seem to disagree. For example, for 𝛼=1
2𝑎𝑛𝑑𝛾=1
2 , we show the graphs of both sides of the
identity at x=t, 𝑠𝑖𝑛𝛾(𝑥𝛾/2) (blue curve) and𝑠𝑖𝑛𝛾(𝑥𝛾)𝑐𝑜𝑠𝛾(−𝑡𝛾
2)+𝑐𝑜𝑠𝛾(𝑥𝛾)𝑠𝑖𝑛𝛾(−𝑡𝛾
2)
(yellow curve). The graphs of both sides of the identity show that the blue and the yellow curves
do not match up. However, we know that when 𝛾 takes integral values, both curves overlap. We
tried different 𝛾 and n values of B-polys, and these curves still did not overlap. It is concluded
that certain traditional trigonometry identities may not be valid in fractional calculus.
77
Figure 15: Two graphs of the identity Eq. (107) are presented to show that both sides of the
identity do not agree for fractional values of 𝛾 . The left side of the identity is 𝑠𝑖𝑛𝛾(𝑥𝛾/2) and
the right side of the identity is 𝑠𝑖𝑛𝛾(𝑥𝛾)𝑐𝑜𝑠𝛾(−𝑥𝛾
2)+𝑐𝑜𝑠𝛾(𝑥𝛾)𝑠𝑖𝑛𝛾(−𝑥𝛾
2) at t = x. The values
for 𝛼=1
2 and 𝛾=1
2 are used. It is shown that the blue curve and the yellow curve do not agree
or overlap. Hence, in general, the identity is invalid when fractional calculus is considered.
Example 4: We consider a final example of the partial fractional-order differential equation with
an initial condition as the generalized fractional-order cosine function,
2𝑑𝛾𝑈(𝑥,𝑡)
𝑑𝑡𝛾+𝑑𝛾𝑈(𝑥,𝑡)
𝑑𝑥𝛾=0.
(108)
A numerical solution is sought out using initial condition 𝑈(𝑥,0)=𝑓(𝑥)=𝑐𝑜𝑠𝛾(𝑥𝛾)
[79], in the intervals 0≤𝑥≤1𝑎𝑛𝑑0≤𝑡≤1. The assumed approximate solution contains an
initial condition that has a generalized cosine function. The efficiency of the numerical result is
based on the number of terms that are kept in the summation of the cosine formula. In this
example, we kept k = 15 terms in the summation of the generalized cosine function.
78
The exact solution of Eq. (108) is 𝑈𝑒𝑥𝑎𝑐𝑡(𝑥,𝑡)=(𝑐𝑜𝑠𝛾(𝑥𝛾)𝑐𝑜𝑠𝛾(−𝑡𝛾2
⁄ )−
𝑠𝑖𝑛𝛾(𝑥𝛾)𝑠𝑖𝑛𝛾(−𝑡𝛾2
⁄)). The assumed solution, 𝑈(𝑥,𝑡)=∑𝑏𝑗𝑖𝐵𝑗(𝛼,𝑡)𝐵𝑖(𝛼,𝑥)
𝑛
𝑖,𝑗=0 +
𝑐𝑜𝑠𝛾(𝑥𝛾), is substituted into Eq. (108), and the result is presented below,
2𝑑𝛾
𝑑𝑡𝛾(∑𝑏𝑗𝑖𝐵𝑗(𝛼,𝑡)𝐵𝑖(𝛼,𝑥)
𝑛
𝑖,𝑗=0 +𝑐𝑜𝑠𝛾(𝑥𝛾))
+𝑑𝛾
𝑑𝑥𝛾(∑𝑏𝑗𝑖𝐵𝑗(𝛼,𝑡)𝐵𝑖(𝛼,𝑥)
𝑛
𝑖,𝑗=0 +𝑐𝑜𝑠𝛾(𝑥𝛾))=0.
(109)
After further simplification and applying Caputo’s derivative, we obtain,
∑𝑏𝑗𝑖[2⟨𝐵𝑖(𝛼,𝑥)|𝐵𝑚(𝛼,𝑥)⟩
𝑛
𝑖,𝑗=0 ⟨𝐷𝑡𝛾𝐵𝑗(𝛼,𝑡)|𝐵𝑛(𝛼,𝑡)⟩
−⟨𝐷𝑥𝛾𝐵𝑖(𝛼,𝑥)|𝐵𝑚(𝛼,𝑥)⟩⟨𝐵𝑗(𝛼,𝑡)|𝐵𝑛(𝛼,𝑡)⟩]
=⟨⟨−𝑓𝛾(𝛼,𝑥)|𝐵𝑚(𝛼,𝑥)⟩|𝐵𝑛(𝛼,𝑡)⟩.
(110)
Where 𝑓𝛾(𝛾,𝑥)=𝑑𝛾
𝑑𝑥𝛾(𝑐𝑜𝑠𝛾(𝑥𝛾))=−𝑠𝑖𝑛𝛾(𝑥𝛾)with 𝛼=𝛾. The current technique
leads to a system of (𝑛+1)× (𝑛+1) equations. This system of equations may be summarized
in the matrix equation of the form 𝑋𝐵=𝑊, where the elements of matrix 𝐵=
{𝑏1
1,𝑏2
1,𝑏3
1,…,𝑏1
2,𝑏2
2,𝑏3
2,…,} are the unknown constants. The matrix elements of the column
matrix 𝑊, and operational matrix 𝑋 in terms of inner products of B-polys are given as,
𝑋𝑚,𝑛 =∑[2⟨𝐵𝑖(𝛼,𝑥)|𝐵𝑚(𝛼,𝑥)⟩⟨𝐷𝑡𝛾𝐵𝑗(𝛼,𝑡)|𝐵𝑛(𝛼,𝑡)⟩
𝑛
𝑖,𝑗=0 −⟨𝐷𝑥𝛾𝐵𝑖(𝛼,𝑥)|𝐵𝑚(𝛼,𝑥)⟩⟨𝐵𝑗(𝛼,𝑡)|𝐵𝑛(𝛼,𝑡)⟩],
(111)
79
𝑊𝑚,𝑛=⟨⟨−𝑓𝛾(𝛾,𝑥)|𝐵𝑚(𝛼,𝑥)⟩|𝐵𝑛(𝛼,𝑡)⟩=∬sin𝛾(𝑥𝛾)𝐵𝑚(𝛼,𝑥)𝐵𝑛(𝛼,𝑡)𝑑𝑥𝑑𝑡.
𝑅,𝑇
0
The partial fractional-order differential equation (108) is now converted into an
operational matrix equation 𝑋𝐵=𝑊. By deleting rows and the corresponding columns of the
equation (111), the initial condition is imposed on the operational matrix equation 𝑋, so that the
result has the correct behavior at t=0 and x=0. The inverse of matrix X is multiplied by the
column matrix W to solve the matrix equation 𝐵=𝑋−1𝑊 to yield the values of the unknown
coefficients 𝑏𝑗𝑖. The resulting approximate solution is composed of the product of the expansion
coefficients and the B-poly basis set as given in Eq. (86). The technique provides an
approximate solution 𝑈(𝑥,𝑡) of Eq. (86) using n=15 B-polys of fractional-order 𝛼=1
2 and
fractional differential-order 𝛾=1
2 ,
𝑈(𝑥,𝑡)= 1.0−1.0𝑥+0.5𝑥2−0.166𝑥3+0.0417𝑥4−0.0083𝑥5+0.00139𝑥6−
0.000198𝑥7+0.0000248𝑥8−0.00000275𝑥9+2.75×10−7𝑥10+𝑡(−0.25+
0.25𝑥+8.1×10−8𝑥3 2
⁄−0.125𝑥2+0.00000574𝑥5 2
⁄+0.0416𝑥3+0.0001𝑥7 2
⁄−
0.0106𝑥4+0.000518𝑥9 2
⁄+0.00132𝑥5+0.00081𝑥11 2
⁄−0.00097𝑥6+
0.000317𝑥13 2
⁄−0.0000371𝑥7)+𝑡2(0.0312+1.52×10−8√𝑥−0.0312𝑥+
0.00000179𝑥3 2
⁄+0.0156𝑥2+0.0000437𝑥5 2
⁄−0.00534𝑥3+0.000291𝑥7 2
⁄+
0.00082𝑥4+0.000561𝑥9 2
⁄−0.000729𝑥5+0.000257𝑥11 2
⁄−0.0000324𝑥6)+
𝑡7 2
⁄(−5.31×10−8−0.000758√𝑥−0.00000312𝑥+0.000518𝑥3 2
⁄−
0.0000364𝑥2−0.000128𝑥5 2
⁄)+𝑡9 2
⁄(−1.87×10−7+0.0000432√𝑥−
0.00000404𝑥−0.0000178𝑥3 2
⁄)+𝑡11 2
⁄(−2.15×10−7−0.00000121√𝑥−
0.00000159𝑥)+𝑡6(2.45×10−7+3.59×10−7√𝑥−9.5×10−7𝑥)+𝑡15 2
⁄(1.83×
(112)
80
10−9√𝑥)+𝑡7(−2.52×10−8+3.24×10−8√𝑥)+𝑡13 2
⁄(−8.×10−8+2.06×
10−7√𝑥)+𝑡5(−0.0000083+9.9×10−7√𝑥+0.00000518𝑥)+𝑡4(0.000162+
8.5×10−7√𝑥−0.000167𝑥+0.0000132𝑥3 2
⁄+0.0000518𝑥2)+𝑡3(−0.0026+
2.24×10−7√𝑥+0.0026𝑥+0.0000091𝑥3 2
⁄−0.00133𝑥2+0.000085𝑥5 2
⁄+
0.000276𝑥3)+𝑡5 2
⁄(0.0106√𝑥−7.18×10−7𝑥−0.00707𝑥3 2
⁄−0.0000218𝑥2+
0.0029𝑥5 2
⁄−0.00017𝑥3−0.000514𝑥7 2
⁄−0.000368𝑥4+0.000503𝑥9 2
⁄)+
𝑡3 2
⁄(−0.106√𝑥−4.06×10−8𝑥+0.0707𝑥3 2
⁄−0.00000359𝑥2−0.0283𝑥5 2
⁄−
0.0000729𝑥3+0.0082𝑥7 2
⁄−0.000425𝑥4−0.00114𝑥9 2
⁄−0.000736𝑥5+
0.00091𝑥11 2
⁄−0.00031𝑥6+0.0000376𝑥13 2
⁄)+√𝑡(0.637√𝑥−0.424𝑥3 2
⁄−
1.22×10−7𝑥2+0.169𝑥5 2
⁄−0.00000718𝑥3−0.0485𝑥7 2
⁄−0.000109𝑥4+
0.011𝑥9 2
⁄−0.00051𝑥5−0.00124𝑥11 2
⁄−0.000736𝑥6+0.00084𝑥13 2
⁄−
0.000265𝑥7+0.00003𝑥15 2
⁄)
From the above result Eq. (112), it is noted that the approximate solution is converged
and accurate. We have experimented with various values of fractional-order γ of the differential
equation while keeping the same fractional-order polynomials basis set, the result remains the
same at the level of the desired accuracy. It is noted that when n = 6 set of B-polys is used, the
absolute error is 10-3, and when n = 15 set of B-polys is used, the absolute error reduces to 10-7.
It is concluded that with the increasing number of n sets of B-polys, higher order of accuracy is
attainable. A 3D plot of the estimated and exact results (112) is presented in Fig. 16 for the
purpose of comparison. Also, note that when t = x is substituted in the Eq. (112), the absolute
error in one dimension also goes to 10−7. The absolute error 3D graph is also presented in Fig.
16 showing error in the converged solution is of the order of 10−7.
81
Figure 16: A plot of the absolute error between approximate (fx) and exact (sol) solutions is
introduced on the left-hand for t = x, Eq. (108). The one-dimensional graph shows the overlap of
both results is pretty good. On the right-hand side, a 3D plot of the absolute error between
approximate and exact results is also presented in the intervals 𝑡∈[0,1] and 𝑥∈[0,1]. The
figure represents the efficacy of the numerical solution is of the order of 10−7.
From the traditional trigonometric rule, we know that the following is a valid identity,
𝑐𝑜𝑠(𝑥+𝑡)=cos(𝑥)cos(𝑡)−sin(𝑥)sin(𝑡).
(113)
However, in Ref. [103] the authors state that this trigonometry identity may not be true in
fractional calculus. In this example, we have computationally proven that the above identity is no
longer valid in fractional calculus, i.e.
𝑐𝑜𝑠𝛾(𝑥𝛾+(−𝑡𝛾
2))≠𝑐𝑜𝑠𝛾(𝑥𝛾)𝑐𝑜𝑠𝛾(−𝑡𝛾
2)−𝑠𝑖𝑛𝛾(𝑥𝛾)𝑠𝑖𝑛𝛾(−𝑡𝛾
2).
(114)
For further verification, we have plotted both sides of the identity Eq. (114), and they
seem to disagree. For example, for 𝛼=1
2𝑎𝑛𝑑𝛾=1
2 , we show the graphs of both sides of the
identity at x=t, 𝑐𝑜𝑠𝛾(𝑥𝛾/2) (blue curve) and (𝑐𝑜𝑠𝛾𝑥𝛾𝑐𝑜𝑠𝛾(−𝑡𝛾
2)−𝑠𝑖𝑛𝛾𝑥𝛾𝑠𝑖𝑛𝛾(−𝑡𝛾
2))
(yellow curve). The graphs of both sides of the identity show that the blue and the yellow curves
82
do not agree. However, we know that when 𝛾 takes integral values, both curves overlap. We tried
different n values for B-polys and fractional values of 𝛾, and these curves still did not overlap. It
is concluded that this trigonometry identity may not be valid in fractional calculus.
Figure 17: Two graphs of the identity are presented to show that both sides of the identity do not
match for fractional-order values of 𝛾 . The left side of the identity is 𝑐𝑜𝑠𝛾(𝑥𝛾/2) and the right
side of the identity is (𝑐𝑜𝑠𝛾(𝑥𝛾)𝑐𝑜𝑠𝛾(−𝑥𝛾
2)−𝑠𝑖𝑛𝛾(𝑥𝛾)𝑠𝑖𝑛𝛾(−𝑥𝛾
2)) at t = x. The value for 𝛼=
1
2 and 𝛾=1
2 are used. It is shown that the blue curve and the yellow curve do not agree. Hence, in
general, the identity does not hold true when fractional calculus is considered.
Error Analysis
We have performed the calculations in the absence of a grid for solving linear fractional
partial differential based on fractional B-polys. The fractional-order B-polys basis sets are
defined on the intervals 𝑥∈[0,1] and 𝑡∈[0,1]. Our approximated results are dependent on the
chosen (n) number of B-polys and the fractional-order modified Bhatti-polynomials. In this
section, we present error analysis based on the increasing number of B-poly basis sets, and it is
noted that the accuracy improves. The absolute error analysis for example 4 is presented for the
83
exact and approximate results. As you may have seen in example 4, in the final calculation, we
have used the number k=15 in the summation of the generalized formula for 𝑐𝑜𝑠𝛾(𝑥𝛾/2)=
∑(−1)𝑘𝑥(2𝑘+1)
Γ(2𝑘𝛾+𝛾+1)
𝑛
𝑘=0 and also used n=15 for the B-poly basis set in both the x and t variables. Here
we want to show that as we set n = 6, the B-poly basis set would have only seven B-polys in it.
We perform the calculations on example 4; it is observed that the absolute error among
solutions is of the order of 10-3. Next, we use n = 10, giving us 11 B-poly sets. The absolute error
among solutions reduces to the level of 10-6. Finally, we use n = 15, comprising 16 B-polys in
the basis set. It is observed the error reduces to 10-7. We note that n=15 leads to a 256 x 256-
dimensional operational matrix which is already a big matrix to invert. We had to increase the
accuracy of the program to handle this big matrix in the Mathematica symbolic program. Beyond
these limits, finding an accurate inversion of the matrix becomes problematic. Please note that
increasing the number of terms in the summation (k-values in the initial conditions) also helps
reduce error in the approximate solutions of the linear partial fractional differential equations.
We can observe from the graphs (Fig. 19, Fig. 18, and Fig. 16) that the absolute error decreases
as we steadily increase the size of the fractional B-poly basis set. The appropriate agreement
between exact solutions begins to converge as n is increased to include additional fractional B-
polys in both variables, see Figs. 12-19. Because of the analytic nature of the fractional B-polys,
all the calculations are carried out without a grid representation on the intervals of integration.
We have also presented the absolute error regarding 3D graphs in Figs. 12-19. Clearly, the error
is systematically decreased as the number of B-polys basis sets is increased in the calculations. It
also shows that the method works well to provide a converged solution that is comparable with
the exact solution. The CPU time for the calculation notably rises as we include a larger set of
fractional B-polys in the computations.
84
Figure 18: The absolute error between exact and approximate results of example 4 with n=6
basis set of B-polys is presented in both variables (x, t); please see the 3D graph. The estimated
result is not converging yet. The values for 𝛼=1
2𝑎𝑛𝑑𝛾=1
2 are used.
Figure 19: The absolute error analysis between exact and approximate results of example 4 with
n = 10 basis set of B-polys is presented in both variables (x, t); please see the 3D graph. The
estimated result is not converged yet. The values for 𝛼=1
2𝑎𝑛𝑑𝛾=1
2 are used.
85
Results and Discussions
In the current study, we have investigated the 2D modified fractional Bhatti-polys basis
set technique for determining solutions of the partial fractional differential equations. Four
examples of linear partial fractional differential equations have been presented with various
initial conditions, and their semi-analytic solutions are also provided. Furthermore, a great
explanation of the 2D fractional algorithm process has been given to calculate approximate
solutions of the linear fractional differential equation. We have estimated the results using the
Galerkin method [71] in both variables (x, t). The graphs of the converged solutions have been
provided in Figures 12-17. For the first example, the estimated solution was precise, which was
equivalent to the exact resolution after ignoring the tiny contributions. As the number of
fractional B-polynomial basis set in the approximate solutions Eq. (86) was increased, the
accuracy of the numerical solutions [53] increased. In our second, third, and fourth examples, we
have used a value of n = 15 for the basis set of fractional B-polys in two variables (x, t). We also
present 3D graphs of the precise and approximated results of the absolute error in Figs. 12-19. In
every case, the accuracy of the solutions has been different because different B-poly basis sets
and different sizes of the operational matrix have been used. Also, the numerical efficiency of
the inverted matrix depends on the size of the matrix. In the last three examples, we have used
the series representation of the generalized sine and cosine functions, and this requires many
terms in the summation to be included. When variable t is set equal to x for 1D error analysis,
the absolute errors among approximate and exact results have been examined. The precision
appears to be identical in 1D and 3D error analyses. It is concluded that the present technique
performed well in resolving linear fractional-order differential equations utilizing an operational
matrix scheme[71], [104] as exhibited by the graphs and data shown in the study. We have
86
analyzed all integrations and performed computations using Wolfram Mathematica symbolic
program version-12 [105] for both x and t variables over closed intervals.
The technique has presented great possibilities for solving linear multidimensional
fractional differential equation problems in chemistry, physics, genetics, and other related
disciplines. Nonlinear partial fractional differential equations will be investigated in another
paper. Many authors [70], [100] have recently constructed operational matrices using B-polys
methods to explain 1D partial differential equations. We have successfully expanded this
technique to solve the 2D linear fractional differential equations. In our study, we also have
shown detailed error investigation for the fourth problem that can be applied to other examples.
The CPU time for computing the first example was performed in less than 1 minute, while for
examples 2-4, it took 5-30 minutes of CPU time since those required a more extensive B-ploys
basis set and higher dimensions of the operational matrix.
This paper presents an expanded form of this technique [53] to determine solutions to
linear partial fractional differential problems using fractional-order basis sets. This technique
works well for resolving the equations connected to a complicated system of linear fractional-
order differential problems where no known solutions exist. In forthcoming publications, we may
explore this method in solving 2D nonlinear partial fractional-order differential equations.
87
CHAPTER VI
A METHOD TO SOLVE ONE-DIMENSIONAL NONLINEAR FRACTIONAL
DIFFERENTIAL EQUATION USING B-POLYNOMIALS
Md. Habibur Rahman, Muhammad I. Bhatti* and Nicholas Dimakis
University of Texas Rio Grande Valley, Edinburg, Texas, 78539
*Corresponding author: [email protected]
Abstract
This article applies the fractional Bhatti-Polynomial bases to solve one-
dimensional nonlinear fractional differential equations (NFDEs). We derive a semi-analytical
solution from a matrix equation using an operational matrix which is constructed from the
terms of the NFDE using Caputo's fractional derivative of fractional B-polynomials (B-polys).
The results obtained using the prescribed method agree well with the analytical and numerical
solutions presented by other authors. The legitimacy of this method is demonstrated by
using it to calculate the approximate solutions to four NFDEs. The estimated solutions to the
differential equations have also been compared with other known numerical and exact
solutions. It is also noted that for solving the NFDEs, the present method provides a higher
order of precision compared to the various finite difference methods. The current
technique could be effortlessly extended to solving complex linear, nonlinear, partial, and
fractional differential equations in multivariable problems.
Keywords: Fractional B-Polynomials, Fractional differential equations, Nonlinear
partial fractional differential equation, fractional B-polynomials in multiple variables
88
Introduction
The formulation of fractional calculus was started over 300 years ago. It can be traced
back to Leibniz's letter to L'Hôpital, in which he first discussed the meaning of the one-half order
derivative [106]. Although fractional calculus is as old as conventional calculus, it was not as
widely used in engineering and science at its conception. Due to rapid advancements in the fields
of mathematical physics, differential equations, interface chaos, and probability [107], [108], as
well as in other fields of science and engineering [109]–[114], fractional differential equations
have become a subject of interest and a rapidly growing area of research. It is also used to
describe a wide range of complex phenomena in different fields, such as anomalous diffusion,
systems identification, wave propagation, continuous-time random walk dynamical systems,
fractional electrical circuits, control theory, sub-diffusive systems, chaos synchronization, signal
processing, viscoelasticity, fluid flow, and more [115]–[121]. Seismic analysis, viscoelastic
materials, and viscous damping have all been successfully modeled in recent years using
fractional differential equations (FDEs) [111], [122]–[125]. The nonlinear oscillation of an
earthquake can be modeled using fractional derivatives, and a fluid-dynamic traffic model using
fractional derivatives can eliminate the deficiency caused by the assumption of continuous traffic
flow [107], [111], [125] as a result, developing robust methods for solving FDEs is essential.
Many fractional-order differential equations have unknown exact solutions; thus, various
numerical methods have been employed to provide approximate solutions. Unfortunately, each
method has its own set of limitations, and while no single method can solve every problem, most
techniques excel at solving specific problems.
There have been numerous approaches proposed to solve fractional differential equations;
the widespread ones are the one with operational method [126], [127], the Fourier transform
89
method [128] the iteration method [106], the iterative Laplace transform method (ILTM) [122],
the Bernoulli wavelet method [129], the Spectral method [130], and the Laplace transform
method [111], [125]. The approaches vary in their strengths and weaknesses, but from the variety
comes a new factor in solving FDEs, namely computation time, with some methods requiring a
significant amount of computational time to accomplish solutions to fractional-order differential
equations.
In this article, we present a novel technique known as the modified fractional Bhatti
Polynomials method. With this method, we have successfully solved linear and nonlinear
differential equations; see references [28], [33], [34], [131]. The method is effective, and the
results obtained thus far are encouraging and reliable. Four examples are provided in this article
to explain the dependability and efficacy of the method. The results of the method are also
compared with existing techniques, and excellent agreement has been found between the results;
in all the cases, the present semi-analytical results are superior in accuracy.
Caputo's Fractional differential operator
The fractional-order derivative of Caputo is explained as follows [111],
𝐷𝛾𝑓(𝑥)=𝐽𝑚−𝛾𝐷𝑚𝑓(𝑥)=1
𝛤(𝑚−𝛽)∫(𝑥−𝑡)𝑚−𝛾−1𝑓(𝑚)(𝑡)𝑑𝑡,
𝑥
0
𝑓𝑜𝑟𝑚−1<𝛾≤𝑚, 𝑐𝑜𝑛𝑡𝑖𝑛𝑜𝑢𝑠𝑤ℎ𝑒𝑟𝑒𝑚∈𝑁,𝑥>0,𝑓∈𝐶−1
𝑚,
(1)
𝐷𝛾 is Caputo's fractional derivative operator. The Caputo's fractional derivative for any
constant, C, is zero such that: 𝐷𝛾𝐶=0, and the fractional derivative 𝐷𝑥𝛾𝑥𝛼 is given by:
𝐷𝛾𝑥𝛼={0𝑓𝑜𝑟𝛼∈𝑁0𝑎𝑛𝑑𝛼<[𝛾]
𝛤(𝛼+1)𝑥𝛼−𝛾
𝛤(𝛼+1−𝛾)𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒.
(2)
90
The fractional order of the function is denoted by α, and the fractional order of the
derivative is given by 𝛾. The unknown function 𝑦(𝑥) is expanded in terms of the generalized
fractional-order B-polynomials 𝐵𝑖,𝑛(𝛼,𝑥), which can be regarded as an approximate solution to
the one-dimensional NFD equation:
𝑦(𝑥)=∑𝑏𝑖𝐵𝑖,𝑛(𝛼,𝑥)+𝑓(𝑥).
𝑛
𝑖=0
(3)
In variable x, 𝐵𝑖,𝑛(𝛼,𝑥) is ith fractional-order B-poly with α as a fractional-order
parameter and 𝑓(𝑥) is the initial condition imposed on the solution. The expansion coefficients
𝑏𝑖represent the expansion coefficients that are determined in the Galerkin scheme of
minimization in Eq. (3). Fractional calculus of differentiation can be accomplished using
Caputo's derivative as a linear operator:
𝐷𝑥𝛾(∑𝑏𝑖𝐵𝑖,𝑛(𝛼,𝑥)+𝑓(𝑥)
𝑛
𝑖=0 )=∑𝑏𝑖(𝐷𝑥𝛾(𝐵𝑖,𝑛(𝛼,𝑥)))
𝑛
𝑖=0 +𝐷𝑥𝛾(𝑓(𝑥)).
(4)
The generalized fractional-order B-poly basis and some of its properties that may be
useful in determining a solution to the nonlinear fractional-order differential equation are briefly
discussed in the following section.
Fractional-order B-Poly basis
The generalization of fractional B-polys 𝐵𝑖,𝑛(𝛼,𝑥) in terms of single variable x over the
interval [0, R] is defined in [33], [34],
𝐵𝑖,𝑛(𝛼,𝑥)=∑𝛽𝑖,𝑘
𝑛
𝑘=0 (𝑥
𝑅)𝛼𝑘.
(5)
91
The fractional-order parameter α represents the fractional B-polys and the Eq. (5)
provides an (n+1) fractional-order B-polynomial basis set. In Eq. (5), the factor 𝛽𝑖,𝑘 is defined as:
𝛽𝑖,𝑘 =(−1)𝑖−𝑘(𝑛
𝑘)(𝑘𝑖),
(6)
and the binomial coefficient is defined as: (𝑛
𝑘)= 𝑛!
𝑘!(𝑛−𝑘)!. Using a simple symbolic code
prewritten with any value of n supported over an interval [0, R], it is possible to produce a
fractional B-poly basis set. The boundary conditions are typically associated with the first and
the last polynomials in the basis set. For example, the fractional basis set for n =3 and 𝛼
=1
2,3
5,2
3,3
4, and 4
5, are given for various values of fractional order in Table 2.
92
Table 2: When 𝛾=𝛼=1
2,3
5,2
3,3
4, and 4
5, the corresponding basis set and their derivative for n = 3
are given in this table. The Gamma function of fractional order is represented by 𝛤[𝑎/𝑏].
𝛼
𝛾
𝑛
Basis Set (n+1)
Caputo’s Derivative of Basis set (Equation (4))
1
2
1
2
3
{1−3√𝑥+3𝑥
−𝑥3 2
⁄,3√𝑥−6𝑥
+3𝑥3 2
⁄,3𝑥
−3𝑥3 2
⁄,𝑥3 2
⁄}
{−3√𝜋
2+6√𝑥
√𝜋−3√𝜋𝑥
4,3√𝜋
2−12√𝑥
√𝜋+9√𝜋𝑥
4,6√𝑥
√𝜋
−9√𝜋𝑥
4,3√𝜋𝑥
4}
3
5
3
5
3
{1−3𝑥3 5
⁄
+3𝑥6 5
⁄
−𝑥9 5
⁄,3𝑥3 5
⁄
−6𝑥6 5
⁄
+3𝑥9 5
⁄,3𝑥6 5
⁄
−3𝑥9 5
⁄,𝑥9 5
⁄}
{−3Γ[8
5]+3𝑥3 5
⁄Γ[11
5]
Γ[8
5]−𝑥6 5
⁄Γ[14
5]
Γ[11
5],3Γ[8
5]−6𝑥3 5
⁄Γ[11
5]
Γ[8
5]
+3𝑥6 5
⁄Γ[14
5]
Γ[11
5],3𝑥3 5
⁄Γ[11
5]
Γ[8
5]
−3𝑥6 5
⁄Γ[14
5]
Γ[11
5],𝑥6 5
⁄Γ[14
5]
Γ[11
5]}
2
3
2
3
3
{1−3𝑥2 3
⁄
+3𝑥4 3
⁄
−𝑥2,3𝑥2 3
⁄
−6𝑥4 3
⁄
+3𝑥2,3𝑥4 3
⁄
−3𝑥2,𝑥2}
{−3Γ[5
3]−2𝑥4 3
⁄
Γ[7
3]+3𝑥2 3
⁄Γ[7
3]
Γ[5
3],3Γ[5
3]+6𝑥4 3
⁄
Γ[7
3]
−6𝑥2 3
⁄Γ[7
3]
Γ[5
3],−6𝑥4 3
⁄
Γ[7
3]+3𝑥2 3
⁄Γ[7
3]
Γ[5
3],2𝑥4 3
⁄
Γ[7
3]}
3
4
3
4
3
{1−3𝑥3 4
⁄
+3𝑥3 2
⁄
−𝑥9 4
⁄,3𝑥3 4
⁄
−6𝑥3 2
⁄
+3𝑥9 4
⁄,3𝑥3 2
⁄
−3𝑥9 4
⁄,𝑥9 4
⁄}
{9√𝜋𝑥3 4
⁄
4Γ[7
4]−3Γ[7
4]−4𝑥3 2
⁄Γ[13
4]
3√𝜋,−9√𝜋𝑥3 4
⁄
2Γ[7
4]+3Γ[7
4]
+4𝑥3 2
⁄Γ[13
4]
√𝜋,9√𝜋𝑥3 4
⁄
4Γ[7
4]
−4𝑥3 2
⁄Γ[13
4]
√𝜋,4𝑥3 2
⁄Γ[13
4]
3√𝜋}
4
5
4
5
3
{1−3𝑥4 5
⁄
+3𝑥8 5
⁄
−𝑥12 5
⁄,3𝑥4 5
⁄
−6𝑥8 5
⁄
+3𝑥12 5
⁄,3𝑥8 5
⁄
−3𝑥12 5
⁄,𝑥12 5
⁄}
{−3Γ[9
5]+3𝑥4 5
⁄Γ[13
5]
Γ[9
5]−𝑥8 5
⁄Γ[17
5]
Γ[13
5],3Γ[9
5]−6𝑥4 5
⁄Γ[13
5]
Γ[9
5]
+3𝑥8 5
⁄Γ[17
5]
Γ[13
5],3𝑥4 5
⁄Γ[13
5]
Γ[9
5]
−3𝑥8 5
⁄Γ[17
5]
Γ[13
5],𝑥8 5
⁄Γ[17
5]
Γ[13
5]}
93
Method for approximating solutions of one-dimensional NFDEs
Using the Galerkin method [34] and the generalized fractional-order B-poly basis set, we
exploit a method to seek practical solutions to nonlinear fractional-order differential equations
(NFDEs). Using the recently developed method [28], [33], [34], [38], [131]–[133], we transform
the fractional-order NFDE into an operational matrix with initial and boundary conditions
imposed on it. To construct the operational matrix, we substitute Eq. (3) into the given NFDE,
then Caputo's derivative operator is applied to the basis set used in the expansion in each term of
the NFDE, and both sides of the NFDE are multiplied by the elements of the fractional B-poly
basis set, 𝐵𝑖(𝛼,𝑥). Finally, the integrations are carried out using the symbolic program
Mathematica [134], [135] over the closed interval [0, R] into interaction matrix. For example, the
integration matrix over the closed interval of the two fractional-order B-polys is given in the
closed symbolic formula:
𝑚𝑖,𝑗 =(𝐵𝑖,𝑛(𝛼,𝑥),𝐵𝑗,𝑛(𝛼,𝑥))=∑𝛽𝑖,𝑘
𝑛
𝑘=𝑖 (𝑥
𝑅)𝛼𝑘∑𝛼𝑗,𝑘
𝑛
𝑙=𝑗 (𝑥
𝑅)𝛼𝑙 𝑅
(𝑘+𝑙)𝛼.
(7)
The Caputo's derivative defined in Eq. (2) is applied to the fractional B-poly basis set,
leading to the following closed results:
𝐷𝑥𝛾(𝐵𝑖,𝑛(𝛼,𝑥))=∑𝛼𝑖,𝑘
𝑛
𝑘=𝑖 𝐷𝑥𝛾(𝑥
𝑅)𝛼𝑘 =∑𝛽𝑖,𝑘
𝑅𝛼𝑘
𝑛
𝑘=𝑖 𝛤(𝛼𝑘+1)
𝛤(𝛼𝑘+1−𝛾)𝑥𝛼𝑘−𝛾,
𝑑𝑖,𝑙
(𝛾)(𝑥)=(𝐷𝑥𝛾𝐵𝑖,𝑛(𝛼,𝑥),𝐵𝑙,𝑛(𝛼,𝑥))=⟨𝐷𝑥𝛾𝐵𝑖,𝑛(𝛼,𝑥)|𝐵𝑙,𝑛(𝛼,𝑥)⟩
= ∑ 𝛽𝑖,𝑘
𝑛
𝑘=𝑖,𝑙=𝑗 𝛽𝑙,𝑘 𝛤(𝛼𝑘+1)
𝛤(𝛼𝑘+1−𝛾)𝑅1−𝛾
((𝑘+𝑙)𝛼+1−𝛾),
(8)
and the integrals of some arbitrary function are given by,
94
𝑊𝑚=∫ 𝑓(𝑥)𝐵𝑚(𝛼,𝑥)𝑑𝑥.
𝑅
0
(9)
With the help of the above analytic formulas, Eqs. (7-9), the operational matrix is
formulated. The inverse of the operational matrix is required to find out the unknown
coefficients 𝑏𝑖 of the linear combination in Eq. (3). In the next section, we will apply our method
and demonstrate how one can obtain a desirable solution to the nonlinear fractional-order
differential equation. The method will be applied to four examples to demonstrate that it works
appropriately for approximating the solutions with greater accuracy. We will also explain how
the inverse of the operational matrix is calculated using the symbolic program Mathematica 13.0
[134], [135]. Plots of the approximate and exact solutions will be presented for the purpose of
making comparisons. Also, the absolute error analysis of the fourth example will be elaborated to
show that by including larger basis set of the fractional B-polys and increasing the number of
iterations used to solve NFDEs, the accuracy of the solution is enhanced considerably. In the
following sections, for the sake of simplicity, we will drop the subscript 𝑛 from the fractional B-
poly basis, so that 𝐵𝑖,𝑛(𝛼,𝑥)=𝐵𝑖(𝛼,𝑥).
Example 1. Consider a one-dimensional nonlinear fractional-order differential equation,
𝐷𝛾𝑦(𝑥)−𝑦(𝑥)+2𝑦2(𝑥)=0,
(10)
where the value of 0≤𝛾≤1 with the initial condition 𝑦(0)=1/3. The exact solution of this
equation when 𝛾=1 is known to be: 𝑦𝑒𝑥𝑎𝑐𝑡(𝑥)=1/(2+𝑒−𝑥). Using fractional B-poly basis, a
solution may be approximated as 𝑦𝑎𝑝𝑝(𝑥)=∑𝑏𝑖𝐵𝑖(𝛼,𝑥)+𝑦0,
𝑛
𝑖=0 with the initial
condition,𝑦0=1/3. By substituting the approximate solution into the Eq. (10), we obtain,
95
𝑑𝛾
𝑑𝑥𝛾(∑𝑏𝑖𝐵𝑖(𝛼,𝑥)+𝑦0
𝑛
𝑖=0 )−(∑𝑏𝑖𝐵𝑖(𝛼,𝑥)+𝑦0)
𝑛
𝑖=0
+2(∑𝑏𝑗𝐵𝑗(𝛼,𝑥)+𝑦0
𝑛
𝑗=0 )(∑𝑏𝑖𝐵𝑖(𝛼,𝑥)+𝑦0
𝑛
𝑖=0 )=0.
(11)
We evaluate Eq. (11) by computing the Caputo fractional derivative, multiplying both
sides of the equation by the elements of the fractional B-polys basis set, 𝐵𝑙(𝛼,𝑥), and carrying
out the integration over the interval [0, R],
∑𝑏𝑖
𝑛
𝑖[⟨𝐷𝛾𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩−⟨𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩+4⟨𝑦0𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩
+2∑𝑏𝑗
𝑛
𝑗⟨𝐵𝑖(𝛼,𝑥)𝐵𝑗(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩]
=⟨(𝑦0−2𝑦0
2)𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩.
(12)
The above equation may be rewritten in the matrix form,
⇒𝐵[𝐴−𝐶+𝐷+𝐸]=𝑊,
(13)
where coefficients of the 4th term, E, exhibit nonlinearity via its coefficients. The matrices of Eq.
(13) are given below:
𝐴=⟨𝐷𝛾𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩=∫ 𝐷𝛾(𝐵𝑖(𝛼,𝑥))𝐵𝑙(𝛼,𝑥)𝑑𝑥,
𝑅
0
𝐶=⟨𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩=∫ 𝐵𝑖(𝛼,𝑥)𝐵𝑙(𝛼,𝑥)𝑑𝑥,
𝑅
0
𝐷=4⟨𝑦0𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩=4∫ 𝑦0𝐵𝑖(𝛼,𝑥)𝐵𝑙(𝛼,𝑥)𝑑𝑥,
𝑅
0
(14)
96
𝐸=2∑𝑏𝑗𝑔𝑖𝑗𝑙
𝑛
𝑗=2∑𝑏𝑗⟨𝐵𝑖(𝛼,𝑥)𝐵𝑗(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩
𝑛
𝑗
=2∑𝑏𝑗∫ 𝐵𝑖(𝛼,𝑥)𝐵𝑗(𝛼,𝑥)𝐵𝑙(𝛼,𝑥)𝑑𝑥
𝑅
0
𝑛
𝑗,
𝑊=⟨(𝑦0−2𝑦02)𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩=∫(𝑦0−2𝑦0
2)𝐵𝑖(𝛼,𝑥)𝐵𝑙(𝛼,𝑥)𝑑𝑥,
𝑅
0.
For the initial estimates of the coefficients, we ignore nonlinear terms in equation (13) to
calculate matrix B for unknown coefficients b and the equation is solved to obtain initial guess,
𝐵[𝐴−𝐶+𝐷]=𝑊.
(15)
The new estimate for the value of matrix𝐵can be calculated using equation (13) and
using initial estimate from equation (15). After a few more iterations, we revise our estimated
solution for comparison. We also solve the nonlinear fractional differential equation for different
fractional values of 𝛾 by repeating the same procedure. The graphs are plotted for various values
of 𝛾 alongside the exact solution to observe the deviation from the exact solution for 𝛾=1 integral
value, Fig. 2.
97
Figure 20 shows that when 𝛾 is equal to 1, the graphs of numerical convergent 𝑓(𝑥)
solution and exact (𝑠𝑜𝑙) solution are shown on the left side. The order of absolute error between
estimated 𝑓(𝑥) and accurate (𝑠𝑜𝑙) solutions is given on the right side which is of the order
Figure 20: The approximate solution 𝑓(𝑥) and the precise solution (sol) are shown in Fig. 1 for
the case 𝛾=1 in equation (11), demonstrating that both solutions overlap. In the picture on the
right, the absolute error between the exact and approximate solutions is displayed. The absolute
error is of the order of 10−9.
Figure 21: Various fractional values of 𝛾=1,4
5,3
4,2
3,1
2 are used in
equation (11) and the plots of approximate solutions are presented in this
figure. It is noted that all the graphs approximately intersect at one point
(x ≅ 1).
98
of10−9. Higher accuracy can be accomplished if the number of fractional B-poly sets is
increased. In the references [136], [137], the absolute error is 10−3 using their numerical
technique. As a result, our method produces a highly accurate solution.
Example 2. Let us consider the following nonlinear fractional differential equation:
𝐷𝛾𝑦(𝑥)=2𝑦(𝑥)−𝑦2(𝑥)+1.
(16)
Where the value of 0≤𝛾≤1 and the boundary conditions 𝑦(0)=𝑦(𝑅)=0 are
imposed on Eq. (16). The exact solution of this equation when 𝛾=1 is known: 𝑦𝑒𝑥𝑎𝑐𝑡(𝑥)=1+
√2tanh(√2𝑥+1
2log(√2−1
√2+1)). Using the fractional B-poly basis, a solution may be
approximated as 𝑦𝑎𝑝𝑝(𝑥)=∑𝑏𝑖𝐵𝑖(𝛼,𝑥)+𝑦0
𝑛
𝑖, with the initial condition, 𝑦0=0. By plugging
this approximate solution into Eq. (16), we get the following expression:
𝑑𝛾
𝑑𝑥𝛾(∑𝑏𝑖𝐵𝑖(𝛼,𝑥)+𝑦0
𝑛
𝑖)
=2(∑𝑏𝑖𝐵𝑖(𝛼,𝑥)+𝑦0)
𝑛
𝑖
−(∑𝑏𝑖𝐵𝑖(𝛼,𝑥)+𝑦0
𝑛
𝑖)(∑𝑏𝑗𝐵𝑗(𝛼,𝑥)+𝑦0
𝑛
𝑗)+1.
(17)
We evaluate Eq. (17) by applying the Caputo fractional derivative on the first term,
multiplying both sides of the equation by the elements of the fractional B-polys basis set,
𝐵𝑙(𝛼,𝑥), and then carrying out the integration over the interval [0, R] on both sides. We get,
99
∑𝑏𝑖
𝑛
𝑖[⟨𝐷𝛾𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩−2⟨𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩+2⟨𝑦0𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩
+∑𝑏𝑗
𝑛
𝑗⟨𝐵𝑖(𝛼,𝑥)𝐵𝑗(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩]
=⟨(2𝑦0−𝑦02+1)𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩.
(18)
The above equation can be represented in the matrix form:
𝐵[𝐴−𝐶+𝐷+𝐸]=𝑊,
(19)
with the elements of each matrix in Eq. (19) are given,
𝐴=⟨𝐷𝛾𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩=∫ 𝐷𝛾(𝐵𝑖(𝛼,𝑥))𝐵𝑙(𝛼,𝑥)𝑑𝑥
𝑅
0,
𝐶=2⟨𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩=2∫ 𝐵𝑖(𝛼,𝑥)𝐵𝑙(𝛼,𝑥)𝑑𝑥
𝑅
0,
𝐷=2⟨𝑦0𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩=2∫ 𝑦0𝐵𝑖(𝛼,𝑥)𝐵𝑙(𝛼,𝑥)𝑑𝑥,
𝑅
0
𝐸=∑𝑏𝑗
𝑛
𝑗𝑔𝑖𝑗𝑙 =𝑏𝑗⟨𝐵𝑖(𝛼,𝑥)𝐵𝑗(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩=𝑏𝑗∫ 𝐵𝑖(𝛼,𝑥)𝐵𝑗(𝛼,𝑥)𝐵𝑙(𝛼,𝑥)𝑑𝑥
𝑅
0,
𝑊=⟨(2𝑦0−𝑦02+1)𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩=∫ (2𝑦0−𝑦02+1)𝐵𝑖(𝛼,𝑥)𝐵𝑙(𝛼,𝑥)𝑑𝑥
𝑅
0.
(20)
For the initial guess of the coefficients, we can solve Eq. (21) by neglecting the nonlinear
term. We calculate matrix elements of column matrix 𝐵 as an initial guess by solving the
equation,
𝐵[𝐴−𝐶+𝐷]=𝑊.
(21)
We substitute the elements of 𝐵into nonlinear Eq. (19) to obtain the revised estimate of
the unknown coefficients. The process of iteration is repeated until a convergent solution is
100
found. We have used the same procedure to solve the nonlinear fractional differential equation
for various fractional values of 𝛾. The graphs of the solutions for several values of 𝛾 are shown
in Fig. 22. In each, the accuracy was desirable as shown in Fig. 23.
Figure 23: Nonlinear Equation (16) is solved using fractional B-polys to produce the plots of various
approximate solutions shown in the figure. To produce these graphs, various fractional-order values for
𝛾=1,4
5,3
4,2
3,1
2,3
5 are used.
Figure 22: A 1-D graph of the approximate solution f(x) and exact (sol) solution is displayed on
the left side, demonstrating how well the two solutions overlap for 𝛾=1. In the picture on the
right side, the absolute error between the precise and approximate solutions is displayed. The
graph shows that the nonlinear solutions is converged and has high accuracy of the solution in the
order of 10-10.
101
Higher-order accuracy can be accomplished if the number of fractional B-Polys and
iterations are increased. In the other references [136], [137], the absolute error is of the order of
10−3 and 10−6, but our method accomplishes higher-order computational accuracy. Figure 22
depicts solutions of example 2 with various fractional order 𝛾=1,4
5,3
4,2
3,1
2,3
5 .
Tables 3, 4, and 5 show that our method produces excellent results when compared with
approximate other methods for various fractional values[136]–[138] shown in these tables. The
accuracy of our method, as pointed out in the earlier works [28], [38], [139], was dependent on
the number of basis sets and the number of iterations to achieve a converged solution. We accept
convergent values after a few iterations are performed. Tables 3, 4, and 5 compare present results
and other results obtained from using various finite difference methods.
102
Table 3: For fractional order derivative 𝛾=1, the comparison of the fractional B-poly approach
to other methods [136]–[138] is shown, where N denotes the number of B-polys in a basis set.
𝛾=1
Our Results
x
N=8
N=10
N=12
N=15
𝑦ℎ𝑎𝑎𝑟[136
]
𝑦ℎ𝑝𝑚[138
]
𝑦𝑐𝑤𝑚
[137]
𝑦𝑒𝑥𝑎𝑐𝑡[138
]
0.
0
0
0
0
0
0
0
0
0
0.
1
0.11028
4
0.11029
5
0.11029
5
0.11029
5
0.110295
0.110294
0.110311
0.110295
0.
2
0.24196
5
0.24197
7
0.24197
7
0.24197
7
0.241977
0.241965
0.241995
0.241977
0.
3
0.39508
7
0.39510
5
0.39510
5
0.39510
5
0.395105
0.395106
0.395123
0.395105
0.
4
0.56779
6
0.56781
3
0.56781
2
0.56781
2
0.567813
0.568115
0.567829
0.567812
0.
5
0.75599
9
0.75601
5
0.75601
4
0.75601
4
0.756015
0.757564
0.756029
0.756014
0.
6
0.95354
6
0.95356
7
0.95356
6
0.95356
6
0.953567
0.958259
0.953576
0.953566
0.
7
1.15293
0
1.15295
0
1.15295
0
1.15295
0
1.152949
1.163459
1.152955
1.152949
0.
8
1.34635
0
1.34636
0
1.34636
0
1.34636
0
1.346364
1.365240
1.346365
1.346364
0.
9
1.52689
0
1.52691
0
1.52691
0
1.52691
0
1.526911
1.554960
1.526909
1.526911
1.
0
1.68948
0
1.68950
0
1.68950
0
1.68950
0
1.689499
1.723810
1.689494
1.689498
103
Table 4: For fractional order derivatives 𝛾=0.5 and 0.75, the comparison of the current
fractional B-poly approach to other methods [137], [138] is presented, where N denotes the
number of basis sets.
𝛾=0.5
𝛾=0.75
Our Results
Our Results
X
N=6
N=7
N=8
𝑦ℎ𝑝𝑚[1
38]
𝑦𝑐𝑤𝑚
[137]
N=6
N=7
N=8
𝑦ℎ𝑝𝑚[1
38]
𝑦𝑐𝑤𝑚
[137]
0.
0
0
0
0
0
0
0
0
0
0
0
0.
1
0.593
69
0.593
86
0.593
08
0.32173
0.59276
0.246
23
0.245
52
0.2453
4
0.21687
0.31073
0.
2
0.932
64
0.933
41
0.933
06
0.62967
0.93318
0.476
13
0.475
36
0.47490
0.42889
0.58431
0.
3
1.173
53
1.173
82
1.174
02
0.94094
1.17398
0.711
54
0.710
40
0.7099
1
0.65461
0.82217
0.
4
1.347
05
1.346
97
1.346
84
1.25074
1.34665
0.939
52
0.938
76
0.9384
0
0.89140
1.02497
0.
5
1.474
28
1.474
33
1.473
88
1.54944
1.47389
1.149
29
1.149
28
1.1488
9
1.13276
1.19862
0.
6
1.570
25
1.570
69
1.570
44
1.82546
1.57057
1.334
21
1.334
75
1.3342
1
1.37024
1.34915
0.
7
1.645
38
1.646
07
1.646
22
2.06652
1.64620
1.491
57
1.492
39
1.49187
1.59428
1.48145
0.
8
1.706
35
1.706
9
1.707
02
2.26063
1.70688
1.622
13
1.623
37
1.6230
4
1.79488
1.59924
0.
9
1.756
69
1.757
02
1.756
59
2.39684
1.75664
1.729
35
1.731
4
1.7311
8
1.96223
1.70530
1.
0
1.797
24
1.798
26
1.798
39
2.46600
1.79822
1.818
8
1.820
87
1.8209
3
2.08738
1.80176
104
Table 5: For fractional order derivatives 𝛾=0.6 and 0.8, the comparison of the current
fractional B-poly approach to other methods [136] is shown, where N denotes the number of B-
polys in the basis set.
𝛾=.6
𝛾=.8
x
N=6
N=7
N=8
𝑦ℎ𝑎𝑎𝑟[136
]
x
N=6
N=7
N=8
𝑦ℎ𝑎𝑎𝑟[136
]
0.0
0
0
0
0
0.0
0
0
0
0
0.1
0.41362
0.4135
3
0.4131
1
0.42697
0.1
0.2089
9
0.2079
6
0.2078
8
0.211942
0.2
0.72356
0.7230
4
0.7222
1
0.73068
0.2
0.4143
0
0.4133
9
0.4132
3
0.41000
0.3
0.98636
0.9857
4
0.9855
7
0.99053
0.3
0.6330
5
0.6315
5
0.6313
6
0.635149
0.4
1.2006
1.2000
1
1.2000
4
1.19393
0.4
0.8545
0
0.8529
9
0.8528
6
0.85002
0.5
1.37101
1.3701
0
1.3699
5
1.37789
0.5
1.0675
2
1.0664
4
1.0663
3
1.06670
0.6
1.50529
1.5033
9
1.5034
2
1.48518
0.6
1.2633
4
1.2625
3
1.2623
9
1.27900
0.7
1.61173
1.6080
8
1.6088
6
1.62822
0.7
1.4363
1.4353
9
1.4352
5
1.43387
0.8
1.6978
1.6915
3
1.6931
4
1.85574
0.8
1.5839
8
1.5829
8
1.5829
1.58714
0.9
1.7693
1.7591
3
1.7612
5
2.16946
0.9
1.707
1.7063
9
1.7063
2
1.70696
1.0
1.8298
1.8134
2
1.8171
9
2.73017
1.0
1.8087
4
1.8081
4
1.8080
8
1.80913
105
Example 3. Consider another nonlinear fractional differential equation of the form,
𝐷𝛾𝑦(𝑥)=−𝑦2(𝑥)+1,
(22)
where the value of 0≤𝛾≤1 and the boundary conditions 𝑦(0)=𝑦(𝑅)=0 are given. The
exact solution of this equation (22) when 𝛾=1 is 𝑦𝑒𝑥𝑎𝑐𝑡(𝑥)=(𝑒2𝑥−1)/(𝑒2𝑥+1). An
approximate solution may be written as 𝑦𝑎𝑝𝑝(𝑥)=∑𝑏𝑖𝐵𝑖(𝛼,𝑥)+𝑦0
𝑛
𝑖, where initial condition,
𝑦0=0 is given. Plugging in the approximate solution into the Eq. (22), we get,
𝑑𝛾
𝑑𝑥𝛾(∑𝑏𝑖𝐵𝑖(𝛼,𝑥)+𝑦0
𝑛
𝑖)=−(∑𝑏𝑖𝐵𝑖(𝛼,𝑥)+𝑦0
𝑛
𝑖)(∑𝑏𝑗𝐵𝑗(𝛼,𝑥)+𝑦0
𝑛
𝑗)+1.
(23)
Evaluating Eq. (23) by computing the Caputo’s fractional derivative in the first term,
multiplying both sides of the equation by the elements of the fractional B- basis set, 𝐵𝑙(𝛼,𝑥), and
integrating over the interval [0, R] on both sides of the Eq. (23), we obtain,
∑𝑏𝑖
𝑛
𝑖[⟨𝐷𝛾𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩+2⟨𝑦0𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩
+∑𝑏𝑗
𝑛
𝑗⟨𝐵𝑖(𝛼,𝑥)𝐵𝑗(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩]=⟨(1−𝑦02)𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩.
(24)
The following is a representation of the above equation in matrix form,
𝐵[𝐴+𝐷+𝐸]=𝑊,
(25)
where elements of each matrix are given as follows:
𝐴=⟨𝐷𝛾𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩=∫ 𝐷𝛾(𝐵𝑖(𝛼,𝑥))𝐵𝑙(𝛼,𝑥)𝑑𝑥,
𝑅
0
(26)
106
𝐷=2⟨𝑦0𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩=2∫ 𝑦0𝐵𝑖(𝛼,𝑥)𝐵𝑙(𝛼,𝑥)𝑑𝑥
𝑅
0,
𝐸=∑𝑏𝑗
𝑛
𝑗𝑔𝑗𝑘𝑙 =𝑏𝑗⟨𝐵𝑖(𝛼,𝑥)𝐵𝑗(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩=𝑏𝑗∫ 𝐵𝑖(𝛼,𝑥)𝐵𝑗(𝛼,𝑥)𝐵𝑙(𝛼,𝑥)𝑑𝑥.
𝑅
0
𝑊=⟨(1−𝑦02)𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩=∫ (1−𝑦02)𝐵𝑖(𝛼,𝑥)𝐵𝑙(𝛼,𝑥)𝑑𝑥
𝑅
0.
By ignoring the nonlinear term, we may solve the equation for the initial guess such that.
𝐵[𝐴+𝐷]=𝑊.
(27)
The initial guess of Eq. (27) is substituted for the nonlinear term to obtain a new guess for matrix
B by solving equation (25). This process is repeated until convergent values of the coefficients of
matrix B are obtained. We have solved the nonlinear fractional differential equation Eq. (25) for
different fractional-order values of 𝛾 using the same procedure. The graphs are also plotted of
various solutions for values of 𝛾 in Figs. 24 and 25. The order of absolute error between
estimated 𝑓(𝑥) and accurate (𝑠𝑜𝑙) solutions on the right side are 10−9. Higher accuracy can be
accomplished if the number of fractional B-Poly basis sets and iterations increases. In the
reference [136], the absolute error is of the order 10−3 for the same problem. Our method
performs exceptionally well in terms of computation accuracy.
107
Tables 6, 7, and 8 show that the fractional B-poly method produces excellent results when
compared to other numerical methods for various fractional-order values [136], [138]. The
Figure 24: For fractional derivative 𝛾=1 in equation (22), a 1-D graphic of our approximate
solution, f(x), and the precise solution, (sol), are displayed on the left side showing that the two
solutions overlap. The image on the right shows a plot of the absolute error between the
approximate and precise solutions on the order of 10-9.
Figure 25: The fractional differential equation (22), which considers
various fractional orders γ=1,4
5,3
4,2
3,1
2, is shown in the picture
along with a plot of several approximations. All the graphs intersect
at one point (x≅1).
108
accuracy of our method [28], [38] was determined by the number of B-polys used in the basis set
and the number of iterations used. The convergent values are obtained after a few iterations.
Table 6: For fractional order 𝛾=1, the comparison of the fractional B-poly approach to other
methods [136], [138] has been presented, where N denotes the number of B-polys in the basis
set.
𝛾=1
Our Results
x
N=6
N=7
N=8
N=9
𝑦ℎ𝑎𝑎𝑟[136]
𝑦ℎ𝑝𝑚[138]
𝑦𝑒𝑥𝑎𝑐𝑡[138]
0.0
0
0
0
0
0
0
0
0.1
0.099652
0.099669
0.099668
0.099668
0.099668
0.099668
0.099668
0.2
0.197366
0.197376
0.197376
0.197375
0.197375
0.197375
0.197375
0.3
0.291303
0.291314
0.291313
0.291313
0.291313
0.291312
0.291313
0.4
0.379935
0.379950
0.379949
0.379949
0.379949
0.379944
0.379949
0.5
0.462105
0.462118
0.462117
0.462117
0.462117
0.462078
0.462117
0.6
0.537042
0.537050
0.537050
0.537050
0.537050
0.536857
0.537050
0.7
0.604362
0.604369
0.604368
0.604368
0.604368
0.603631
0.604368
0.8
0.664028
0.664038
0.664037
0.664037
0.664037
0.661706
0.664037
0.9
0.716290
0.716298
0.716298
0.716298
0.716298
0.709919
0.716298
1.0
0.761589
0.761595
0.761594
0.761594
0.761594
0.746032
0.761594
109
Table 7: For fractional order 𝛾=0.75, the comparison of the fractional B-poly method to other
numerical methods [138] is presented, where N denotes the number of B-polys in the basis set.
𝛾=0.75
Our Results
X
N=6
N=7
N=8
N=9
𝑦ℎ𝑝𝑚[138]
0.0
0
0
0
0
0
0.1
0.190088
0.190101
0.190101
0.190097
0.184795
0.2
0.309960
0.309973
0.309976
0.309960
0.313795
0.3
0.404581
0.404612
0.404614
0.404612
0.414562
0.4
0.481611
0.481632
0.481632
0.481630
0.492889
0.5
0.545088
0.545090
0.545090
0.545081
0.462117
0.6
0.597783
0.597781
0.597783
0.597777
0.597393
0.7
0.641807
0.641819
0.641820
0.641821
0.631772
0.8
0.678832
0.678850
0.678850
0.678847
0.660412
0.9
0.710173
0.710175
0.710175
0.710170
0.687960
1.0
0.736827
0.736837
0.736837
0.736834
0.718260
110
Table 8: For fractional order 𝛾=0.5, the comparison of the fractional B-poly method to other
numerical methods [136], [138] is presented, where N denotes the number of B-polys in the basis
set.
𝛾=0.5
Our Results
X
N=6
N=7
N=8
N=9
𝑦ℎ𝑎𝑎𝑟[136]
𝑦ℎ𝑝𝑚[138]
0.0
0
0
0
0
0
0
0.1
0.330098
0.330097
0.330107
0.330096
0.324691
0.273875
0.2
0.436815
0.436838
0.43684
0.436846
0.432214
0.454125
0.3
0.504894
0.504896
0.504889
0.504888
0.504115
0.573932
0.4
0.553794
0.553780
0.553781
0.553775
0.553825
0.644422
0.5
0.591194
0.591187
0.591194
0.591194
0.590729
0.674137
0.6
0.621003
0.621012
0.621015
0.621019
0.622213
0.671987
0.7
0.645476
0.645491
0.645487
0.645487
0.643153
0.648003
0.8
0.666022
0.666024
0.666021
0.666016
0.667030
0.613306
0.9
0.683566
0.683553
0.683558
0.683560
0.680422
0.579641
1.0
0.698737
0.698755
0.698748
0.698743
0.695251
0.558557
Example 4. Consider a final example of the nonlinear fractional differential equation,
𝐷𝛾𝑦(𝑥)−𝑦(𝑥)+𝑦2(𝑥)+𝑦′2(𝑥)−1=0,
(28)
where the value of 0≤𝛾≤1 and the boundary conditions are specified as 𝑦(0)=𝑦(𝑅)=0.
The exact solution of the fractional differential equation for integral order 𝛾=2
111
is 𝑦𝑒𝑥𝑎𝑐𝑡(𝑥)=1+cos(𝑥). To determine the approximate solution, we consider 𝑦𝑎𝑝𝑝(𝑥)=
∑𝑏𝑖𝐵𝑖(𝛼,𝑥)+𝑦0
𝑛
𝑖, with initial conditions, 𝑦0=2&𝑦0
′=0. Putting the approximate solution
into Eq. (28), we get
𝑑𝛾
𝑑𝑥𝛾(∑𝑏𝑖𝐵𝑖(𝛼,𝑥)+𝑦0
𝑛
𝑖)−(∑𝑏𝑖𝐵𝑖(𝛼,𝑥)+𝑦0
𝑛
𝑖)
+(∑𝑏𝑖𝐵𝑖(𝛼,𝑥)+𝑦0
𝑛
𝑖)(∑𝑏𝑗𝐵𝑗(𝛼,𝑥)+𝑦0
𝑛
𝑗)
+(∑𝑏𝑖𝐵𝑖′(𝛼,𝑥)+𝑦0
′
𝑛
𝑖)(∑𝑏𝑗𝐵𝑗′(𝛼,𝑥)+𝑦0
′
𝑛
𝑗)−1=0.
(29)
Evaluating Eq. (29) by applying the Caputo fractional derivative on the first term,
multiplying both sides of the equation by the elements of the fractional B-polys basis set,
𝐵𝑙(𝛼,𝑥), and then integrating both sides of the equation over the interval [0, R]. The expression
in Eq. (29) is transformed into,
∑𝑏𝑖
𝑛
𝑖[⟨𝐷𝛾𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩−⟨𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩
+∑𝑏𝑗
𝑛
𝑗⟨𝐵𝑖(𝛼,𝑥)𝐵𝑗(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩+2⟨𝑦0𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩
+2⟨𝑦0
′𝐵𝑖′(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩+∑𝑏𝑗
𝑛
𝑗⟨𝐵𝑖′(𝛼,𝑥)𝐵𝑗′(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩]
=⟨(1−𝑦02−𝑦0
′2+𝑦0)𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩.
(30)
The following is a nonlinear matrix representation of the above equation,
112
𝐵[𝐴−𝐶+𝐷+𝐸+𝐹+𝐺]=𝑊,
(31)
where elements of each matrix in Eq. (31) are given as follows:
𝐴=⟨𝐷𝛾𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩=∫ 𝐷𝛾(𝐵𝑖(𝛼,𝑥))𝐵𝑙(𝛼,𝑥)𝑑𝑥,
𝑅
0
𝐶=⟨𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩=∫ 𝐵𝑖(𝛼,𝑥)𝐵𝑙(𝛼,𝑥)𝑑𝑥
𝑅
0,
𝐷=2⟨𝑦0𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩=2∫ 𝑦0𝐵𝑖(𝛼,𝑥)𝐵𝑙(𝛼,𝑥)𝑑𝑥
𝑅
0,
𝐸=2⟨𝑦0
′𝐵𝑖′(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩=2∫ 𝑦0
′𝐵𝑖′(𝛼,𝑥)𝐵𝑙(𝛼,𝑥)𝑑𝑥
𝑅
0,
𝐹=∑𝑏𝑗
𝑛
𝑗𝑔𝑖𝑗𝑙 =𝑏𝑗⟨𝐵𝑖(𝛼,𝑥)𝐵𝑗(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩=𝑏𝑗∫ 𝐵𝑖(𝛼,𝑥)𝐵𝑗(𝛼,𝑥)𝐵𝑙(𝛼,𝑥)𝑑𝑥,
𝑅
0
𝐺=∑𝑏𝑗
𝑛
𝑗𝑔𝑖𝑗𝑙
′=𝑏𝑗⟨𝐵𝑖′(𝛼,𝑥)𝐵𝑗′(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩=𝑏𝑗∫ 𝐵𝑖′(𝛼,𝑥)𝐵𝑗′(𝛼,𝑥)𝐵𝑙(𝛼,𝑥)𝑑𝑥,
𝑅
0
𝑊=⟨(1−𝑦02−𝑦0
′2+𝑦0)𝐵𝑖(𝛼,𝑥)|𝐵𝑙(𝛼,𝑥)⟩
=∫ (1−𝑦02−𝑦0
′2+𝑦0)𝐵𝑖(𝛼,𝑥)𝐵𝑙(𝛼,𝑥)𝑑𝑥
𝑅
0.
(32)
To obtain an initial guess for the nonlinear terms, we ignore nonlinear matrices F and G
and solve the following linear matrix equation,
𝐵[𝐴−𝐶+𝐷+𝐸]=𝑊.
(33)
Using the initial guess from Eq. (33) and substituting it into Eq. (31), we get a new
approximation for the unknown coefficients, which are utilized in the linear combination to
construct the approximate solution. This process is repeated to obtain a converged solution to the
fractional differential equation (28). To solve the nonlinear fractional differential equation for
113
different fractional values of 𝛾, an iterative scheme has been used. The convergent solutions are
plotted for various fractional order of 𝛾.
Figure 26 shows the absolute error between estimated 𝑓(𝑥) and exact (𝑠𝑜𝑙) solutions is
of the order of 10−14. This is the highest accuracy achieved after choosing a set of 8 fractional
polynomials. Further accuracy is attainable by increasing the set of fractional polys, but it
requires a larger number of iterations to obtain convergent solutions. In the references [136],
[140], the absolute error is of the order 10−3. As a result, our technique performs better in terms
of computational accuracy. Figure 27 depicts graphs of various fractional orders 𝛾=
2,1.5,1.7,1.9of nonlinear differential equations. The plots of various approximate solutions to
various nonlinear fractional-order differential equations are presented to show the smoothness of
the curves in the graphs.
114
Tables 8, 9, and 10 show that our method predicts excellent solutions compared to other
approximate solutions for various fractional order values[136], [140]. The accuracy of our
method depends on the number of fractional B-polys and iterations used. We receive convergent
values after a few iterations, and only converged solutions are reported. Comparisons between
Figure 26: For integral order 𝛾=2 in nonlinear differential equation (28), a 1-D graph of approximation
𝑓(𝑥), and the precise solution (sol) is displayed on the left side. It shows that the two solutions essentially
overlap each other. The plot on the right shows the absolute error between the precise and approximate
results. The accuracy of the results is much higher as compared to other results in the literature.
Figure 27: This graph displays the plots of several approximate solutions for various
values of the fractional order 𝛾=2,1.9,1.7,1.5 used in Equation (28). All the results
converged after 10 number of iterations.
115
the available results are given in these tables. In some cases present results are superior in
accuracy as well as less computational time is involved. The error analysis is also carried out
using N = 6, 7, and 8 number of fractional B-polys basis set. The error in the solution sets
decreased as we increased the number of fractional B-polys for solving fractional order
differential equations. The results of this observation are shown in the Tables presented in this
paper.
Table 9: For an integral value of 𝛾=2, the comparison of the fractional B-poly method to other
methods[136], [140] is presented, where N denotes the number of fractional B-polys in the basis
set.
𝛾=2
Our results
x
N=6
N=7
N=8
𝑦ℎ𝑎𝑎𝑟[136]
𝑦𝑣𝑚[140]
𝑦ℎ𝑝𝑚[140]
𝑦𝑒𝑥𝑎𝑐𝑡[140]
0.0
2.000000
2.000000
2.000000
2.000000
2.000000
2.000000
2.000000
0.1
1.995004
1.995004
1.995004
1.995004
1.994996
1.995013
1.995004
0.2
1.980067
1.980067
1.980067
1.980067
1.979933
1.980200
1.980067
0.3
1.955336
1.955336
1.955336
1.955336
1.954661
1.956013
1.955336
0.4
1.921061
1.921061
1.921061
1.921061
1.918928
1.923200
1.921061
0.5
1.877583
1.877583
1.877583
1.877583
1.872374
1.882813
1.877583
0.6
1.825336
1.825336
1.825336
1.825336
1.814535
1.836200
1.825336
0.7
1.764842
1.764842
1.764842
1.764842
1.744831
1.785013
1.764842
0.8
1.696707
1.696707
1.696707
1.696707
1.662565
1.731200
1.696707
0.9
1.621610
1.621610
1.621610
1.621610
1.566914
1.677013
1.621610
1.0
1.540302
1.540302
1.540302
1.540302
1.456919
1.6250000
1.540302
116
Table 10: For fractional order of 𝛾=1.3&1.5, the comparison of the fractional B-poly method
to other methods [136] is shown, where N denotes the number of B-polys in the basis sets.
𝛾=1.3
𝛾=1.5
Our Result
Our Result
x
N=6
N=7
N=8
𝑦ℎ𝑎𝑎𝑟[136]
N=6
N=7
N=8
𝑦ℎ𝑎𝑎𝑟[136]
0.0
2.00000
2.00000
2.00000
2.00000
2.00000
2.00000
2.00000
2.00000
0.1
1.92197
1.96543
1.94497
1.96361
1.97465
1.97496
1.97511
1.97607
0.2
1.84223
1.87730
1.85567
1.90615
1.92798
1.92855
1.92880
1.93585
0.3
1.78373
1.74742
1.73137
1.84812
1.86749
1.86815
1.86843
1.88705
0.4
1.74370
1.60752
1.59675
1.79521
1.79691
1.79762
1.79794
1.83225
0.5
1.71384
1.49123
1.48813
1.75255
1.71985
1.72063
1.72099
1.77573
0.6
1.68703
1.41712
1.42255
1.70324
1.64033
1.64116
1.64153
1.72092
0.7
1.66028
1.38228
1.39164
1.65990
1.56273
1.56353
1.56387
1.66355
0.8
1.63486
1.37065
1.38017
1.61351
1.49145
1.49217
1.49248
1.61749
0.9
1.61419
1.37032
1.38129
1.57941
1.43030
1.43092
1.43119
1.56404
1.0
1.59976
1.37954
1.38997
1.55214
1.38178
1.38228
1.38250
1.51286
117
Table 11: For fractional order of 𝛾=1.7&1.9, the comparison of the fractional B-poly method
to other methods [136] is given, where N denotes the number of B-polys in the basis sets.
𝛾=1.7
𝛾=1.9
Our Result
Our Result
x
N=6
N=7
N=8
𝑦ℎ𝑎𝑎𝑟[13
6]
N=6
N=7
N=8
𝑦ℎ𝑎𝑎𝑟[13
6]
0.
0
2.000000
2.000000
2.000000
2.000000
2.000000
2.000000
2.000000
2.000000
0.
1
1.986884
1.986925
1.986948
1.987162
1.993098
1.993103
1.993107
1.988393
0.
2
1.957458
1.957558
1.957610
1.958162
1.974312
1.974326
1.974338
1.969996
0.
3
1.915545
1.915685
1.915754
1.917257
1.944746
1.944769
1.944788
1.948187
0.
4
1.863133
1.863297
1.863378
1.873808
1.905153
1.905184
1.905208
1.907802
0.
5
1.801916
1.802104
1.80220
1.820153
1.856259
1.856260
1.856326
1.862716
0.
6
1.733598
1.733812
1.73392
1.768749
1.798833
1.798876
1.798912
1.811744
0.
7
1.659979
1.660211
1.660326
1.705186
1.733716
1.733766
1.733806
1.741617
0.
8
1.582999
1.583239
1.583359
1.640477
1.661831
1.661886
1.661930
1.680558
0.
9
1.504768
1.505014
1.505137
1.570663
1.584186
1.584246
1.584294
1.597211
1.
0
1.427572
1.427817
1.427939
1.506804
1.501885
1.501949
1.502000
1.526272
118
Conclusion
This study solves the one-dimensional (1-D) fractional-order nonlinear differential
equations using the fractional Bhatti polynomial bases set. It is shown that this method of
fractional B-polys works well to solve such types of equations, and the accuracy of the results
can be adjusted by increasing the size of the basis set. The method is applied to four 1-D
nonlinear fractional differential equations (NFDEs) with initial and boundary conditions
imposition. The method can predict highly accurate solutions. In addition to its efficacy, the
predicted results are compared to other methods [136]–[138], [140]. The method's versatility is
shown to solve various types of differential equations with various values of the fractional order.
The solution of the fractional-order differential equation is expanded in terms of the linear
combination of coefficients which are determined using the Galerkin scheme [141]. In all 4
examples considered, the linear matrix equation is attained by setting the nonlinear terms of the
equation equal to zero. We solve the linear part of the differential equation to obtain an initial
guess for the unknown coefficients 𝑏𝑖 of column matrix B. For example, see Eq. (33), the
nonlinear matrix equation is inverted to obtain a new guess in the iteration process, Eq. (31). The
iteration procedure is continued until convergent values of the expansion coefficients are
obtained. Final expansion coefficients are used to approximate the solution (Eq. 3) of the 1-D
fractional-order nonlinear differential equations.
Figures 20 through 27 show the graphs of the convergent solutions of the NFDEs. The
error analysis is also carried out for both integral and fractional order differential equations and
is presented in the Tables with basis sets. It is clear as the number of iterations increases, the
results converge quickly after about 10 iterations, and by increasing the fractional B-poly basis
set, the accuracy increases rapidly. We compared our solutions to the exact solutions for
119
nonfractional cases and found that the agreement is also excellent between fractional-order cases.
Our absolute errors between approximate and exact solutions range between 10−8 and 10−14 in
all the examples. The results for fractional-order cases are significantly better as compared to
other references [136]–[138], [140]. We revealed that our state-of-the-art method could solve a
variety of examples with greater precision, as reported in Tables 3-11 and Figures 20-27. In this
article, we provide the graphs for different fractional orders in Figures 21, 23, 25, and 27. All
calculations, including integrations, fractional differentiations, iterations of results, and matrix
inversions, are carried out using Wolfram Mathematica symbolic program version-13 [135].
Our method revealed a strong potential for solving nonlinear fractional-order differential
equations with a higher degree of precision, easy to use, and ease of implementation. The
computing time for solving examples of integral-order differential equations is less than a
minute. For fractional orders examples, the computing time is about 10-30 minutes for
accomplishing converged results.
120
CHAPTER VII
CONCLUSION
In conclusion, this thesis presents a new method for solving one-dimensional (1-D)
and two-dimensional (2-D) fractional-order nonlinear and linear partial differential equations
using the fractional Bhatti polynomial bases set. The method proves to be effective in
accurately solving various types of equations, including nonlinear and linear fractional
differential equations and nonlinear and linear partial differential equations, with different
initial and boundary conditions.
For nonlinear and linear partial differential equations, the method demonstrates its
efficacy by providing highly accurate solutions. The accuracy of the results can be improved by
increasing the size of the fractional B-polynomial basis set. The method is compared to other
existing methods, and it outperforms them in terms of accuracy and precision. The
absolute errors between the approximate and exact solutions are found to be remarkably
small, further emphasizing the superiority of the proposed method.
In the case of 2-D linear partial fractional differential equations, the method
is successfully extended to handle such problems. The Galerkin method is employed, and
the convergence of the solutions is achieved within a reasonable number of iterations. The
accuracy of the solutions depends on the chosen basis set and the size of the operational matrix
used. The
121
method showcases its potential in solving complex systems of linear and nonlinear differential
equations in multiple scientific disciplines.
The computational implementation of the method is facilitated by utilizing Wolfram
Mathematica symbolic program, ensuring efficient and reliable computations. The computational
time required for solving integral-order differential equations is minimal, while for fractional
orders, it ranges from 10 to 30 minutes, depending on the complexity of the problem and the size
of the basis set.
Overall, this thesis highlights the robustness and versatility of the proposed method in
solving fractional differential equations. It provides a valuable tool for researchers in various
fields, including chemistry, physics, genetics, and other related disciplines, enabling them to
obtain precise and reliable solutions. The presented technique opens doors for further research in
nonlinear partial fractional differential equations and offers great potential for solving linear and
nonlinear multidimensional fractional differential equation problems.
