Next Article in Journal
A Novel Treatment of Fuzzy Fractional Swift–Hohenberg Equation for a Hybrid Transform within the Fractional Derivative Operator
Next Article in Special Issue
Chirped Periodic and Solitary Waves for Improved Perturbed Nonlinear Schrödinger Equation with Cubic Quadratic Nonlinearity
Previous Article in Journal
Some Dynamic Inequalities via Diamond Integrals for Function of Several Variables
Previous Article in Special Issue
Approximate Solutions of Nonlinear Partial Differential Equations Using B-Polynomial Bases
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Technique to Solve Linear Fractional Differential Equations Using B-Polynomials Bases

by
Muhammad I. Bhatti
* and
Md. Habibur Rahman
Department of Physics and Astronomy, University of Texas Rio Grande Valley, Edinburg, TX 78539, USA
*
Author to whom correspondence should be addressed.
Fractal Fract. 2021, 5(4), 208; https://doi.org/10.3390/fractalfract5040208
Submission received: 9 October 2021 / Revised: 25 October 2021 / Accepted: 2 November 2021 / Published: 11 November 2021
(This article belongs to the Special Issue Frontiers in Fractional Schrödinger Equation)

Abstract

:
A multidimensional, modified, fractional-order B-polys technique was implemented for finding solutions of 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 was employed. Moreover, minimization of error in the coefficients was found by employing the Galerkin method. Before the Galerkin method was applied, the linear FPDE was transformed into an operational matrix equation that was inverted to provide the values of the unknown coefficients in the approximate solution. A valid multidimensional solution was determined when an appropriate number of basis sets and fractional-order of B-polys were chosen. In addition, initial conditions were applied to the operational matrix to seek proper solutions in multidimensions. The technique was applied to four examples of linear FPDEs and the agreements between exact and approximate solutions were found to be excellent. The current technique can be expanded to find multidimensional fractional partial differential equations in other areas, such as physics and engineering fields.

1. Introduction

In real-world scientific phenomena, most problems follow either linearity or nonlinearity in their systems. In different fields, for example, engineering, computer science, and chemistry [1,2,3], fractional-order differential equations emerge more often. The most physical phenomena are described by the differential systems, which are integral-order systems. Many systems could be expressed with the help of the fractional differential equation [4,5,6,7,8,9]. Due to the materials and chemical properties of the real-world problems, and their memory and genetic characteristics, many physical problems follow fractional dynamical behavior [9,10,11]. The partial fractional-order differential equations are becoming a useful tool to model the physical phenomena [9]. For this reason, there has been an urgent need to find a solution to the fractional-order problems. However, there has been difficulty in obtaining the accurate analytical or numerical results of the most fractional-order differential model equations. There is a need for a suitable technique to find solutions to the fractional-order differential problems, linear and nonlinear. In our current paper, we aim to apply a technique to resolve multivariable linear fractional-order differential equations. Nonlinear fractional-order partial differential equations would be considered in 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 [12,13,14,15] the variational iteration procedure [16], Adams–Bashfort–Mowlton Method [17], the Lagrange characteristic approach [18,19], Adomian decomposition method [20], the finite difference procedure [21], the differential transformation method [22], the finite element technique [23], the fractional sub equation procedure [24], the ( G /G)-expansion method [25], first integral approach [26], Jacobi–Gauss collocation method [27], the spectral collocation method (SCM) [28,29], and the fractional complex transform technique [30]. 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 [31,32,33,34,35,36,37] that is significantly able to solve a variety of multivariable linear fractional-order differential equations. We chose the fractional-order B-poly due to its well-defined basis set and precision [33]. With these basis sets, it can be demonstrated that an arbitrary function can be represented to the desired accuracy and is directly differentiable over a closed interval. In several papers [31,32,33,34,35,36], using the B-poly basis of fractional-order and a generalized Galerkin method, the authors were able to find solutions to the fractional-order partial differential equations.
In an earlier work [31,32,33,34,35,36], the authors used a similar technique to find the solution to ordinary nonlinear and linear multidimensional differential equations. The current study focuses on linear fractional-order differential equations using the generalized Galerkin method [36] and the B-poly basis of fractional-order. This technique has the special advantage of the unitary partition property and the continuity of the generalized fractional-order B-polys over an interval [0, R], which are seamlessly differentiated. With the help of fractional-order B-polys, a fractional-order differential equation is transformed into an operational matrix using a matrix formalism that provides greater flexibility to the application of 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. In the following sections, we present analytical formulism to employ Caputo’s fractional-order derivative on the polynomials, present the process used to create fractional-order basis sets, and develop an algorithm to resolve various linear fractional-order partial differential equations. We apply this technique to four examples. Finally, we shall present an error analysis of one of the fourth considered examples.

2. Caputo’s Fractional Differential-Order Operator

