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

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 orders 𝑂(𝜏^{min(1+𝛼,2−𝛼)}+ℎ⁴) and 𝑂(𝜏^2−𝛼+ℎ⁴) 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.

1. Introduction

In this paper, we consider numerical approximations of the following 3D nonlinear time-fractional convection-diffusion equation: 𝜕^𝛼𝑢𝜕𝑡^𝛼_Ω×(0,𝑇]=Δ𝑢−𝑝𝑢_𝑥−𝑞𝑢_𝑦−𝑟𝑢_𝑧+𝑠𝑢+𝜙(𝑢)+𝑓_Ω×(0,𝑇], (1) with the initial and boundary conditions 𝑢𝑥,𝑦,𝑧,0_Ω=𝑢₀𝑥,𝑦,𝑧_Ω,𝑢𝑥,𝑦,𝑧,𝑡_{𝜕Ω×(0,𝑇]}=0, (2) 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 𝜕^𝛼𝑢𝜕𝑡^𝛼=𝑡^−𝛼Γ()𝑢1−𝛼⁰+1Γ(2−𝛼)^𝑡₀𝜕𝑢𝑥,𝑦,𝑧,𝑠𝑑𝜕𝑠(𝑡−^𝑠𝑠)^𝛼. (3) It should be noted that Grünward-Letnikov fractional derivative will become Caputo fractional derivative when 𝑢⁰=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 𝑂(𝜏+ℎ⁴) in [11], and a higher order 𝑂(𝜏^2−𝛼+ℎ⁴) 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 𝑂(𝜏^2−𝛼+ℎ⁴) 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 Section 2, 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]³, 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: 𝜙𝑢^𝑚+1=2𝜙𝑢^𝑚−𝜙𝑢^𝑚−1+𝑂𝜏²,1≤𝑚≤𝑀−1. (4)

Next we discuss the linear equation 𝜕^𝛼𝑢^𝑚+1𝜕𝑡^𝛼=Δ𝑢^𝑚+1−𝑝𝑢_𝑥^𝑚+1−𝑞𝑢_𝑦^𝑚+1−𝑟𝑢_𝑧^𝑚+1+𝑠𝑢^𝑚+1+2𝜙𝑢^𝑚−𝜙𝑢^𝑚−1+𝑓^𝑚+1 (5) 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 𝑢=𝜆₁(𝑥)𝜆₂𝑦𝜆₃(𝑧)𝑣, (6) where 𝑣=𝑣(𝑥,𝑦,𝑧,𝑡) and 𝜆₁(𝑥),𝜆₂(𝑦), and 𝜆₃(𝑧) are functions to be determined. Substituting (6) into the first formula of (1), we have 𝜆₁𝜆₂𝜆₃𝜕^𝛼𝑣𝑡^𝛼=𝜆−Δ𝑣_1𝑥𝑥𝜆₂𝜆₃−𝑝𝜆_1𝑥𝜆₂𝜆₃+𝜆₁𝜆_2𝑦𝑦𝜆₃−𝑞𝜆₁𝜆_2𝑦𝜆₃+𝜆₁𝜆₂𝜆_3𝑧𝑧−𝑟𝜆₁𝜆₂𝜆_3𝑧𝑣+2𝜆_1𝑥−𝑝𝜆₁𝜆₂𝜆₃𝑣_𝑥+2𝜆_2𝑦−𝑞𝜆₂𝜆₁𝜆₃𝑣_𝑦+2𝜆_3𝑧−𝑟𝜆₃𝜆₁𝜆₂𝑣_𝑧. (7) Setting the coefficients of 𝑣_𝑥,𝑣_𝑦, and 𝑣_𝑧 to zero leads to 2𝜆_1𝑥−𝑝𝜆₁=0,2𝜆_2𝑦−𝑞𝜆₂=0,2𝜆_3𝑧−𝑟𝜆₃=0. (8) Then, we have 𝜆₁(𝑥)=𝑒^(𝑝/2)𝑥,𝜆₂𝑦=𝑒^(𝑞/2)𝑦,𝜆₃(𝑧)=𝑒^(𝑟/2)𝑧. (9)

Set 𝜆𝑥,𝑦,𝑧=𝜆₁(𝑥)𝜆₂𝑦𝜆₃(𝑧)=𝑒^{(𝑝/2)𝑥+(𝑞/2)𝑦+(𝑟/2)𝑧}. (10)

