Crank-Nicolson ADI Galerkin Finite Element Methods for Two Classes of Riesz Space Fractional Partial Differential Equations

: In this paper, two classes of Riesz space fractional partial differential equations including space-fractional and space-time-fractional ones are considered. These two models can be regarded as the generalization of the classical wave equation in two space dimensions. Combining with the Crank-Nicolson method in temporal direction, efﬁcient alternating direction implicit Galerkin ﬁnite element methods for solving these two fractional models are developed, respectively. The corresponding stability and convergence analysis of the numerical methods are discussed. Numerical results are provided to verify the theoretical analysis.

The fractional operator ∂ 2α 2 ∂|y| 2α 2 is also in the Riesz sense, and can be given similarly. Replacing the second order partial derivative in time in (1) with Caputo derivative operator C D β 0,t with order β ∈ (1, 2), we obtain the following space-time-fractional partial differential equation in two space dimensions, in which the Caputo derivative operator C D β 0,t is defined by ∂s 2 ds.
(1) and (2) reduce to the classical partial differential wave equations in two spatial dimensions. So, these two fractional models (1) and (2) can be regarded as the generalization of the classical wave equations. In recent years, fractional approach has been received many attentions as a powerful modeling methodology, and it is widely applied in materials and mechanics, anomalous diffusion, turbulence, wave propagation, etc. In particular, one of the resulting equations may be the space-fractional model (1) or the space-time-fractional model (2). Recently, the fractional model which involves the Riesz derivative in space or the Caputo derivative in time has attracted many researchers' attention. Li et al. [Li, Yi and Kurths (2018)] introduced the fractional convection operator with Riesz derivative. By employing a continuous time random walks (CTRWs) scheme on one-dimensional lattice and choosing certain kind of two-sided power-law jump length distribution, Li et al. [Li and Yi (2019)] obtained a fractional convention equation containing Riesz derivative. Cai et al. [Cai and Li (2019a)] further studied the properties of Riesz derivative and clarified the relationship between Riesz derivative and fractional Laplacian. In addiction, they also investigated the regularity of the solution to a corresponding fractional differential equation [Cai and Li (2019b)]. Aman et al. [Aman, Al-Mdallal and Khan (2018)] investigated the effect of second order slip on magnetohydrodynamic flow of a fractional Maxwell fluid by considering the Caputo derivative in time. Al-Mdallal et al. [Al-Mdallal, Abro and Khan (2018)] also considered a kind of non-Newtonian models in a porous medium, namely fractional Walter's liquid model by fractional approach in . Bira et al. [Bira, Raja Sekhar and Zeidan (2018)] studied some nonlinear time-fractional evolution equations arising in some important physical phenomena.
It is well known that the analytical solutions of fractional partial differential equations are not always readily obtained, so seeking efficient and reliable numerical treatments is of great importance in real applications. Many numerical methods have been proposed for solving fractional partial differential equations in one and several space variables [Li and Zeng (2015); Li and Cai (2019)]. Moreover, one can refer to the review paper [Li and Chen (2018)] for more details. For the Galerkin finite element methods, recently Li et al. [Li and Wang (2020)] established the discontinuous Galerkin finite element method for a Caputo-type nonlinear conservation law. However, the numerical discretized strategy can not be directly extended to deal with the space-fractional models. About 2006, Ervin et al. developed a rigorous theoretical framework for the wellposedness of a Galerkin weak formulation to fractional elliptic differential equations with constant diffusivity coefficient Roop (2006, 2007)]. Subsequently, many researchers extended the theoretical framework to other numerical methods for numerically solving different space-fractional models, see Zeng et al. [Zeng, Liu, Li et al. (2014)] and the references therein. It seems that the numerical study for fractional models (1) or (2) is scarce, especially using the finite element methods. Besides, solving a space-fractional partial differential equation always leads to expensive computational cost due to the nonlocal properties of Riesz derivative, that is, one needs to calculate the inverse of the dense coefficient matrix with computational cost of N 3 x and storage of N 2 x , which N x is the stepsize number. The computational cost will increase more rapidly if we solve the space-time one, let alone solve the high-dimensional one. Although some algorithms have been proposed to deal with such problem, such as high-order algorithms and fast algorithms (cf. [Wang and Basu (2012)]), more proper and efficient methods with rigorous analysis need to be further developed. It is known that ADI methods have proved to be valuable techniques for solving classical partial differential equations in serval space dimensions. Their effectiveness is to rely on the fact that they reduce the solution of a multidimensional problem to the solution of sets of independent one-dimensional problems. In this paper, we aim to apply the ADI technique with finite element methods to solving the multidimensional fractional problems (1) and (2). Especially, we propose a Crank-Nicolson ADI Galerkin finite element scheme, that is, we use the finite element method in space, and Crank-Nicolson method in time, then combine the alternating direction implicit method to obtain a numerical scheme for (1). Based on the same discretized technique in space, we further apply the modified L1 method for Caputo derivative in (2), then obtain an efficient discretized scheme for (2). It is worthy to mention that the modified L1 method is different from the classical L1 one. For the classical L1 method, the discretization of the Caputo derivative is considered on the grid points of the form {t 0 , t 1 , t 2 · · · }, while the modified L1 method is on {t 0 , t 1/2 , t 3/2 , · · · }. So the two proposed methods naturally lead to the Crank-Nicolson schemes since they are derived on the direct discretization for the first-order/Caputo derivative at the half grid point in time. The rest of this paper is constructed as follows. We first introduce the preliminaries in Section 2, then we derive the Galerkin ADI finite element schemes in Section 3. The corresponding theoretical analyses, including stability analysis and error estimate of the derived schemes, are presented in Section 4. In Section 5, numerical results are demonstrated to verify the effectiveness of the numerical schemes. Some conclusions and further discussions are given in Section 6. As a notation convention, we use the notation C to denote a generic positive constant, which may change at different occurrences, but is always independent of the mesh spacing.

