Non-polynomial B-spline and shifted Jacobi spectral collocation techniques to solve time-fractional nonlinear coupled Burgers’ equations numerically

This paper proposes two numerical approaches for solving the coupled nonlinear time-fractional Burgers’ equations with initial or boundary conditions on the interval [0, L]. The first method is the non-polynomial B-spline method based on L1-approximation and the finite difference approximations for spatial derivatives. The method has been shown to be unconditionally stable by using the Von-Neumann technique. The second method is the shifted Jacobi spectral collocation method based on an operational matrix of fractional derivatives. The proposed algorithms’ main feature is that when solving the original problem it is converted into a nonlinear system of algebraic equations. The efficiency of these methods is demonstrated by applying several examples in time-fractional coupled Burgers equations. The error norms and figures show the effectiveness and reasonable accuracy of the proposed methods.


Introduction
Many phenomena in engineering and applied sciences can be represented successfully using fractional calculus models such as anomalous diffusion, materials, and mechanics, signal processing, biological systems, finance, etc. (see, for instance, [1][2][3][4][5][6][7]). There is a tremendous interest in fractional differential equations, as the theory of fractional derivatives itself and its applications have been intensively developed. The theory of fractional differential equations describes the reality of life more powerfully and systematically. In recent years, several researchers have studied differential equations of fractional order through diverse techniques [8][9][10][11].
The time-fractional Burgers' equation is important since it is a kind of sub-diffusion convection equation. Several different methods for solving the equation have been developed such as the local fractional Riccati differential equation method in [12], the homotopy analysis transform method [13], the finite difference method [14], the variational iteration method [15]. The study of coupled Burgers equations is significant for t-dimensional. The system is a simple sedimentation model or the evolution of scaled volume concentrations of two types of fluid suspension or colloid particles under the influence of gravity. Many powerful methods had been developed to find analytic or numerical solutions of coupled Burgers' equations such as homotopy perturbation method [16], differential transformation method [17], non-polynomial spline method [18], septic b-spline collocation method [19], Galerkin quadratic b-spline method [20], Adomian decomposition method [21], meshless radial point interpolation method [22].
Several vital analytical and numerical techniques have been proposed to solve the coupled nonlinear time-fractional Burger equations NLTFBEs. Prakash et al. [23] suggested an analytical algorithm based on the homotopy perturbation Sumudu transform method to investigate the coupled NLFBEs. Hoda et al. [24] introduced the Laplace-Adomian decomposition method, the Laplace variational iteration method, and the reduced differential transformation method for solving the one-dimensional and two-dimensional fractional coupled Burgers' equations. In [25] Liu and Hou explicitly applied the generalized two-dimensional differential transform method to solve the coupled space-and time-fractional Burgers equations (STFBEs). Heydari and Avazzadeh [26] proposed an effective numerical method based on Hahn polynomials to solve the nonsingular variableorder time-fractional coupled Burgers' equations. The authors in [27] suggested a hybrid spectral exponential Chebyshev approach based on a spectral collection method to solve the coupled TFBEs. Veeresha and Prakasha [28] and Singh et al. [29] presented the qhomotopy analysis transform method to solve the coupled TFBEs and STFBEs, respectively. The coupled STFBEs have also been solved using the Adomian decomposition method by Chan and An [30]. Islam and Akbar [31] obtained exact general solutions of the coupled STFBEs by using the generalized (G /G)-expansion method with the assistance of the complex fractional transformation. Prakash et al. [32] proposed the fractional variational iteration method to solve the same problem. Hussein [33] proposed a continuous and discrete-time weak Galerkin finite element approach for solving two-dimensional time-fractional coupled Burgers' equations. Hussain et al. [34] obtained the numerical solutions of the coupled TFBEs using the meshfree spectral method.
In comparison with local methods, spectral methods are often efficient and highly accurate systems. Convergence speed is one of the spectral methods' great advantages. Furthermore, spectral methods have a high level of precision (often so-called "exponential convergence"). The primary idea of all spectral method versions is to express the solution to the problem as a finite sum of certain foundation functions and then to choose the coefficients, to minimize the difference between precise and approximate solutions. Recently, the classical spectral methods have been developed to obtain accurate solutions for linear and nonlinear FDEs. Spectral approaches with the assistance of operational matrices of orthogonal polynomials have been considered to approximate the solution of FDEs (see, for example, [35][36][37][38][39]).
One of the disadvantages of the non-polynomial method is that the time step size must be small enough. The main advantage of the proposed methods is that they are easy to implement. Also, the solutions can be obtained with high accuracy with relatively fewer spatial grid nodes compared with other numerical techniques. For this reason and because the current methods can be directly applied to other applications, we are motivated to apply these techniques for coupled Burgers equations.
In this paper we develop two accurate numerical methods to approximate the numerical solutions of the coupled TFBEs. The first method is the non-polynomial B-spline method [8,[40][41][42] based on the L1-approximation and finite difference approximations for spatial derivatives. The second method is the shifted Jacobi spectral collocation method [43][44][45] with the assistance of the operational matrix of fractional and integer-order derivatives. The collocation approach proposed in this paper is somewhat different from those collocation methods commonly discussed in the literature. Now, we consider the time-fractional coupled Burgers' equations of the form subject to the initial and boundary conditions where α 1 and α 2 are parameters describing fractional derivatives, x and t are the space and time variables, respectively. u and v are the velocity components to be determined. f (x, t) and g(x, t) are continuous functions on their domains. The functions p(x), q(x), f 1 (t), f 2 (t), g 1 (t), g 2 (t) are sufficient smooth functions. The fractional derivatives of order α 1 and α 2 in Eqs. (1) and (2) are treated in the sense of Liouville-Caputo defined by Jerome and Oldham [6]. In the case of α 1 = α 2 = 1, Eqs. (1) and (2) are reduced to the classical coupled Burgers equations.

Definition 1 ([1])
A real function u(t), t > 0, is said to be in the space C Ω , Ω ∈ R if there exists a real number p > Ω such that u(t) = t p u 1 (t), where u 1 (t) ∈ C(0, ∞), and it is said to be in the space C n Ω if and only if u (n) ∈ C Ω , n ∈ N.

Non-polynomial B-spline method
In this section, we take a spline function of the form: where ω is the frequency of the hyperbolic part of spline functions which will be used to raise the accuracy of the method.

Convergence analysis
According to Remark 2, we have chosen 2γ + β = h 2 , where γ = h 2 12 and β = 5h 2 6 . Let us rewrite Eqs. (21) and (22) as follows: T and a matrix Q is given as a block matrix where and square matrices Q 0 , Q 1 , and Z x , x = δ, η, ξ , ζ are given by (1) and (2) at nodal points x i , i = 1, 2, . . . , N , and then we have where T = [T 1i , T 2i ] T is the local truncation error described in Remark 2. From Eqs. (25) and (27), we can write the error equation as follows: For sufficiently small step h, the diagonal blocks Q 11 and Q 22 are invertible and the following condition holds: According to [48], matrix Q is invertible. Moreover, From Eq. (28) and norm inequalities, we have Since and from classifications of the matrices Q 11 , Q 12 , Q 21 , and Q 22 defined in Eq. (26), we have This shows that Eq. (25) is a second-order convergence method in the case 2γ + β = h 2 , γ = h 2 /12, and a fourth-order convergence method in the case 2γ + β = h 2 , γ = h 2 /12.

Stability analysis of the method
The stability analysis of the difference schemes listed in Eqs. (19) to (22) is discussed by assuming the nonlinear terms δ n r and η n r , r = i -1, i, i + 1, as local constants D and E respectively.
LetŨ n i andṼ n i be the approximate solutions of Eqs. (19) to (22) and define With the above definition and regarding Eqs. (19) and (21), we can get the following roundoff error equations: where , and a 8 = -βσ 1 ϕ α 1 n-1 . The Von Neumann method assumes that where I = √ -1.
The analytic form of shifted Jacobi polynomial P (μ,η) L,j (x) is given by The values of the shifted Jacobi polynomials at the boundary points are given by Suppose that f (x) is a square-integrable function with respect to the shifted Jacobi weight function ω as follows: where The shifted Jacobi-Gauss quadrature is used to approximate the previous integral as follows: where x Now, if we approximate f (x) in Eq. (39) by the first (N + 1)-terms, then which can be approximated by the first (M + 1) × (N + 1) terms with the truncation error with the coefficient matrix A given by a 00 a 01 · · · a 0N a 10 a 11 · · · a 1N . . .

Theorem 1 ([39])
Let τ ,M (t) be shifted Jacobi vector defined in Eq. (42), and let α > 0. Then where D α is the (M + 1) × (M + 1) operational matrix of derivatives of order α in the Liouville-Caputo sense and is defined by where α is the floor function and Similarly, the fractional derivative α of L,N (x) can be expressed as in Eq. (45): where D α is the (N + 1) × (N + 1) operational matrix of derivatives of order α.

Time-fractional coupled Burgers' equation
We are going to consider the time-fractional coupled Burgers' Eqs. (1) and (2), which may be written as follows: with the initial-boundary conditions Using Eq. (43), we approximate u(x, t), v(x, t), f (x, t), and g(x, t) as follows: where U and V are unknown coefficients (M + 1) × (N + 1) matrices, while A and B are defined by Eq. (43), where Using Theorem 1, we have

∂v(x, t) ∂x
Substituting Eqs. (49) and (50) in Eqs. (47) and (48) and collocating at (M + 1) × (N -1) points, we have where x j , j = 0, 1, . . . , N -2, are the roots of shifted Jacobi polynomial P System (51) can be combined with Eq. (52) to form the system of 2(M + 1)(N + 1) nonlinear algebraic equations in 2(M + 1)(N + 1) undefined coefficients, which can be resolved using Newton's iterative approach. Consequently, u M,N (x, t) and v M,N (x, t) can be calculated by the formulae given in Eq. (49). The convergence and error analysis of the shifted Jacobi polynomial has been considered in [50,51]. The Caputo fractional derivative of the shifted Jacobi polynomials satisfies the following estimate: where C is a positive generic constant and q = max{μ, η, -1/2}. [50,Theorem 3] proved that the expansion coefficient a i,j in Eq. (43) satisfies the following estimate: Finally, by [51,Theorem 3], we find that the truncation error of solutions of Eqs. (47) and (48) obtained by shifted Jacobi polynomial satisfies the following estimates:

Conclusion
In this paper, we solved the coupled TFBEs by two different methods. Firstly, we developed the non-polynomial B-spline method based on L1-formula to approximate the Liouville-Caputo time-fractional derivative. The study of stability using the Von Neumann method showed that the scheme is unconditionally stable. Secondly, we applied the shifted Jacobi spectral collocation method based on the operational matrix of fractional derivatives in the Liouville-Caputo sense with the aid of Jacobi-Gauss quadrature. From the tables and figures introduced in Sect. 4, it is clear that the SJSCM is more accurate and stable than the non-polynomial B-spline method for all different values α 1 , α 2 , and N . We also note that the accuracy of the non-polynomial B-spline method increases whenever the value of α 1 and α 2 decreases. The validity is tested by solving three problems of the presented methods. The elicited results confirm the high precision of the methods presented.