Fractional Bernstein Series Solution of Fractional Diffusion Equations with Error Estimate

In the present paper, we introduce the fractional Bernstein series solution (FBSS) to solve the fractional diffusion equation, which is a generalization of the classical diffusion equation. The Bernstein polynomial method is a promising one and can be generalized to more complicated problems in fractional partial differential equations. To get the FBSS, we first convert all terms in the problem to matrix forms. Then, the fundamental matrix equation is obtained and thus, the solution is obtained. Two error estimation methods based on a residual correction procedure and the consecutive approximations are incorporated to find the estimate and bound of the absolute error. The perturbation and stability analysis of the method is given. We apply the method to some illustrative examples. The numerical results are compared with the exact solutions and known second-order methods. The outcomes of the numerical examples are very encouraging and show that the FBSS is highly useful in solving fractional partial problems. The results show the accuracy and effectiveness of the method.


Introduction
Fractional derivatives have been used to model many problems in science, e.g., physics [1][2][3], medicine [4], hydrology [5], biomedical problems [6], dynamics of particles [7] and applied sciences [8]. The linear fractional diffusion equation is considered for scientists and engineers [9]. Fractional space derivatives are used in the modeling of anomalous diffusion. In a diffusion model, replacement of the second derivative by a fractional derivative causes enhanced diffusion, also called superdiffusion [10].
In the present study, we consider the following one-dimensional fractional diffusion equation (FDE), on a finite domain L < x < R, 1 < α ≤ 2 and d(x) > 0. We assume the initial condition c(x, t = 0) = F(x) for L < x < R and the boundary conditions c(x = L, t) = 0 and y(x = R, t) = b R (t).
The fractional derivative in Equation (1) is the Caputo fractional derivative of order α [11]. The basic property of the Caputo derivative is as follows.
One effective method to solve the problems is the Bernstein polynomials method (BPM), also called the Bernstein series solution method [20][21][22][23][24][25]. Positive linear operators of Bernstein with sequence of these operators are established to solve to solve two-variables equations [26]. The operational matrices of the BPM have been used to get the numerical solutions of a class of third-order ordinary differential equations [27]. The multi-stage BPM, a generalization of the standard BPM, was applied to fractional-order stiff systems [28]. The BPM with new modifications was employed to solve fractional differential equations [29]. Moreover, the same method was used to solve some types of ordinary differential equations [30][31][32]. A special type of the singular Emden-Fowler problems was solved by BPM [33]. Recently, some linear and non-linear systems of ordinary differential equations have been solved by BPM and the accuracy has been improved by a residual correction procedure [34]. A novel error estimation method for the parametric non-intrusive reduced order model based on machine learning is presented in [35].
In this paper, the fractional Bernstein series solution (FBSS) method is introduced and applied to solve Equation (1) numerically. The method comprises the fractional Bernstein polynomials and collocation method. We approximate the exact solution of Equation (1) by f α p n,n such that f α p n,n satisfies Equation (1) on the collocation nodes.
This paper is organized as follows. Section 2 presents necessary definitions and theorems. Section 3, the main section, discusses the matrix forms of the solution f α p n,n with its derivatives. Then, we get the fundamental matrix equation of the FDE. Employing the conditions and applying the Gauss elimination procedure yields the unknown matrix. Thus, we obtain the FBSS. In Section 4, two different error analysis methods are presented to get the upper bound of the absolute error with the corrected FBSS. Moreover, an upper bound obtained by the generalized Taylor theorem is presented. To decide the stability, the perturbation and stability analysis of the method is done in Section 5. Section 6 provides the numerical results to illustrate the FBSS method for different n values. We compare the method with other methods used to solve the problem. The results show that the present method gives accurate solutions and is also efficient. We test the effect of small perturbations to the approximate solutions both theoretically and practically for various values of n. Section 7 summarizes the results.

Preliminaries and Notations
Bernstein polynomials of n-th degree are given by the following equation [29].
Similarly, the fractional Bernstein polynomials B α k,n (x) are constituted by using x → x α , so Equation (2) becomes We will use the definition of Caputo's derivative, which is the modification of the definition of Riemann-Liouville. It has some advantages for solvinginitial value problems [36].
Definition 1 ([37,38]). The Riemann-Liouville integral operator of order α > 0 for a ≥ 0 is defined as follows Definition 2 ([37,38]). The Caputo fractional dervative of f of order α > 0 for a ≥ 0 is given as follows, We will use the following multivariate fractional Taylor's theorem [39] to bound the absolute error.

Numerical Method
Let be the FBSS of Equation (1). Let us find the matrix forms of f α p α n,n = ∂ α f α p n,n ∂x α , 1 < α ≤ 2 and ∂ f α p n,n ∂t .
First note that f α p n,n can be written as follows: We write f α B n (x) as Therefore, we can write f α p n,n and f α p α n,n in the forms respectively. As f α B α n (x) can be written as is obtained where Axioms 2021, 10, 6

