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

Authors: Shuying Zhai; Xinlong Feng; Zhifeng Weng
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

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 differentiation method for time derivative, two new high-order compact ADI algorithms with ordersO(τmin(1+α,2−α)+h4) and O(τ2−α+h4) are presented. The resulting schemes in each ADI solution step corresponding to a tridiagonal matrix equation can be solved by the Thomas algorithm which makes the computation cost effective. Numerical experiments are shown to demonstrate the high accuracy and robustness of two new schemes.


---

## Body

## 1. Introduction

In this paper, we consider numerical approximations of the following 3D nonlinear time-fractional convection-diffusion equation:(1)∂αu∂tα|Ω×(0,T]=(Δu-pux-quy-ruz+su+ϕ(u)+f)|Ω×(0,T],
with the initial and boundary conditions
(2)u(x,y,z,0)|Ω=u0(x,y,z)|Ω,u(x,y,z,t)|∂Ω×(0,T]=0,
where 0<α<1, p, q, r, and s are constants, and ϕ(u) is some reasonable nonlinear function of u=u(x,y,z,t). Ω is a continuous convex domain in 3D space with suitable boundary conditions prescribed on its boundary ∂Ω. The solution u and the source terms ϕ(u) and f=f(x,y,z,t) are assumed to be sufficiently smooth and have the necessary continuous partial derivatives up to certain orders. The fractional derivative  ∂αu/∂tα is Grünward-Letnikov fractional derivative defined by
(3)∂αu∂tα=t-αΓ(1-α)u0+1Γ(2-α)∫0t‍∂u(x,y,z,s)∂sds(t-s)α.
It should be noted that Grünward-Letnikov fractional derivative will become Caputo fractional derivative when u0=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–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–10] are convergent with order not higher than two for the space variable, and it is interesting to seek high-order numerical methods for FDEs. Recently, such problem was solved by the compact finite difference scheme with convergence order O(τ+h4) in [11], and a higher order O(τ2-α+h4) one can be found in [12, 13]. Lin and Xu [14], Jiang and Ma [15], and Wei et al. [16–18] proposed high order method based on the finite spectral method, the finite element method (FEM), and local discontinuous Galerkin (LDG) method for space discretization, respectively.For the two-dimensional (2D) FDEs, Zhuang and Liu [19] considered a time-fractional diffusion equation on a finite domain, and a first-order implicit difference approximation was proposed to solve the equation. Additional references [20–22] treat explicit or/and implicit finite difference methods (FDMs) for 2D time FDEs, as well as FEM [23], LDG method [24], and spectral method [25] for such problems. Recently, Zhang and Sun [26] presented new numerical techniques for time FDEs and ADI schemes were considered, but both methods have only second-order accuracy in spatial direction. Then, Cui [27] proposed a high-order ADI compact finite difference scheme with order O(τ2-α+h4) for time-fractional diffusion equation. References [28–34] address FDM and FEM for space FDEs, respectively.In this paper, we focus on the 3D nonlinear time-fractional convection-diffusion equation. To the authors’ knowledge, the research for the 3D case is very lacking. The reason, of necessity, is that 3D problems present many computational challenges because of the prohibitive memory and time requirement often necessary to solve them on grids of sufficient resolution. Traditional numerical schemes have low accuracy and thus require fine discretization. Moreover, being opposite to differential equations of integer order, in which derivatives depend only on the local behaviour of the function, fractional differential equations accumulate the whole information of the function in a weighted form. This makes realistic 3D simulations of FDEs prohibitively expensive. Therefore, it is interesting to discuss high-order numerical methods for 3D time-fractional convection-diffusion equation.It is well known that Padé approximation has the advantages of the fourth-order accuracy to approximate the second-order derivatives and of keeping the desirable tridiagonal nature of the finite difference equations. So, the main purpose of this paper is to solve the nonlinear time-fractional convection-diffusion problem using the Padé approximation combined with ADI method and the one-dimensional tridiagonal Thomas algorithm. However, the convection terms will impede implementation of tridiagonal Thomas algorithm. So, we propose a transformation based on the compact finite difference method [35]. In this method, the original time-fractional convection-diffusion equation is transformed to a time-fractional diffusion equation; the resulting time-fractional diffusion problem is then solved by an efficient high-order method. Another advantage of our method is that the nonlinear term is transformed into a linear term which makes the computation cost effective and the approximation solutions more accurate.This paper is organized as follows. In Section2, two new high-order ADI algorithms are presented for solving 3D nonlinear time-fractional convection-diffusion equation. In Section 3, numerical examples are carried out to verify the high efficiency of our method. Finally, conclusions are drawn in Section 4.


