Maximum Norm Error Estimates of ADI Methods for a Two-Dimensional Fractional Subdiffusion Equation

Abstract

This paper is concerned with two alternating direction implicit (ADI) finite difference methods for solving a two-dimensional fractional subdiffusion equation. An explicit error estimate for each of the two methods is provided in the discrete maximum norm. It is shown that the methods have the same order as their truncation errors with respect to the discrete maximum norm. Numerical results are given to confirm the theoretical analysis results.

1. Introduction

Fractional differential equations and fractional calculus arise in various application problems in science and engineering [1–16]. Various numerical methods have been developed for the computation of fractional differential equations [17–34]. Fractional subdiffusion equations describe a special type of anomalous diffusion [35], and it is a more difficult task to solve this kind of equation numerically.

Numerical works for fractional subdiffusion equations are mostly focused on one-dimensional problems due to the memory effect in fractional derivatives; see, for example, [19, 20, 20–26, 31–33, 36–42]. A two-dimensional anomalous subdiffusion equation was numerically treated in [43, 44], where explicit and implicit finite difference schemes were proposed. Chen et al. [28] extended their work in [43] to a variable-order subdiffusion equation. Liu et al. [45] developed an implicit meshless approach based on the radial basis function for the numerical simulation of a two-dimensional subdiffusion problem. Chen and Liu [18] considered an implicit difference scheme for a three-dimensional fractional advection-diffusion equation, and a Richardson extrapolation was applied to improve the accuracy.

The complexity of the fractional differential equations comes from the involving fractional derivatives that are nonlocal and have the character of history dependence and universal mutuality. This means that the computations would be costly if the implicit schemes were applied, especially for solving multidimensional problems [43, 44]. Some researchers have explored some techniques for reducing this cost. These techniques include the adaptive technique [46] and the matrix transfer technique [47, 48]. It is well known that alternating direction implicit (ADI) methods are unconditionally stable as the traditional implicit methods. On the other hand, they reduce a multidimensional problem to a series of independent one-dimensional problems, and thus the computational complexities and the computational cost can be greatly reduced. Therefore, ADI methods for fractional differential equations have the potential to significantly reduce the computational cost, while maintaining the stability of the numerical methods. The works in [29, 49–53] treated ADI finite difference methods for space fractional diffusion equations, and the work in [54] discussed ADI finite difference methods for fractional diffusion wave equations. Recently, Cui [55] derived an ADI compact finite difference scheme for a two-dimensional fractional subdiffusion equation, where the Grünwald formula is used to approximate the temporal Riemann-Liouville fractional derivative, and the spatial derivatives are approximated by a compact finite difference scheme. Another way to treat fractional subdiffusion problem is to transform the original subdiffusion equation into an equivalent equation by replacing the temporal Riemann-Liouville fractional derivative by the temporal Caputo fractional derivative; see, for example, [40, 42, 56]. An advantage of this approach is that the 𝐿₁ approximation (see [1, 24, 25, 40, 57, 58]) can be used to deal with the temporal fractional derivative, and so the resulting scheme has the better temporal accuracy than the first order without the Crank-Nicolson technique (see [25, 40–42]). Based on the previous approach, Zhang and Sun [59] constructed two ADI finite difference schemes, called 𝐿₁-ADI and BD-ADI schemes, for the following two-dimensional problem of subdiffusion equation with the temporal Caputo fractional derivative: ^𝐶₀𝒟^𝛼_𝑡𝑢,(],𝑢,(],𝑢,∈𝑥,𝑦,𝑡=Δ𝑢𝑥,𝑦,𝑡+𝑓𝑥,𝑦,𝑡𝑥,𝑦∈Ω,𝑡∈0,𝑇𝑥,𝑦,𝑡=𝜑𝑥,𝑦,𝑡𝑥,𝑦∈𝜕Ω,𝑡∈0,𝑇𝑥,𝑦,0=𝜓𝑥,𝑦𝑥,𝑦Ω, (1) where Ω=(0,𝑙₁)×(0,𝑙₂), Ω=Ω∪𝜕Ω, 𝜕Ω is the boundary of Ω, 𝑇 is a positive constant, Δ is the two-dimensional Laplacian, and ^𝐶₀𝒟^𝛼_𝑡 denotes the temporal Caputo fractional derivative operator defined as ^𝐶₀𝒟^𝛼_𝑡𝑢=1𝑥,𝑦,𝑡Γ(1−𝛼)^𝑡₀𝜕𝑢𝑥,𝑦,𝑠()𝜕𝑠𝑡−𝑠^−𝛼d𝑠,0<𝛼<1. (2)

The main concern in that paper is the construction of the schemes and error estimates in the discrete 𝐻¹-norm. Since 𝐻¹-norm error estimates do not provide immediate insight on the phase error occurring during the time evolution, it is more preferable to give error estimates in the discrete maximum norm when we measure computation errors in practice. In this paper, we continue the investigation of the paper [59], by establishing a maximum norm error estimate for the ADI discretizations. It is known that an 𝐻¹-norm error estimate does not imply a maximum norm error estimate for two-dimensional problems. We here present a technique of 𝐻² discrete energy analysis in order to obtain an explicit maximum norm error estimate.

The outline of the paper is as follows. In Section 2, we derive ADI finite difference schemes for (1) and present our main results of the maximum norm error estimates. The proof of the main error results is given in Section 3. In Section 4, we give some numerical results demonstrating the accuracy of the schemes in the discrete maximum norm. Section 5 contains some concluding remarks.