Combining (1) and (5)–(10), we obtain the following time-fractional diffusion equation satisfied by 𝑣: 𝜕^𝛼𝑣^𝑚+1𝜕𝑡^𝛼_Ω=Δ𝑣^𝑚+1+𝑝𝑠−²+𝑞²+𝑟²4𝑣^𝑚+1+𝜆⁻¹2𝜙𝜆𝑣^𝑚−𝜙𝜆𝑣^𝑚−1+𝑓^𝑚+1_Ω,𝑣𝑚≥1,𝑥,𝑦,𝑧,0_Ω=𝜆⁻¹𝑢₀𝑥,𝑦,𝑧_Ω,𝑣𝑥,𝑦,𝑧,𝑡^𝑚_𝜕Ω=0,𝑚≥0. (11)

In the following, we will introduce a high-order compact difference method for (11).

The time-fractional derivative 𝜕^𝛼𝑣/𝜕𝑡^𝛼 at 𝑡_𝑚+1 is estimated by 𝜕^𝛼𝑣^𝑚+1𝜕𝑡^𝛼=𝑡^−𝛼_𝑚+1Γ()𝑣1−𝛼⁰+1Γ()×1−𝛼_{𝑚𝑙=0𝑡𝑙+1𝑡𝑙}𝜕𝑣𝑥,𝑦,𝑧,𝜉𝑑𝜕𝜉𝑡^𝜉_𝑚+1−𝜉^𝛼=𝑡^−𝛼_𝑚+1Γ()𝑣1−𝛼⁰+1Γ()×1−𝛼_𝑚𝑙=0𝑣^𝑙+1−𝑣^𝑙𝜏^𝑡𝑙+1_𝑡𝑙𝑑_𝜉𝑡_𝑚+1−𝜉^𝛼+𝐸_𝜏^𝑚+1=𝑡^−𝛼_𝑚+1Γ()𝑣1−𝛼⁰+1Γ(2−𝛼)^𝑚_𝑙=0𝑏_𝑙𝑣^{𝑚+1−𝑙}−𝑣^𝑚−𝑙𝜏^𝛼+𝐸_𝜏^𝑚+1, (12) where 𝑏_𝑙=(𝑙+1)^1−𝛼−𝑙^1−𝛼 and the truncation error 𝐸_𝜏^𝑚+1 has the following form: 𝐸_𝜏^𝑚+1=1Γ()×1−𝛼_{𝑚𝑙=0𝑡𝑙+1𝑡𝑙}𝜕𝑣𝑥,𝑦,𝑧,𝜉−𝑣𝜕𝜉^𝑙+1−𝑣^𝑙𝜏×𝑑_𝜉𝑡_𝑚+1−𝜉^𝛼. (13) If we expend 𝑣(𝑥,𝑦,𝑧,𝜉) in Taylor series about 𝜉 at 𝑡_𝑙+1/2 in (13) and use the analysis in [14], we will get 𝐸_𝜏^𝑚+1≤𝐶_𝑣1Γ(1−𝛼)^𝑚_𝑙=0^𝑡𝑙+1_𝑡𝑙𝑡_𝑙+1+𝑡_𝑙−2𝜉𝑡_𝑚+1−𝜉^𝛼𝑑_𝜉+𝑂𝜏²≤𝐶_𝑣𝜏^2−𝛼, (14) where 𝐶_𝑣 is a positive constant only depending on 𝑣.