The explanation of the fractional-order derivative of Caputo is provided as [3]
D γ f ( x ) = J m γ D m f ( x ) = 1 Γ ( m β ) 0 x ( x t ) m γ 1 f ( m ) ( t ) d t ,     f o r   m 1 < γ m ,   m N ,   x > 0 ,   f C 1 m ,
where D γ are Caputo’s fractional operator and fractional derivative in Caputo’s sense is D γ f ( x ) , Equation (1). Caputo’s derivative of a constant is zero, i.e., D γ C = 0 and a fractional derivative of the polynomial D γ x α is given by
D γ x α = { 0                                                       f o r   α   N 0   a n d   α < [ γ ] Γ ( α + 1 )     x α γ Γ ( α + 1 γ )                                         o t h e r w i s e .
Here, α denotes the order of the fractional function. The unknown two-variable dependent function U ( x , t ) is expanded as a product of two generalized fractional-order B-polynomials, B j , m ( α , t ) B i , n ( α , x ) , which may be considered as an approximate outcome to the FPD equation represented by
U ( x , t ) = i ,   j = 0 n b j i   B j , m ( α , t ) B i , n ( α , x ) ,
where B j , m ( α , t ) is a j-th and m-degree fractional-order B-poly in variable t   o r   x , with α as a fractional-order parameter over a given interval. The expansion coefficients b j i in Equation (3) 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
D x γ ( i , j = 0 n b j i   B j , m ( α , t ) B i , n ( α , x ) ) = i , j = 0 n b j i ( B j , m ( α , t )   D x γ ( B i , n ( α , x ) ) ) .
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 of the linear fractional-order partial differential equation.

3. Fractional-Order B-Poly Basis

The generalized form of fractional B-polys B i , n ( α , x ) in terms of variable x or t over an interval [0, R] or [0, T] are defined in Refs. [36,38]:
B i , n ( α , x ) = i = 0 n β i , k ( x R ) α k .
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 Equation (5). The factor β i , k in Equation (5) is defined as
β i , k = ( 1 ) i k ( n k ) ( k i ) ,
where this binomial coefficient is defined as, ( n k ) = n ! k ! ( n k ) ! . For convenience, if i < 0 or i > n we can set B i , k ( α ,   x ) = 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 the first and last polynomial of the basis set. As an example, when n = 10 and fractional order α = 1 2 , 5 3 , 9 4 are chosen in Equation (5), the corresponding basis sets of B-polys are plotted in Figure 1. Graphs of these fractional-order B-polys show how these B-polys add up to 1 at any given point, x. Such B-polys sets may be used to represent an arbitrary function with higher accuracy.

4. Technique for Approximating Solutions

Using the Galerkin method [36] and the generalized fractional-order B-poly basis set, we exploit a technique to seek practical solutions to fractional-order partial differential equations. We intend to apply Caputo’s fractional derivative to the fractional-order B-ploys. The 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, Equation (3) is substituted into the fractional-order differential equation and by-products are separated in terms of integral products in both variables x and t. Finally, both sides of the fractional equation are multiplied with fractional B-polys basis elements, B m ( α , x )   B n ( α , t ) , and the integrations are carried out using symbolic program, Mathematica, over the closed intervals [0, R] and [0, T], respectively. For example, the integration of the two fractional B-polys is given in the closed symbolic formula
m i , j = ( B i , n ( α , x ) , B j , n ( α , x ) ) = k = i n β i , k ( x R ) α k l = j n α i , k ( x R ) α l R ( k + l ) α
Caputo’s derivative defined in Equation (2) is applied to the fractional B-ploy basis set, leading to the following closed results:
D x γ ( B i , n ( α , x ) ) = k = i n α i , k D x γ ( x R ) α   k = k = i n β i , k R α   k   Γ ( α k + 1 ) Γ ( α k + 1 γ )   x α   k γ , d i , j ( γ ) ( x ) = ( D x γ B i , n ( α , x ) , B j , n ( α , x ) ) =   D x γ B i , n ( α , x ) | B j , n ( α , x ) = k = i , l = j n β i , k β j , k Γ ( α k + 1 ) Γ ( α k + 1 γ )   R 1 γ ( ( k + l ) α + 1 γ ) ,
and the integrals of some arbitrary functions are given
F ( x , t ) = (   f ( x , t ) , B i , n ( α , x ) ) = k = i n β i , k R α k 0 R f ( x , t )   x α k d x , W m , n =   0 R , T f ( x , t )   B m ( α , x ) B n ( α , t ) d x   d t .
With the help of these analytic formulas Equations (7)–(9), the operational matrix is constructed. The inverse of the operational matrix is required to find out the unknown coefficients b j i of the linear combination in Equation (3). In the next section, we will describe our technique, as well as how to obtain a desirable result for the linear fractional-order partial differential equation. The technique will be employed in four examples to demonstrate that it works appropriately for approximating the accurate solutions. Plots of the approximate as well as exact solutions will be presented for comparison. An absolute error analysis of the fourth example will be introduced to show that, when the basis set of the fractional B-polys is enlarged, the accuracy of the solution is increased. 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 subscript n and m from the fractional B-polys, i . e . , B i , n ( α , x ) = B i   ( α , x )     and   B j , m ( α , t ) = B j   ( α , t ) .
Example 1: 
Let us introduce a linear partial fractional-order differential equation of the form
2 d γ U ( x , t ) d t γ + d γ U ( x , t ) d x γ = 0 .
The ideal solution to Equation (10) is U e x a c t ( x , t ) = ( x γ t γ / 2 ) . A numerical solution is sought out in the intervals 0   x   1   &   0   t   1   using   initial   condition   U ( x ,   0 ) =   f   ( x ) = x γ . The assumed solution, U a p p ( x , t ) = i , j = 0 n b j i   B j ( α , t )   B i ( α , x ) + x γ , is substituted into the Equation (10) and the result is presented below:
2 d γ d t γ ( i , j = 0 n b j i   B j ( α , t )   B i ( α , x ) + x γ ) + d γ d x γ ( i , j = 0 n b j i   B i ( α , t )   B i ( α , x ) + x γ ) = 0
Caputo’s derivative operator is applied to Equation (11). The product of fractional B-polys B m ( α , x )   B n ( α , t ) from the basis set is multiplied on both sides of the Equation (11) and the integration on both variables is calculated over the intervals using a symbolic program. This operation provides the following equation
i , j = 0 n b j i [ 2 B i ( α , x ) | B m ( α , x ) D t γ B j ( α , t ) | B n ( α , t ) D x γ   B i ( α , x ) | B m ( α , x )   B j ( α , t ) | B n ( α , t ) ] = f γ ( α , x ) | B m ( α , x ) | B n ( α , t ) ,
where f γ ( x ) = D x γ ( x γ ) = Γ ( γ + 1 ) with α = γ . The current technique leads to a system of ( n + 1 ) × ( n + 1 ) equations. The elements of matrix B = { b 1 1 , b 2 1 ,   b 3 1 ,   , b 1 2 , b 2 2 ,   b 3 2 ,   ,   } are the unknown constants that are involved in those equations. After further simplification, the right-hand side column matrix W and the matrix elements of operational matrix X in terms of inner products of B-polys are given
  X m , n = i , j = 0 n [ 2   B i ( α , x ) | B m ( α , x ) D t γ B j ( α , t ) | B n ( α , t ) D x γ   B i ( α , x ) | B m ( α , x )   B j ( α , t ) | B n ( α , t ) ] , W m , n = f γ ( x ) | B m ( α , x ) | B n ( α , t ) = 0 R , T Γ ( γ + 1 )   B m ( α , x ) B n ( α , t ) d x   d t .
The partial fractional-order differential Equation (10) is now transformed into a matrix equation X   B = W . By deleting the rows and corresponding columns of the equation (13), the initial conditions are imposed on the operational matrix equation X and the corresponding matrix W, so that the solution vanishes at t = 0 and x = 0. The operational matrix X was coded in 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 b j i . The emerging estimated result is composed of the linear combination of the B-poly basis set via Equation (3). The process provides a valid approximate solution U a p p ( x , t ) of the Equation (10) using B-polys of fractional-order α = 1 2 and fractional differential-order of γ = 1 2 in Equation (10) is given below:
U a p p ( x , t ) = x γ + t γ ( 0.5 + 0 . × 10 30 x γ )   x γ t γ / 2 .
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 γ   and   α . To solve the fractional-order partial differential Equation (10), we choose n = 1 and α = 1 2 order B-poly basis set { 1 t , t } and { 1 x , x } in variables t and x, respectively. The corresponding coefficient values 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 } , as seen in the last column of Table 1. A 3D plot of the estimated and the exact results of Equation (10) are presented in Figure 1 for comparison, and an excellent agreement can be seen between both results at the level of machine accuracy. Note that when t = x is substituted into Equation (14), 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 results, both exact and approximate, shows that both results have excellent reliability. The absolute error in the 3D graph is also presented on the right-hand side in Figure 2. The 3D graph shows that the absolute error in the converged solution is of the order of 10 17 .
Example 2: 
Consider another example of fractional-order linear partial differential equation with different initial condition U ( x ,   0 ) =   f   ( x ) = E α , 1 ( x α )
2 d γ U ( x , t ) d t γ + d γ U ( x , t ) d x γ = 0 .  
The ideal solution of the Equation (15) is U e x a c t ( x , t ) = E α ,   1 ( x α t α / 2 ) . The function E α ,   β ( z ) , is called the Mittag–Leffler function [39] and is described as E α ,   β ( z ) = k = 0 Z k Γ ( k   α + β ) . In the summation of Mittag–Leffler function, we only kept k = 15 in the summation of terms. Therefore, the accuracy of the numerical solution will likely depend on the number of terms that we would keep in the summation of the Mittag–Leffler function. According to Equation (3), an estimated solution of Equation (15) using the initial condition may be assumed as U a p p ( x , t ) = i = 0 n a i ( α , t )   B i ( α , x ) + E α , 1 ( x α ) . After substituting this expression into the Equation (15). The Galerkin method, [29] and [32], is also applied to the presumed solution to obtain
2 d γ d t γ ( i , j = 0 n b j i   B j ( α , t )   B i ( α , x ) + E α , l ( x α ) ) + d γ d x γ ( i , j = 0 n b j i   B j ( α , t )   B i ( α , x ) + E α , l ( x α ) ) = 0 .
Caputo’s fractional derivative is applied to Equation (16), and the product of fractional B-polys B m ( α , x )   B n ( α , t ) from the basis set is multiplied on both sides of the Equation (16). The resulting integration of both variables (t and x) is calculated over the intervals 0 x 1 and 0 t 1 , respectively. After further simplification of the Equation (16), we obtain
i , j = 0 n b j i [ 2 B i ( α , x ) | B m ( α , x ) D t γ B j ( α , t ) | B n ( α , t ) D x γ B i ( α , x ) | B m ( α , x )   B j ( α , t ) | B n ( α , t ) ] = f γ ( α , x ) | B m ( α , x ) | B n ( α , t ) ,
where the fractional-order derivative of the Mittag-Leffler function f γ ( α , x ) = d γ d x γ ( E α , 1 ( x α ) ) = E α , 1 ( x α ) , with α = γ is used. The current technique leads to a system of ( n + 1 ) ×   ( n + 1 ) equations. This system of equations may be summarized in the matrix equation X   B = W , where the elements of matrix B = { b 1 1 , b 2 1 ,   b 3 1 ,   , b 1 2 , b 2 2 ,   b 3 2 ,   ,   } are the unknown constants. The right-hand side column matrix elements of W and the matrix elements of operational matrix X are given as
X m , n =   i , j = 0 n [ 2   B i ( α , x ) | B m ( α , x ) D t γ B j ( α , t ) | B n ( α , t ) D x γ B i ( α , x ) | B m ( α , x ) B j ( α , t ) | B n ( α , t ) ] , W m , n =   f γ ( γ , x ) | B m ( α , x ) | B n ( α , t ) = 0 R , T f γ ( γ , x )   B m ( α , x ) B n ( α , t ) d x   d t .
By deleting the rows and corresponding columns of Equation (18), the initial condition was imposed on the operational matrix, to make sure the solution vanishes at x = 0 and t = 0. The operational matrix X is inverted using Mathematica symbolic program and multiplied with the column matrix W to solve equation B = X 1 W and yield values of the unknown coefficients b j i . The emerging estimated result is composed of the B-poly basis set and the coefficients via Equation (3). The process provides an approximate solution U a p p ( x , t ) to Equation (15) using B-polys of fractional-order α = 1 2 and a fractional differential-order of γ = 1 2 . The final approximate solution is provided below:
U a p p ( x , t ) =   1.0 + t 11 / 2 ( 1.654 × 10 6 1.914 × 10 6 x ) + t 6 ( 2.998 × 10 7 + 3.826 × 10 7 x ) + t 5 ( 8.106 × 10 6 + 9.183 × 10 6 x + 8.138 × 10 6 x ) + 1.128 x + 1.0 x + 0.7522 x 3 / 2 + 0.5 x 2 + 0.3009 x 5 / 2 + 0.1667 x 3 + 0.0860 x 7 / 2 + 0.0417 x 4 + 0.0191 x 9 / 2 + 8.333 × 10 3 x 5 + 3.473 × 10 3 x 11 / 2 + 0.00139 x 6 + 5.344 × 10 4 x 13 / 2 + 1.984 × 10 4 x 7 + 7.125 × 10 5 x 15 / 2 + t 9 / 2 ( 3.729 × 10 5 4.210 × 10 5 x 3.731 × 10 5 x 2.807 × 10 5 x 3 / 2 ) + t 4 ( 1.628 × 10 4 + 1.837 × 10 4 x + 1.628 × 10 4 x + 1.224 × 10 4 x 3 / 2 ) + t 7 / 2 ( 6.717 × 10 4 7.579 × 10 4 x 6.717 × 10 4 x 5.053 × 10 4 x 3 / 2 3.358 × 10 4 x 2 ) + t 3 ( 2.604 × 10 3 + 2.939 × 10 3 x + 2.604 × 10 3 x + 1.959 × 10 3 x 3 / 2 + 1.302 × 10 3 x 2 + 7.836 × 10 4 x 5 / 2 + 4.34 × 10 4 x 3 ) + t 5 / 2 ( 9.403 × 10 3 0.0106 x 9.403 × 10 3 x 0.0071 x 3 / 2 0.0047 x 2 0.0028 x 5 / 2 0.0016 x 3 8.084 × 10 4 x 7 / 2 ) + t 2 ( 0.03125 + 0.0353 x + 0.0313 x + 0.0235 x 3 / 2 + 0.0156 x 2 + 9.403 × 10 3 x 5 / 2 + 5.208 × 10 3 x 3 + 2.687 × 10 3 x 7 / 2 + 1.302 × 10 3 x 4 + 5.970 × 10 4 x 9 / 2 + 2.604 × 10 4 x 5 ) + t 3 / 2 ( 0.0940 0.1061 x 0.0940 x 0.0707 x 3 / 2 0.047 x 2 0.0283 x 5 / 2 0.0157 x 3 8.084 × 10 3 x 7 / 2 3.918 × 10 3 x 4 1.796 × 10 3 x 9 / 2 7.836 × 10 4 x 5 3.266 × 10 4 x 11 / 2 1.306 × 10 4 x 6 5.025 × 10 5 x 13 / 2 ) + t ( 0.25 + 0.2821 x + 0.25 x + 0.188 x 3 / 2 + 0.125 x 2 + 0.0752 x 5 / 2 + 0.0417 x 3 + 0.0215 x 7 / 2 + 0.0104 x 4 + 4.776 × 10 3 x 9 / 2 + 2.083 × 10 3 x 5 + 8.684 × 10 4 x 11 / 2 + 3.472 × 10 4 x 6 + 1.336 × 10 4 x 13 / 2 + 4.960 × 10 5 x 7 ) + t ( 0.5642 0.6366 x 0.5642 x 0.4244 x 3 / 2 0.2821 x 2 0.1698 x 5 / 2 0.094 x 3 0.0485 x 7 / 2 0.0235 x 4 0.0108 x 9 / 2 4.702 × 10 3 x 5 1.96 × 10 3 x 11 / 2 7.836 × 10 4 x 6 3.015 × 10 4 x 13 / 2 1.119 × 10 4 x 7 4.02 × 10 5 x 15 / 2 )
To solve the fractional order partial differential Equation (15), we chose 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 x + 105 x 455 x 3 / 2 + 1365 x 2 3003 x 5 / 2 + 5005 x 3 6435 x 7 / 2 + 6435 x 4 5005 x 9 / 2 + 3003 x 5 1365 x 11 / 2 + 455 x 6 105 x 13 / 2 + 15 x 7 x 15 / 2 , 15 x 210 x + 1365 x 3 / 2 5460 x 2 + 15,015 x 5 / 2 30,030 x 3 + 45,045 x 7 / 2 51,480 x 4 + 45,045 x 9 / 2 30,030 x 5 + 15,015 x 11 / 2 5460 x 6 + 1365 x 13 / 2 210 x 7 + 15 x 15 / 2 , 105 x 1365 x 3 / 2 + 8190 x 2 30,030 x 5 / 2 + 75,075 x 3 135,135 x 7 / 2 + 180,180 x 4 180,180 x 9 / 2 + 135,135 x 5 75,075 x 11 / 2 + 30,030 x 6 8190 x 13 / 2 + 1365 x 7 105 x 15 / 2 , 455 x 3 / 2 5460 x 2 + 30,030 x 5 / 2 100,100 x 3 + 225,225 x 7 / 2 360360 x 4 + 420,420 x 9 / 2 360,360 x 5 + 225,225 x 11 / 2 100,100 x 6 + 30,030 x 13 / 2 5460 x 7 + 455 x 15 / 2 , 1365 x 2 15,015 x 5 / 2 + 75,075 x 3 225,225 x 7 / 2 + 450,450 x 4 630,630 x 9 / 2 + 630,630 x 5 450,450 x 11 / 2 + 225,225 x 6 75,075 x 13 / 2 + 15,015 x 7 1365 x 15 / 2 , 3003 x 5 / 2 30,030 x 3 + 135,135 x 7 / 2 360,360 x 4 + 630,630 x 9 / 2 756,756 x 5 + 630,630 x 11 / 2 360,360 x 6 + 135,135 x 13 / 2 30,030 x 7 + 3003 x 15 / 2 , 5005 x 3 45,045 x 7 / 2 + 180,180 x 4 420,420 x 9 / 2 + 630,630 x 5 630,630 x 11 / 2 + 420,420 x 6 180,180 x 13 / 2 + 45,045 x 7 5005 x 15 / 2 , 6435 x 7 / 2 51,480 x 4 + 180,180 x 9 / 2 360,360 x 5 + 450,450 x 11 / 2 360,360 x 6 + 180,180 x 13 / 2 51,480 x 7 + 6435 x 15 / 2 , 6435 x 4 45,045 x 9 / 2 + 135,135 x 5 225,225 x 11 / 2 + 225,225 x 6 135,135 x 13 / 2 + 45,045 x 7 6435 x 15 / 2 , 5005 x 9 / 2 30,030 x 5 + 75,075 x 11 / 2 100,100 x 6 + 75,075 x 13 / 2 30,030 x 7 + 5005 x 15 / 2 , 3003 x 5 15,015 x 11 / 2 + 30,030 x 6 30,030 x 13 / 2 + 15,015 x 7 3003 x 15 / 2 , 1365 x 11 / 2 5460 x 6 + 8190 x 13 / 2 5460 x 7 + 1365 x 15 / 2 , 455 x 6 1365 x 13 / 2 + 1365 x 7 455 x 15 / 2 , 105 x 13 / 2 210 x 7 + 105 x 15 / 2 , 15 x 7 15 x 15 / 2 , x 15 / 2 } .
To obtain 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 set, the accuracy of the numerical solution increases. A 3D plot of the absolute error, between the estimated and exact solution of Equation (15), is presented for comparison in Figure 3, which presents the reliability between both results at the level of 10 6 . Note that when t = x is substituted in the solution Equation (19), the absolute error can be observed of the order of 10 6 in one-dimension. The absolute error between the solutions, both approximate and exact, is found to be in good agreement.
Example 3: 
Consider another partial fractional differential equation with an initial condition depending on the generalized fractional-order sine function
2 d γ U ( x , t ) d t γ + d γ U ( x , t ) d x γ = 0 .  
A numerical solution is pursued using the initial condition U ( x , 0 ) = f   ( x ) = s i n γ ( x γ ) in intervals 0 x 1   &   0 t 1 . The generalized definitions of sine and cosine function [3] are given
s i n γ ( x γ ) = k = 0 ( 1 ) k   x 2 k Γ ( 2 k γ + 1 ) ,   and   c o s γ ( x γ ) = k = 0 ( 1 ) k   x ( 2 k + 1 ) Γ ( 2 k γ + γ + 1 ) .
In this example, the assumed approximate solution contains an initial condition that has a generalized sine function as 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 Equation (20) is U e x a c t ( x , t ) = s i n γ ( x γ )   c o s γ ( t γ / 2 ) + c o s γ ( x γ )   s i n γ ( t γ / 2 ) . The approximated solution U a p p ( x , t ) = i , j = 0 n b j i   B j ( α , t )   B i ( α , x ) + s i n γ ( x γ ) is substituted into the Equation (20) and the result is given below:
2 d γ d t γ ( i , j = 0 n b j i   B j ( α , t )   B i ( α , x ) + s i n γ ( x γ ) ) + d γ d x γ ( i , j = 0 n b j i   B j ( α , t )   B i ( α , x ) + s i n γ ( x γ ) ) = 0 .
The Caputo’s derivative is applied to the above expression and the product of fractional B-polys B m ( α , x )   B n ( α , t ) from the basis, sets are multiplied on both sides of the Equation (21). The integration on both variables (x and t) is carried out over the intervals, respectively, and after further simplification, Equation (21) may be written in the form
i , j = 0 n b j i [ 2 B i ( α , x ) | B m ( α , x ) D t γ B j ( α , t ) | B n ( α , t ) D x γ B i ( α , x ) | B m ( α , x )   B j ( α , t ) | B n ( α , t ) ] = f γ ( α , x ) | B m ( α , x ) | B n ( α , t ) .    
The fractional-order derivative of the sine function is f γ ( γ , x ) = d γ d x γ ( s i n γ ( x γ ) ) = c o s γ ( x γ ) with α = γ . The technique leads to a system of ( n + 1 )   ×   ( n + 1 ) operational matrix. This system of equations may be summarized in the matrix equation of the form X   B = W , where the elements of matrix B = { b 1 1 , b 2 1 ,   b 3 1 ,   , b 1 2 , b 2 2 ,   b 3 2 ,   ,   } are the unknown constants. After further simplification, the right-hand side column matrix W and the matrix elements of operational matrix X in terms of inner products of B-polys are given as
X m , n =   i , j = 0 n [ 2   B i ( α , x ) | B m ( α , x ) D t γ B j ( α , t ) | B n ( α , t ) D x γ B i ( α , x ) | B m ( α , x )   B j ( α , t ) | B n ( α , t ) ] , W m , n =   f γ ( γ , x ) | B m ( α , x ) | B n ( α , t ) = 0 R , T c o s γ ( x γ )   B m ( α , x ) B n ( α , t ) d x   d t .
To construct an appropriate solution to Equation (20), the partial fractional order differential Equation (20) is transformed into a matrix equation X   B = W . By deleting rows and the corresponding columns of the Equation (23), the initial conditions are imposed on the operational matrix equation, X 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 inverse of the matrix X is multiplied by the column matrix W and, by solving the matrix equation B = X 1 W , the values of the unknown coefficients b j i are determined. The resulting approximate solution is composed of the product of the B-poly basis set and the expansion coefficients, as in Equation (3). The technique provides the approximate solution U a p p ( x , t ) to Equation (20) using 16 B-polys with n = 15, fractional-order α = 1 2 and fractional differential-order of γ = 1 2 . The approximate solution is provided,
U a p p ( x , t ) = 2.06 × 10 8 t 7 + t 13 / 2 ( 1.277 × 10 7 + 1.09 × 10 7 x ) + t 6 ( 1.783 × 10 7 + 7.839 × 10 7 x ) + t 11 / 2 ( 1.411 × 10 6 + 9.554 × 10 7 x 3.475 × 10 6 x ) + t 5 ( 3.806 × 10 7 7.588 × 10 6 x 4.062 × 10 6 x + 1.254 × 10 5 x 3 / 2 ) + t 9 / 2 ( 3.766 × 10 5 + 1.931 × 10 6 x + 3.083 × 10 5 x + 1.401 × 10 5 x 3 / 2 ) + t 4 ( 2.423 × 10 7 + 1.854 × 10 4 x 7.464 × 10 6 x 1.012 × 10 4 x 3 / 2 4.062 × 10 5 x 2 ) + t 7 / 2 ( 6.715 × 10 4 + 1.12 × 10 6 x 6.779 × 10 4 x + 2.317 × 10 5 x 3 / 2 + 2.775 × 10 4 x 2 + 1.009 × 10 4 x 5 / 2 ) + t 3 ( 4.38 × 10 8 2.938 × 10 3 x 3.848 × 10 6 x + 1.977 × 10 3 x 3 / 2 5.971 × 10 5 x 2 6.475 × 10 4 x 5 / 2 2.167 × 10 4 x 3 + 4.587 × 10 4 x 7 / 2 ) + t 5 / 2 ( 9.403 × 10 3 + 1.779 × 10 7 x + 9.402 × 10 3 x + 1.045 × 10 5 x 3 / 2 4.746 × 10 3 x 2 + 1.298 × 10 4 x 5 / 2 + 1.295 × 10 3 x 3 + 4.035 × 10 4 x 7 / 2 8.027 × 10 4 x 4 ) + t 2 ( 0.0353 x 5.2405 × 10 7 x 0.0235 x 3 / 2 2.309 × 10 5 x 2 + 9.491 × 10 3 x 5 / 2 2.388 × 10 4 x 3 2.22 × 10 3 x 7 / 2 6.5 × 10 4 x 4 + 1.223 × 10 3 x 9 / 2 3.858 × 10 4 x 5 ) + t 3 / 2 ( 0.0940 0.0940 x + 1.186 × 10 6 x 3 / 2 + 0.047 x 2 + 4.181 × 10 5 x 5 / 2 0.0158 x 3 + 3.707 × 10 4 x 7 / 2 + 0.0032 x 4 + 8.968 × 10 4 x 9 / 2 1.605 × 10 3 x 5 + 4.839 × 10 4 x 11 / 2 1.148 × 10 5 x 6 1.151 × 10 5 x 13 / 2 ) + t ( 0.2821 x 1.417 × 10 8 x + 0.1881 x 3 / 2 2.096 × 10 6 x 2 0.0752 x 5 / 2 6.157 × 10 5 x 3 + 0.0217 x 7 / 2 4.777 × 10 4 x 4 3.947 × 10 3 x 9 / 2 1.04 × 10 3 x 5 + 1.779 × 10 3 x 11 / 2 5.144 × 10 4 x 6 + 1.174 × 10 5 x 13 / 2 + 1.136 × 10 5 x 7 ) + t ( 0.5642 + 0.5642 x + 2.406 × 10 8 x 3 / 2 0.2821 x 2 + 2.847 × 10 6 x 5 / 2 + 0.0940 x 3 + 7.168 × 10 5 x 7 / 2 0.0237 x 4 + 4.943 × × 10 4 x 9 / 2 + 3.885 × 10 3 x 5 + 9.783 × 10 4 x 11 / 2 1.605 × 10 3 x 6 + 4.466 × 10 4 x 13 / 2 9.84 × 10 6 x 7 9.205 × 10 6 x 15 / 2 ) + x ( 1.128 × 10 6 0.752 x + 0.301 x 2 0.0860 x 3 + 0.0191 x 4 3.474 × 10 3 x 5 + 5.344 × 10 4 x 6 7.125 × 10 5 x 7 + 8.383 × 10 6 x 8 ) .
From the above result, it is noted that the desired approximate solution is converged and reached the desired accuracy. To find the solution to fractional-order partial differential Equation (20), we used the same fractional-order B-poly basis set as in Example 2. With the higher number (n) of fractional B-polys, a higher order of accuracy is achievable at the expense of computer CPU time. A 3D plot of the estimated and exact results of Equation (20) is presented in Figure 4 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 the Equation (24), the absolute error can be observed at the same level as 10 7 exhibiting the great aspect of constancy in one-dimension x.
It is further noted that, from the traditional trigonometric identity, we know that
s i n ( x + t ) = s i n ( x ) c o s ( t ) + c o s ( x ) s i n ( t ) .
However, in Ref. [39] the authors state that, in fractional calculus, this kind of trigonometry identity does not hold. In this example, we have computationally proven that the above identity is no longer valid in fractional calculus, i.e.,
s i n γ ( x γ + ( t γ 2 ) ) s i n γ ( x γ )   c o s γ ( t γ 2 ) + c o s γ ( x γ )   s i n γ ( t γ 2 ) .
For further verification, we have plotted both sides of the identity Equation (26) and Figure 5 and they seem to disagree. For example, for α = 1 2   and   γ = 1 2 , we show the graphs of both sides of the identity at x = t, sin γ ( x γ / 2 ) (blue curve) and s i n γ ( x γ )   c o s γ ( t γ 2 ) + c o s γ ( x γ )   s i n γ ( t γ 2 ) (yellow curve). From the graphs of both sides of the identity, we see that the blue and the yellow curves do not match. However, we know that when γ takes integral values, both curves overlap. We tried different γ and n values of B-polys; these curves still did not overlap. It is concluded that, in fractional calculus, certain traditional trigonometry identities may not be valid.
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 d γ U ( x , t ) d t γ + d γ U ( x , t ) d x γ = 0 .  
A numerical solution is sought using initial condition U ( x , 0 ) = f   ( x ) = c o s γ ( x γ ) [3], in the intervals 0 x 1   and   0 t 1 . The assumed approximate solution contains an initial condition that has a generalized cosine function. The exact solution of Equation (27) is U e x a c t ( x , t ) = ( c o s γ ( x γ )   c o s γ ( t γ / 2 ) s i n γ ( x γ )   s i n γ ( t γ / 2 ) ) . The assumed solution, U a p p ( x , t ) = i , j = 0 n b j i   B j ( α , t )   B i ( α , x ) + c o s γ ( x γ ) , is substituted into Equation (27) and the result is presented below:
2 d γ d t γ ( i , j = 0 n b j i   B j ( α , t )   B i ( α , x ) + c o s γ ( x γ ) ) + d γ d x γ ( i , j = 0 n b j i   B j ( α , t )   B i ( α , x ) + c o s γ ( x γ ) ) = 0 .
After further simplification and applying the Caputo’s derivative, we obtain
i , j = 0 n b j i [ 2 B i ( α , x ) | B m ( α , x ) D t γ B j ( α , t ) | B n ( α , t ) D x γ B i ( α , x ) | B m ( α , x )   B j ( α , t ) | B n ( α , t ) ] = f γ ( α , x ) | B m ( α , x ) | B n ( α , t ) ,
where f γ ( γ , x ) = d γ d x γ ( c o s γ ( x γ ) ) = s i n γ ( x γ ) with α = γ . The current technique leads to a system of ( n + 1 ) ×   ( n + 1 ) equations. This system of equations may be summarized in the matrix equation of the form X   B = W , where the elements of matrix B = { b 1 1 , b 2 1 ,   b 3 1 ,   , b 1 2 , b 2 2 ,   b 3 2 ,   ,   } are the unknown constants. The matrix elements of the column matrix W , and operational matrix X are given as
X m , n = i , j = 0 n [ 2   B i ( α , x ) | B m ( α , x ) D t γ B j ( α , t ) | B n ( α , t ) D x γ B i ( α , x ) | B m ( α , x )   B j ( α , t ) | B n ( α , t ) ] , W m , n =   f γ ( γ , x ) | B m ( α , x ) | B n ( α , t ) = 0 R , T s i n γ ( x γ ) B m ( α , x ) B n ( α , t ) d x   d t .
The partial fractional-order differential Equation (27) is now converted into an operational matrix equation X   B = W . By deleting the rows and corresponding columns of Equation (30), the initial condition is imposed on the operational matrix equation X , 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 B = X 1 W to yield the values of the unknown coefficients b j i . The resulting approximate solution is composed of the product of the expansion coefficients and the B-poly basis set, as given in Equation (3). The technique provides the approximate solution U a p p ( x , t ) of Equation (27) using n = 15 B-polys of fractional-order α = 1 2   and fractional differential-order γ = 1 2 .
U a p p ( x , t ) = 1.0 1.0 x + 0.5 x 2 0.166 x 3 + 0.0417 x 4 8.3 × 10 3 x 5 + 1.39 × 10 3 x 6 1.98 × 10 4 x 7 + 2.48 × 10 5 x 8 2.75 × 10 6 x 9 + 2.75 × 10 7 x 10 + t ( 0.25 + 0.25 x + 8.1 × 10 8 x 3 / 2 0.125 x 2 + 5.74 × 10 6 x 5 / 2 + 0.0416 x 3 + 1 × 10 4 x 7 / 2 0.0106 x 4 + 5.18 × 10 4 x 9 / 2 + 1.32 × 10 3 x 5 + 8.1 × 10 4 x 11 / 2 9.7 × 10 4 x 6 + 3.17 × 10 4 x 13 / 2 3.71 × 10 5 x 7 ) + t 2 ( 0.0312 + 1.52 × 10 8 x 0.0312 x + 1.79 × 10 6 x 3 / 2 + 0.0156 x 2 + 4.37 × 10 5 x 5 / 2 5.34 × 10 3 x 3 + 2.91 × 10 4 x 7 / 2 + 8.2 × 10 4 x 4 + 5.61 × 10 4 x 9 / 2 7.29 × 10 4 x 5 + 2.57 × 10 4 x 11 / 2 3.24 × 10 5 x 6 ) + t 7 / 2 ( 5.31 × 10 8 7.58 × 10 4 x 3.12 × 10 6 x + 5.18 × 10 4 x 3 / 2 3.64 × 10 5 x 2 1.28 × 10 4 x 5 / 2 ) + t 9 / 2 ( 1.87 × 10 7 + 4.32 × 10 5 x 4.04 × 10 6 x 1.78 × 10 5 x 3 / 2 ) + t 11 / 2 ( 2.15 × 10 7 1.21 × 10 6 x 1.59 × 10 6 x ) + t 6 ( 2.45 × 10 7 + 3.59 × 10 7 x 9.5 × 10 7 x ) + t 15 / 2 ( 1.83 × 10 9 x ) + t 7 ( 2.52 × 10 8 + 3.24 × 10 8 x ) + t 13 / 2 ( 8 × 10 8 + 2.06 × 10 7 x ) + t 5 ( 8.3 × 10 6 + 9.9 × 10 7 x + 5.18 × 10 6 x ) + t 4 ( 1.62 × 10 4 + 8.5 × 10 7 x 1.67 × 10 4 x + 1.32 × 10 5 x 3 / 2 + 5.18 × 10 5 x 2 ) + t 3 ( 0.0026 + 2.24 × 10 7 x + 2.6 × 10 3 x + 9.1 × 10 6 x 3 / 2 1.33 × 10 3 x 2 + 8.5 × 10 5 x 5 / 2 + 2.76 × 10 4 x 3 ) + t 5 / 2 ( 0.0106 x 7.18 × 10 7 x 7.07 × 10 3 x 3 / 2 2.18 × 10 5 x 2 + 2.9 × 10 3 x 5 / 2 1.7 × 10 4 x 3 5.14 × 10 4 x 7 / 2 3.68 × 10 4 x 4 + 5.03 × 10 4 x 9 / 2 ) + t 3 / 2 ( 0.106 x 4.06 × 10 8 x + 0.0707 x 3 / 2 3.59 × 10 6 x 2 0.0283 x 5 / 2 7.29 × 10 5 x 3 + 8.2 × 10 3 x 7 / 2 4.25 × 10 4 x 4 1.14 × 10 3 x 9 / 2 7.36 × 10 4 x 5 + 9.1 × 10 4 x 11 / 2 3.1 × 10 4 x 6 + 3.76 × 10 5 x 13 / 2 ) + t ( 0.637 x 0.424 x 3 / 2 1.22 × 10 7 x 2 + 0.169 x 5 / 2 7.18 × 10 6 x 3 0.0485 x 7 / 2 1.09 × 10 4 x 4 + 0.011 x 9 / 2 5.1 × 10 4 x 5 1.24 × 10 3 x 11 / 2 7.36 × 10 4 x 6 + 8.4 × 10 4 x 13 / 2 2.65 × 10 4 x 7 + 3.0 × 10 5 x 15 / 2 ) .
From the above result of Equation (31), 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 remained the same at the desired level of 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, a higher order of accuracy is attainable. In Table 2, we compare our calculated U a p p ( x , t ) values with the exact values U e x a c t ( x , t ) at various points for x and t. The absolute differences between the solutions are provided in the last column of Table 2, which shows excellent agreement. A 3D plot of the estimated and exact results of Equation (31) is presented in Figure 6 for the purpose of comparison. Note that when t = x is substituted in Equation (31), the absolute error in one dimension also goes to 10 7 . The absolute error in the 3D graph is also presented in Figure 6 showing that the error in the converged solution is of the order of 10 7 .
From the traditional trigonometric rule, we know that the following is a valid identity
c o s ( x + t ) = c o s ( x )   c o s ( t ) s i n ( x )   s i n ( t )
However, in Ref. [39], the authors state that, in fractional calculus, this trigonometry identity may not be true. In this example, we have computationally proven that the above identity is no longer valid in fractional calculus, i.e.,
c o s γ ( x γ + ( t γ 2 ) ) c o s γ ( x γ )   c o s γ ( t γ 2 ) s i n γ ( x γ )   s i n γ ( t γ 2 ) .
For further verification, we have plotted both sides of the identity Equation (33) and they seem to disagree as shown in Figure 7. For example, for α = 1 2   and   γ = 1 2 , we show the graphs of both sides of the identity at x = t,  c o s γ ( x γ / 2 ) (blue curve) and ( c o s γ x γ   c o s γ ( t γ 2 ) s i n γ x γ   s i n γ ( t γ 2 ) ) (yellow curve). The graphs of both sides of the identity show that the blue and the yellow curves 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 γ , these curves still did not overlap. It is concluded that, in fractional calculus, this trigonometry identity may not be valid.