2. ADI Schemes and Maximum Norm Error Estimates

We partition Ω with nonisotropic uniform mesh sizes ℎ₁ and ℎ₂ in the 𝑥 and 𝑦 directions, respectively. The integers 𝑀₁=𝑙₁/ℎ₁ and 𝑀₂=𝑙₂/ℎ₂. The mesh points (𝑥_𝑖,𝑦_𝑗)=(𝑖ℎ₁,𝑗ℎ₂) (0≤𝑖≤𝑀₁,0≤𝑗≤𝑀₂). Let Ω_ℎ and 𝜕Ω_ℎ be the sets of mesh points lying in Ω and on 𝜕Ω, respectively, and let Ω_ℎ=Ω_ℎ∪𝜕Ω_ℎ. For any grid function 𝑣={𝑣_𝑖,𝑗∣0≤𝑖≤𝑀₁,0≤𝑗≤𝑀₂}, we denote 𝛿_𝑥𝑣_𝑖,𝑗=𝑣_𝑖+1,𝑗−𝑣_𝑖,𝑗ℎ₁,𝛿_𝑦𝑣_𝑖,𝑗=𝑣_𝑖,𝑗+1−𝑣_𝑖,𝑗ℎ₂,𝛿²_𝑥𝑣_𝑖,𝑗=𝛿_𝑥𝑣_𝑖,𝑗−𝛿_𝑥𝑣_{𝑖−1,𝑗}ℎ₁,𝛿²_𝑦𝑣_𝑖,𝑗=𝛿_𝑦𝑣_𝑖,𝑗−𝛿_𝑦𝑣_{𝑖,𝑗−1}ℎ₂,Δ_ℎ𝑣_𝑖,𝑗=𝛿²_𝑥𝑣_𝑖,𝑗+𝛿²_𝑦𝑣_𝑖,𝑗. (3)

For a positive integer 𝑁, we let 𝜏=𝑇/𝑁 be the time step. Define 𝑡_𝑛=𝑛𝜏 and 𝑈^𝑛_𝑖,𝑗=𝑢𝑥_𝑖,𝑦_𝑗,𝑡_𝑛,𝑓^𝑛_𝑖,𝑗=𝑓𝑥_𝑖,𝑦_𝑗,𝑡_𝑛,𝜑^𝑛_𝑖,𝑗=𝜑𝑥_𝑖,𝑦_𝑗,𝑡_𝑛,𝜓_𝑖,𝑗=𝜓𝑥_𝑖,𝑦_𝑗. (4)

For the temporal approximation, we introduce the operators 𝛿_𝑡𝑈^𝑛_𝑖,𝑗=1𝜏𝑈^𝑛+1_𝑖,𝑗−𝑈^𝑛_𝑖,𝑗,𝐷^𝛼_𝑡𝑈^𝑛_𝑖,𝑗=1𝜇𝑈^𝑛_𝑖,𝑗−^𝑛−1_𝑘=1𝑎_{𝑛−𝑘−1}−𝑎_𝑛−𝑘𝑈^𝑘_𝑖,𝑗−𝑎_𝑛−1𝑈⁰_𝑖,𝑗, (5) where 𝜇=𝜏^𝛼Γ(2−𝛼) and 𝑎_𝑘=(𝑘+1)^1−𝛼−𝑘^1−𝛼.

Using the Taylor expansion and the 𝐿₁ approximation of ^𝐶₀𝒟^𝛼_𝑡𝑢(𝑥,𝑦,𝑡) (see [25, 58, 59]), we have ^𝐶₀𝒟^𝛼_𝑡𝑢𝑥_𝑖,𝑦_𝑗,𝑡_𝑛=𝐷^𝛼_𝑡𝑈^𝑛_𝑖,𝑗+𝑅_𝑡^𝑛_𝑖,𝑗,Δ𝑢𝑥_𝑖,𝑦_𝑗,𝑡_𝑛=Δ_ℎ𝑈^𝑛_𝑖,𝑗−𝑅_𝑥^𝑛_𝑖,𝑗−𝑅_𝑦^𝑛_𝑖,𝑗,𝑥_𝑖,𝑦_𝑗∈Ω_ℎ,1≤𝑛≤𝑁, (6) where 𝑅_𝑥^𝑛_𝑖,𝑗=ℎ²₁6¹₀𝜕⁴𝑢𝜕𝑥⁴𝑥_𝑖+𝑠ℎ₁,𝑦_𝑗,𝑡_𝑛+𝜕⁴𝑢𝜕𝑥⁴𝑥_𝑖−𝑠ℎ₁,𝑦_𝑗,𝑡_𝑛(1−𝑠)³𝑅d𝑠,_{𝑦𝑛𝑖,𝑗}=ℎ²₂6¹₀𝜕⁴𝑢𝜕𝑦⁴𝑥_𝑖,𝑦_𝑗+𝑠ℎ₂,𝑡_𝑛+𝜕⁴𝑢𝜕𝑦⁴𝑥_𝑖,𝑦_𝑗−𝑠ℎ₂,𝑡_𝑛(1−𝑠)³d𝑠. (7) There exists a positive constant 𝐶(𝛼) independent of 𝜏, ℎ₁, ℎ₂ and the time level 𝑛 such that 𝑅_𝑡^𝑛_𝑖,𝑗≤𝐶(𝛼)max_{0≤𝑡≤𝑡𝑛}𝜕²𝑢𝜕𝑡²𝑥_𝑖,𝑦_𝑗,𝑡𝜏^2−𝛼. (8) Substituting (6) into (1), we obtain 𝐷^𝛼_𝑡𝑈^𝑛_𝑖,𝑗=Δ_ℎ𝑈^𝑛_𝑖,𝑗+𝑓^𝑛_𝑖,𝑗−𝑅_𝑡^𝑛_𝑖,𝑗−𝑅_𝑥^𝑛_𝑖,𝑗−𝑅_𝑦^𝑛_𝑖,𝑗,𝑥_𝑖,𝑦_𝑗∈Ω_ℎ,1≤𝑛≤𝑁. (9)

