Spectral methods for two-dimensional space and time fractional Bloch-Torrey equations

In this paper, we consider the numerical approximation of the space and time fractional Bloch-Torrey equations. A fully discrete spectral scheme based on a finite difference method in the time direction and a Galerkin-Legendre spectral method in the space direction is developed. In order to reduce the amount of computation, an alternating direction implicit (ADI) spectral scheme is proposed. Then the stability and convergence analysis are rigorously established. Finally, numerical results are presented to support our theoretical analysis.


1.
Introduction. Although fractional calculus has a history almost as long as the integer calculus, only in the past decades does the fractional calculus undergo a rapid development. Recently, fractional differential equations can be used to describe some phenomena or process with hereditary and memory in physics and chemical processes, materials, control theory, biology, finance, and so on (see [7,11,13]). In physics, fractional derivatives are used to model anomalous diffusion, where diffusing particles are not locally homogeneous and the disorder is not well-approximated by assuming a unified change in diffusion constant (see [6,10]).
In this paper, we consider the following two-dimensional space and time fractional Bloch-Torrey equations ( [8,9]): with the initial and boundary conditions given by u(x, y, 0) = u 0 (x, y), (x, y) ∈ Ω, (2) u(x, y, t)| (x,y)∈∂Ω = 0, t ∈ (0, T ], where 0 < α < 1, 1/2 < β < 1, K x , K y > 0, Ω = (a, b) × (c, d). The functions f (x, y, t) and u 0 (x, y) are known smooth functions. The Caputo fractional derivative C 0 D α t u(x, y, t) is defined by This space and time fractional Bloch-Torrey equations are a class of new fractional diffusion models. This model is more appropriate than the classical diffusion model to be applied to analysis the diffusion images of human brain tissues and provide new insights for further investigations of issue structures.
For the space and time fractional Bloch-Torrey equations, Yu et al. [4,17] established the fundamental solutions and exact solutions. Yu et al. [18,19] proposed a finite difference method for 3D and 2D space and time fractional Bloch-Torrey equations. Bu et al. [3,20] proposed a finite element method for two-dimensional space and time fractional Bloch-Torrey equations and analysed its stability and convergence. In this paper, we aim to solve the problem (1)-(3) using a fully discrete numerical method. The method is based on a finite difference scheme in time and spectral approximation using Legendre functions in space. We give a detailed analysis for the stability and convergence of the fully discrete scheme. The scheme is unconditionally stable and convergent with order O(τ 2 + N β−r ).
The remainder of the paper is organized as follows. In Section 2, some preliminaries and notations are shown. In section 3, we propose a fully discrete scheme for the problem (1)-(3). In section 4, the stability and convergence of the fully discrete scheme are given. In section 5, we propose an ADI spectral method, and analyse its stability and convergence. We do some numerical experiments in Section 6. Finally, some conclusions are given in Section 7.
2. Preliminaries and notations. In this section, we show some notation and lemmas that are need in the following sections.
Let Ω be a finite domain satisfying Ω = I x × I y = (a, b) × (c, d), and denote by (·, ·) the inner product on the space L 2 (Ω) with the L 2 -norm · L 2 (Ω) and the maximum norm · L ∞ (Ω) . If it does not cause confusion, we also define (·, ·) as the inner product on the interval I x or I y . Let µ be a nonnegative real number, and denote by H µ (Ω) and H µ 0 (Ω) the usual Sobolev spaces with the norm · H µ (Ω) and the seminorm | · | H µ (Ω) .
Firstly, we introduce some spaces that are used in the formulation of the numerical algorithms. We first introduce the space J µ L (Ω), J µ R (Ω) and J µ S (Ω) in R 2 (see [5]). , and denote J µ L (Ω) (or J µ L,0 (Ω)) as the closure of C ∞ (Ω) (or C ∞ 0 (Ω)) with respect to · J µ L (Ω) , where C ∞ 0 (Ω) is the space of smooth functions with compact support in Ω.
, whereû is the extension of u by zero outside Ω. Furthermore, if µ = n−1/2, n ∈ N, and u ∈ J µ L,0 (Ω), then there exists a positive constant C independent of u such that |û| J µ L (R 2 ) ≤ C|u| J µ L (Ω) . Similarly to Lemma 3.1.4 in [14], we can obtain the following lemma.
whereû is the extension of u by zero outside Ω.
Then there exist positive constants C 1 and C 2 independent of u such that and H µ (R 2 ) are equivalence with equivalent norms and seminorms. If µ > 0, and H µ 0 (R 2 ) are equivalence with equivalent norms and seminorms.
c,y u , where C 1 , C 2 , C 3 and C 4 are positive constants independent of u.
In the following, C or C j denotes a generic positive constant independent τ, N , and the constant C will not be the same. Lemma 2.10. Let 1 < β < 2. Then for any u ∈ H β 0 (Ω) and u ∈ H Then, we introduce some function spaces. As in [15], the function spaces V x,0 N and V y,0 N can be expressed as is the Legendre polynomial defined by the following recurrence relation in [16]: The explicit formula of the Jacobi polynomial is stated as follows (see [16]): The Jacobi polynomials can also be generated by the three-term recurrence formula (see [16] for more details). Therefore, the function space At the end of this section, we introduce a interpolation operator and the projector Π β,0 N . The Legendre-Gauss-Lobatto (LGL) interpolation operator I N : , k, l = 0, 1, · · · , N, where x k and y l are LGL points on the intervalsĪ x andĪ y , respectively. Denote Then the orthogonal projection operator Π β,0 N : From [2], we have the following properties of the projector Π 1,0 N . Lemma 2.11. Let s and r be real numbers satisfying 0 < s < r. Then there exist a projector Π 1,0 N and a positive constant C depending only on r such that, for any function u ∈ H s 0 (Ω) ∩ H r 0 (Ω), the following estimate holds: Similarly, we obtain the following lemma about the projector Π β,0 N . Lemma 2.12. Let β and r be arbitrary real numbers satisfying 0 < β < 1, β < r, β = 1/2. Then there exist a projector Π 1,0 N and a positive constant C independent of N such that, for any function u ∈ H β 0 (Ω) ∩ H r 0 (Ω), the following estimate holds: 3. The fully discrete scheme. In this section, we present the fully discrete scheme for (1) applying the Legendre spectral method. Let Suppose that u(x, y, t) is sufficiently smooth with respect to time. At each time level n, the temporal derivative of (6) is discretized by the L2-1 σ formula (see [1] where c where a (α,σ) 0 For the the coefficients {c (k,α,σ) j } k j=0 , we have the following lemma (see [1]). Let then we have δ α t u k+σ = k j=0 g k+1 j (u j+1 − u j ). Using the same arguments as in proof of Lemma 1 in [1], we obtain the following lemma.
4. Stability and convergence of the fully discrete scheme. In this section, we study the stability and convergence of the scheme (12). For the stability, we have the following theorem. Proof. Taking v N = uk N in (12) gives To estimate the second term on the left-hand side of (14), applying Lemma 2.10 and Lemma 2.5, we firstly deduce that Using Young's inequality and Lemma 2.9, we find Applying Corollary 3.3, (15) and (16), from (14), we obtain Let us rewrite the inequality (17) in the form

HONG LU, JI LI AND MINGJI ZHANG
Noticing that Denote Then, (19) can be rewritten as When k = 0, the estimate (13) holds immediately from the above inequality (20). When k = 1, 2, · · · , M − 1, we can obtain the inequality (13) applying the mathematical induction method. Suppose that we have proven that For l = k + 1, using the inequality (20) and the hypothesis (21), we deduce The proof is completed.
For the error analysis of the fully discrete scheme (12), we have Proof.
From the initial equation (1) and the fully discrete scheme (12), for all v N ∈ V 0 N , we have the following error equation, where r k+σ u,τ = C 0 D α t u(t k+σ ) − δ α t u(t k+σ ), and

SPECTRAL METHODS FOR FRACTIONAL BLOCH-TORREY EQUATIONS 3365
We estimate every term on the right-hand side of (23). For the first term, applying Hölder's inequality, Younger's inequality and Lemma 2.9, we infer where C 1,u is a positive constant depending on u. For the term R k+σ 1 , using Taylor's expansion at t = t k+σ , uk − u k+σ = O(τ 2 ), then we have For the term R k+σ 2 , applying Lemma 2.12, we deduce Using Lemma 2.9, Lemma 2.12 and (7), one has can be estimated as Substituting (24)-(27) into (23), and taking v N =ẽk N , we obtain where C u is a positive constant dependent on u. Let us rewrite (28) in the form Following the same lines as in the proof of Theorem 4.1, we obtain Then applying the triangular inequality e k+1 N ≤ ẽ k+1 N + ê k+1 N , we deduce the desired result.

5.
Alternating direction implicit spectral method. Since the problem is defined in the two-dimensional spatial domain, alternating direction implicit (ADI) method can significantly reduce the computation time and storage requirements. This leads us to propose an ADI spectral scheme. If we add the term B(uk N , v N ) to the left-hand side of (12), we obtain the ADI spectral scheme where v N ∈ V 0 N and k = 0, 1, · · · , M − 1. For the stability of the scheme, we have the following theorem.
Theorem 5.1. The scheme (32) is unconditionally stable in the sense that for all τ > 0, it holds Proof. The proof is almost the same as in the proof of Theorem 4.1, one just notices that Theorem 5.2. Let u be the exact solution, {u k N } M k=0 be the solution of the problem with the initial condition u 0 In comparison with the proof of Theorem 4.2 for (12), one need only notice the additional term, Then following the same lines as in the proof of Theorem 4.2, we get the desired result.
6. Numerical experiment. In this section, we give the detailed implementation of the proposed method and present two numerical examples to show the efficiency of our method. The algorithms are implemented by Matlab.

Numerical examples.
In this subsection, we present numerical examples to verify our theoretical analysis.
To confirm the temporal accuracy, we choose N big enough to eliminate the error caused by spatial discretization. We take N = 125. Tables 1-2 show the errors and temporal convergence rate of u(T )−u k N (T = 1) in L 2 discrete norm for different α and β. From which, we can see the temporal accuracy is second-order, independent of α, which is consistent with our theoretical analysis. Next, we check the spatial accuracy with respect to the polynomial degree N . By fixing the time step small enough to avoid the contamination of the temporal error. We take τ = 1.0 × 10 −3 . Fig. 1 show the errors with respect to polynomial degree N in semi-log scale with different values of α and β for Example 6.1. From which, we can see that the errors decay exponentially. Fig. 2 shows the errors with respect to polynomial degree N in log-log scale for Example 6.2.
7. Summary and discussion. In this work, we have developed a fully discrete scheme for two-dimensional space and time fractional Bloch-Torrey equations. The  scheme employs the spectral approximation using Legendre functions in space and the so-called L2-1 σ formula for the discretization of the time Caputo fractional derivative. The proposed scheme has been proved to be unconditionally stable and convergent with order O(τ 2 + N β−r ). We have presented some numerical experiments to confirm the theoretical analysis. For the problems defined in unbounded spatial domain, one could use the Hermite functions or rational functions as basis functions to approximate the exact solutions. The fully discrete scheme proposed in this paper can be extended to solve the space and time fractional diffusion equations, but the computation work may be huge. In the future, we will try to solve this problem.