Abstract
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 . Here, and are the fractional orders which both range from 0 to 1, and and 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.
1. Introduction
In this paper, we focus on the numerical method for the time-fractional modified subdiffusion equation [1]:with the initial condition and the homogeneous Dirichlet boundary condition. Here, , is the rectangle domain, , and is the Laplacian defined by . The parameters and are some fixed positive constants, and are two given functions, , and 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 one-dimensional 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 in the new defined norm which is equivalent with -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–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 ill-conditioned 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 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.
2. The Compact Difference Scheme with Fast Solver
In this section, we derive the fast compact difference scheme for (1).
We first introduce the temporal discretization. The time stepsize is given by with the positive integer . The grid point is denoted by for . Let . For , the modified L1 method for the approximation of Riemann-Liouville derivative with at is described as:where ([22], Lemma 1). The operator in (2) is given bywhere
So, applying the difference discretization with , we derive thatwhere the local truncation error and .
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 unless otherwise noted. Denote the domain . Let be a positive integer. The spatial stepsize is then denoted as and for The fully discrete grids in space are denoted by . We further denote that and the boundary . The space of grid function is denoted as . For the grid function with the index vector at th position, we denote the compact difference operator as , with the difference operator . Here, is the identity operator, and So, the fourth-order spatial approximation of for is given by .
Combining the compact difference approximation in space with (4), we havewhere the local truncation error is given by . Here, . Omitting the small term , we obtain the following Crank-Nicolson compact difference scheme for (1): find of for , such thatwhere and .
If we solve the discretized system (6) directly, the computational cost will be 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 [11], which greatly improves the computational performance. Since the discrete sine transform of the grid function at the th position is provided by , it follows from the definition of the compact difference operator thatwhere and . One can refer to [11] or [10] for the derivation. Denote the index set . Therefore, the scheme (6) is equivalent to
The computational procedure is described as follows:(a)For , we first computed and from and by means of DST(b)And then we solve equation (10) from which the numerical solution is obtained from by the inverse of DST.
3. 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 , the discrete -norm is given by with the discrete inner product . The discrete seminorm and norm are denoted as and . Here, . One can readily have the equivalence of and for any in view of the embedding theorem.
We shall first need the following lemma.
Lemma 1. The operator given by (3) satisfies the inequality:where .
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 , we get
Notice that the difference operator is bounded in discrete inner product ([10], Theorem 2):with the notation . By the identity , the Lemma 1 yieldswhere , and . With , we write the above inequality as:
We sum up from 1 to and replace with to get
By the Cauchy-Schwarz inequality and the inequality with the equivalence of and , we obtain the estimate:from which we derive that
Next, we consider the case for the scheme (6) for the estimate of . By a similar procedure, we take the discrete inner product for (6) with when , we havefrom which we havewhere and . 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 and and the equivalence of the and , we further getwhich implies that the 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 , , , and . 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 , then we have the discrete -norm error estimate: For ,
Proof. The error equation can be obtained by subtracting (6) from (5), that is, by letting the error for , we haveIt follows from Theorem 2 thatwhich leads to the desired convergence result.
4. 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 and the first-order time derivative at as follows:
Here, and are the starting weights which are chosen such that the above schemes are exact for some power functions with and , 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 ,
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 and , one may refer to [19].
5. Numerical Examples
In this part, we present two numerical examples to verify the accuracy and effectiveness of the scheme (9). The -norm error at is obtained by , and the convergence orders in time and in space are calculated by and , respectively. For simplicity, we set the parameters and in (1) to be one and restrict the computational domain to be . 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:whereThe exact solution is with the two given nonnegative parameters and .
We verify the accuracy of the proposed scheme (9) using two cases: the smooth and nonsmooth solutions. We first let and . The numerical results are obtained at 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 for different fractional orders and , which is in agreement with the theoretical analysis.
Next, for the nonsmooth case, we let and . One can see that the first-order partial derivative of with respect to is , which is unbounded at when . By using the fast Crank-Nicolson compact difference scheme with correction terms (9), we compute the -norm errors at and present the results in Tables 1 and 4. We can see from Table 1 that when , 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 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 and varying the stepsize in each spatial direction simultaneously, we obtain the CPU execution time at 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.
6. 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 solution 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.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
This research was funded by Guangxi Natural Science Foundation (grant numbers 2018GXNSFBA281020 and 2018GXNSFAA138121).