2.1. Construction of 𝐿₁-ADI and BD-ADI Schemes

In order to construct an ADI scheme, we add the term 𝜇²𝐷^𝛼_𝑡𝛿²_𝑥𝛿²_𝑦𝑈^𝑛_𝑖,𝑗 to (9). This yields 𝐷^𝛼_𝑡𝑈^𝑛_𝑖,𝑗+𝜇²𝛿²_𝑥𝛿²_𝑦𝑈^𝑛_𝑖,𝑗=Δ_ℎ𝑈^𝑛_𝑖,𝑗+𝑓^𝑛_𝑖,𝑗+𝑅₁^𝑛_𝑖,𝑗,𝑥_𝑖,𝑦_𝑗∈Ω_ℎ,1≤𝑛≤𝑁, (10) where 𝑅₁^𝑛_𝑖,𝑗=𝜇²𝐷^𝛼_𝑡𝛿²_𝑥𝛿²_𝑦𝑈^𝑛_𝑖,𝑗−𝑅_𝑡^𝑛_𝑖,𝑗−𝑅_𝑥^𝑛_𝑖,𝑗−𝑅_𝑦^𝑛_𝑖,𝑗. (11)

It was shown in [59] that |𝜇²𝐷^𝛼_𝑡𝛿²_𝑥𝛿²_𝑦𝑈^𝑛_𝑖,𝑗|≤𝐶_𝐿(𝛼)𝜏^2𝛼 for a positive constant 𝐶_𝐿(𝛼) independent of 𝜏, ℎ₁, ℎ₂ and the time level 𝑛. Thus, there exists a positive constant 𝐶₁ independent of 𝜏, ℎ₁, ℎ₂ and the time level 𝑛 such that 𝑅₁^𝑛_𝑖,𝑗≤𝐶₁𝜏^{min{2𝛼,2−𝛼}}+ℎ²₁+ℎ²₂. (12)

Denote by 𝑢^𝑛_𝑖,𝑗 the finite difference approximation to 𝑈^𝑛_𝑖,𝑗, and let 𝑏_𝑛,𝑘=𝑎_{𝑛−𝑘−1}−𝑎_𝑛−𝑘. After multiplying (10) by 𝜇 and then dropping the term 𝜇(𝑅₁)^𝑛_𝑖,𝑗, we derive a finite difference scheme as follows: 𝐼−𝜇𝛿²_𝑥𝐼−𝜇𝛿²_𝑦𝑢^𝑛_𝑖,𝑗=^𝑛−1_𝑘=1𝑏_𝑛,𝑘𝑢^𝑘_𝑖,𝑗+𝜇²𝛿²_𝑥𝛿²_𝑦𝑢^𝑘_𝑖,𝑗+𝑎_𝑛−1𝑢⁰_𝑖,𝑗+𝜇²𝛿²_𝑥𝛿²_𝑦𝑢⁰_𝑖,𝑗+𝜇𝑓^𝑛_𝑖,𝑗,𝑥_𝑖,𝑦_𝑗∈Ω_ℎ𝑢,1≤𝑛≤𝑁,_{𝑛𝑖,𝑗}=𝜑^𝑛_𝑖,𝑗,𝑥_𝑖,𝑦_𝑗∈𝜕Ω_ℎ𝑢,1≤𝑛≤𝑁,_0𝑖,𝑗=𝜓_𝑖,𝑗,𝑥_𝑖,𝑦_𝑗∈Ω_ℎ. (13)By introducing the intermediate variable 𝑣^𝑛_𝑖,𝑗, we obtain the following 𝐿₁-ADI scheme (see [59]): 𝐼−𝜇𝛿²_𝑥𝑣^𝑛_𝑖,𝑗=^𝑛−1_𝑘=1𝑏_𝑛,𝑘𝑢^𝑘_𝑖,𝑗+𝜇²𝛿²_𝑥𝛿²_𝑦𝑢^𝑘_𝑖,𝑗+𝑎_𝑛−1𝑢⁰_𝑖,𝑗+𝜇²𝛿²_𝑥𝛿²_𝑦𝑢⁰_𝑖,𝑗+𝜇𝑓^𝑛_𝑖,𝑗,𝑥_𝑖,𝑦_𝑗∈Ω_ℎ𝑣,1≤𝑛≤𝑁,_𝑛0,𝑗=𝐼−𝜇𝛿²_𝑦𝜑^𝑛_0,𝑗,𝑣^𝑛_𝑀1,𝑗=𝐼−𝜇𝛿²_𝑦𝜑^𝑛_𝑀1,𝑗,0≤𝑗≤𝑀₂𝑢,1≤𝑛≤𝑁,_0𝑖,𝑗=𝜓_𝑖,𝑗,𝑥_𝑖,𝑦_𝑗∈Ω_ℎ,𝐼−𝜇𝛿²_𝑦𝑢^𝑛_𝑖,𝑗=𝑣^𝑛_𝑖,𝑗,𝑥_𝑖,𝑦_𝑗∈Ω_ℎ𝑢,1≤𝑛≤𝑁,_𝑛𝑖,0=𝜑^𝑛_𝑖,0,𝑢^𝑛_𝑖,𝑀2=𝜑^𝑛_𝑖,𝑀2,0≤𝑖≤𝑀₁,1≤𝑛≤𝑁. (14)