Preliminaries
In this section, we introduce some notations and lemmas that are useful in the following numerical study.
Definition 2.4. For µ > 0, define the semi-norm |u| H µ (Ω) = |ω| µ F(u)(ω) L 2 (Ω) and the norm , and denote H µ (Ω) (or H µ 0 (Ω)) as the closure of C ∞ (Ω) (or C ∞ 0 (Ω)) with respect to · H µ (Ω) . Here, F(u)(ω) is the Fourier transformation of function u(x, y). The above mentioned spaces are equivalent which are shown in the lemma below. Lemma 2.1. For µ > 0, µ = n−1/2, n ∈ N. Then the spaces J µ L,0 (Ω), J µ R,0 (Ω), J µ S,0 (Ω), and H µ 0 (Ω) are equivalent, with equivalent seminorms and norms. Lemma 2.2. Let β ∈ (1, 2), Then ∀u ∈ H β 0 (Ω) and ∀v ∈ H β/2 0 (Ω), we have Similarly, the results for ( RL D β c,y u, v) and ( RL D β y,d u, v) also hold. Lemma 2.3. Let µ > 0 and denoteû as the extension of u by zero outside Ω . Then Similarly, the results for ( RL D µ c,y u, RL D µ y,d u) also hold. Lemma 2.4. Let µ 1 , µ 2 > 0 and denoteû as the extension of u by zero outside Ω . Then 3 The ADI Galerkin finite element schemes 3.1 Space-fractional partial differential equation We first rewrite the initial-boundary problem (1) as a first-order one by setting φ = ∂u/∂t. That is, with L x u = K x ∂ 2α 1 u ∂|x| 2α 1 and L y u = K y ∂ 2α 2 u ∂|y| 2α 2 . For the finite difference of temporal direction, let t n = nτ for n = 0, 1, · · · , n T , where τ = T /n T is the time step size and n T is given positive integer. Denote t n+1/2 = (t n + t n+1 )/2 for n = 0, 1, · · · , n T − 1. Suppose that u(x, y, t) is suitably smooth with respect to the variable t, then we readily obtain and where u n+1/2 = u n+1 +u n 2 and δ t u n+1/2 = u n+1 −u n τ . Now we consider (3) on the time level (x, y, t n+1/2 ), namely, Using the approximations (4) and (5), we have Note that u n+1/2 = τ 2 φ n+1/2 + u n , then (6) can be rewritten as   Adding a small term τ 4 16 L x L y δ t φ n+1/2 = O(τ 2 ) to the left-hand side of the first equation in (7), one has Define the finite element space X r h as the set of piecewise polynomials with degree at most r(r ≥ 1). And denote the interpolation operator by I h in X r h . Here, h denote the maximal length of the sides of each finite element on Ω, so h is a parameter which decreases as the element is made finer. X r h is obvious the finite dimensional subspace of the space V with Suppose that u n h ∈ X r h is the numerical solution of u n and the initial values are approximated by u 0 h denotes an appropriate projection of φ 0 and φ 1 onto X r h . Thus, we obtain the ADI Galerkin finite element scheme as follows: from which we obtain the following approximation to u: Here the bilinear form A is defined by and the bilinear form B in the perturbation term is defined by Eq. (8) present the ADI Galerkin method in inner product form. We can rewrite them as the matrix form. To this end, suppose X r h = X r h,x ⊗ X r h,y , where X r h,x and X r h,y are finitedimensional subspaces of H α1 0 (I x ) and H α2 0 (I y ), respectively. And their corresponding bases are denoted as {ϕ i (x)} Nx i=1 and {ϑ j (y)} Ny j=1 . Here, N x , N y are the positive integers. The uniform partitions of I x and I y are then given by and denote c,y ϑ k (y), RL D α2 y,d ϑ m (y)).
Here, Φ n+1 , RHS n , U n , F n+ 1 2 ∈ R Nx×Ny . Therefore, denoting the intermediate variable Φ * n+1 ∈ R Nx×Ny , the alternating direction Galerkin finite element scheme can be executed as follows. First we solve from which we obtain the approximation to U n+1 : 3.2 Space-time-fractional partial differential equation We firstly present the modified L1 scheme of the Caputo derivative C D β 0,t with β ∈ (1, 2) below, where the truncation error The coefficients b n and B n are defined by One can readily observe that the coefficients (14) have the following properties, Now considering the Eq. (2) on the time level (x, y, t n+1/2 ), and using the approximations (12) and (4), we get h − u n h , then multiplying by µτ on both sides of (15), one has and Adding two small terms (τ µ) 2 4 L x L yũ n+1 and (τ µ) 2 4B 2 0 L x L yũ 1 to the left-hand side of (16) and (17), respectively, we obtain the following ADI Galerkin scheme for (2) in inner product form: Findũ n+1 h ∈ X r h for n = 0, 1, 2, · · · , n T − 1 such that with initial values u 0 h = Π 1,0 h φ 0 and φ 0 h = Π 1,0 h φ 1 . From (18) we can obtain the following approximation to u: Here the bilinear forms A and B are defined by (10) and (11), respectively. The matrix form of (18) can be formulated similarly to (8), we omit the details here.