of 19
Substituting Equation (8) into Equation (7) we get Note that Q n (t) can be given as whereȲ The term ∂Q n (t) ∂t can be written as: where Substituting Equations (9) and (13) into Equation (5) yields the matrix forms of f α p n,n (x, t), f α p α n,n (x, t) and ∂ f α p n,n ∂t as respectively. The use of (15) and (16) in (1) gives the fundamental matrix equation Inserting the collocation points For the conditions, we first put t = 0, x = L, and x = R into Equation (14), respectively. Then, we obtain the following matrix relations by substituting the collocation nodes, Applying the Gauss elimination method to the augmented matrix [W,G] and deleting the zero rows gives the system [W,Ḡ], whereW is square matrix and has full rank. Then, the unknown matrix A can be obtained as

Error Analysis
In this section, first an upper bound for the absolute error is given. Then, we give two error estimation methods that can be applied easily and are useful practically. The first one is the residual correction procedure. The second one is specifying the consecutive approximations which is similar to the error analysis of the RK4 method.
Theorem 2. Let f α p n,n and y be the FBSS and the exact solution of Equation (1), respectively. By using the above notations, the absolute error is bounded as follows, provided that D kα f ∈ C(D), where D is a rectangular domain contains the collocation nodes.
Proof. Adding and subtracting the term P α n (x, t), the two-dimensional truncated generalized Taylor polynomial defined in (4), into the left-hand side and applying triangle inequality yields the desired result.
Let R be the function defined as Adding the term R into both sides of Equation (1) yields the following fractional differential equation for the absolute error where e n,n = y − f α p n,n . The initial and boundary conditions for the problem are converted to e n,n (L, t) = e n,n (R, t) = 0.
We get an approximate solution, e m,m n,n , for the absolute error by applying the method to Equation (21) with the conditions (22) on the nodes (x i , t j ) : 0 ≤ i, j ≤ m ⊂ Ω. Then, the absolute error e n,n can be estimated by e m,m n,n provided that e n,n − e m,m n,n < ε is small.
Corollary 1. Let f α p n,n be the FBSS. Then, f α p n,n + e m,m n,n is another approximate solution, corrected FBSS, of Equation (1) and its error function is e n,n − e m,m n,n . Moreover, if e n,n − e m,m n,n < y − e n,n , then f α p n,n + e m,m n,n is a better approximation than p n,n in any given norm · .
Let f α p n,n and f α p s,s be any two FBSS of Equation (1), y be the exact solution of Equation (1). Then, by using the triangle inequality, we find the following inequality, If e m,m < e n,n , then we can write as follows e n,n = C e m,m , C > 1.
Therefore, we can bound the error as • Therefore, we can bound the error by Then, we can bound the error e m,m well in case of C ≥ 2, even when the exact solution is not known. A similar argument may be proposed when the error sequence is decreasing (or increasing). Thus, in case of a decreasing (or increasing) error sequence, one of the following bounds is satisfied

Numerical Stability
In this section, we will study the perturbation analysis with stability estimation of the linear systems obtained by the FBSS method, which is similar to that in [40], for a given problem. Perturbing the initial or boundary conditions yields the following perturbed solutions, f α p per n,n (x, t). There occur two cases: Case 1: Only the vectorḠ on the right hand side of (20) is perturbed, i.e., Let us show the perturbed solution of (24) as A p = A+∆A, where ∆A is the perturbation of the solution resulting from the perturbations in the initial and boundary conditions. Then, the change in the solution caused by the initial change is bounded as [41] ∆A A ≤ cond(W) ∆Ḡ Ḡ .

Case 2:
The perturbations might be occurred bothW andḠ, As the same notation in Case 1, the change in A caused by perturbing the initials is bounded above as [41] ∆A Thus, for Case 1, and therefor we can specify the effect of the little changes in the initial and boundary conditions on the FBSS by measuring cond(W). We omit the result for Case 2 by simplicity.

Numerical Results and Discussion
In this section, three examples are provided to illustrate the properties and effectiveness of the technique.