Adding the term 𝜏𝜇𝛿²_𝑥𝛿²_𝑦𝛿_𝑡𝑈^𝑛_𝑖,𝑗 to (9), we have 𝐷^𝛼_𝑡𝑈^𝑛_𝑖,𝑗+𝜏𝜇𝛿²_𝑥𝛿²_𝑦𝛿_𝑡𝑈^𝑛_𝑖,𝑗=Δ_ℎ𝑈^𝑛_𝑖,𝑗+𝑓^𝑛_𝑖,𝑗+𝑅₂^𝑛_𝑖,𝑗, (15) where 𝑅₂^𝑛_𝑖,𝑗=𝜏𝜇𝛿²_𝑥𝛿²_𝑦𝛿_𝑡𝑈^𝑛_𝑖,𝑗−𝑅_𝑡^𝑛_𝑖,𝑗−𝑅_𝑥^𝑛_𝑖,𝑗−𝑅_𝑦^𝑛_𝑖,𝑗. (16)

Since |𝜏𝜇𝛿²_𝑥𝛿²_𝑦𝛿_𝑡𝑈^𝑛_𝑖,𝑗|≤𝐶_𝐵𝐷(𝛼)𝜏^1+𝛼 for a positive constant 𝐶_𝐵𝐷(𝛼) independent of 𝜏, ℎ₁, ℎ₂ and the time level 𝑛 (see [59]), there exists a positive constant 𝐶₂ independent of 𝜏, ℎ₁, ℎ₂ and the time level 𝑛 such that 𝑅₂^𝑛_𝑖,𝑗≤𝐶₂𝜏^{min{1+𝛼,2−𝛼}}+ℎ²₁+ℎ²₂. (17)

By (15), we obtain the following finite difference scheme: 𝐼−𝜇𝛿²_𝑥𝐼−𝜇𝛿²_𝑦𝑢^𝑛_𝑖,𝑗=𝜇²𝛿²_𝑥𝛿²_𝑦𝑢^𝑛−1_𝑖,𝑗+^𝑛−1_𝑘=1𝑏_𝑛,𝑘𝑢^𝑘_𝑖,𝑗+𝑎_𝑛−1𝑢⁰_𝑖,𝑗+𝜇𝑓^𝑛_𝑖,𝑗,𝑥_𝑖,𝑦_𝑗∈Ω_ℎ𝑢,1≤𝑛≤𝑁,_{𝑛𝑖,𝑗}=𝜑^𝑛_𝑖,𝑗,𝑥_𝑖,𝑦_𝑗∈𝜕Ω_ℎ𝑢,1≤𝑛≤𝑁,_0𝑖,𝑗=𝜓_𝑖,𝑗,𝑥_𝑖,𝑦_𝑗∈Ω_ℎ. (18)

It is equivalent to the following BD-ADI scheme (see [59]): 𝐼−𝜇𝛿²_𝑥𝑣^𝑛_𝑖,𝑗=𝜇²𝛿²_𝑥𝛿²_𝑦𝑢^𝑛−1_𝑖,𝑗+^𝑛−1_𝑘=1𝑏_𝑛,𝑘𝑢^𝑘_𝑖,𝑗+𝑎_𝑛−1𝑢⁰_𝑖,𝑗+𝜇𝑓^𝑛_𝑖,𝑗,𝑥_𝑖,𝑦_𝑗∈Ω_ℎ𝑣,1≤𝑛≤𝑁,_𝑛0,𝑗=𝐼−𝜇𝛿²_𝑦𝜑^𝑛_0,𝑗,𝑣^𝑛_𝑀1,𝑗=𝐼−𝜇𝛿²_𝑦𝜑^𝑛_𝑀1,𝑗,0≤𝑗≤𝑀₂𝑢,1≤𝑛≤𝑁,_0𝑖,𝑗=𝜓_𝑖,𝑗,𝑥_𝑖,𝑦_𝑗∈Ω_ℎ,𝐼−𝜇𝛿²_𝑦𝑢^𝑛_𝑖,𝑗=𝑣^𝑛_𝑖,𝑗,𝑥_𝑖,𝑦_𝑗∈Ω_ℎ𝑢,1≤𝑛≤𝑁,_𝑛𝑖,0=𝜑^𝑛_𝑖,0,𝑢^𝑛_𝑖,𝑀2=𝜑^𝑛_𝑖,𝑀2,0≤𝑖≤𝑀₁,1≤𝑛≤𝑁. (19)

2.2. Maximum Norm Error Estimates

For any grid function 𝑣={𝑣_𝑖,𝑗∣0≤𝑖≤𝑀₁,0≤𝑗≤𝑀₂}, we define its maximum norm ‖𝑣‖_∞ by ‖𝑣‖_∞=max_{0≤𝑖≤𝑀1,0≤𝑗≤𝑀2}𝑣_𝑖,𝑗. (20)