Stability and convergence
In this section, we first study the stability and convergence of the fully discrete scheme (8).
The seminorm | · | α and norm · α are equivalent if u ∈ H α1 0 ∩ H α2 0 , see [Zeng, Liu, Li et al. (2014)]. This result is demonstrated by the following lemma. Lemma 4.1. If u ∈ H α1 0 ∩ H α2 0 , then the seminorm | · | α and norm · α are equivalent, that is, there exists positive constants C 1 < 1 and C 2 independent of u, such that Define the operator Π h : V → X r h . Suppose that u ∈ H l (Ω), 0 < l < r + 1, 0 ≤ s ≤ l, then the following approximation property holds, where the constant C is only dependent on Ω.

Proof. By Cauchy-Schwarz inequality and Lemma 2.3 for the definition of A, we obtain
That is, Taking u h as Π h u, and utilizing Lemma 4.1 and the estimate (20), we obtain the desired result.
Lemma 4.3. If w ∈ H m ∩ H α1 0 ∩ H α2 0 and 1 ≤ m ≤ r + 1, then we have Proof. One can readily derive the above results from the symmetry of these bilinear forms A and B defined by (10) and (11). Thus we complete the proof. Combining Lemma 2.4 with the above proof for equality (22), we readily obtain the following result. Lemma 4.4. If w ∈ H m ∩ H α1 0 ∩ H α2 0 and 1 ≤ m ≤ r + 1, then we have B(w, w) ≥ 0.
We shall need the following version of Gronwall's inequality (cf. [Zeng, Cao and Li (2013)]). Lemma 4.5. Assume that {κ n } and {p n } are nonnegative sequences, and the sequence {ρ n } satisfies where g 0 ≥ 0. Then the sequence {ρ n } satisfies Now, we are in the position to present the stability for the fully discrete scheme (8).
Theorem 4.1. The fully discrete scheme (8) associated with (9) is unconditionally stability in the sense that where C is a positive constant independent of mesh spacing, B is a bilinear form defined by (11), and ρ n+1/2 is the perturbation term for (9), that is δ t u n+ 1 2 = φ n+ 1 2 + ρ n+ 1 2 .
Proof. Let v = φ n+ 1 2 h or, equivalently, v = δ t u n+1/2 h − ρ n+1/2 in the equivalent form (19) of (8), then we obtain Note that So, it follows from Lemma 4.3 that (24) can be rewritten as Using Cauchy-Schwarz inequality and the equivalence of | · | α and · α , one has Similarly, we have For sufficiently small τ , we can derive from (25) that Then summing up n from 0 to m − 1 for (26) gives From Lemma 4.4, the third term on the left-hand side of (27) is nonnegative. Removing the nonnegative term, we get Comparing each term in the above inequality with that in (23), then applying Lemma 4.5, we can conclude that This completes the proof. Next, we study the convergence analysis for (8).
Theorem 4.2. Suppose that r ≥ 1, u and u n+1 h (0 ≤ n ≤ n T − 1) are the solutions of (1) and the fully discrete scheme (8) ) , m ≥ r + 1, then there exists a positive constant C independent of τ and h such that Proof. It is convenient to use the equivalent forms of (8) and (9) below.
By Theorem 4.1, we can derive that Next, we estimate the right-hand side of the above inequality. For the initial errors e 0 andê 0 , we have The third term on the right-hand side of (31) has the following estimate.

