Fast High-Order Difference Scheme for the Modified Anomalous Subdiffusion Equation Based on Fast Discrete Sine Transform

The modified anomalous subdiffusion equation plays an important role in the modeling of the processes that become less anomalous as time evolves. In this paper, we consider the efficient difference scheme for solving such time-fractional equation in two space dimensions. By using the modified L1 method and the compact difference operator with fast discrete sine transform technique, we develop a fast Crank-Nicolson compact difference scheme which is proved to be stable with the accuracy of Oðτmin ð1+α,1+βÞ + h4Þ. Here, α and β are the fractional orders which both range from 0 to 1, and τ and h are, respectively, the temporal and spatial stepsizes. We also consider the method of adding correction terms to efficiently deal with the nonsmooth problems. Numerical examples are provided to verify the effectiveness of the proposed scheme.


Introduction
In this paper, we focus on the numerical method for the timefractional modified subdiffusion equation [1]: with the initial condition uðx, 0Þ = vðxÞ and the homogeneous Dirichlet boundary condition. Here, x = ðx 1 , x 2 Þ, Ω is the rectangle domain, T > 0, and Δ is the Laplacian defined by Δ = ∂ 2 /∂x 2 1 + ∂ 2 /∂x 2 2 . The parameters κ 1 and κ 2 are some fixed positive constants, f and v are two given functions, 0 < α, β < 1, and RL D γ 0,t is the Riemann-Liouville derivative of order γ given by: where Γð·Þ is the Gamma function.
Anomalous diffusion is ubiquitous in nature and it can be characterized by the method of mean square particle displacement at the microscopic level. When the mean square displacement (MSD) is linear with time, the particle is precisely in classical Brownian motion. If the MSD grows either sublinearly or superlinearly with time, then this phenomenon is regarded as the subdiffusion and superdiffusion, respectively. Numerous experimental studies have shown that the anomalous diffusion may adequately describe a number of physical processes in recent decades [2,3]. Equation (1) is an important class of anomalous diffusion models in which the physical processes are observed to become less anomalous as time evolves [4].
In [4], the author presented the solution of a onedimensional modified anomalous subdiffusion equation on an infinite domain. The analytical solution the author obtained contains the infinite series of Fox special functions, which is of complex form that makes it difficult to apply to practical numerical simulations. So, one needs to resort to the numerical methods for efficiently solving equation (1). Many efficient numerical methods for solving fractional models have emerged in recent years, see the book [5] and the two review papers [2,6]. For equation (1), some numerical schemes have been developed. Ding and Li applied two kinds of high-order compact finite difference methods to construct efficient numerical schemes. The stability and convergence analysis are proved by the Fourier method [7]. In [8], the authors developed the compact difference scheme based on the second-order compact approximation formula of the first-order derivative. The two papers mentioned above both focus on the one-dimensional case.
For the two-dimensional case, Chen and Li employed the modified L1 method and compact difference method and proposed a compact alternating direction implicit scheme. By utilizing the energy method, they proved that their scheme is stable with an accuracy of Oðτ 2 min ðα,βÞ + h 4 Þ in the new defined norm which is equivalent with H 1 -norm, under the assumption that the solution is sufficiently smooth [1]. Such assumption may be too restrictive to limit the scope of application of their scheme. To address this issue, Chen proposed two robust fully discrete finite element methods by convolution quadrature in time. He proved that the schemes are convergent under data regularity without relying on the assumption of the solution regularity. In addition, he also proposed a Crank-Nicolson finite element scheme to numerically compare and verify the robustness of the convolution-based schemes in solving nonsmooth solution problems, but no detailed theoretical analysis of the scheme was given [9]. It seems that the numerical methods for equation (1) have not been sufficiently studied. This motivates us to design efficient numerical schemes for (1), especially for high-dimensional problems where the solutions are not sufficiently smooth.
As the further work on the high-dimensional equation (1), we focus on designing numerical schemes that are computationally efficient and can handle the nonsmooth solution case. In [10], Li et al. implemented the fourth-order compact difference operator by a fast discrete sine transform (DST) via the FFT algorithm, which greatly reduces the computational cost and storage requirement. Notice that the DST technology can avoid solving matrix inversion directly and has been successfully applied in the discretization of classical models, such as Poisson equation [11] and general order semilinear evolution equations [12], just to name a few. On the other hand, the weak singularity of the fractional model has gradually attracted the attention of scholars in the fractional community, and some kinds of methods have been proposed to resolve this issue, such as nonuniform meshes [13][14][15][16] and convolution quadrature [9,17,18]. The method of adding correction terms is also an efficient way of dealing with nonsmooth solutions problems. However, this method is generally not very stable as the starting weights need to be obtained through a linear system which involves the illconditioned exponential Vandermonde matrix. To resolve this issue, Zeng et al. theoretically and numerically shown that the accuracy of numerical solution can be efficiently improved with only a few correction terms [19]. Since then, a variety of numerical schemes based on the addition of correction terms have emerged for fractional problems with nonsmooth solutions, see [3,20,21]. To the best of our knowledge, it seems that the method of adding correction terms with DST for solving equation (1) has not been considered in the existing literatures yet.
The contributions of this paper are as follows. First, we apply the modified L1 method to discrete the Riemann-Liouville derivative and compact difference operator with DST to discrete the Laplacian and then naturally obtain the fast Crank-Nicolson compact difference scheme for the two-dimensional problem (1), see (7). Second, the stability and error estimate are rigorously proved by the energy method, see Theorems 2 and 3. Specially, we improve the convergence results in [1] and guarantee computational efficiency but without sacrificing the accuracy of the scheme. Note that the small term added during the construction of the ADI scheme in [1] destroys the accuracy of their original scheme. Third, by using the method of adding correction terms, we successfully improve the accuracy of the proposed scheme in solving the nonsmooth problem with no impact on the stability of the original scheme, see (9). Finally, we provide concrete numerical tests to show the effectiveness of the scheme in solving the high-dimensional problem with nonsmooth solution, see Table 1 and Figures 1-3. The rest of the paper is organized as follows. In Section 2, we derive the fast Crank-Nicolson compact difference scheme by the modified L1 method and the compact difference operator with DST. In Section 3, we prove that the proposed numerical scheme is stable with an accuracy of Oðτ min ð1+α,1+βÞ + h 4 Þ under the assumption that the solution is sufficiently smooth. To weaken the assumption and make the scheme more robust in solving nonsmooth solution problems, we present the improved version in Section 4 with the method of adding correction terms. Numerical examples are given in Section 5 to confirm the effectiveness of the proposed scheme. Finally, we present the conclusions of this paper in Section 6.
Throughout this paper, we shall let the symbol c (with or without subscript) be a positive constant which may vary at different locations but is always independent of the temporal and spatial stepsizes.
Next, we consider the spatial discretization for (4). In order to make our discussion more general, we follow the notations presented in [10] and always set the symbol d = 2 unless otherwise noted. Denote the domain The fully discrete grids in space are denoted by where the local truncation error is given by R n+1/2 Omitting the small term R n+1/2 xt , we obtain the following Crank-Nicolson compact difference scheme for (1): find u n h of uðx h , t n Þ for n ≥ 1, such that  Journal of Function Spaces If we solve the discretized system (6) directly, the computational cost will be OððM 1 M 2 ⋯ M d Þ 2 Þ on each time level due to the calculation of matrix inversion. Next, we employ the fast discrete sine transform based on FFT to reduce the computational cost to OððM 1 M 2 ⋯ M d Þ log ðM 1 M 2 ⋯ M d Þ Þ [11], which greatly improves the computational performance. Since the discrete sine transform of the grid function v h at the kth position is provided by where s j k = cos ðj k π/M k Þ and 1 ≤ j k ≤ M k − 1. One can refer to [11] or [10] for the derivation. Denote the index set ν The computational procedure is described as follows: (a) For n ≥ 0, we first computedû n ν and f ∧ n+1/2 from u n h and f n+1/2 by means of DST (b) And then we solve equation (10) from which the numerical solution u n h is obtained fromû n ν by the inverse of DST.