Let 𝑈^𝑛_𝑖,𝑗 be the value of the solution of (1) at the mesh point (𝑥_𝑖,𝑦_𝑗,𝑡_𝑛), and let 𝑢^𝑛_𝑖,𝑗 be the solution of the 𝐿₁-ADI scheme (14) or the BD-ADI scheme (19). We now present our main results of the maximum norm estimate for the error 𝑒^𝑛_𝑖,𝑗=𝑈^𝑛_𝑖,𝑗−𝑢^𝑛_𝑖,𝑗 in the following two theorems. Their proofs will be given in the next section.

Theorem 1. Assume that the solution 𝑢(𝑥,𝑦,𝑡) of (1) is sufficiently smooth, and let 𝑢^𝑛_𝑖,𝑗 be the solution of the 𝐿₁-ADI scheme (14). Then 𝜏^𝑛_𝑙=1𝑈^𝑙−𝑢^𝑙_∞^1/2≤𝐶^∗₁𝜏^{min{2𝛼,2−𝛼}}+ℎ²₁+ℎ²₂,1≤𝑛≤𝑁, (21) where 𝐶^∗₁=(𝑙²₀/2)√𝑇𝐶₁ and 𝑙₀=max{𝑙₁,𝑙₂}.

Theorem 2. Assume that the solution 𝑢(𝑥,𝑦,𝑡) of (1) is sufficiently smooth, and let 𝑢^𝑛_𝑖,𝑗 be the solution of the BD-ADI scheme (19). Then 𝜏^𝑛_𝑙=1𝑈^𝑙−𝑢^𝑙_∞^1/2≤𝐶^∗₂𝜏^{min{1+𝛼,2−𝛼}}+ℎ²₁+ℎ²₂,1≤𝑛≤𝑁, (22) where 𝐶^∗₂=(𝑙²₀/2)√𝑇𝐶₂ and 𝑙₀=max{𝑙₁,𝑙₂}.

Theorems 1 and 2 show that the ADI difference solution 𝑢^𝑛_𝑖,𝑗 from (14) or (19) converges to the analytical solution 𝑈^𝑛_𝑖,𝑗 of (1) in the discrete maximum norm. We also see from the estimates (12), (17), (21), and (22) that the 𝐿₁-ADI scheme (14) and the BD-ADI scheme (19) have the same order as their truncation error with respect to the discrete maximum norm.

3. Proof of the Main Results

Let 𝒱_ℎ be the set of all grid functions 𝑣={𝑣_𝑖,𝑗∣0≤𝑖≤𝑀₁,0≤𝑗≤𝑀₂} defined in Ω_ℎ and vanishing on 𝜕Ω_ℎ. For arbitrary 𝑣,𝑤∈𝒱_ℎ, we define the following inner products: (𝑣,𝑤)=ℎ₁ℎ₂^{𝑀1−1𝑀}_𝑖=1²⁻¹_𝑗=1𝑣_𝑖,𝑗𝑤_𝑖,𝑗,𝛿_𝑥𝑣,𝛿_𝑥𝑤_𝑥=ℎ₁ℎ₂^{𝑀1−1𝑀}_𝑖=0²⁻¹_𝑗=1𝛿_𝑥𝑣_𝑖,𝑗𝛿_𝑥𝑤_𝑖,𝑗,𝛿_𝑦𝑣,𝛿_𝑦𝑤_𝑦=ℎ₁ℎ₂^{𝑀1−1𝑀}_𝑖=1²⁻¹_𝑗=0𝛿_𝑦𝑣_𝑖,𝑗𝛿_𝑦𝑤_𝑖,𝑗. (23)

For any 𝑣∈𝒱_ℎ, we introduce the following norms: ‖𝑣‖=(𝑣,𝑣)^1/2,𝛿_𝜈𝑣_𝜈=𝛿_𝜈𝑣,𝛿_𝜈𝑣_𝜈^1/2,|𝑣|𝜈=𝑥,𝑦₁=𝛿_𝑥𝑣²_𝑥+𝛿_𝑦𝑣²_𝑦^1/2,𝛿_𝑥𝛿²_𝑦𝑣_𝑥=𝛿_𝑥𝛿²_𝑦𝑣,𝛿_𝑥𝛿²_𝑦𝑣_𝑥^1/2,𝛿_𝑦𝛿²_𝑥𝑣_𝑦=𝛿_𝑦𝛿²_𝑥𝑣,𝛿_𝑦𝛿²_𝑥𝑣_𝑦^1/2. (24)

Using a simple calculation, we have that for arbitrary 𝑣,𝑤∈𝒱_ℎ, 𝛿²_𝜈𝑣,𝑤=−𝛿_𝜈𝑣,𝛿_𝜈𝑤_𝜈,𝜈=𝑥,𝑦. (25)

Before proving Theorems 1 and 2, we first introduce the following embedding theorem from [60, page 281].

Lemma 3. For any 𝑣∈𝒱_ℎ, one has ‖𝑣‖_∞≤𝑙²₀2𝑙₁𝑙₂⁻¹Δ_ℎ𝑣_∞,𝑙₀=max𝑙₁,𝑙₂. (26)