As for the spatial derivatives, we can use the fourth-order Padé scheme. That is, the spatial derivatives can be approximated by 𝛿Δ𝑣=^𝑥2ℎ²1+(1/12)𝛿_𝑥2+𝛿_𝑦2ℎ²1+(1/12)𝛿_𝑦2+𝛿_𝑧2ℎ²1+(1/12)𝛿_𝑧2×𝑣+𝑂ℎ⁴, (15) where 𝛿²_𝑥, 𝛿²_𝑦, and 𝛿²_𝑧 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: ((𝑚+1)𝜏)^−𝛼Γ()𝑣1−𝛼^{0𝑖𝑗𝑘}+1Γ(2−𝛼)^𝑚_𝑙=0𝑏_𝑙𝑣^{𝑚+1−𝑙}_𝑖𝑗𝑘−𝑣^𝑚−𝑙_𝑖𝑗𝑘𝜏^𝛼=𝛿_𝑥2ℎ²1+(1/12)𝛿_𝑥2+𝛿_𝑦2ℎ²1+(1/12)𝛿_𝑦2+𝛿_𝑧2ℎ²1+(1/12)𝛿_𝑧2𝑣^𝑚+1_𝑖𝑗𝑘+𝑝𝑠−²+𝑞²+𝑟²4𝑣^𝑚+1_𝑖𝑗𝑘+𝜆⁻¹2𝜙^𝑚_𝑖𝑗𝑘−𝜙^𝑚−1_𝑖𝑗𝑘+𝑓^𝑚+1_𝑖𝑗𝑘,1≤𝑖,𝑗,𝑘≤𝑁−1,1≤𝑚≤𝑀−1, (16) with the initial discretization 𝑣⁰_𝑖𝑗𝑘=𝜆𝑥_𝑖,𝑦_𝑗,𝑧_𝑘𝑢₀𝑥_𝑖,𝑦_𝑗,𝑧_𝑘,0≤𝑖,𝑗,𝑘≤𝑁, (17) and the boundary discretization 𝑣^𝑚_0𝑗𝑘=𝑣^𝑚_𝑁𝑗𝑘=𝑣^𝑚_𝑖0𝑘=𝑣^𝑚_𝑖𝑁𝑘=𝑣^𝑚_𝑖𝑗0=𝑣^𝑚_𝑖𝑗𝑁=0,0≤𝑖,𝑗,𝑘≤𝑁,0≤𝑚≤𝑀, (18) where 𝑣^𝑚_𝑖𝑗𝑘 denotes 𝑣(𝑖ℎ,𝑗ℎ,𝑘ℎ,𝑚𝜏). In the following, 𝑣^𝑚 will be written in short for 𝑣^𝑚_𝑖𝑗𝑘 if there is no confusion about the notation.

