New High-Order Compact ADI Algorithms for 3 D Nonlinear Time-Fractional Convection-Diffusion Equation

Numerical approximations of the three-dimensional (3D) nonlinear time-fractional convection-diffusion equation is studied, which is firstly transformed to a time-fractional diffusion equation and then is solved by linearization method combined with alternating direction implicit (ADI) method. By using fourth-order Padé approximation for spatial derivatives and classical backward differentiationmethod for time derivative, two new high-order compact ADI algorithms with ordersO(τmin(1+α,2−α) +h) andO(τ+h) are presented.The resulting schemes in each ADI solution step corresponding to a tridiagonal matrix equation can be solved by theThomas algorithm which makes the computation cost effective. Numerical experiments are shown to demonstrate the high accuracy and robustness of two new schemes.


Introduction
In this paper, we consider numerical approximations of the following 3D nonlinear time-fractional convection-diffusion equation: Ω×(0,] = (Δ −   −   −   +  +  () + )     Ω×(0,] , with the initial and boundary conditions  (, , , 0)    Ω =  0 (, , )    Ω ,  (, , , )    Ω×(0,] = 0, where 0 <  < 1, , , , and  are constants, and () is some reasonable nonlinear function of  = (, , , ).Ω is a continuous convex domain in 3D space with suitable boundary conditions prescribed on its boundary Ω.The solution  and the source terms () and  = (, , , ) are assumed to be sufficiently smooth and have the necessary continuous partial derivatives up to certain orders.The fractional derivative   /  is Grünward-Letnikov fractional derivative defined by (, , , )    ( − )  . ( It should be noted that Grünward-Letnikov fractional derivative will become Caputo fractional derivative when  0 = 0. Fractional differential equations (FDEs) have become popular in recent years both in mathematics and in applications.They have been increasingly used in physical and chemical processes [1][2][3], relaxation in polymer systems, dynamics of protein molecules, and the diffusion of contaminants in complex geological formations [4].Because of the absence or appropriate mathematical methods, most of the analytical solutions for FDEs are difficult to obtain.However, several numerical methods have been proposed to solve them.
For the one-dimensional (1D) time FDEs, the schemes in [5][6][7][8][9][10] are convergent with order not higher than two for the space variable, and it is interesting to seek highorder numerical methods for FDEs.Recently, such problem was solved by the compact finite difference scheme with convergence order ( + ℎ 4 ) in [11], and a higher order

Two New High-Order Compact ADI Algorithms
For simplicity of the presentation, we consider a rectangular domain Ω = [0, 1] 3 , which is divided into an ×× mesh with the spatial step size ℎ  = ℎ  = ℎ  = ℎ = 1/, while the time step size is chosen as  = / and  is the number of time coordinate direction.
To obtain an efficient scheme, it is important to approximate the nonlinear term () in (1) explicitly.Otherwise, one would need to perform a nonlinear iteration in each time step which is quite time consuming.Let   = (, , ,   ) and   = (, , ,   ).The following linear method is used to transform () into linear term: Next we discuss the linear equation for 1 ≤  ≤  − 1.
In order to keep the fourth-order accuracy in space and the tridiagonal nature of the scheme by utilizing the Padé approximation, we first introduce a transformation that is similar to that of Liao [35] to eliminate the convection terms in (1). Let where V = V(, , , ) and  1 (),  2 (), and  3 () are functions to be determined.Substituting (6) into the first formula of (1), we have Setting the coefficients of V  , V  , and V  to zero leads to Then, we have Set  (, , ) =  1 ()  2 ()  3 () =  (/2)+(/2)+(/2) .(10) Combining ( 1) and ( 5)-( 10), we obtain the following timefractional diffusion equation satisfied by V: In the following, we will introduce a high-order compact difference method for (11).
The time-fractional derivative   V/  at  +1 is estimated by where   = ( + 1) 1− −  1− and the truncation error  +1  has the following form: If we expend V(, , , ) in Taylor series about  at  +1/2 in (13) and use the analysis in [14], we will get where  V is a positive constant only depending on V.
As for the spatial derivatives, we can use the fourthorder Padé scheme.That is, the spatial derivatives can be approximated by where  2  ,  2  , and  2  are the second-order central difference operators along and and -directions, respectively.
Combining with ( 12) and ( 15), then the implicit compact finite difference scheme for (11) is given as follows: with the initial discretization and the boundary discretization where V   denotes V(ℎ, ℎ, ℎ, ).In the following, V  will be written in short for V   if there is no confusion about the notation.By reformulating ( 16), we obtain an equivalent form where Due to the nature of three-time-level scheme, we need an additional method to compute the numerical solution at the first time step.
For the special case  = 0, the scheme simply reads with It is clear that ( 21) must be solved using some relaxation procedure because  1 is nonlinear.
For simplicity of the presentation, we write ( 19) and ( 21) in a unified form Multiplying ) on both sides of (23), we have In the following, two new high-order algorithms are proposed for (24) that will coincide with the ADI methods employed.First, some finite difference operations are introduced that will be used hereafter as follows: Algorithm 1.
(1) For  ≥ 1, add splitting term to the left-hand side of ( 24), and rewrite it into an equivalent form (2) Introducing two intermediate variables V * and V * * , scheme SI can be solved in three steps as Note that the intermediate values of V * * and V * on the boundary are easily obtained from (28c) and (28b).