Stability and Error Estimates
In this part, we demonstrate the stability and error estimates for the compact difference scheme (6).
We first introduce some useful notations. For any grid function v ∈ V h , the discrete L 2 -norm is given by ∥v∥ = ffiffiffiffiffiffiffiffiffiffiffiffi ffi ðv, vÞ h p with the discrete inner product ðu, . Here, ∇ h = ðδ 1 , δ 2 ,⋯,δ d Þ. One can readily have the equivalence of jvj 1 and ∥v∥ 1 for any v ∈ V h in view of the embedding theorem.
We shall first need the following lemma.
where v n ∈ V h , n ≥ 0.
Proof. The proof of the lemma can be obtained in view of Lemma 4.2 in [22] or Lemma 4.4 in [1], thus, the details are omitted here.
We are ready to present the stability of the scheme (6).

Theorem 2.
The Crank-Nicolson compact difference scheme (6) is stable in the sense that Proof. By taking the discrete inner product on both sides of (6) with 2τu n+1/2 h , we get Notice that the difference operator Δ h is bounded in discrete inner product ([10], Theorem 2): where write the above inequality as: We sum up n from 1 to m and replace m with n to get By the Cauchy-Schwarz inequality and the inequality kv h k ≤ ckv h k 1 with the equivalence of kv h k 1 and k∇ h v h k, we obtain the estimate: from which we derive that Next, we consider the case n = 0 for the scheme (6) for the estimate of G 1 . By a similar procedure, we take the discrete inner product for (6) with 2τu 1/2 h when n = 0, we have from which we have where . Utilizing the Cauchy-Schwarz inequality again, we arrive at the estimate for the last two terms on the right-hand of the above inequality: By letting the constants ε 1 = B 0 /ð4A 0 Þ and ε 2 = B 0 /4 and the equivalence of the ∥v h ∥ 1 and jv h j 1 , we further get which implies that the G 1 has the following estimate: Therefore, the inequality (8) yields By the mean value theorem, one can readily check that the coefficients appearing in the above inequality are all bounded, that is, we formally have Thus, the proof is completed.
By means of the error equation and the stability conclusion, we have the following convergence result.
Theorem 3. Suppose that u ∈ C 2 ð0, T ; C 6 ðΩÞÞ, then we have the discrete L 2 -norm error estimate: For n ≥ 1, Proof. The error equation can be obtained by subtracting (6) from (5), that is, by letting the error e n h = uðx h , t n Þ − u n h for It follows from Theorem 2 that which leads to the desired convergence result.