5. Error Analysis

We performed the calculations in the absence of a grid to solve linear fractional partial differential based on fractional B-polys. The fractional-order B-polys basis sets are defined on the intervals x [ 0 ,   1 ] and t [ 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 an error analysis based on the increasing number of B-poly basis sets; it is noted that the accuracy improves. The absolute error analysis for Example 4 is presented for the exact and approximate results. As you may have seen in Example 4, in the final calculation, we used the number k = 15 in the summation of the generalized formula for c o s γ ( x γ / 2 ) = k = 0 n ( 1 ) k   x ( 2 k + 1 ) Γ ( 2 k γ + γ + 1 ) and also used n = 15 for the B-poly basis set in both 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 performed the calculations in Example 4; it is observed that the absolute error among solutions is of the order of 10−3. Next, we used n = 10, which would give us 11 B-poly sets. The absolute error among solutions reduces to the level of 10−6. Finally, we use n = 15, which would comprise 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 × 256-dimensional operational matrix, which is already a large matrix to invert. We had to increase the accuracy of the program to handle this matrix in the Mathematica symbolic program. Beyond these limits, it becomes problematic to find an accurate inversion of the matrix. 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 (Figure 8 and Figure 9) that the absolute error decreases as we steadily increase the size of the fractional B-poly basis set. Due to the analytic nature of the fractional B-polys, all the calculations are carried out without a grid representation on the intervals of integration. We also presented the absolute error in terms of 3D graphs in Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9. Clearly, the error is systematically decreased as the number of B-polys basis sets is increased in the calculations. The method provides 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.

6. Results and Discussions

In the current study, we investigated the 2D modified fractional Bhatti-polys basis set technique to determine the solutions to the partial fractional differential equations. Four examples of the linear partial fractional differential equations have been presented with various initial conditions and their semi-analytic solutions are also provided. Furthermore, an explanation of the 2D fractional algorithm process has been provided to calculate approximate solutions of the linear fractional differential equation. We estimated the results using the Galerkin method [40] in both variables (x, t). The graphs of the converged solutions are provided in Figure 2, Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7. In the first example, the estimated solution was precise, which was equivalent to the exact solution after ignoring the tiny contributions. As the number of fractional B-polynomial basis set in the approximate solutions to Equation (3) was increased, the accuracy of the numerical solutions [36] 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 the approximated results of the absolute error in Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9. In every case, the accuracy of the solutions was different because different B-poly basis sets and different operational matrix sizes were used. The numerical efficiency of the inverted matrix depends on the size of the matrix. In the last three examples, we used the series representation of the generalized sine and cosine functions; this requires the inclusion of many terms in the summation. When variable t was equal to x for 1D error analysis, the absolute errors among approximate and exact results were examined. The precision appears to be the same in both 1D and 3D error analyses. It is concluded that the present technique performed well in resolving linear fractional-order differential equations utilizing operational matrix scheme [40,41], as exhibited by the graphs and data shown in the study. We performed all integrations analytically and performed computations using Wolfram Mathematica symbolic program version-12 [42] for both x and t variables over the 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. Recently, many authors [43,44] have 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 showed a detailed error investigation for the fourth problem that can be applied to other examples. The CPU time for computing the first example was less than 1 min, while examples 2–4 took 5–30 min of CPU time since those required a larger B-ploys basis set and higher dimensions of the operational matrix.
In this paper, we presented an expanded form of this technique [36] 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 there are no known solutions. We may explore this method’s potential to solve 2D nonlinear partial fractional-order differential equations in forthcoming publications.

Author Contributions

Conceptualization, M.I.B.; software, M.I.B.; methodology, M.I.B.; formal analysis, M.I.B. and M.H.R.; validation, M.I.B., and M.H.R.; data curation, M.I.B. and M.H.R.; resources, M.I.B.; project administration, M.I.B.; investigation, M.I.B. and M.H.R.; supervision, M.I.B.; writing—original draft preparation, M.I.B.; writing—review and editing, M.I.B. and M.H.R. We confirm that the authors listed in the manuscript have been approved by all of us. All authors have read and agreed to the published version of the manuscript.

Funding

No external funding has been received for this research.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data that supports the findings are available within the article.

Acknowledgments

The authors would like to thank the Department for helping us to use the departmental computational facilities which have been utilized to execute the results of this study paper.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Miller, K.S.; Ross, B. An Introduction to the Fractional Calculus and Fractional Differential Equations; John Wiley: Hoboken, NJ, USA, 1993. [Google Scholar]
  2. Kilbas, A.; Srivastava, H.M.; Trujillo, J.J. Theory and Applications of Fractional Differential Equations; North-Holland Mathematics Studies, 204, Fractional Calculus and Applied Analysis; Elsevier: Amsterdam, The Netherlands, 2006. [Google Scholar]
  3. Igor, P. An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution. In Fractional Differential Equations; Mathematics in Science and Engineering Series; Academic Press: Cambridge, MA, USA, 1998; Volume 198. [Google Scholar]
  4. Bagley, R.L.; Calico, R.A. Fractional order state equations for the control of viscoelasticallydamped structures. J. Guid. Control Dyn. 1991, 14, 304–311. [Google Scholar] [CrossRef]
  5. Sun, H.; Abdelwahab, A.; Onaral, B. Linear approximation of transfer function with a pole of fractional power. IEEE Trans. Autom. Control 1984, 29, 441–444. [Google Scholar] [CrossRef]
  6. Ichise, M.; Nagayanagi, Y.; Kojima, T. An analog simulation of non-integer order transfer functions for analysis of electrode processes. J. Electroanal. Chem. 1971, 33, 253–265. [Google Scholar] [CrossRef]
  7. Laskin, N. Fractional market dynamics. Phys. A Stat. Mech. Appl. 2000, 287, 482–492. [Google Scholar] [CrossRef]
  8. Kusnezov, D.; Bulgac, A.; Dang, G.D. Quantum Lévy Processes and Fractional Kinetics. Phys. Rev. Lett. 1999, 82, 1136–1139. [Google Scholar] [CrossRef] [Green Version]
  9. Li, C.; Zhang, F. A survey on the stability of fractional differential equations. Eur. Phys. J. Spec. Top. 2011, 193, 27–47. [Google Scholar] [CrossRef]
  10. Tavazoei, M.S.; Haeri, M. A note on the stability of fractional order systems. Math. Comput. Simul. 2009, 79, 1566–1576. [Google Scholar] [CrossRef]
  11. Torvik, P.J.; Bagley, R.L. On the Appearance of the Fractional Derivative in the Behavior of Real Materials. J. Appl. Mech. Trans. ASME 1984, 51, 294–298. [Google Scholar] [CrossRef]
  12. Jawad, A.J.M.; Petković, M.D.; Biswas, A. Modified simple equation method for nonlinear evolution equations. Appl. Math. Comput. 2010, 217, 869–877. [Google Scholar]
  13. Zayed, E.M.E. A note on the modified simple equation method applied to Sharma–Tasso–Olver equation. Appl. Math. Comput. 2011, 218, 3962–3964. [Google Scholar] [CrossRef]
  14. Wei, Y.; Guo, Y.; Li, Y. A new numerical method for solving semilinear fractional differential equation. J. Appl. Math. Comput. 2021, 1–23. [Google Scholar] [CrossRef]
  15. Mahor, T.C.; Mishra, R.; Jain, R. Analytical solutions of linear fractional partial differential equations using fractional Fourier transform. J. Comput. Appl. Math. 2021, 385, 113202. [Google Scholar] [CrossRef]
  16. Inc, M. The approximate and exact solutions of the space- and time-fractional Burgers equations with initial conditions by variational iteration method. J. Math. Anal. Appl. 2008, 345, 476–484. [Google Scholar] [CrossRef] [Green Version]
  17. Roman, P. Application of the Adams-Bashfort-Mowlton Method to the Numerical Study of Linear Fractional Oscillators Models. In AIP Conference Proceedings; AIP Publishing LLC.: Melville, NY, USA, 2021. [Google Scholar]
  18. Srivastava, H.M.; Alomari, A.-K.N.; Saad, K.M.; Hamanah, W.M. Some Dynamical Models Involving Fractional-Order Derivatives with the Mittag-Leffler Type Kernels and Their Applications Based upon the Legendre Spectral Collocation Method. Fractal Fract. 2021, 5, 131. [Google Scholar] [CrossRef]
  19. Guy, J. Lagrange characteristic method for solving a class of nonlinear partial differential equations of fractional order. Appl. Math. Lett. 2006, 19, 873–880. [Google Scholar] [CrossRef] [Green Version]
  20. Safari, M.; Ganji, D.D.; Moslemi, M. Application of He’s Variational Iteration Method and Adomian’s Decomposition Method to the Fractional KdV-Burgers-Kuramoto Equation. Comput. Math. Appl. 2009, 58, 2091–2097. [Google Scholar] [CrossRef] [Green Version]
  21. Cui, M. Compact finite difference method for the fractional diffusion equation. J. Comput. Phys. 2009, 228, 7792–7804. [Google Scholar] [CrossRef]
  22. Momani, S.; Odibat, Z.M. Fractional Green Function for Linear Time-Fractional Inhomogeneous Partial Differential Equations in Fluid Mechanics. J. Appl. Math. Comput. 2007, 24, 167. [Google Scholar] [CrossRef]
  23. Huang, Q.; Huang, G.; Zhan, H. A Finite Element Solution for the Fractional Advection-Dispersion Equation. Adv. Water Resour. 2008, 31, 1578–1589. [Google Scholar] [CrossRef]
  24. Guo, S.; Mei, L.; Li, Y.; Sun, Y. The Improved Fractional Sub-Equation Method and Its Applications to the Space-Time Fractional Differential Equations in Fluid Mechanics. Phys. Lett. Sect. A Gen. At. Solid State Phys. 2012, 376, 407–411. [Google Scholar] [CrossRef]
  25. Zheng, B. (G′/G)-Expansion Method for Solving Fractional Partial Differential Equations in the Theory of Mathematical Physics. Commun. Theor. Phys. 2012, 58, 623–630. [Google Scholar] [CrossRef]
  26. Lu, B. The first integral method for some time fractional differential equations. J. Math. Anal. Appl. 2012, 395, 684–693. [Google Scholar] [CrossRef] [Green Version]
  27. Peykrayegan, N.; Ghovatmand, M.; Skandari, M.H.N. On the convergence of Jacobi-Gauss collocation method for linear fractional delay differential equations. Math. Methods Appl. Sci. 2021, 44, 2237–2253. [Google Scholar] [CrossRef]
  28. Toh, Y.T.; Phang, C. Collocation Method for Fractional Differential Equation via Laguerre Polynomials with Eigenvalue Degree. In AIP Conference Proceedings; AIP Publishing LLC.: Melville, NY, USA, 2021. [Google Scholar]
  29. Khader, M.M.; Saad, K.M.; Baleanu, D.; Kumar, S. A spectral collocation method for fractional chemical clock reactions. Comput. Appl. Math. 2020, 39, 394. [Google Scholar] [CrossRef]
  30. Su, W.-H.; Yang, X.-J.; Jafari, H.; Baleanu, D. Fractional complex transform method for wave equations on Cantor sets within local fractional differential operator. Adv. Differ. Equ. 2013, 2013, 97. [Google Scholar] [CrossRef] [Green Version]
  31. Bhatti, M.I. Solution of Fractional Harmonic Oscillator in a Fractional B-poly Basis. Phys. Tech. Sci. 2014, 2, 8. [Google Scholar] [CrossRef]
  32. Bhatti, M.; Bracken, P.; Dimakis, N.; Herrera, A. Solution of mathematical model for gas solubility using fractional-order Bhatti polynomials. J. Phys. Commun. 2018, 2, 085013. [Google Scholar] [CrossRef]
  33. Bhatti, M.I. Analytic Matrix Elements of the Schrödinger Equation. Adv. Stud. Theor. Phys. 2013, 7, 11. [Google Scholar] [CrossRef]
  34. I Bhatti, M.; Hinojosa, E. 26 Results of hyperbolic partial differential equations in B-poly basis. J. Phys. Commun. 2020, 4, 095010. [Google Scholar] [CrossRef]
  35. Bhatti, M.I.; Rahman, H.; Dimakis, N. Approximate Solutions of Nonlinear Partial Differential Equations Using B-Polynomial Bases. Fractal Fract. 2021, 5, 106. [Google Scholar] [CrossRef]
  36. Bhatti, M.I.; Bracken, P. Solutions of differential equations in a Bernstein polynomial basis. J. Comput. Appl. Math. 2007, 205, 272–280. [Google Scholar] [CrossRef] [Green Version]
  37. Jong, S.-G.; Choe, H.-C.; Ri, Y.-D. A New Approach for an Analytical Solution for a System of Multi-term Linear Fractional Differential Equations. Iran. J. Sci. Technol. Trans. A Sci. 2021, 45, 955–996. [Google Scholar] [CrossRef]
  38. Bhatti, M.I. Solutions of the Harmonic Oscillator Equation in a B-Polynomial Basis. Adv. Stud. Theor. Phys. 2009, 3, 451–460. [Google Scholar]
  39. Massopust, P.R.; Zayed, A.I. On the Invalidity of Fourier Series Expansions of Fractional Order. Fract. Calc. Appl. Anal. 2015, 18, 1507–1517. [Google Scholar] [CrossRef] [Green Version]
  40. Fletcher, C.A.J. Computational Galerkin Methods; Springer: Berlin/Heidelberg, Germany, 1984. [Google Scholar]
  41. Alinhac, S. Hyperbolic Partial Differential Equations; Springer: Berlin/Heidelberg, Germany, 2009. [Google Scholar]
  42. Wolfram Research, Inc. Mathematica 12.0; Wolfram Research, Inc.: Champaign, IL, USA, 2019. [Google Scholar]
  43. Singh, A.K.; Singh, V.K.; Singh, O.P. The Bernstein Operational Matrix of Integration. Appl. Math. Sci. 2009, 3, 2427–2436. [Google Scholar]
  44. Farouki, R.T. Legendre–Bernstein basis transformations. J. Comput. Appl. Math. 2000, 119, 145–160. [Google Scholar] [CrossRef] [Green Version]
Figure 1. (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 B i ( 1 2 , x ) ,   B i ( 5 3 , x ) ,   B i ( 9 4 , x ) , and x are dimensionless quantities.
Figure 1. (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 B i ( 1 2 , x ) ,   B i ( 5 3 , x ) ,   B i ( 9 4 , x ) , and x are dimensionless quantities.
Fractalfract 05 00208 g001
Figure 2. 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, Equation (14). The 1D plot of the absolute error between approximate and exact results is also presented in the intervals t [ 0 ,   1 ] and x [ 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.
Figure 2. 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, Equation (14). The 1D plot of the absolute error between approximate and exact results is also presented in the intervals t [ 0 ,   1 ] and x [ 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.
Fractalfract 05 00208 g002
Figure 3. 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 t = x is substituted in Equation (19). On the right-hand side, a 3D plot of the absolute error between exact and estimated solutions is also presented in the intervals x [ 0 ,   1 ] and t [ 0 ,   1 ] . Both graphs show that the efficiency of the numerical solutions is of the order of 10 6 . This kind of accuracy was observed when n = 15 and α = 1 2 order polynomial basis set was used.
Figure 3. 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 t = x is substituted in Equation (19). On the right-hand side, a 3D plot of the absolute error between exact and estimated solutions is also presented in the intervals x [ 0 ,   1 ] and t [ 0 ,   1 ] . Both graphs show that the efficiency of the numerical solutions is of the order of 10 6 . This kind of accuracy was observed when n = 15 and α = 1 2 order polynomial basis set was used.
Fractalfract 05 00208 g003
Figure 4. 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 Equation (24). The plot shows that 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 t [ 0 ,   1 ] and x [ 0 ,   1 ] . The figure represents the efficacy of the numerical solutions is of the order of 10 7 . This accuracy occurred with n = 15 number of fractional order B-poly basis set in x variable and the same basis set was used in t variable.
Figure 4. 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 Equation (24). The plot shows that 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 t [ 0 ,   1 ] and x [ 0 ,   1 ] . The figure represents the efficacy of the numerical solutions is of the order of 10 7 . This accuracy occurred with n = 15 number of fractional order B-poly basis set in x variable and the same basis set was used in t variable.
Fractalfract 05 00208 g004aFractalfract 05 00208 g004b
Figure 5. Two graphs of the identity Equation (26) are presented to show that both sides of the identity do not agree for fractional values of γ . The left side of the identity is sin γ ( x γ / 2 ) and the right side of the identity is s i n γ ( x γ )   c o s γ ( x γ 2 ) + c o s γ ( x γ )   s i n γ ( x γ 2 ) at t = x. The values for α = 1 2   and γ = 1 2 are used. The blue curve and the yellow curve do not agree or overlap. Hence, in general, the identity is invalid when fractional calculus is considered.
Figure 5. Two graphs of the identity Equation (26) are presented to show that both sides of the identity do not agree for fractional values of γ . The left side of the identity is sin γ ( x γ / 2 ) and the right side of the identity is s i n γ ( x γ )   c o s γ ( x γ 2 ) + c o s γ ( x γ )   s i n γ ( x γ 2 ) at t = x. The values for α = 1 2   and γ = 1 2 are used. The blue curve and the yellow curve do not agree or overlap. Hence, in general, the identity is invalid when fractional calculus is considered.
Fractalfract 05 00208 g005
Figure 6. A plot of the absolute error between approximate (fx) and exact (sol) solutions is introduced on the left-hand for t = x, Equation (31). The one-dimensional graph shows that 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 t [ 0 ,   1 ] and x [ 0 ,   1 ] . The figure represents the effectiveness of the numerical solution is of the order of 10 7 .
Figure 6. A plot of the absolute error between approximate (fx) and exact (sol) solutions is introduced on the left-hand for t = x, Equation (31). The one-dimensional graph shows that 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 t [ 0 ,   1 ] and x [ 0 ,   1 ] . The figure represents the effectiveness of the numerical solution is of the order of 10 7 .
Fractalfract 05 00208 g006
Figure 7. 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 c o s γ ( x γ / 2 ) and the right side of the identity is ( c o s γ ( x γ )   c o s γ ( x γ 2 ) s i n γ ( x γ )   s i n γ ( x γ 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.
Figure 7. 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 c o s γ ( x γ / 2 ) and the right side of the identity is ( c o s γ ( x γ )   c o s γ ( x γ 2 ) s i n γ ( x γ )   s i n γ ( x γ 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.
Fractalfract 05 00208 g007
Figure 8. 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 3D graph. The estimated result is not converged. The values for α = 1 2   and   γ = 1 2 are used.
Figure 8. 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 3D graph. The estimated result is not converged. The values for α = 1 2   and   γ = 1 2 are used.
Fractalfract 05 00208 g008
Figure 9. 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 3D graph. The estimated result is not converged. The values for α = 1 2   and   γ = 1 2 are used.
Figure 9. 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 3D graph. The estimated result is not converged. The values for α = 1 2   and   γ = 1 2 are used.
Fractalfract 05 00208 g009
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.
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.
α γ nBasis SetCaputo’s Derivative of Basis Set (Equation (2))
1/21/21 { 1 x , x } { π 2 , π 2 }
3/43/41 { 1 x 3 / 4 , x 3 / 4 } { Γ ( 7 4 ) ,   Γ ( 7 4 ) }
5/35/31 { 1 x 5 / 3 , x 5 / 3 } { Γ ( 8 3 ) ,   Γ ( 8 3 ) }
5/45/41 { 1 x 5 / 4 , x 5 / 4 } { Γ ( 9 4 ) ,   Γ ( 9 4 ) }
9/49/41 { 1 x 9 / 4 , x 9 / 4 } { Γ ( 13 4 ) ,   Γ ( 13 4 ) }
9/59/51 { 1 x 9 / 5 , x 9 / 5 } { Γ ( 14 5 ) ,   Γ ( 14 5 ) }
Table 2. For different values of x and t, we compare our calculated approximated U a p p ( x , t ) value with the exact value U e x a c t ( x , t ) . We also provide their absolute difference.
Table 2. For different values of x and t, we compare our calculated approximated U a p p ( x , t ) value with the exact value U e x a c t ( x , t ) . We also provide their absolute difference.
xt U e x a c t ( x , t ) U a p p ( x , t ) Absolute Difference
| U e x a c t ( x , t ) U a p p ( x , t ) |
0.10.10.9410970.9410971.036 × 10−11
0.20.20.8867840.8867841.091 × 10−10
0.30.30.8366740.8366745.120 × 10−10
0.40.40.7904200.7904201.660 × 10−9
0.50.50.7476680.7476684.333 × 10−9
0.60.60.7081520.7081529.785 × 10−9
0.70.70.6715930.6715931.991 × 10−8
0.80.80.6377420.6377423.744 × 10−8
0.90.90.6063760.6063766.617 × 10−8
1.01.00.5772880.5772881.112 × 10−7
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Bhatti, M.I.; Rahman, M.H. Technique to Solve Linear Fractional Differential Equations Using B-Polynomials Bases. Fractal Fract. 2021, 5, 208. https://doi.org/10.3390/fractalfract5040208

AMA Style

Bhatti MI, Rahman MH. Technique to Solve Linear Fractional Differential Equations Using B-Polynomials Bases. Fractal and Fractional. 2021; 5(4):208. https://doi.org/10.3390/fractalfract5040208

Chicago/Turabian Style

Bhatti, Muhammad I., and Md. Habibur Rahman. 2021. "Technique to Solve Linear Fractional Differential Equations Using B-Polynomials Bases" Fractal and Fractional 5, no. 4: 208. https://doi.org/10.3390/fractalfract5040208

Article Metrics

Back to TopTop