Proof of Theorem 1. Let 𝑒^𝑛_𝑖,𝑗=𝑈^𝑛_𝑖,𝑗−𝑢^𝑛_𝑖,𝑗. Then by (10) and (13), 𝑒^𝑛_𝑖,𝑗−𝜇Δ_ℎ𝑒^𝑛_𝑖,𝑗+𝜇²𝛿²_𝑥𝛿²_𝑦𝑒^𝑛_𝑖,𝑗=𝑤^𝑛−1_𝑖,𝑗+𝜇²𝛿²_𝑥𝛿²_𝑦𝑤^𝑛−1_𝑖,𝑗+𝜇𝑅₁^𝑛_𝑖,𝑗,𝑥_𝑖,𝑦_𝑗∈Ω_ℎ𝑒,1≤𝑛≤𝑁,_{𝑛𝑖,𝑗}=0,𝑥_𝑖,𝑦_𝑗∈𝜕Ω_ℎ𝑒,1≤𝑛≤𝑁,_0𝑖,𝑗=0,𝑥_𝑖,𝑦_𝑗∈Ω_ℎ, (27) where 𝑤^𝑛−1_𝑖,𝑗=^𝑛−1_𝑘=1𝑏_𝑛,𝑘𝑒^𝑘_𝑖,𝑗+𝑎_𝑛−1𝑒⁰_𝑖,𝑗. (28)
This implies that 𝑒^𝑛={𝑒^𝑛_𝑖,𝑗∣0≤𝑖≤𝑀₁,0≤𝑗≤𝑀₂}∈𝒱_ℎ for each 𝑛. Taking the inner product of the first equation in (27) and −Δ_ℎ𝑒^𝑛, we get 𝑒^𝑛+𝜇²𝛿²_𝑥𝛿²_𝑦𝑒^𝑛,−Δ_ℎ𝑒^𝑛+−𝜇Δ_ℎ𝑒^𝑛,−Δ_ℎ𝑒^𝑛=𝑤^𝑛−1+𝜇²𝛿²_𝑥𝛿²_𝑦𝑤^𝑛−1,−Δ_ℎ𝑒^𝑛+𝜇𝑅₁^𝑛,−Δ_ℎ𝑒^𝑛,1≤𝑛≤𝑁. (29)
It follows from (25) that 𝑒^𝑛+𝜇²𝛿²_𝑥𝛿²_𝑦𝑒^𝑛,−Δ_ℎ𝑒^𝑛=𝑒^𝑛2₁+𝜇²𝛿_𝑦𝛿²_𝑥𝑒^𝑛2_𝑦+𝛿_𝑥𝛿²_𝑦𝑒^𝑛2_𝑥. (30)
Similarly, 𝑤^𝑛−1+𝜇²𝛿²_𝑥𝛿²_𝑦𝑤^𝑛−1,−Δ_ℎ𝑒^𝑛=^𝑛−1_𝑘=1𝑏_𝑛,𝑘𝛿_𝑥𝑒^𝑘,𝛿_𝑥𝑒^𝑛_𝑥+𝛿_𝑦𝑒^𝑘,𝛿_𝑦𝑒^𝑛_𝑦+𝜇²𝛿_𝑦𝛿²_𝑥𝑒^𝑘,𝛿_𝑦𝛿²_𝑥𝑒^𝑛_𝑦+𝜇²𝛿_𝑥𝛿²_𝑦𝑒^𝑘,𝛿_𝑥𝛿²_𝑦𝑒^𝑛_𝑥. (31)
Since 𝑏_𝑛,𝑘=𝑎_{𝑛−𝑘−1}−𝑎_𝑛−𝑘>0,^𝑛−1_𝑘=1𝑏_𝑛,𝑘≤1, (32) we have from Cauchy-Schwarz inequality that 𝑤^𝑛−1+𝜇²𝛿²_𝑥𝛿²_𝑦𝑤^𝑛−1,−Δ_ℎ𝑒^𝑛≤12^𝑛−1_𝑘=1𝑏_𝑛,𝑘𝑒^𝑘2₁+𝜇²𝛿_𝑦𝛿²_𝑥𝑒^𝑘2_𝑦+𝛿_𝑥𝛿²_𝑦𝑒^𝑘2_𝑥+12𝑒^𝑛2₁+𝜇²𝛿_𝑦𝛿²_𝑥𝑒^𝑛2_𝑦+𝛿_𝑥𝛿²_𝑦𝑒^𝑛2_𝑥,𝜇𝑅₁^𝑛,−Δ_ℎ𝑒^𝑛≤𝜇2𝑅₁^𝑛2+Δ_ℎ𝑒^𝑛2. (33)
Substituting (30) and (33) into (29) gives the following: 𝜇Δ_ℎ𝑒^𝑛2+𝐸^𝑛≤𝐸^𝑛−1+𝜇𝑅₁^𝑛2,1≤𝑛≤𝑁, (34) where 𝐸^𝑛=^𝑛_𝑘=1𝑎_𝑛−𝑘𝑒^𝑘2₁+𝜇²𝛿_𝑦𝛿²_𝑥𝑒^𝑘2_𝑦+𝛿_𝑥𝛿²_𝑦𝑒^𝑘2_𝑥,𝐸⁰=0. (35)
This implies that 𝜇^𝑛_𝑙=1Δ_ℎ𝑒^𝑙2+𝐸^𝑛≤𝜇^𝑛_𝑙=1𝑅₁^𝑙2,1≤𝑛≤𝑁. (36)
Since 𝐸^𝑛≥0 and by (12), ‖(𝑅₁)^𝑙‖²≤𝑙₁𝑙₂𝐶²₁(𝜏^{min{2𝛼,2−𝛼}}+ℎ²₁+ℎ²₂)², the estimate (21) follows from (36) and (26) immediately.