## 2. 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 N×N×N mesh with the spatial step size hx=hy=hz=h=1/N, while the time step size is chosen as τ=T/M and M is the number of time coordinate direction.To obtain an efficient scheme, it is important to approximate the nonlinear termϕ(u) in (1) explicitly. Otherwise, one would need to perform a nonlinear iteration in each time step which is quite time consuming. Let um=u(x,y,z,tm) and fm=f(x,y,z,tm). The following linear method is used to transform ϕ(u) into linear term:
(4)ϕ(um+1)=2ϕ(um)-ϕ(um-1)+O(τ2),1≤m≤M-1.Next we discuss the linear equation(5)∂αum+1∂tα=Δum+1-puxm+1-quym+1-ruzm+1+sum+1+2ϕ(um)-ϕ(um-1)+fm+1
for 1≤m≤M-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(6)u=λ1(x)λ2(y)λ3(z)v,
where v=v(x,y,z,t) and λ1(x),λ2(y), and λ3(z) are functions to be determined. Substituting (6) into the first formula of (1), we have
(7)λ1λ2λ3(∂αvtα-Δv)=(λ1xxλ2λ3-pλ1xλ2λ3+λ1λ2yyλ3-qλ1λ2yλ3+λ1λ2λ3zz-rλ1λ2λ3z)v+(2λ1x-pλ1)λ2λ3vx+(2λ2y-qλ2)λ1λ3vy+(2λ3z-rλ3)λ1λ2vz.
Setting the coefficients of vx,vy, and vz to zero leads to
(8)2λ1x-pλ1=0,2λ2y-qλ2=0,2λ3z-rλ3=0.
Then, we have
(9)λ1(x)=e(p/2)x,λ2(y)=e(q/2)y,λ3(z)=e(r/2)z.Set(10)λ(x,y,z)=λ1(x)λ2(y)λ3(z)=e(p/2)x+(q/2)y+(r/2)z.Combining (1) and (5)–(10), we obtain the following time-fractional diffusion equation satisfied by v:
(11)∂αvm+1∂tα|Ω=Δvm+1+(s-p2+q2+r24)vm+1+λ-1(2ϕ(λvm)-ϕ(λvm-1)+fm+1)|Ω,+λ-1(2ϕ(λvm)-ϕ(λvm-1)+fm1)m≥1,v(x,y,z,0)|Ω=λ-1u0(x,y,z)|Ω,v(x,y,z,tm)|∂Ω=0,m≥0.In the following, we will introduce a high-order compact difference method for (11).The time-fractional derivative∂αv/∂tα at tm+1 is estimated by
(12)∂αvm+1∂tα=tm+1-αΓ(1-α)v0+1Γ(1-α)×∑l=0m‍∫tltl+1‍∂v(x,y,z,ξ)∂ξdξ(tm+1-ξ)α=tm+1-αΓ(1-α)v0+1Γ(1-α)×∑l=0m‍vl+1-vlτ∫tltl+1‍dξ(tm+1-ξ)α+Eτm+1=tm+1-αΓ(1-α)v0+1Γ(2-α)∑l=0m‍blvm+1-l-vm-lτα+Eτm+1,
where bl=(l+1)1-α-l1-α and the truncation error Eτm+1 has the following form:
(13)Eτm+1=1Γ(1-α)×∑l=0m‍∫tltl+1‍(∂v(x,y,z,ξ)∂ξ-vl+1-vlτ)×dξ(tm+1-ξ)α.
If we expend v(x,y,z,ξ) in Taylor series about ξ at tl+1/2 in (13) and use the analysis in [14], we will get
(14)Eτm+1≤Cv(1Γ(1-α)∑l=0m‍∫tltl+1‍tl+1+tl-2ξ(tm+1-ξ)αdξ+O(τ2))≤Cvτ2-α,
where Cv is a positive constant only depending on v.As for the spatial derivatives, we can use the fourth-order Padé scheme. That is, the spatial derivatives can be approximated by(15)Δv=(δx2h2(1+(1/12)δx2)+δy2h2(1+(1/12)δy2)+δz2h2(1+(1/12)δz2))×v+O(h4),
where δx2, δy2, and δz2 are the second-order central difference operators along x- and y- and z-directions, respectively.Combining with (12) and (15), then the implicit compact finite difference scheme for (11) is given as follows:
(16)((m+1)τ)-αΓ(1-α)vijk0+1Γ(2-α)∑l=0m‍blvijkm+1-l-vijkm-lτα=[δx2h2(1+(1/12)δx2)+δy2h2(1+(1/12)δy2)+δz2h2(1+(1/12)δz2)δx2h2(1+(1/12)δx2)+δy2h2(1+(1/12)δy2)]vijkm+1+(s-p2+q2+r24)vijkm+1+λ-1(2ϕijkm-ϕijkm-1+fijkm+1),1≤i,j,k≤N-1,1≤m≤M-1,
with the initial discretization
(17)vijk0=λ(xi,yj,zk)u0(xi,yj,zk),0≤i,j,k≤N,
and the boundary discretization
(18)v0jkm=vNjkm=vi0km=viNkm=vij0m=vijNm=0,0≤i,j,k≤N,0≤m≤M,
where vijkm denotes v(ih,jh,kh,mτ). In the following, vm will be written in short for vijkm if there is no confusion about the notation.By reformulating (16), we obtain an equivalent form
(19)vm+1=[(κ0h2)-1δx2(1+(1/12)δx2)+(κ0h2)-1δy2(1+(1/12)δy2)+(κ0h2)-1δz2(1+(1/12)δz2)]×vm+1+Fm+1,1≤m≤M-1,
where
(20)κ0=τ-αΓ(2-α)+p2+q2+r24-s,Fm+1=1κ0[τ-αvmΓ(2-α)-τ-α∑l=1mbl(vm+1-l-vm-l)Γ(2-α)-((m+1)τ)-αv0Γ(1-α)+λ-1(2ϕm-ϕm-1+fm+1)].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 casem=0, the scheme simply reads
(21)v1=[(κ0h2)-1δx2(1+(1/12)δx2)+(κ0h2)-1δy2(1+(1/12)δy2)+(κ0h2)-1δz2(1+(1/12)δz2)]×v1+F1
with
(22)F1=1κ0[ατ-αΓ(2-α)v0+λ-1(ϕ1+f1)].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
(23)vm=[(κ0h2)-1δx2(1+(1/12)δx2)+(κ0h2)-1δy2(1+(1/12)δy2)+(κ0h2)-1δz2(1+(1/12)δz2)]×vm+Fm,1≤m≤M.Multiplying(1+(1/12)δx2)(1+(1/12)δy2)(1+(1/12)δz2) on both sides of (23), we have
(24)(1+112δx2)(1+112δy2)(1+112δz2)vm=(κ0h2)-1δx2(1+112δy2)(1+112δz2)vm+(κ0h2)-1δy2(1+112δx2)(1+112δz2)vm+(κ0h2)-1δz2(1+112δx2)(1+112δy2)vm+(1+112δx2)(1+112δy2)(1+112δz2)Fm,1≤m≤M.
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:
(25)Lx=1+112δx2-1κ0h2δx2,Ly=1+112δy2-1κ0h2δy2,Lz=1+112δz2-1κ0h2δz2,LF=(1+112δx2)(1+112δy2)(1+112δz2),LH=((1+112δx2)δy2δz2(κ0h2)2+(1+112δy2)δx2δz2(κ0h2)2+(1+112δz2)δx2δy2(κ0h2)2)-δx2δy2δz2(κ0h2)3.Algorithm 1.
(1) Form≥1, add splitting term
(26)I:LH(vm-vm-1)=O(τ1+2α)
to the left-hand side of (24), and rewrite it into an equivalent form
(27)SI:LxLyLzvm=LFFm+LHvm-1.
(2) Introducing two intermediate variables v* and v**, scheme SI can be solved in three steps as(28a)LxΔv*=LFFm+LHvm-1,(28b)LyΔv**=v*,(28c)Lzv=v**.
Note that the intermediate values ofv** and v* on the boundary are easily obtained from (28c) and (28b).Algorithm 2.
(1) For the first time step,v1 can be obtained from Algorithm 1.
(2) Form≥2, add splitting term
(29)II:LH(vm-2vm-1+vm-2)=O(τ2+2α)
to the left-hand side of (24), and rewrite it into an equivalent form
(30)SII:LxLyLzvm=LFFm+LH(2vm-1-vm-2).
(3) Use the same way as the second step of Algorithm1.Remark 3.
It is clear that the splitting error ofII 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:
(31)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 O(τmin(1+α,2-α)+h4) and O(τ(2-α)+h4), respectively. The resulting ADI schemes in each ADI solution step corresponding to a strictly diagonally dominant matrix equation which can be solved using the one-dimensional tridiagonal Thomas algorithm with a considerable saving in computing time.Remark 4.
The previous analysis shows that splitting termII 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.


