Compact alternating direction implicit method for two-dimensional time fractional diffusion equation

Abstract

High-order compact finite difference scheme with operator splitting technique for solving two-dimensional time fractional diffusion equation is considered in this paper. A Grünwald–Letnikov approximation is used for the Riemann–Liouville time derivative, and the second order spatial derivatives are approximated by the compact finite differences to obtain a fully discrete implicit scheme. Alternating direction implicit (ADI) method is used to split the original problem into two separate one-dimensional problems. The local truncation error is analyzed and the stability is discussed by the Fourier method. The proposed scheme is suitable when the order of the time fractional derivative γ lies in the interval $\left[\frac{1}{2}\text{,}1\right)$. A correction term is added to maintain high accuracy when $\gamma \in \left(0\text{,}\frac{1}{2}\right)$. Numerical results are provided to verify the accuracy and efficiency of the proposed algorithm.

Introduction

Fractional differential equations (FDEs) have been intensely studied in the recent years, as their applications can be found in physical, biological, geological and financial systems. Fractional diffusion equation is used to describe phenomena of anomalous diffusion in transport processes, and detailed discussion can be found in the review article [1] and monograph [2].

Several numerical methods have been proposed for solving the space and/or time FDEs up to now, see, e.g., [3], [4], [5], [6], [7], [8], [9], [10], [11], [12]. These schemes are convergent with order of two for the space variable(s), and it is interesting to seek high-order numerical methods for FDEs. Compact finite difference schemes have the advantage of high-order of accuracy and the desirable tridiagonal nature of the finite-difference equations (see papers [13], [14]), and recently, one-dimensional fractional sub-diffusion equation was solved by the compact finite difference scheme with convergence order $\mathcal{O}(\tau +h^4)$ in [15], and a higher order $\mathcal{O}({\tau }^{2-\gamma }+h^4)$ one can be found in [16].

For solving fractional equations numerically, many research papers are focused on one-dimensional problems. For the two-dimensional cases, there are papers [17], [18] where the fractional derivatives were for the space variables, and for the time fractional diffusion problem, a kernel-basis representation for the spatial variables with both adaptive time stepping and adaptive spatial basis selection approach was considered recently in [19].

The main purpose of this paper is to solve the multidimensional fractional diffusion problem using the compact difference scheme with the operator-splitting techniques, that is, using the compact ADI scheme. ADI methods are approximate factorization methods which replace the solutions of multidimensional problems by sequences of one-dimensional problems. As the resulting algebraic systems are tridiagonal so they can be solved easily by the Thomas algorithm, thus these methods can reduce the computational cost greatly. They were proposed by Peaceman, Rachford and Douglas in 1950’s, see papers [20], [21], [22], [23], and the book [24] for details. Compact ADI schemes for partial differential equations of integer order were considered in [25], [26], and most recently in [27], [28], [29]. For fractional-order differential equations, two ADI schemes were given in [30] with the spatial fractional derivatives, and a fast alternating-direction implicit finite difference method for space-fractional diffusion equations was given in [31]. The two-dimensional case is discussed in this paper and extension to 3D problems can be given in a similar way. As stated above, the proposed algorithm has high accuracy and good efficiency.

The model problem considered here is the two-dimensional time fractional diffusion equation describing subdiffusive phenomena with a non-homogeneous term,

$$\frac{\partial u(x\text{,}y\text{,}t)}{\partial t}={}_0{D}_t^{1-\gamma }\left[K_1\frac{{\partial }^{2}u}{\partial x^2}+K_2\frac{{\partial }^{2}u}{\partial y^2}\right]+f(x\text{,}y\text{,}t)\text{,}(x\text{,}y)\in (L_1\text{,}{L}_2)\times (L_3\text{,}{L}_4)\text{,}0<t<T\text{,}$$

where K₁ and K₂ are positive constants, and ${}_0{D}_t^{1-\gamma }v(0<\gamma <1)$ denotes the Riemann–Liouville fractional derivative of the function v(x, y, t), defined by [2], i.e.,

$${}_0{D}_t^{1-\gamma }v=\frac{1}{\Gamma (\gamma )}\frac{\partial }{\partial t}{\int }_0^t\frac{v(x\text{,}y\text{,}\tau )}{(t-\tau {)}^{1-\gamma }}d\tau \text{,}$$

and in the limit γ → 1, Eq. (1) corresponds to the ordinary (or Browning) diffusion equation.

The initial condition for (1) is

$$u(x\text{,}y\text{,}0)=w(x\text{,}y)\text{,}(x\text{,}y)\in (L_1\text{,}{L}_2)\times (L_3\text{,}{L}_4)$$

and the Dirichlet boundary condition for (1) is given by

$$u(L_1\text{,}y\text{,}t)={\phi }_1(y\text{,}t)\text{,}u(L_2\text{,}y\text{,}t)={\phi }_2(y\text{,}t)\text{,}u(x\text{,}{L}_3\text{,}t)={\psi }_1(x\text{,}t)\text{,}u(x\text{,}{L}_4\text{,}t)={\psi }_2(x\text{,}t)\text{,}t⩾0\text{.}$$

The paper is organized as follows. In Section 2, after approximating the second order spatial derivatives by the compact finite differences and using the Grünwald–Letnikov discretization for the time fractional derivative, we give an implicit compact ADI difference scheme. Challenging and interesting, we find that we must consider some terms that can be discarded for differential equations of integer order. This makes both the discussion and computation difficult. The local truncation error is discussed and the stability is analyzed using the Fourier analysis in Section 3. The local truncation error shows the loss of order for $\gamma \in \left(0\text{,}\frac{1}{2}\right)$ and the scheme is proved unconditionally stable for any γ ∈ (0, 1). We find it advantageous to add a correction term in our original compact ADI scheme for $\gamma \in \left(0\text{,}\frac{1}{2}\right)$. Finally, some numerical examples are given in Section 4 to verify the theoretical conclusions. This paper closes with a summary in Section 5. The correction for a mistake in a previous paper [15] is given in the appendix.