Algorithm 2.
(1) For the first time step, V 1 can be obtained from Algorithm 1.
(2) For  ≥ 2, add splitting term to the left-hand side of ( 24), and rewrite it into an equivalent form (3) Use the same way as the second step of Algorithm 1.
Remark 3. It is clear that the splitting error of II is higher order in  than the local error for the time-fractional derivative approximation, so, for the time accuracy, scheme SII can ensure the local truncation error to be 2 −  order.However, as for the splitting term I, it should be 1+ ≥ 2− in order to maintain the time accuracy.That is, the order of local error for scheme SI is determined by the following equation: min (1 + , 2 − ) = { 1 + ,  ∈ (0, 0.5] , 2 − ,  ∈ (0.5, 1) .
Based on the previous analysis, we propose two linear ADI schemes, and the truncation errors are ( min(1+,2−) + ℎ 4 ) and ( (2−) +ℎ 4 ), respectively.The resulting ADI schemes in each ADI solution step corresponding to a strictly diagonally dominant matrix equation which can be solved using the onedimensional tridiagonal Thomas algorithm with a considerable saving in computing time.
Remark 4. The previous analysis shows that splitting term II is better than I when  ∈ (0, 0.5].However, when  ∈ (0.5, 1), the splitting error is essentially negligible if the time step size is small enough, and so both schemes can obtain the same accuracy in this case.Nevertheless, numerical experiments in Section 3 show that scheme SII is much better than schem SI, regardless of the value of .We would like to comment the fact that if the splitting error is much larger than the truncation error, as pointed out by Douglas and Kim [36], a good splitting term will play an important role in the accuracy of the solution, especially for the analytical solutions that decrease exponentially in time.So, we think that scheme SII is a better choice even with an additional computational cost than scheme SI.

Numerical Experiments
In this section, two numerical examples are presented to demonstrate the high efficiency and accuracy of the new method.The numerical results got from Algorithm 2 will be compared with Algorithm 1.The  ∞ -norm error (denoted as  ∞ -error) and the  2 -norm error (denoted as  2 -error) of the numerical solutions are shown.Note that the treatment of space is the same for two schemes all of them can achieve the same results when  is small enough to eliminate the influence of the temporal approximation.So, in the following experiments, we only choose Algorithm 2 to investigate the spatial convergence rate.Moreover, we compute the convergent order of accuracy for both schemes.Let us consider two mesh sizes Δ and Δℎ on Ω  and Ω ℎ , respectively.The  2 -norm errors of these two grids are denoted as ‖ ⋅ ‖  0 and ‖ ⋅ ‖ ℎ 0 .If we set the convergent order of accuracy to be Rate, then we have the following form: So, the convergent order of accuracy Rate can be computed as The order of accuracy is formally defined when the mesh size approaches to zero.Therefore, when the mesh size is relatively large, the discretization scheme may not achieve its formal order of accuracy.Problem 5. We firstly consider a pure diffusion problem.The exact analytical solution and the corresponding forcing term and initial condition are given by      = Δ +  2 +  (, , , ) ,  (, , , ) =   sin () sin () sin () ,  0 (, , ) = sin () sin () sin () .(34) Since it is difficult to accurately solve   /  based on  except for  = 0.5, in this problem, we only discuss the case with  = 0.5 and the corresponding We first investigate the time accuracy with  = 1 for both algorithms.To this end,  = 40 is chosen such that the errors stemming from the spatial approximation are negligible.Numerical results with  = 0.5 are presented in Table 1.It is clearly show that Algorithm 2 is much better than Algorithm 1.For instance, in Table 1, the maximum absolute error of the computed solution by Algorithm 2 is around 2.2764 − 3 when  = 40,  = 50.Whereas for Algorithm 1 with the same mesh size, the maximum absolute error is around 9.9234 − 2. Similar comparison can be made to reach similar conclusions.Meanwhile, we can see that the convergence rate for Algorithm 1 is around 1.5, which agrees with the theoretical estimates min(1 + , 2 − ).Howevre, for Algorithm 2, the  2 -rate is approximately 2.6, which is superconvergent.For better comparison of two algorithms, we plot the  ∞ errors (as shown in Figure 1) at  = 1 and  = 0.5 with  = 40 and  = 40 for  = 0.5.Comparing both subplots, we can see that Algorithm 2 has substantially higher accuracy than Algorithm 1.
For the spatial convergence rate, we choose  = 0.01 and  = 100 to eliminate the influence of the temporal approximation.Numerical results using Algorithm 2 are presented in Table 2.The convergence order is indeed 4.
Numerical results for time with  = 0.2, 0.5, 0.8 are presented in Tables 3-5, respectively.Again, as expected, it is shown that Algorithm 2 has higher accuracy in comparison with Algorithm 1.For instance, in Table 5, the maximum absolute error of the computed solution by Algorithm 2 is around 8.1067 − 4 when  = 20,  = 50.Whereas for Algorithm 1 with the same mesh size, the maximum absolute error is around 6.0692 − 2.
However, we note that the convergence orders for Algorithm 1 with  = 0.2 and 0.5 are unsatisfying, as shown in Tables 3 and 4, respectively.We would like to comment on the fact that the order of accuracy is formally defined when the mesh size approaches to zero.Therefore, when the mesh size is relatively large, the discretization scheme may not achieve its formal order of accuracy.Moreover, the data show that the order is increasing now; hence, it will achieve its formal order 1 +  with the decrease of the time step .
Figure 3 shows the  ∞ errors at  = 2 and  = 0.5 with  = 50 and  = 40 for different .Again, as expected, it indicates that Algorithm 2 has substantially higher accuracy than Algorithm 1 regardless of .
Table 6 presents space results with  = 0.01 and  = 100 for different  and , from which we can find that Algorithm 2 is stable and convergent and it can obtain fourthorder accuracy when mesh becomes finer, regardless of the value of .
Figure 4 shows the contour plot of the exact solution in the - plane and numerical solutions using Algorithm 2 for different  at  = 0.5 and  = 0.01 with  = 0.5 and  = 100.We can see that when  = 20, the numerical solution is a good approximation of the exact solution.

Conclusions
In this paper, we have proposed two high-order compact ADI algorithms with order ( min(1+,2−) + ℎ 4 ) and ( 2− + ℎ 4 ), respectively, for solving the 3D nonlinear time-fractional convection-diffusion equation.Since it only involves threepoint stencil for each one-dimensional operator, it can be solved by applying the one-dimensional tridiagonal Thomas algorithm with a considerable saving in computing time.Numerical experiments have shown that Algorithm 2 is very significantly reduced in the splitting perturbation error and is superconvergent, which is much better than Algorithm 1.
Compared with existing deducing methods, our approach is advantageous because (1) the original time-fractional convection-diffusion equation is first transformed to a time-fractional diffusion equation, and the resulting time-fractional diffusion problem is then solved by ADI method, this is novel and much more simpler than other traditional methods; (2) the nonlinear term is transformed into a linear term which makes the computation cost effective and the approximation solutions more accurate; (3) we show that the critical modification in the ADI procedures can virtually eliminate the perturbation errors and produce identical approximations of the solution of the differential problem.

Table 6 :
= 0.01 and  = 100 with different  and  using Algorithm 2 for Problem 6.