Proof of Theorem 2. The proof follows from the similar argument as that in the proof of Theorem 1 and we give a sketch. Let 𝑒^𝑛_𝑖,𝑗=𝑈^𝑛_𝑖,𝑗−𝑢^𝑛_𝑖,𝑗. By (15) and (18), 𝑒^𝑛_𝑖,𝑗−𝜇Δ_ℎ𝑒^𝑛_𝑖,𝑗+𝜇²𝛿²_𝑥𝛿²_𝑦𝑒^𝑛_𝑖,𝑗=𝜇²𝛿²_𝑥𝛿²_𝑦𝑒^𝑛−1_𝑖,𝑗+𝑤^𝑛−1_𝑖,𝑗+𝜇𝑅₂^𝑛_𝑖,𝑗,𝑥_𝑖,𝑦_𝑗∈Ω_ℎ𝑒,1≤𝑛≤𝑁,_{𝑛𝑖,𝑗}=0,𝑥_𝑖,𝑦_𝑗∈𝜕Ω_ℎ𝑒,1≤𝑛≤𝑁,_0𝑖,𝑗=0,𝑥_𝑖,𝑦_𝑗∈Ω_ℎ, (37) where 𝑤^𝑛−1_𝑖,𝑗 is defined by (28). Taking the inner product of the first equation in (37) and −Δ_ℎ𝑒^𝑛, we get 𝑒^𝑛,−Δ_ℎ𝑒^𝑛+−𝜇Δ_ℎ𝑒^𝑛,−Δ_ℎ𝑒^𝑛+𝜇²𝛿²_𝑥𝛿²_𝑦𝑒^𝑛−𝑒^𝑛−1,−Δ_ℎ𝑒^𝑛=𝑤^𝑛−1,−Δ_ℎ𝑒^𝑛+𝜇𝑅₂^𝑛,−Δ_ℎ𝑒^𝑛,1≤𝑛≤𝑁. (38)
Since 𝑒^𝑛,−Δ_ℎ𝑒^𝑛=𝑒^𝑛2₁,−𝜇Δ_ℎ𝑒^𝑛,−Δ_ℎ𝑒^𝑛=𝜇Δ_ℎ𝑒^𝑛2,𝜇²𝛿²_𝑥𝛿²_𝑦𝑒^𝑛−𝑒^𝑛−1,−Δ_ℎ𝑒^𝑛≥𝜇²2𝛿_𝑦𝛿²_𝑥𝑒^𝑛2_𝑦+𝛿_𝑥𝛿²_𝑦𝑒^𝑛2_𝑥−𝛿_𝑦𝛿²_𝑥𝑒^𝑛−12_𝑦−𝛿_𝑥𝛿²_𝑦𝑒^𝑛−12_𝑥,𝑤^𝑛−1,−Δ_ℎ𝑒^𝑛≤12^𝑛−1_𝑘=1𝑏_𝑛,𝑘𝑒^𝑘2₁+𝑒^𝑛2₁,𝜇𝑅₂^𝑛,−Δ_ℎ𝑒^𝑛≤𝜇2𝑅₂^𝑛2+Δ_ℎ𝑒^𝑛2, (39) we obtain that 𝜇Δ_ℎ𝑒^𝑛2+𝐹^𝑛≤𝐹^𝑛−1+𝜇𝑅₂^𝑛2,1≤𝑛≤𝑁, (40) where 𝐹^𝑛=^𝑛_𝑘=1𝑎_𝑛−𝑘𝑒^𝑘2₁+𝜇²𝛿_𝑦𝛿²_𝑥𝑒^𝑛2_𝑦+𝛿_𝑥𝛿²_𝑦𝑒^𝑛2_𝑥,𝐹⁰=0. (41)
Therefore 𝜇^𝑛_𝑙=1Δ_ℎ𝑒^𝑙2+𝐹^𝑛≤𝜇^𝑛_𝑙=1𝑅₂^𝑙2,1≤𝑛≤𝑁. (42)
The estimate (22) follows from the previous inequality, (17), and (26) immediately.

4. Numerical Results

In this section, we give some numerical results to demonstrate the accuracy of the 𝐿₁-ADI scheme (14) and the BD-ADI scheme (19) with respect to the discrete maximum norm. Some detailed numerical comparisons of these two ADI schemes with the implicit scheme proposed in [43] can be found in [59]. The dependence of the solution of (1) on the anomalous diffusion exponent 𝛼 was also exhibited in [59] through some numerical results.

In our numerical computations, we take an equal mesh size in each of the space directions; that is, ℎ₁=ℎ₂=ℎ. We compute the discrete maximum norm error of the numerical solution 𝑢^𝑛_𝑖,𝑗 by 𝖾𝗋𝗋𝗈𝗋(𝜏,ℎ)=𝜏^𝑁_𝑙=1𝑈^𝑙−𝑢^𝑙_∞^1/2, (43) and its convergence orders by 𝗈𝗋𝖽𝖾𝗋(𝜏)=log₂𝖾𝗋𝗋𝗈𝗋(2𝜏,ℎ)(),(ℎ)𝖾𝗋𝗋𝗈𝗋𝜏,ℎ𝗈𝗋𝖽𝖾𝗋=log₂𝖾𝗋𝗋𝗈𝗋(𝜏,2ℎ)(),𝖾𝗋𝗋𝗈𝗋𝜏,ℎ (44) where 𝑈^𝑛_𝑖,𝑗 represents the value of the exact analytic solution of (1) at (𝑥_𝑖,𝑦_𝑗,𝑡_𝑛).

