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

Authors: Lijuan Nong; An Chen
Journal: Journal of Function Spaces
 (2021)

---

## 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 ofOτmin1+α,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.


---

## Body

## 1. Introduction

In this paper, we focus on the numerical method for the time-fractional modified subdiffusion equation [1]:(1)∂tux,t=κ1RLD0,t1−α+κ2RLD0,t1−βΔu+fx,t,x,t∈Ω×0,T,with the initial condition ux,0=vx and the homogeneous Dirichlet boundary condition. Here, x=x1,x2, Ω is the rectangle domain, T>0, and Δ is the Laplacian defined by Δ=∂2/∂x12+∂2/∂x22. The parameters κ1 and κ2 are some fixed positive constants, f and v are two given functions, 0<α,β<1, and RLD0,tγ is the Riemann-Liouville derivative of order γ given by:(2)RLD0,tγu·,t=1Γn−γ∂n∂tn∫0tt−sn−γ−1u·,sds,n−1<γ<n,n∈ℕ,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 ofOτ2minα,β+h4 in the new defined norm which is equivalent with H1-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.Table 1
TheL2-norm errors in time for nonsmooth case in Example 1 with h=1/64.
α,βnTm=0m=1m=3L2 errorRateL2 errorRateL2 errorRate0.3,0.8803.58E-03—1.71E-03—3.53E-04—1602.24E-030.671.03E-030.732.07E-040.773201.42E-030.666.28E-040.721.22E-040.766409.06E-040.653.83E-040.717.31E-050.740.5,0.6801.45E-03—7.86E-04—2.93E-04—1608.52E-040.764.11E-040.941.54E-040.933205.24E-040.702.15E-040.938.05E-050.936403.36E-040.641.13E-040.934.21E-050.940.7,0.4802.04E-03—1.05E-03—3.14E-04—1601.22E-030.735.86E-040.851.72E-040.873207.54E-040.703.28E-040.849.43E-050.876404.76E-040.661.84E-040.835.20E-050.86Figure 1
Comparison of CPU execution time between original and fast schemes with fixedα=0.3 and β=0.8.Figure 2
Comparison of CPU execution time between original and fast schemes with fixedα=0.5 and β=0.6Figure 3
Comparison of CPU execution time between original and fast schemes with fixedα=0.7 and β=0.4The rest of the paper is organized as follows. In Section2, 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τmin1+α,1+β+h4 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 τ=T/nT with the positive integer nT. The grid point is denoted by tn=nτ for n≥0. Let tn+1/2=tn+tn+1/2. For gt∈C20,T, the modified L1 method for the approximation of Riemann-Liouville derivative RLD0,tγgt with γ∈0,1 at t=tn+1/2 is described as:(3)RLD0,tγgtn+12=D¯τγgn+1/2+Rn+1/2,where ∣Rn+1/2∣≤cτ2−γmaxt∈0,T∣g′′t∣ ([22], Lemma 1). The operator D¯τγ in (2) is given by(4)D¯τγgn+1/2=b0γgtn+12−∑k=1nbn−kγ−bn−k+1γgtk−12−bnγ−Bnγgt12−Anγgt0,where(5)bkγ=τ−γΓ2−γk+11−γ−k1−γ,Bkγ=τ−γΓ2−γk+121−γ−k1−γ,Anγ=Bnγ−1−γτ−γΓ2−γn+12−γ.So, applying the difference discretization∂tux,tn+1/2=δtun+1/2+Oτ2 with δtun+1/2=un+1−un/τ, we derive that(6)δtun+1/2=κ1D¯τ1−α+κ2D¯τ1−βΔun+1/2+fn+1/2+Rxn+1/2,where the local truncation error Rxn+1/2=Oτmin1+α,1+β and un+1/2=un+1+un/2.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 Ω=x1L,x1R×x2L,x2R×⋯×xdL,xdR. Let Mk1≤k≤d be a positive integer. The spatial stepsize is then denoted as hk=xkR−xkL/Mk and xk,jk=xkL+jkhk for jk=0,1,⋯,Mk. The fully discrete grids in space are denoted by Ω¯h=xh=x1,j1,x2,j2,⋯,xd,jd∣0≤jk≤Mk,1≤k≤d. We further denote that Ωh=Ω¯h∩Ω and the boundary ∂Ωh=Ω¯h∩∂Ω. The space of grid function is denoted as Vh=v∣v=vhxhandvh=0forxh∈∂Ωh. For the grid function vh=vxh with the index vector h=i1,i2,⋯,id at kth position, we denote the compact difference operator as Δ¯kvik=δk2/Hkvik, with the difference operator Hkvik≔I+hk2/12δk2vik. Here, I is the identity operator, δk2vik=δkvik+1/2−δkvik−1/2/hk and δkvik+1/2≔vik+1−vik/hk. So, the fourth-order spatial approximation of Δvxh for xh∈Ωh is given by Δ¯hvh≔∑kΔ¯kvh.Combining the compact difference approximation in space with (4), we have(7)δtuxhn+1/2=κ1D¯τ1−α+κ2D¯τ1−βΔ¯huxhn+1/2+fn+1/2+Rxtn+1/2,where the local truncation error is given by Rxtn+1/2=Oτmin1+α,1+β+h4. Here, h4=h14+h24+⋯+hd4. Omitting the small term Rxtn+1/2, we obtain the following Crank-Nicolson compact difference scheme for (1): find uhn of uxh,tn for n≥1, such that(8)δtuhn+1/2=κ1D¯τ1−α+κ2D¯τ1−βΔ¯huhn+1/2+fn+1/2,where uh0=vxh and uxhxh∈∂Ωh=0.If we solve the discretized system (6) directly, the computational cost will be OM1M2⋯Md2 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 OM1M2⋯MdlogM1M2⋯Md [11], which greatly improves the computational performance. Since the discrete sine transform of the grid function vh at the kth position is provided by vik=∑jk=1Mk−1vjk^sinikjkπ/Mk, it follows from the definition of the compact difference operator Δ¯k that(9)vjk′′^≈vjk^12hk2·sjk−1sjk+5=vjk^λjk,Mk,where sjk=cosjkπ/Mk and 1≤jk≤Mk−1. One can refer to [11] or [10] for the derivation. Denote the index set ν=j1,j2,⋯,jd∣1≤jk≤Mk−1,1≤k≤d. Therefore, the scheme (6) is equivalent to(10)δtu^νn+1/2=κ1D¯τ1−α+κ2D¯τ1−β∑k=1dλjk,Mku^νn+1/2+f∧n+1/2.The computational procedure is described as follows:(a)
Forn≥0, we first computed u^νn and f∧n+1/2 from uhn and fn+1/2 by means of DST(b)
And then we solve equation (10) from which the numerical solution uhn is obtained from u^νn 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 functionv∈Vh, the discrete L2-norm is given by ∥v∥=v,vh with the discrete inner product u,vh=∏k=1dhk∑xh∈Ωhuhvh. The discrete H1 seminorm and H1 norm are denoted as v1=∇hvh2=∑k=1dδkvh2 and ∥v∥1=∥v∥2+v12. Here, ∇h=δ1,δ2,⋯,δd. One can readily have the equivalence of v1 and ∥v∥1 for any v∈Vh in view of the embedding theorem.We shall first need the following lemma.Lemma 1.
The operatorD¯τα given by (3) satisfies the inequality:(11)−2D¯ταvn+1/2,vn+1/2h≤∑k=1nbn−kγvk−1/22−∑k=1n+1bn+1−kγvk−1/22+Anγv02,where vn∈Vh,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(12)uhn+12≤cuh02+c1τα+c2τβ∇huh02+τ∑k=1n+1fk−1/22.Proof.
By taking the discrete inner product on both sides of (6) with 2τuhn+1/2, we get(13)2τδtuhn+1/2,uhn+1/2h=2τκ1D¯τ1−α+κ2D¯τ1−βΔ¯huhn+1/2,uhn+1/2h+2τfn+1/2,uhn+1/2h.Notice that the difference operatorΔ¯h is bounded in discrete inner product ([10], Theorem 2):(14)32Δhuhn+1/2,uhn+1/2h<Δ¯huhn+1/2,uhn+1/2h<Δhuhn+1/2,uhn+1/2h,with the notation Δhuhn=∑k=1dδk2uhn. By the identity Δhuhn+1/2,uhn+1/2h=−∇huhn+1/2,∇huhn+1/2h, the Lemma 1 yields(15)uhn+12−uhn2≤τ∑k=1nbn−k∇huhk−1/22−∑k=1n+1bn+1−k∇huhk−1/22+An∇huh02+2τfn+1/2,uhn+1/2h,where bk=bk1−α+bk1−β,Bk=Bk1−α+Bk1−β, and An=An1−α+An1−β. With Gn=∥uhn∥2+τ∑k=1nbn−k∇huhk−1/22, we write the above inequality as:(16)Gn+1≤Gn+τAn∇huh02+2τfn+1/2,uhn+1/2h.We sum upn from 1 to m and replace m with n to get(17)Gn+1≤G1+τ∑k=1nAk∇huh02+2τ∑k=1nfk+1/2,uhk+1/2h.By the Cauchy-Schwarz inequality and the inequalityvh≤cvh1 with the equivalence of vh1 and ∇hvh, we obtain the estimate:(18)2τ∑k=1nfk+1/2,uhk+1/2h≤τ∑k=1n+1bn+1−k∇huhk−1/22+τ∑k=2n+11bn+1−kfk−1/22,from which we derive that(19)∥uhn+1∥2≤G1+τ∑k=1nAk∇uh02+τ∑k=2n+11bn+1−kfk−1/22.Next, we consider the casen=0 for the scheme (6) for the estimate of G1. By a similar procedure, we take the discrete inner product for (6) with 2τuh1/2 when n=0, we have(20)2τδtuh1/2,uh1/2h=2τκ1D¯τ1−α+κ2D¯τ1−βΔ¯huh1/2,uh1/2h+2τf1/2,uh1/2h,from which we have(21)uh12+2τB0∇huh1/22=uh02+2A0τ∇huh0,∇huh1/2h+2τf1/2,uh1/2h,where B0=κ1B01−α+κ2B01−β and A0=κ1A01−α+κ2A01−β. Utilizing the Cauchy-Schwarz inequality again, we arrive at the estimate for the last two terms on the right-hand of the above inequality:(22)2A0τ∇huh0,∇huh1/2h+2τf1/2,uh1/2h≤2A0τ14ε1∇huh02+ε1∇huh1/22+2τ14ε2f1/22+ε2uh1/22.By letting the constantsε1=B0/4A0 and ε2=B0/4 and the equivalence of the ∥vh∥1 and vh1, we further get(23)2A0τ∇huh0,∇huh1/2h+2τf1/2,uh1/2h≤2τA02B0∇huh02+B0τ∇huh1/22+2τ1B0f1/22,which implies that the G1 has the following estimate:(24)G1=uh12+τb0uh1/22≤uh12+τB0uh1/22≤uh02+2τA02B0∇huh02+2τ1B0f1/22.Therefore, the inequality (8) yields(25)uhn+12≤uh02+2τA02B0∇huh02+2τ1B0f1/22+τ∑k=1nAk∇uh02+τ∑k=2n+11bn+1−kf1/22.By the mean value theorem, one can readily check that the coefficients appearing in the above inequality are all bounded, that is, we formally haveA02/B0≤B0=c1τα−1+c2τβ−1, 1/B0≤c3τmin1−α,1−β, τ∑k=1nAk≤c4τα+c5τβ, and ∑k=2n+11/bn+1−k≤c5maxT1−α,T1−β. Thus, the proof is completed.By means of the error equation and the stability conclusion, we have the following convergence result.Theorem 3.
Suppose thatu∈C20,T;C6Ω, then we have the discrete L2-norm error estimate: For n≥1,(26)utn−uhn≤cτmin1+α,1+β+h4.Proof.
The error equation can be obtained by subtracting (6) from (5), that is, by letting the error ehn=uxh,tn−uhn for xh∈Ωh, we have(27)δtehn+1/2=κ1D¯τ1−α+κ2D¯τ1−βΔ¯hehn+1/2+Rxtn+1/2.
It follows from Theorem2 that(28)ehn2≤ceh02+c1τα+c2τβ∇heh02+τ∑k=1nRxtk−1/22≤cτmin1+α,1+β+h42,which 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 RLD0,tγgt and the first-order time derivative dgt/dt at t=tn+1/2 as follows:(29)RLD0,tγgtt=tn+1/2≈D¯τγgtn+1/2red+∑k=1mwn,kγgtk−g0,dgtdtt=tn+1/2≈δtgtn+1/2red+∑k=1mwn,k1gtk−g0.Here,wn,kγ and wn,k1 are the starting weights which are chosen such that the above schemes are exact for some power functions gt=tζj with 0<ζj<ζj+1 and 0≤j≤m, that is, they can be determined by the two linear systems:(30)∑k=1mwn,kγtkζj=Γ1+ζjΓ1+ζj−γtn+1/2ζj−γ−D¯τγtn+1/2ζj,∑k=1mwn,k1tkζj=ζjtn+12ζj−1−tn+1ζj−tnζjτ,respectively. So, we have the following fast Crank-Nicolson compact difference scheme with correction terms: for n≥0,(31)δtu^νn+1/2+∑k=1mwn,k1u^νk=κ1D¯τ1−α+κ2D¯τ1−β∑k=1dλjk,Mku^νn+1/2+∑k=1mwn,k1−α+wn,k1−β∑k=1dλjk,Mku^νk+f∧n+1/2.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].