Numerical Implementation for Nonsmooth Problems
In general, the solution of equation (1) may not have the regularity required in Theorem 3. If the nonsmooth solution problems are directly solved by the fast Crank-Nicolson compact difference scheme (7), unsatisfactory accuracy may be obtained. In this part, we apply the method of adding suitable correction terms when dealing with such nonsmooth issue. Following the idea presented in [3], we, respectively, take the numerical approximations of RL D γ 0,t gðtÞ and the firstorder time derivative dgðtÞ/dt at t = t n+1/2 as follows: 6 Journal of Function Spaces Here, w ðγÞ n,k and w ð1Þ n,k are the starting weights which are chosen such that the above schemes are exact for some power functions gðtÞ = t ζ j with 0 < ζ j < ζ j+1 and 0 ≤ j ≤ m, that is, they can be determined by the two linear systems: respectively. So, we have the following fast Crank-Nicolson compact difference scheme with correction terms: for n ≥ 0, The execution procedure of the above scheme is similar to that of (7). We can observe that the scheme (9) is stable and effective in solving nonsmooth problems, which will be verified by numerical examples in the next section. We remark that the method of adding correction terms is based on the assumption that the problem solution can be divided into two terms: low regularity and high regularity terms (with respect to time). Such assumption is valid for equation (1) in view of the solution formulation discussed in [4]. By using the starting weights in the correction terms, one can improve the accuracy of the proposed scheme for dealing with the nonsmooth solution problem. For further details about the parameters m and ζ j , one may refer to [19].