Example 4. We consider the problem (1) in the domain Ω=(0,𝜋)×(0,𝜋), and let 𝑇=1/2. Assume that the solution 𝑢(𝑥,𝑦,𝑡) of this problem is 𝑢(𝑥,𝑦,𝑡)=𝑡²sin(𝑥)sin(𝑦). It can be checked that the corresponding known functions are given by 𝑓(𝑥)𝑦𝑡𝑥,𝑦,𝑡=2sinsin^2−𝛼Γ(3−𝛼)+𝑡²,𝜑𝑥,𝑦,𝑡=𝑡²sin(𝑥)sin𝑦,𝜓𝑥,𝑦=0. (45)

In Table 1, we present the maximum norm error 𝖾𝗋𝗋𝗈𝗋(𝜏,ℎ) and the temporal convergence order 𝗈𝗋𝖽𝖾𝗋(𝜏) of the numerical solution 𝑢^𝑛_𝑖,𝑗 by the 𝐿₁-ADI scheme (14) and the BD-ADI scheme (19). We see that the 𝐿₁-ADI scheme has the temporal accuracy of min{2𝛼,2−𝛼} in the discrete maximum norm, and the best temporal accuracy is attained at 𝛼=2/3. It is also seen that the BD-ADI scheme generates the temporal accuracy of min{1+𝛼,2−𝛼} in the discrete maximum norm, and it gets the best temporal accuracy when 𝛼=1/2. These observations coincide well with the theoretical analysis.

Table 1 

The maximum norm error and the temporal accuracy of 𝑢^𝑛_𝑖,𝑗 for Example 4.

Table 2 gives the maximum norm error 𝖾𝗋𝗋𝗈𝗋(𝜏,ℎ) and the spatial convergence order 𝗈𝗋𝖽𝖾𝗋(ℎ) of the numerical solution 𝑢^𝑛_𝑖,𝑗 by the 𝐿₁-ADI scheme (14) and the BD-ADI scheme (19). As expected from the theoretical analysis, these two schemes have the second-order spatial accuracy.

Table 2 

The maximum norm error and the spatial accuracy of 𝑢^𝑛_𝑖,𝑗 for Example 4.

Example 5. We consider the subdiffusion equation 𝑢_𝑡=𝑥,𝑦,𝑡₀𝒟_𝑡^1−𝛼Δ𝑢𝑥,𝑦,𝑡+𝑒^𝑥+𝑦(1+𝛼)𝑡^𝛼−2Γ(2+𝛼)Γ()𝑡1+2𝛼^2𝛼,𝑥,𝑦∈Ω,𝑡∈[0,1], (46) with boundary and initial conditions 𝑢𝑥,𝑦,𝑡=𝑒^𝑥+𝑦𝑡^1+𝛼,(],𝑢∈𝑥,𝑦∈𝜕Ω,𝑡∈0,1𝑥,𝑦,0=0,𝑥,𝑦Ω, (47) where Ω=(0,1)×(0,1) and ₀𝒟_𝑡^1−𝛼 denotes the Riemann-Liouville fractional derivative operator defined as ₀𝒟_𝑡^1−𝛼𝑦(𝑡)=1Γ(𝛼)dd𝑡^𝑡₀𝑦(𝑠)(𝑡−𝑠)^1−𝛼d𝑠,0<𝛼<1. (48)

Operating Riemann-Liouville fractional derivative operator ₀𝒟_𝑡^𝛼−1 on both sides of (46), we obtain the equivalent problem of the subdiffusion equation with the temporal Caputo fractional derivative [40, 42, 59] as follows: ^𝐶₀𝒟^𝛼_𝑡𝑢𝑥,𝑦,𝑡=Δ𝑢𝑥,𝑦,𝑡+𝑒^𝑥+𝑦Γ(2+𝛼)𝑡−2𝑡^1+𝛼,[],𝑢𝑥,𝑦∈Ω,𝑡∈0,1𝑥,𝑦,𝑡=𝑒^𝑥+𝑦𝑡^1+𝛼,(],𝑢∈𝑥,𝑦∈𝜕Ω,𝑡∈0,1𝑥,𝑦,0=0,𝑥,𝑦Ω. (49)

We now solve the aforementioned problem by the 𝐿₁-ADI scheme (14) and the BD-ADI scheme (19). Tables 3 and 4 give the maximum norm error 𝖾𝗋𝗋𝗈𝗋(𝜏,ℎ) and the convergence orders 𝗈𝗋𝖽𝖾𝗋(𝜏) and 𝗈𝗋𝖽𝖾𝗋(ℎ) of the numerical solution 𝑢^𝑛_𝑖,𝑗. It is easily seen that the numerical results confirm the theoretical analysis results.

Table 3 

The maximum norm error and the temporal accuracy of 𝑢^𝑛_𝑖,𝑗 for Example 5.

Table 4 

The maximum norm error and the spatial accuracy of 𝑢^𝑛_𝑖,𝑗 for Example 5.

5. Conclusions

We have studied two ADI finite difference methods for a two-dimensional fractional subdiffusion equation. An explicit error estimate for each of the two methods has been provided in the discrete maximum norm. It has been shown that the methods have the same order as their truncation errors with respect to the discrete maximum norm. The maximum norm error estimates presented here are more preferable for measuring computation errors in practice, compared to the 𝐻¹-norm error estimates in [59]. Numerical results have confirmed the theoretical analysis results.

Acknowledgments

This work was supported in part by E-Institutes of Shanghai Municipal Education Commission no. E03004 and Shanghai Leading Academic Discipline Project no. B407.