## 3. 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 Algorithm2 will be compared with Algorithm 1. The L∞-norm error (denoted as L∞-error) and the L2-norm error (denoted as L2-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ΔH and Δh on ΩH and Ωh, respectively. The L2-norm errors of these two grids are denoted as ∥·∥0H and ∥·∥0h. If we set the convergent order of accuracy to be Rate, then we have the following form:
(32)(ΔHΔh)Rate=∥·∥0H∥·∥0h.
So, the convergent order of accuracy Rate can be computed as
(33)Rate=log∥·∥0H/∥·∥0hlog(ΔH/Δh).
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(34)∂αu∂tα=Δu+u2+f(x,y,z,t),u(x,y,z,t)=etsin(πx)sin(πy)sin(πz),u0(x,y,z)=sin(πx)sin(πy)sin(πz).
Since it is difficult to accurately solve ∂αu/∂tα based on u except for α=0.5, in this problem, we only discuss the case with α=0.5 and the corresponding
(35)f(x,y,z,t)=(t-αΓ(1-α)+πeterf(t)+3π2(t2+1))×sin(πx)sin(πy)sin(πz)-e2tsin2(πx)sin2(πy)sin2(πz).
We first investigate the time accuracy withT=1 for both algorithms. To this end, N=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.2764e-3 when N=40,M=50. Whereas for Algorithm 1 with the same mesh size, the maximum absolute error is around 9.9234e-2.  Similar comparison can be made to reach similar conclusions.
Meanwhile, we can see that the convergence rate for Algorithm1 is around 1.5, which agrees with the theoretical estimates min(1+α, 2-α). Howevre, for Algorithm 2, the L2-rate is approximately 2.6, which is superconvergent. For better comparison of two algorithms, we plot the L∞ errors (as shown in Figure 1) at T=1 and z=0.5 with N=40 and M=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 chooseT=0.01 and M=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.
Figure2 shows the contour plot of the exact solution in the x-y plane and numerical solutions using Algorithm 2 for different N at z=0.5 and T=0.01 with α=0.5 and M=100. We can see that when N=20, the numerical solution is a good approximation of the exact solution.Table 1
T
=
1
,
N
=
40, and α=0.5 with different M for Problem 5.
M
10
20
30
40
50
Algorithm1
L
∞-error
8.5329
e
-
1
3.9114
e
-
1
2.1763
e
-
1
1.4039
e
-
1
e
9.9234
e
-
2
L
∞-rate
—
1.1253
1.4459
1.5239
1.5547
L
2-error
2.7920
e
-
1
1.2721
e
-
1
7.0539
e
-
2
4.5420
e
-
2
3.2070
e
-
2
L
2-rate
—
1.1341
1.4544
1.5302
1.5596
Algorithm2
L
∞-error
1.6081
e
-
1
2.4544
e
-
2
8.4991
e
-
3
4.0461
e
-
3
2.2764
e
-
3
L
∞-rate
—
2.7119
2.6156
2.5800
2.5776
L
2-error
5.2043
e
-
2
7.8663
e
-
3
2.7077
e
-
3
1.2840
e
-
3
7.2062
e
-
3
L
2-rate
—
2.7259
2.6303
2.5936
2.5884Table 2
T
=
0.01,  M=100, and α=0.5 with different N using Algorithm 2 for Problem 5.
N
5
10
20
30
40
L
∞-error
4.8628
e
-
4
3.5004
e
-
5
2.1860
e
-
6
4.3559
e
-
7
1.4123
e
-
7
L
∞-rate
—
3.7962
4.0012
3.9784
3.9152
L
2-error
1.5138
e
-
4
1.0654
e
-
5
7.1348
e
-
7
1.4566
e
-
7
4.7836
e
-
8
L
2-rate
—
3.8287
3.9004
3.9186
3.8706The plot of theL∞-error using both algorithms at T=1 and z=0.5 with N=40,M=40, and α=0.5 for Problem 5.
(a)
TheL∞-error using Algorithm 1
(b)
TheL∞-error using Algorithm 2The contour plot of the exact solution and numerical solutions for differentN at T=0.01 and  z=0.5 with α=0.5 andM=100 for Problem 5. (a) Exact solution; (b) N=5; (c) N=10; (d) N=20.
(a)
(b)
(c)
(d)Problem 6.
In order to illustrate the effectiveness of our method, the following time-fractional convection-diffusion problem is considered:(36)∂αu∂tα=Δu-2ux-4uy-6uz+u+e-3x-6y-9z(u2-e2x+4y+6z)2+ex+2y+3zf(x,y,z,t),(37)u(x,y,z,t)=(t2+1)x2(1-x)3sin(πy)sin(πz)ex+2y+3z,(38)u0(x,y,z)=x2(1-x)3sin(πy)sin(πz)ex+2y+3z.
The right-hand side function f is specified to satisfy the given equation and the exact solution.
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.1067e-4 when N=20,M=50. Whereas for Algorithm 1 with the same mesh size, the maximum absolute error is around 6.0692e-2.
However, we note that the convergence orders for Algorithm1 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 τ.
Figure3 shows the L∞ errors at T=2 and z=0.5 with N=50 and M=40 for different α. Again, as expected, it indicates that Algorithm 2 has substantially higher accuracy than Algorithm 1 regardless of α.
Table6 presents space results with T=0.01 and M=100 for different α and N, from which we can find that Algorithm 2 is stable and convergent and it can obtain fourth-order accuracy when mesh becomes finer, regardless of the value of α.
Figure4 shows the contour plot of the exact solution in the x-y plane and numerical solutions using Algorithm 2 for different N at z=0.5 and T=0.01 with α=0.5 and M=100. We can see that when N=20, the numerical solution is a good approximation of the exact solution.Table 3
T
=
2,   N=20, and α=0.2 with different M for Problem 6.
M
10
20
30
40
50
Algorithm1
L
∞-error
5.5947
e
-
1
2.9608
e
-
1
1.9977
e
-
1
1.5025
e
-
1
e
1.2017
e
-
1
L
∞-rate
—
0.9181
0.9703
0.9903
1.0012
L
2-error
1.2173
e
-
1
6.4253
e
-
2
4.3325
e
-
2
3.2575
e
-
2
2.6049
e
-
2
L
2-rate
—
0.9219
0.9720
0.9913
1.0019
Algorithm2
L
∞-error
6.1711
e
-
2
1.5422
e
-
2
6.7626
e
-
3
3.7638
e
-
3
2.3877
e
-
3
L
∞-rate
—
2.0006
2.0331
2.0369
2.0396
L
2-error
1.3403
e
-
2
3.3369
e
-
3
1.4631
e
-
3
8.1426
e
-
4
5.1651
e
-
4
L
2-rate
—
2.0060
2.0335
2.0370
2.0398Table 4
T
=
2
,
N
=
20, and α=0.5 with different M for Problem 6.
M
10
20
30
40
50
Algorithm1
L
∞-error
5.1863
e
-
1
2.5840
e
-
1
1.6670
e
-
1
1.2088
e
-
1
e
9.3705
e
-
2
L
∞-rate
—
1.0051
1.0810
1.1173
1.1411
L
2-error
1.1282
e
-
1
5.6080
e
-
2
3.6162
e
-
2
2.6216
e
-
2
2.0321
e
-
2
L
2-rate
—
1.0084
1.0822
1.1180
1.1414
Algorithm2
L
∞-error
5.6623
e
-
2
1.3141
e
-
2
5.4618
e
-
3
2.9080
e
-
3
1.7745
e
-
3
L
∞-rate
—
2.1073
2.1653
2.1910
2.2137
L
2-error
1.2280
e
-
2
2.8416
e
-
3
1.1808
e
-
3
6.2852
e
-
4
3.8345
e
-
4
L
2-rate
—
2.1116
2.1659
2.1918
2.2146Table 5
T
=
2,  N=20, and α=0.8 with different M for Problem 6.
M
10
20
30
40
50
Algorithm1
L
∞-error
4.7068
e
-
1
2.1003
e
-
1
1.2397
e
-
1
8.3453
e
-
2
e
6.0692
e
-
2
L
∞-rate
—
1.1641
1.3002
1.3757
1.4272
L
2-error
1.0234
e
-
1
4.5581
e
-
2
2.6896
e
-
2
1.8105
e
-
2
1.3168
e
-
2
L
2-rate
—
1.1669
1.3010
1.3758
1.4271
Algorithm2
L
∞-error
4.9685
e
-
2
9.7522
e
-
3
3.4660
e
-
3
1.5740
e
-
3
8.1067
e
-
4
L
∞-rate
—
2.3490
2.5514
2.7440
2.9734
L
2-error
1.0752
e
-
2
2.1064
e
-
3
7.5044
e
-
4
3.4373
e
-
4
1.8083
e
-
4
L
2-rate
—
2.3517
2.5455
2.7141
2.8783Table 6
T
=
0.01 and M=100 with different α and N using Algorithm 2 for Problem 6.
N
5
10
20
30
40
α
=
0.2
L
∞-error
2.7154
e
-
4
1.8711
e
-
5
1.1750
e
-
6
2.2446
e
-
7
6.4536
e
-
8
L
∞-rate
—
3.8592
3.9932
4.0825
4.3329
L
2-error
6.3987
e
-
5
4.5423
e
-
6
3.0190
e
-
7
5.9489
e
-
8
1.7738
e
-
8
L
2-rate
—
3.8163
3.9113
4.0060
4.2064
α
=
0.5
L
∞-error
2.4981
e
-
4
1.7177
e
-
5
1.0951
e
-
6
2.2343
e
-
7
7.2771
e
-
8
L
∞-rate
—
3.8623
3.9713
3.9202
3.8994
L
2-error
5.8471
e
-
5
4.1544
e
-
6
2.7991
e
-
7
5.8491
e
-
8
1.9056
e
-
8
L
2-rate
—
3.8150
3.8916
3.8612
3.8983
α
=
0.8
L
∞-error
2.1863
e
-
4
1.5586
e
-
5
1.0322
e
-
6
2.1189
e
-
7
6.8378
e
-
8
L
∞-rate
—
3.8102
3.9165
3.9051
3.9315
L
2-error
5.2327
e
-
5
3.7485
e
-
6
2.4905
e
-
7
5.0023
e
-
8
1.6021
e
-
8
L
2-rate
—
3.8023
3.9118
3.9588
3.9578The plot of theL∞-error using both algorithms at T=2 and  z=0.5 with N=40,M=40, and different α for Problem 6.
(a)
Algorithm1, α=0.2
(b)
Algorithm1, α=0.5
(c)
Algorithm1, α=0.8
(d)
Algorithm2, α=0.2
(e)
Algorithm2, α=0.5
(f)
Algorithm2, α=0.8The contour plot of the exact solution and numerical solutions for differentN at T=0.01 andz=0.5 with α=0.5 and  M=100 for Problem 6. (a) Exact solution; (b) N=5; (c) N=10; (d) N=20.
(a)
(b)
(c)
(d)


## 4. Conclusions

In this paper, we have proposed two high-order compact ADI algorithms with orderO(τmin(1+α,2-α)+h4) and O(τ2-α+h4), respectively, for solving the 3D nonlinear time-fractional convection-diffusion equation. Since it only involves three-point 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.



---

*Source: 246025-2013-05-15.xml*