By reformulating (16), we obtain an equivalent form 𝑣^𝑚+1=𝜅₀ℎ²⁻¹𝛿_𝑥21+(1/12)𝛿_𝑥2+𝜅₀ℎ²⁻¹𝛿_𝑦21+(1/12)𝛿_𝑦2+𝜅₀ℎ²⁻¹𝛿_𝑧21+(1/12)𝛿_𝑧2×𝑣^𝑚+1+𝐹^𝑚+1,1≤𝑚≤𝑀−1, (19) where 𝜅₀=𝜏^−𝛼Γ()+𝑝2−𝛼²+𝑞²+𝑟²4𝐹−𝑠,_𝑚+1=1𝜅₀𝜏^−𝛼𝑣^𝑚Γ()−𝜏2−𝛼^−𝛼∑^𝑚_𝑙=1𝑏_𝑙𝑣^{𝑚+1−𝑙}−𝑣^𝑚−𝑙Γ()−)𝜏)2−𝛼((𝑚+1^−𝛼𝑣⁰Γ(1−𝛼)+𝜆⁻¹2𝜙^𝑚−𝜙^𝑚−1+𝑓^𝑚+1. (20)

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 𝑣¹=𝜅₀ℎ²⁻¹𝛿_𝑥21+(1/12)𝛿_𝑥2+𝜅₀ℎ²⁻¹𝛿_𝑦21+(1/12)𝛿_𝑦2+𝜅₀ℎ²⁻¹𝛿_𝑧21+(1/12)𝛿_𝑧2×𝑣¹+𝐹¹ (21) with 𝐹¹=1𝜅₀𝛼𝜏^−𝛼Γ()𝑣2−𝛼⁰+𝜆⁻¹𝜙¹+𝑓¹. (22)

It is clear that (21) must be solved using some relaxation procedure because 𝜙¹ is nonlinear.

For simplicity of the presentation, we write (19) and (21) in a unified form 𝑣^𝑚=𝜅₀ℎ²⁻¹𝛿_𝑥21+(1/12)𝛿_𝑥2+𝜅₀ℎ²⁻¹𝛿_𝑦21+(1/12)𝛿_𝑦2+𝜅₀ℎ²⁻¹𝛿_𝑧21+(1/12)𝛿_𝑧2×𝑣^𝑚+𝐹^𝑚,1≤𝑚≤𝑀. (23)

Multiplying (1+(1/12)𝛿_𝑥2)(1+(1/12)𝛿_𝑦2)(1+(1/12)𝛿_𝑧2) on both sides of (23), we have 11+𝛿12^𝑥211+𝛿12^𝑦211+𝛿12^𝑧2𝑣^𝑚=𝜅₀ℎ²⁻¹𝛿_𝑥211+𝛿12^𝑦211+𝛿12^𝑧2𝑣^𝑚+𝜅₀ℎ²⁻¹𝛿_𝑦211+𝛿12^𝑥211+𝛿12^𝑧2𝑣^𝑚+𝜅₀ℎ²⁻¹𝛿_𝑧211+𝛿12^𝑥211+𝛿12^𝑦2𝑣^𝑚+11+𝛿12^𝑥211+𝛿12^𝑦211+𝛿12^𝑧2𝐹^𝑚,1≤𝑚≤𝑀. (24) 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: 𝐿_𝑥1=1+𝛿12^2𝑥−1𝜅₀ℎ²𝛿²_𝑥,𝐿_𝑦1=1+𝛿12^2𝑦−1𝜅₀ℎ²𝛿²_𝑦,𝐿_𝑧1=1+𝛿12^2𝑧−1𝜅₀ℎ²𝛿²_𝑧,𝐿_𝐹=11+𝛿12^2𝑥11+𝛿12^2𝑦11+𝛿12^2𝑧,𝐿_𝐻=11+𝛿12^2𝑥𝛿²_𝑦𝛿²_𝑧𝜅₀ℎ²²+11+𝛿12^2𝑦𝛿²_𝑥𝛿²_𝑧𝜅₀ℎ²²+11+𝛿12^2𝑧𝛿²_𝑥𝛿²_𝑦𝜅₀ℎ²²−𝛿²_𝑥𝛿²_𝑦𝛿²_𝑧𝜅₀ℎ²³. (25)

Algorithm 1. (1) For 𝑚≥1, add splitting term 𝐈∶𝐿_𝐻𝑣^𝑚−𝑣^𝑚−1=𝑂𝜏^1+2𝛼 (26) to the left-hand side of (24), and rewrite it into an equivalent form 𝐒𝐈∶𝐿_𝑥𝐿_𝑦𝐿_𝑧𝑣^𝑚=𝐿_𝐹𝐹^𝑚+𝐿_𝐻𝑣^𝑚−1. (27) (2) Introducing two intermediate variables 𝑣^∗ and 𝑣^∗∗, scheme 𝐒𝐈 can be solved in three steps as𝐿_𝑥Δ𝑣^∗=𝐿_𝐹𝐹^𝑚+𝐿_𝐻𝑣^𝑚−1,𝐿_𝑦Δ𝑣^∗∗=𝑣^∗,𝐿_𝑧𝑣=𝑣^∗∗. (28a)(28b)(28c)
Note that the intermediate values of 𝑣^∗∗ and 𝑣^∗ on the boundary are easily obtained from (28c) and (28b).

Algorithm 2. (1) For the first time step, 𝑣¹ can be obtained from Algorithm 1.
(2) For 𝑚≥2, add splitting term 𝐈𝐈∶𝐿_𝐻𝑣^𝑚−2𝑣^𝑚−1+𝑣^𝑚−2=𝑂𝜏^2+2𝛼 (29) to the left-hand side of (24), and rewrite it into an equivalent form 𝐒𝐈𝐈∶𝐿_𝑥𝐿_𝑦𝐿_𝑧𝑣^𝑚=𝐿_𝐹𝐹^𝑚+𝐿_𝐻2𝑣^𝑚−1−𝑣^𝑚−2. (30)
(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: ()=(],().min1+𝛼,2−𝛼1+𝛼,𝛼∈0,0.52−𝛼,𝛼∈0.5,1 (31) Based on the previous analysis, we propose two linear ADI schemes, and the truncation errors are 𝑂(𝜏^{min(1+𝛼,2−𝛼)}+ℎ⁴) and 𝑂(𝜏^(2−𝛼)+ℎ⁴), 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 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.

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 Algorithm 2 will be compared with Algorithm 1. The 𝐿^∞-norm error (denoted as 𝐿^∞-error) and the 𝐿²-norm error (denoted as 𝐿²-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 𝐿²-norm errors of these two grids are denoted as ‖⋅‖^𝐻₀ and ‖⋅‖^ℎ₀. If we set the convergent order of accuracy to be Rate, then we have the following form: Δ𝐻Δℎ^Rate=‖⋅‖^𝐻₀‖⋅‖^ℎ₀. (32) So, the convergent order of accuracy Rate can be computed as ‖⋅‖Rate=log^𝐻₀/‖⋅‖^ℎ₀().logΔ𝐻/Δℎ (33) 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 𝜕^𝛼𝑢𝜕𝑡^𝛼=Δ𝑢+𝑢²,𝑢+𝑓𝑥,𝑦,𝑧,𝑡𝑥,𝑦,𝑧,𝑡=𝑒^𝑡()(),𝑢sin𝜋𝑥sin𝜋𝑦sin𝜋𝑧₀𝑥,𝑦,𝑧=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 𝑓=𝑡𝑥,𝑦,𝑧,𝑡^−𝛼Γ()+√1−𝛼𝜋𝑒^𝑡erf√𝑡+3𝜋²𝑡²()()+1×sin𝜋𝑥sin𝜋𝑦sin𝜋𝑧−𝑒^2𝑡sin²(𝜋𝑥)sin²𝜋𝑦sin²(𝜋𝑧). (35)
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 𝐿²-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.
Figure 2 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.

Table 1 

𝑇=1,𝑁=40, and 𝛼=0.5 with different 𝑀 for Problem 5.

Table 2 

𝑇=0.01,  𝑀=100, and 𝛼=0.5 with different 𝑁 using Algorithm 2 for Problem 5.

(a) The 𝐿∞-error using Algorithm 1

(b) The 𝐿∞-error using Algorithm 2

Figure 1 

The plot of the 𝐿∞-error using both algorithms at 𝑇=1 and 𝑧=0.5 with 𝑁=40,𝑀=40, and 𝛼=0.5 for Problem 5.

(a)

(b)

(c)

(d)

Figure 2 

The contour plot of the exact solution and numerical solutions for different 𝑁 at 𝑇=0.01 and  𝑧=0.5 with 𝛼=0.5 and𝑀=100 for Problem 5. (a) Exact solution; (b) 𝑁=5; (c) 𝑁=10; (d) 𝑁=20.

Problem 6. In order to illustrate the effectiveness of our method, the following time-fractional convection-diffusion problem is considered: 𝜕^𝛼𝑢𝜕𝑡^𝛼=Δ𝑢−2𝑢_𝑥−4𝑢_𝑦−6𝑢_𝑧+𝑢+𝑒^{−3𝑥−6𝑦−9𝑧}𝑢²−𝑒^{2𝑥+4𝑦+6𝑧2}+𝑒^{𝑥+2𝑦+3𝑧}𝑓,𝑢=𝑡𝑥,𝑦,𝑧,𝑡𝑥,𝑦,𝑧,𝑡²+1𝑥²(1−𝑥)³sin𝜋𝑦sin(𝜋𝑧)𝑒^{𝑥+2𝑦+3𝑧}𝑢₀,𝑥,𝑦,𝑧=𝑥²(1−𝑥)³sin𝜋𝑦sin(𝜋𝑧)𝑒^{𝑥+2𝑦+3𝑧}. (36)(37)(38) The right-hand side function 𝑓 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.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 fourth-order 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.

Table 3 

𝑇=2,   𝑁=20, and 𝛼=0.2 with different 𝑀 for Problem 6.

Table 4 

𝑇=2,𝑁=20, and 𝛼=0.5 with different 𝑀 for Problem 6.

Table 5 

𝑇=2,  𝑁=20, and 𝛼=0.8 with different 𝑀 for Problem 6.

Table 6 

𝑇=0.01 and 𝑀=100 with different 𝛼 and 𝑁 using Algorithm 2 for Problem 6.

(a) Algorithm 1, 𝛼=0.2

(b) Algorithm 1, 𝛼=0.5

(c) Algorithm 1, 𝛼=0.8

(d) Algorithm 2, 𝛼=0.2

(e) Algorithm 2, 𝛼=0.5

(f) Algorithm 2, 𝛼=0.8

Figure 3 

The plot of the 𝐿∞-error using both algorithms at 𝑇=2 and  𝑧=0.5 with 𝑁=40,𝑀=40, and different 𝛼 for Problem 6.

(a)

(b)

(c)

(d)

Figure 4 

The contour plot of the exact solution and numerical solutions for different 𝑁 at 𝑇=0.01 and𝑧=0.5 with 𝛼=0.5 and  𝑀=100 for Problem 6. (a) Exact solution; (b) 𝑁=5; (c) 𝑁=10; (d) 𝑁=20.

4. Conclusions

In this paper, we have proposed two high-order compact ADI algorithms with order 𝑂(𝜏^{min(1+𝛼,2−𝛼)}+ℎ⁴) and 𝑂(𝜏^2−𝛼+ℎ⁴), 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.

Acknowledgments

The authors would like to thank the editor and referees for their criticism, valuable comments, and suggestions which helped to improve the results of this paper. This work is in part supported by the NSF of China (nos. 11271313 and 61163027), the China Postdoctoral Science Foundation (nos. 201104702 and 2012M512056), the Key Project of Chinese Ministry of Education (no. 212197), the NSF of Xinjiang Province (no. 2013211B01), and the Excellent Doctor Innovation Program of Xinjiang University (no. XJUBSCX-2012003).