Numerical Examples
In this part, we present two numerical examples to verify the accuracy and effectiveness of the scheme (9). The L 2 -norm error at t = t n is obtained by eðn, hÞ = kuðx h , t n Þ − u n h k, and the convergence orders in time and in space are calculated by log ðeðn, hÞ/eð2n, hÞÞ and log ðeðn, hÞ/eðn, h/2ÞÞ, respectively. For simplicity, we set the parameters κ 1 and κ 2 in (1) to be one and restrict the computational domain to be Ω = ð0, 1Þ 2 . We remark that the numerical tests in this paper are implemented by MATLAB software (R2020a) on an Apple OS platform with a quad-core 2.3 GHz processor and 8 GB of memory. Example 1. (Accuracy). Consider the following problem with zero Dirichlet boundary conditions: where The exact solution is u = sin ðπxÞ sin ðπyÞðc + t γ Þ with the two given nonnegative parameters c and γ.
We verify the accuracy of the proposed scheme (9) using two cases: the smooth and nonsmooth solutions. We first let c = 1 and γ = 2:1. The numerical results are obtained at T = 1 by fast Crank-Nicolson compact difference scheme (9) with no correction terms and demonstrated in Tables 2 and 3. One can observe that accuracy of the scheme is Oðτ min ð1+α,1+βÞ + h 4 Þ for different fractional orders α and β, which is in agreement with the theoretical analysis.
Next, for the nonsmooth case, we let c = 0 and γ = 0:4. One can see that the first-order partial derivative of u with respect to t is ∂ t uðx, y, tÞ = γt γ−1 sin ðπxÞ sin ðπyÞ, which is unbounded at t = 0 when γ = 0:4. By using the fast Crank-Nicolson compact difference scheme with correction terms (9), we compute the L 2 -norm errors at T = 0:5 and present the results in Tables 1 and 4. We can see from Table 1 that when m = 0, that is, no correction term is added to the scheme, the accuracy of the numerical solution suffers from the low regularity of the analytic solution. In contrast, when m is greater than 0, the accuracy of the numerical solution seems to be improved to some extents. Similar phenomenon is also observed in Table 4. This suggests that adding a small number of correction terms does improve the accuracy of the numerical solution in nonsmooth problems. Thus, the fast Crank-Nicolson compact difference scheme with correction terms (9) is valid for solving non-smooth solution problems.
Example 2. (Computational efficiency). In this example, we investigate the computational efficiency of the fast Crank-Nicolson compact difference scheme (7). So, we consider the comparison between results from the schemes with fast solver and the direct solver, that is, fast scheme (7) and original scheme (6). We separately solve the smooth solution case in Example 1 with the two numerical schemes and report the numerical results obtained in Figures 1-3. For the given fractional orders α and β, by fixing the time stepsize τ = 1/4 and 7 Journal of Function Spaces varying the stepsize in each spatial direction simultaneously, we obtain the CPU execution time at T = 1 for both schemes. The comparison shows that the execution time spent using the direct solver in numerical scheme is more expensive than that using the DST technique, especially when the spatial stepsize is getting smaller. This is due to the fact that the direct solver requires to solve matrix inversion on each time level, and such operation would be extremely inefficient when the size of the matrix is large. It is clear that the DST technique can speed up the computational efficiency, thus, the proposed scheme (7) has more potential than the direct solver (6) in high-dimensional problems.

Conclusions
In this paper, we propose the efficient compact difference scheme for solving the modified anomalous subdiffusion equation based on the modified L1 method in time and compact difference operator in space. By combining the DST technology, we improve the effectiveness of the scheme for the two-dimensional problem. The stability and error estimate of the scheme are provided rigorously. We also improve the accuracy of the scheme for the nonsmooth solu-tion problems with the method of adding correction terms. Numerical examples illustrate the effectiveness and accuracy of the proposed scheme.
The results of this paper can be readily generalized to three-dimensional problems. In addition, for inhomogeneous boundary conditions, one can convert them into homogeneous boundary condition problems by variable substitution. For other types of boundary condition problems, such as Neumann, Robin, or other combinations of boundary conditions, we do not discuss them in this paper. In [23], the authors introduced the augmented matched interface and boundary (AMIB) method to efficiently solving the Poisson equation via the FFT. The authors also pointed out that the AMIB method can easily handle different types of boundary conditions. So, it may be possible to combine this method with the correction terms to rapidly solve high-dimensional problems with complex boundary conditions, and this is the possible one of the future research directions.

Data Availability
The data of numerical simulation used to support the findings of this study are included within the article.