It follows that
All this completes the proof. In the following, we demonstrate the stability for the ADI Galerkin scheme (18). Theorem 4.3. Suppose that u n h is the solution of the ADI Galerkin scheme (18) and f ∈ C([0, T ]; L 2 (Ω)). Then it holds that where C is a positive constant and independent of τ and h.
Proof. Notice that the properties of the coefficients B n and b n are nice, one can readily derive the above stability result (32) by using the idea in the proof of Theorem 4.1. Thus we omit the details here. Finally, we present the error estimate for the ADI Galerkin scheme (18). Theorem 4.4. Suppose that r ≥ 1, u and u n+1 h (0 ≤ n ≤ n T − 1) are the solutions of (2) and the fully discrete scheme (18) ∂|y| α 2 ∂ tt u ∈ L ∞ (0, T ; H m (Ω)) , m ≥ r + 1, then there exists a positive constant C independent of τ and h such that Proof. We first consider the case n ≥ 1 in (18). We use its equivalent form which is shown below, The weak form of the Eq. (2) on the time level t n+ 1 which can be written as Here R n+ 1 2 is the form as denoted in (13), with ∂ 3 u/∂t 3 instead of g (3) (t). Multiplying by µ and adding the small term on both side of (35), we obtain Let u * = Π α,0 h u, e = u * − u h , η = u − u * ,ê = ∂u * ∂t − φ h ,η = φ − ∂u * ∂t , and v = χ in (36).
For the last term on the right-hand side of equation (37), we have the following estimate.
So using (38)-(41) gives (33). All this ends the proof. Remark 4.1. It seems that the convergence order in time of the ADI Galerkin scheme (18) will be destroyed when β < 5 3 . However we observe that the convergence order in time always be 3 − β in the numerical tests. The disagreement probably means that the theoretical analysis in (33) may be improved. One possible way to remove the term O(τ 1+β 2 ) is by introducing a much smaller disturbance term to construct the ADI Galerkin scheme, however it seems that the resulting scheme does not improve the numerical results in practical computation but only cause more complicated computations.