Example 1
Let us consider the FDE [17] ∂u(x, t) where 0 < x < 1 and α = 1.8, the diffusion coefficient is 2) x 2.8 6 and the source function The initial condition is and the boundary conditions are The exact solution to this problem is By applying the procedure in Section 3, the matrix equation for Equation (27) is found as The collocation points that will be used are the Chebyshev interpolation nodes By substituting the collocation points in Equation (28), we will obtain W matrix. The conditions matrices for (x, 0) = x 3 , u(0, t) = 0, u(1, t) = e −t are obtained as Then, [W,G] is obtained by combining these matrices. Thus, we obtain the coefficient matrix A. As we perform the method for different n values, we obtain n different FBSS. All calculations are done using the Maple program. The results for n = 5, 10, n = 15, and n = 20 are given in Table 1. The absolute errors are graphed in Figure 1. Table 2 show the absolute errors between consecutive approximations. The exact and approximate solutions with n = 15 are graphed in Figure 2. The values of error functions and the exact solution, at time t = 1, are given in Table 3 and graphed in Figure 3. The upper bounds of absolute error at time t = 1 are given in Figure 4. A comparison of the method with the numerical methods in [16,18] is made for t = 1. The absolute error for n = 5 with the errors obtained by the corrected FBSSs are given in Table 4. We can say that increasing n yields better approximation results and the method provides better results. From Table 2, we deduce that the absolute error e n,n can be estimated approximately by the difference between e n,n and e n+1,n+1 . From Table 3, we can see the method gives more accurate results than the methods in [16,18] for t = 1. From Table 4, increasing m yields more accurate corrected approximate solutions. Table 5 shows the error norms (L 2 , L ∞ ) resulting with CPU time (in seconds) used in the Maple program to find the numerical solutions for different n values. We provide the stability results for the method for some n values in Table 6. Increasing n gives the condition numbers which are increasing. It can be said that the approximate solutions will be more stable around n = 10 by using both the theoretical upper bounds and the numerical results. The method produces approximately 10 5 as a condition number. For n 10, the solutions will be more sensitive to small variations since the number of conditions is large. Thus, around n = 10, we can say that the method works well for this problem.  Table 2. Upper bounds for the absolute errors for n = 9 and n = 10 of Example 1.

Example 2
Let us consider the FDE [17] ∂u(x, t) where 0 < x < 1 and α = 1.4. The diffusion coefficient is The initial and boundary conditions are as follows, respectively, The exact solution of this problem is given by We perform the method for n = 7, 10 and n = 15. The results are given in Table 7. The absolute errors are graphed in Figure 5. A comparison of the method with the numerical method in [17] is given in Table 7. We can say that increasing n yields better approximation results and the method provides better results. The method gives more accurate results than the method in [17]. From Table 8, we conclude that the absolute error e n,n can be bounded approximately by |e n,n − e n+1,n+1 |. The absolute error for n = 5 with the errors obtained by the corrected FBSSs are given in Table 9. The absolute errors for n = 11 and n = 16 at time t = 1 are given in Figure 6. The errors with their upper bounds at t = 1 are given in Figure 7. Table 10 shows the error norms (L 2 , L ∞ ), resulting with CPU time (in seconds) used in the Maple program to find FBSSs for different n values. The method again gives approximately 10 5 as a condition number. Therefore, around n = 10, we can say from Table 11 that the method produces more stable approximations for this problem.

Example 3
Let us consider the fractional convection diffusion equation [42] ∂u(x, t) ∂t where 0 < x < 1 and. The diffusion coefficient is the source function g(x, t) = x 3 cosh(t).
The initial and boundary conditions are as follows, respectively, The exact solution of this problem is given by      By applying the technique in Section 3, for different values of n, the FBSSs for the problem are founded. From Tables 12 and 13, we can see obviously that the results obtained by FBSS and the exact solution of the problem are agreement with each other. We compare our results with [42], can see our results are more accurate. The condition numbers for the systems that related to the method for various n are given in Table 14. We can say from Table 14 that the method is insensitive to small variations for n = 10.    Table 12. Comparisons of the absolute errors for n = 3, 5, 7 of Example 3.

Conclusions
In this paper, we proposed the FBSS to solve the FDE numerically. The method can be easily implemented and is effective. It comprises fractional Bernstein polynomials and the collocation method. First, the fundamental matrix equation is obtained and then it is solved by the Gauss elimination procedure. Applying the initial and boundary conditions, we get the solution for the given n value. Two error estimations, residual correction procedure and estimations are obtained by the difference of the consecutive approximations. We applied the method to some examples. Generally, the results demonstrate that the proposed method achieves better approximation accuracy than some other well-known methods.
From the examples, we can bound approximately the absolute error e n,n by |e n,n − e n+1,n+1 |. On the other hand, one can obtain more accurate results by using the residual correction procedure. Increasing m gives more accurate corrected FBSS. For the stability of the method, we can decide which size n yields more stable approximations by specifying at the number of conditions of the related system. For the examples, the method produces condition numbers approximately as 10 5 for n ≈ 10 approximately. Thus, around n = 10, we can say that the method is suitable for these examples. The subjects of our future works can be exemplified by applying fractional Bernstein series solution for solving fractional integral differential equations and chaotic fractional order systems.