On high-order compact schemes for the multidimensional time-fractional Schrödinger equation

In this article, some high-order compact finite difference schemes are presented and analyzed to numerically solve one- and two-dimensional time fractional Schrödinger equations. The time Caputo fractional derivative is evaluated by the L1 and L1-2 approximation. The space discretization is based on the fourth-order compact finite difference method. For the one-dimensional problem, the rates of the presented schemes are of order O(τ2−α+h4)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$O(\tau ^{2-\alpha }+h^{4})$\end{document} and O(τ3−α+h4)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$O(\tau ^{3-\alpha }+h^{4})$\end{document}, respectively, with the temporal step size τ and the spatial step size h, and α∈(0,1)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\alpha \in (0,1)$\end{document}. For the two-dimensional problem, the high-order compact alternating direction implicit method is used. Moreover, unconditional stability of the proposed schemes is discussed by using the Fourier analysis method. Numerical tests are performed to support the theoretical results, and these show the accuracy and efficiency of the proposed schemes.


Introduction
Nowadays fractional differential equations have been widely studied in many fields, owing to their diverse applications in physics, biology, chemistry, mechanics, and finance theory [1][2][3][4][5][6][7][8][9]. These applications have contributed to the emergence of various fractional differential equations in the mathematical and physical world. In [10][11][12][13][14], the authors have published some new results of fractional operators and their applications. In fact, it is difficult to gain analytic solutions of fractional differential equations. Therefore, it is important to obtain highly accurate numerical methods for solving these fractional differential equations.
For time-fractional partial differential equations, many different numerical methods have been introduced, including theoretical analysis and numerical computing. Authors in [15,16] have proposed various spatial second-order finite difference methods for a onedimensional problem. For two-dimensional time-fractional partial differential equations, several numerical methods have been proposed [17,18]. To improve the numerical accuracy, some fourth-order compact finite difference methods have been proposed for oneand two-dimensional time-fractional partial differential equations. Cui [19] constructed high-order compact alternating direction implicit schemes for two-dimensional timefractional diffusion equations. Gao and Sun [20] focused on the study of spatial sixth-order accurate combined compact alternating direction implicit difference schemes for solving two-dimensional time-fractional advection-diffusion equations. Zhai and Feng [21] investigated several compact alternating direction implicit methods for a two-dimensional time-fractional diffusion equation.
In this paper, we consider the following multidimensional time-fractional Schrödinger equation (TFSE): with the initial and boundary conditions where is a rectangular domain in R d (d = 1, 2), ∂ is the boundary of , i = √ -1, and u 0 and f are known smooth functions. The fractional derivative ∂ α u(x, t)/∂t α is the αth order Caputo time fractional derivative defined by (1.4) with α ∈ (0, 1). In recent years, some researchers have presented numerical solutions for the timefractional Schrödinger equation. Mohebbi et al. [22] investigated meshless method based on collocation method for the numerical solutions of time fractional nonlinear Schrödinger equation. Wei et al. [23] proposed an implicit fully discrete local discontinuous Galerkin method for the TFSE. Li et al. [24] proposed L1-Galerkin finite element method for the numerical and stability analysis of multidimensional TFSEs. In [25], the space-time Jacobi spectral collocation method is used to solve the time-fractional nonlinear Schrödinger equations subject to the appropriate initial and boundary conditions. Chen et al. [26] also proposed linearized compact alternating direction implicit schemes for nonlinear TFSEs.
The main purpose of this paper is to construct some efficient high-order compact difference schemes for solving TFSE (1.1) with (1.2) and (1.3). We apply the L1 and L1-2 formulas to approximate the time-fractional derivative and use compact operators to approximate spatial second-order derivatives. We treat the two-dimensional problem using the compact alternating direction implicit (ADI) scheme. The computational complexity is reduced to some extent. The L1 formula is the main approximation formula for approximating the time-fractional derivative, in which the truncation error is O(τ 2-α ) [27,28]. In [29], Gao et al. proposed a new difference analog of the Caputo fractional derivative with convergence order O(τ 3-α ), called the L1-2 formula.
The remainder of our paper is organized as follows. In Sect. 2, we introduce some notations and useful results, and then two high-order compact finite difference schemes are constructed for the one-dimensional TFSE. We use the Fourier analysis method to investigate the stability of the two schemes. In Sect. 3, we extend our methods to the twodimensional case and propose two compact ADI difference schemes by adding the perturbation term. Furthermore, a stability analysis is presented. In Sect. 4, we present numerical examples and detailed numerical results to confirm our theoretical analysis. Finally, a brief conclusion is provided in Sect. 5.