Numerical results
In this section, we present numerical results to verify the theoretical analysis. For sake of computational convenience, in our numerical tests, we fix the coefficients K x = K y = 1/2, and let Ω = (0, 1) × (0, 1) and T = 1. All the tests are done by using the linear element. Example 5.1. Consider the following initial conditions in the two-dimensional fractional model (1): The source term f is given by The corresponding exact solution is From the ADI Galerkin scheme (8), we obtain the numerical results of Tabs. 1-2 for Eq.
(1). All the numerical errors are computed at t = 1 for different values of α 1 and α 2 . From Tabs. 1-2, we can see that convergent orders in space and in time are (2 − α max )-order and second-order accuracy in seminorm | · | α respectively, which is in line with the theoretical analysis.  u(x, y, 0) = 0, ∂u ∂t t=0 = 0.
The source term f is given by Then we obtain the exact solution below Using the ADI Galerkin scheme (18), we obtain the numerical results of Tabs. 3-4. Here we let τ = 1/128 when test the convergence order in space (see Tab. 3), and let h = τ (3−β)/(2−αmax) when test the convergence order in time (see Tab. 4). The numerical errors are also obtained at t = 1. From Tabs. 3-4, we observe that the convergence orders are 2 − α max and 3 − β in space and in time, respectively, which show better results than that in Theorem 4.4.  (2) in which the exact solution is unknown: f = 1. The corresponding initial conditions are provided by u(x, y, 0) = 0 and ∂u ∂t | t=0 = 0. Since the exact solution is unknown, we determine the convergence order by the ratios of differences between different numerical solutions which are computed by the consecutive halved stepsizes. That is, we compute the spatial convergence order by adopting log 2 (e(h)/e(h/2)) with the error e(h) = |u h − u h/2 | α . By applying the numerical scheme (18), we obtain the numerical results which are shown on Tab. 5. From the numerical results, it can be seen that the | · | α errors will gradually decrease as step sizes shrink at the rate of 1/2. However, the convergence orders are both less than the theoretical ones. Similar numerical results are also observed for the temporal direction which are not presented here. One possible reason is that the right-hand side function f is constant, which leads to the issue that the regularity of analytical solution in (2) does not meet the requirements in Theorem 4.4. The corresponding theoretical analyses of the less regularity issue are not available in this paper and need to be further studied in the future. However from the above numerical study in Example 5.3, one still can get a glimpse of the behavior of the numerical solution numerically with both the orders α 1 and α 2 tend to one. For simplicity of presentation, we focus on the case β = 2 here, that is the fractional model (1), and compute the numerical solution with various α 1 and α 2 by the numerical scheme (8). The numerical results are presented in Fig. 1. From Fig. 1, we can observe that the numerical solution is more closer to that of the integer case as the orders α 1 and α 2 tend to one. We also observe similar phenomenon for the case β ∈ (1, 2) in which numerical results are not demonstrated here. From this perspective, we may conclude that these two fractional models considered in this paper can be regarded as the generalization of the corresponding classical integer-order model. In this work, we studied two classes of Riesz space-fractional partial differential equations including space-fractional and space-time fractional derivatives. Our main contributions in this work are to develop the efficient ADI Galerkin finite element schemes for these two considered models and establish their stability and convergence. Numerical tests are presented to verify the convergence theories. For sake of convenience, the current study relies crucially on the high regularity of the solution. It would be interesting to extend the argument to the less regularity issue, for example the weak singularity of the solution at t = 0. However it seems that the extension is not an easy job and requires further study, since most of the existed results only deal with the time-fractional problems, see Li et al. [Li and Chen (2018)] and the references therein. Conflicts of Interest: The author declares that he has no conflicts of interest to report regarding the present study.