## 5. Numerical Examples

In this part, we present two numerical examples to verify the accuracy and effectiveness of the scheme (9). The L2-norm error at t=tn is obtained by en,h=uxh,tn−uhn, and the convergence orders in time and in space are calculated by logen,h/e2n,h and logen,h/en,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,12. 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:(32)∂tux,y,t=RLD0,t1−α+RLD0,t1−βΔux,y,t+fx,y,t,x,y∈Ω,ux,y,0=csinπxsinπy,where(33)fx,y,t=sinπxsinπyγtγ−1+2π2ctα−1Γα+2π2ctβ−1Γβ+2π2Γγ+1Γγ+αtγ+α−1+2π2Γγ+1Γγ+βtγ+β−1.
The exact solution isu=sinπxsinπyc+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τmin1+α,1+β+h4 for different fractional orders α and β, which is in agreement with the theoretical analysis.Table 2
TheL2-norm errors in time for smooth case in Example 1 with h=1/64.
nT0.3,0.80.5,0.60.7,0.4L2 errorRateL2 errorRateL2 errorRate203.87E-03—1.95E-03—2.55E-03—401.63E-031.257.32E-041.411.00E-031.34806.74E-041.272.67E-041.453.87E-041.371602.76E-041.299.59E-051.481.48E-041.39Table 3
TheL2-norm errors in space for smooth case in Example 1 with τ=T/8000.
M0.3,0.80.5,0.60.7,0.4L2 errorRateL2 errorRateL2 errorRate41.57E-03—1.58E-03—1.58E-03—89.82E-054.009.73E-054.029.75E-054.02167.72E-063.676.30E-063.956.63E-063.88Next, for the nonsmooth case, we letc=0 and γ=0.4. One can see that the first-order partial derivative of u with respect to t is ∂tux,y,t=γtγ−1sinπxsinπ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 L2-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 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.Table 4
TheL2-norm errors in space for nonsmooth case in Example 1 with τ=T/8000.
αMm=0m=1m=2L2 errorRateL2 errorRateL2 errorRate0.3,0.844.01E-04—5.32E-04—5.79E-04—81.59E-041.342.76E-054.271.91E-054.92161.93E-04-0.286.19E-05-1.161.52E-050.330.5,0.645.08E-04—5.87E-04—5.93E-04—85.28E-053.272.59E-054.503.16E-054.23168.72E-05-0.728.49E-061.612.80E-063.500.7,0.444.86E-04—5.74E-04—5.89E-04—87.49E-052.701.34E-055.422.81E-054.39161.09E-04-0.542.09E-05-0.646.22E-062.18


## 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.



---

*Source: 9918955-2021-04-27.xml*