High-order compact difference schemes for the one-dimensional TFSE
In this section, we consider the following one-dimensional (1D) TFSE: where L and T are positive constants. u 0 (x), ϕ(x, t), and f (x, t) are given smooth functions. We first present some notations and useful results, which provide the basis for the theoretical analysis of our numerical methods. Let N x and N be two positive integers, h = L N x be the spatial step size, and τ = T N be the time step size. (2.4) Next we introduce some useful results and design several compact finite difference schemes.
We now investigate the stability of the above schemes by using the Fourier analysis method. For the simplicity of illustration, we consider the case that f (x, t) = 0.
Let the numerical solution be represented by where i = √ (-1), v n is the amplitude at time t n , and θ is the wave number in the x direction. We define the discrete L 2 norm by u n 2 L 2 = h N x j=1 |u n j | 2 . For the fully discrete schemes CS1DI and CS1DII, we have the following stability results. Theorem 1 Schemes CS1DI and CS1DII, defined by (2.9) and (2.10) respectively, are unconditionally stable for α ∈ (0, 1).
Proof We will complete the proof using mathematical induction. For convenience, we only give the proof for scheme CS1DI, and the case for scheme CS1DII can be completed using a similar idea.

High-order compact ADI schemes for the two-dimensional TFSE
In this section, we develop two high-order compact ADI schemes for two-dimensional (2D) TFSE with the following initial and boundary conditions: where L and T are positive constants, = (0, L) × (0, L). u 0 (x, y), ϕ(x, y, t), and f (x, y, t) are given smooth functions.
Let M, N be positive integers, h = L M be the spatial step size, and τ = T N be the time step The notations δ 2 y u n jk and L y u n jk can be defined similarly. Applying the fourth-order Padé scheme for the second-order derivatives, we have Using (3.5) and (3.6), we can obtain the following semi-discrete fourth-order approximation for the 2D problem (3.1): Thus, by using Lemma 2.1 for (3.7) and neglecting the small truncation errors, we obtain the following scheme: (a n-l-1a n-l )u l jka n-1 u 0 Act with the operator (-iτ α (2α))L x L y on both sides of (3.8). By noting that a 0 = 1, after rearranging the terms we have the equivalent forms (a n-l-1a n-l )u l jk + a n-1 u 0 jkiτ α (2α)L x L y f n jk . (3.9) Next we consider the derivation of the ADI difference scheme so that some small splitting terms should be added onto the left-hand side of (3.9) in order to achieve the operator splitting scheme. By adding the splitting term to the left-hand side of (3.9) and rearranging the terms of the resulting scheme, we obtain an approximate scheme for n > 1 as follows: (a n-l-1a n-l )u l jk + a n-1 u 0 Introducing the intermediate variables u 1( * ) and u n( * ) , we can obtain the following highorder compact ADI scheme (CS2DI) for (3.1)-(3.3): CS2DI : = L x L y [ n-1 l=1 (a n-l-1a n-l )u l jk + a n- (3.11) To solve (3.11), we need to give the values of u 1( * ) and u n( * ) on the boundary, and these are obtained from the second and fourth equation in (3.11), respectively. Now, we present the second difference scheme for (3.1)-(3.3). Using Lemma 2.2 for (3.7) and omitting the small terms, the second difference scheme can be given in the form of By acting on both sides of (3.12) with the operator ( -iτ α (2-α) c 0 )L x L y and rearranging the terms, we obtain By adding the splitting term ( iτ α (2-α) c 0 ) 2 δ 2 x δ 2 y (u n jku n-1 jk ) to the left-hand side of (3.13) and rearranging the terms of the resulting scheme, we obtain the following scheme for n > 1: (3.14) Upon introducing the intermediate variables u 1( * ) and u n( * ) and noticing that c 0 = 1 for n = 1, we obtain the following additional high-order compact ADI scheme (CS2DII) for (3.1)-(3.3): CS2DII : To solve (3.15), we need to give the values of u 1( * ) and u n( * ) on the boundary, and these are obtained from the second and fourth equations in (3.15), respectively. From (3.7), it is obvious that schemes CS2DI and CS2DII can maintain fourth-order accuracy in space, and from Lemmas 2.1 and 2.2 we know that schemes (3.9) and (3.13) have temporal convergence rates of O(τ 2-α ) and O(τ 3-α ), respectively. We note the fact that if the splitting errors are much larger than the truncation error, as pointed out by Douglas and Kim [32], a splitting term will play an important role in the accuracy of the solution. This phenomenon has also been pointed out by Zhai et al. [33]. From the splitting term [iτ α (2α)] 2 δ 2 x δ 2 y (u n jku n-1 jk ), we find that the splitting errors are O(τ 2α+1 ) in schemes CS2DI and CS2DII, respectively. This means that the temporal convergence rates for these schemes will tend to α + 1. The numerical results in Sect. 4 further confirm this judgment.
In the following, we investigate the stability of the high-order compact ADI schemes (CS2DI and CS2DII) using the Fourier analysis method. Here, in order to simplify the notations and without loss of generality, we consider the case that f (x, y, t) = 0.
Suppose that the numerical solution is represented by u n jk = v n e i(θ x jh+θ y kh) , (3.16) where i = √ (-1), v n is the amplitude at time t n , and θ x , θ y are the wave numbers in the x and y directions, respectively. We define the discrete L 2 norm by u n 2 For the fully discrete schemes CS2DI and CS2DII, we have the following stability results. CS2DI and CS2DII, defined by (3.11) and (3.15) respectively, are unconditionally stable for α ∈ (0, 1).

Theorem 2 Schemes
Proof We apply mathematical induction to complete the proof. For convenience, we only give the proof for scheme CS2DI, and the case for scheme CS2DII can be completed using a similar idea.

Numerical results
In this section, three numerical examples, which verify the efficiency and accuracy of the proposed schemes, are presented.
Example 1 Let us consider the following 1D time-fractional Schrödinger equation: with the exact solution u(x, t) = (1 + i)t 2 sin πx. The initial and boundary conditions can be obtained from the exact solution. Now, we investigate the spatial convergence rates for schemes CS1DI and CS1DII. We choose N = 1000 to avoid the influence of the temporal approximation. Table 1 gives the L ∞ -norm errors (denoted by L ∞ -error) and convergence orders at time T = 1 for α = 0.25 (or α = 0.5) and various values of N x . From Table 1, we conclude that schemes CS1DI and CS1DII are verified to have fourth-order accuracy in space. Regarding the time accuracy, by fixing N x to eliminate the contamination of the spatial error, the numerical results at T = 1 for different α and N values are presented in Table 2. In Fig. 1, we also present the errors in the L ∞ -norm and L 2 -norm as a function of the time step sizes for α = 0.25, where N x = 50 in scheme CS1DI and N x = 250 in scheme CS1DII. From these results, we find that the convergence order in time for schemes CS1DI and CS1DII is close to 2α and 3α, respectively.   Example 2 Here, we consider the alternative 1D time-fractional Schrödinger equation The exact solution to the problem is u(x, t) = t 2 (cos x + i sin x). The initial and boundary conditions can be obtained from the exact solution.
In this example, we first use schemes CS1DI and CS1DII to verify the spatial convergence rates. Table 3 presents the L ∞ -errors and convergence orders with N = 1000 and various values of N x and α at time T = 1, from which we find that schemes CS1DI and CS1DII both achieve fourth-order accuracy regardless of the value of α. Next, we test the temporal convergence rate for the case that T = 1 and N x = 2000. The results for different N and α values are listed in Table 4. In Fig. 2, we plot the errors in the L ∞ -norm and L 2norm as a function of the time step sizes for N x = 50 in scheme CS1DI and N x = 250 in scheme CS1DII, where α = 0.75. It shows that the slopes of the error curves obtained for   scheme CS1DI is 1.25, which accords with the theoretical estimates 2-α. However, scheme CS1DII yields a temporal approximation order 3-α, i.e., the slopes of the error curves is 2.25. Again, as expected it is clearly shown that scheme CS1DII has higher accuracy than scheme CS1DI. Example 3 We consider the following 2D time-fractional Schrödinger equation: Its exact solution is u(x, y, t) = (1+i)(1-x)(1-y)xyt 2 . The initial and boundary conditions of this problem can be extracted from the exact solution.
First, the spatial accuracy is numerically examined. We decrease the mesh size of h to h/2 and τ to τ /2 4/1+α , and the errors and orders for schemes CS2DI and CS2DII are shown in Table 5 for different α values. We can see that the experimental convergence order is approximately four. Second, the numerical accuracy in time is verified. The numerical results at T = 1 for α = 0.1, 0.5, and 0.9 are listed in Table 6. Meanwhile, the numerical results are plotted in Fig. 3. From these results, we find that schemes CS2DI and CS2DII yield a temporal approximation order close to 1 + α, i.e., the slopes of the error curves are 1.25 and 1.75, respectively, for α = 0.25 and α = 0.75. We also plot the contour graphs of the numerical errors using schemes CS2DI and CS2DII with N = M = 100 for α = 0.1, 0.5, 0.9   in Fig. 4. From these results, we find that two methods have the same temporal convergence rates, and the convergence rate of scheme CS2DII is lower than 3α because of the splitting term.

Conclusions
In this paper, we have proposed two kinds of high-order compact finite difference schemes with order O(τ 2-α + h 4 ) and O(τ 3-α + h 4 ), respectively, for solving the one-dimensional time-fractional Schrödinger equation. Here, the time discretization is replaced by the L1 and L1-2 formulas, and the space discretization is derived using the high-order compact finite difference method. Then, the extension to the two-dimensional problem is considered. We have designed two high-order compact ADI difference schemes with accuracy O(τ 1+α + h 4 ). In the two-dimensional case, it is worth noting that the splitting term can affect the local truncation error for the time accuracy. Moreover, both theoretical analysis and numerical examples show that all schemes are unconditionally stable for α ∈ (0, 1) and have the high accuracy. Evidently, the proposed schemes are easy to be implemented and extended to solve other fractional PDEs.