Fast Compact Difference Scheme for Solving the Two-Dimensional Time-Fractional Cattaneo Equation

: The time-fractional Cattaneo equation is an equation where the fractional order α ∈ ( 1,2 ) has the capacity to model the anomalous dynamics of physical diffusion processes. In this paper, we consider an efﬁcient scheme for solving such an equation in two space dimensions. First, we obtain the space’s semi-discrete numerical scheme by using the compact difference operator in the spatial direction. Then, the semi-discrete scheme is converted to a low-order system by means of order reduction, and the fully discrete compact difference scheme is presented by applying the L2-1 σ formula. To improve the computational efﬁciency, we adopt the fast discrete Sine transform and sum-of-exponentials techniques for the compact difference operator and L2-1 σ difference operator, respectively, and derive the improved scheme with fast computations in both time and space. That aside, we also consider the graded meshes in the time direction to efﬁciently handle the weak singularity of the solution at the initial time. The stability and convergence of the numerical scheme under the uniform meshes are rigorously proven, and it is shown that the scheme has second-order and fourth-order accuracy in time and in space, respectively. Finally, numerical examples with high-dimensional problems are demonstrated to verify the accuracy and computational efﬁciency of the derived scheme.


Introduction
Let Ω ⊂ R 2 be a bounded rectangle domain and T be the fixed time.In this paper, we study the efficient numerical scheme for solving the following time-fractional Cattaneo equation: where the boundary and initial conditions are u(x, t) = 0 on ∂Ω × (0, T] and u(x, 0) = u 0 (x), ∂ t u(x, 0) = u 1 (x) in Ω.Here, x = (x 1 , x 2 ).The parameters κ and µ are the positive constants related to the relaxation time and diffusion properties, respectively.The source term f and the initial conditions u 0 and u 1 are given suitable functions.The Laplacian ∆ is defined by The Caputo derivative C D α 0,t with α ∈ (1, 2] is given by where Γ(•) is the Gamma function.
The time-fractional Cattaneo Equation (1) provides a flexible and efficient way to model the anomalous dynamics of physical diffusion processes, especially the general dynamics crossover behaviors [1,2].In recent years, many efforts have been made in the theoretical and numerical study of such a time-fractional equation or its variants.See [2,3] for the theoretical case and [1,[4][5][6] for the numerical case, just to name a few.One may refer to these two review papers [7,8] or this book [9] for further details.
In [4], Zhao and Sun proposed a compact Crank-Nicolson scheme for solving the time-fractional Cattaneo Equation (1).Ren and Gao also constructed an alternating direction implicit (ADI) compact difference scheme by adding a small term to improve the computational efficiency for the two-dimensional problem [5].The temporal convergence order of both of the above two methods is the order 3 − α, where 1 < α < 2. Later on, Chen and Li proposed a temporal second-order ADI Galerkin scheme by applying the Galerkin finite element method in space and the L2-1 σ method in time [6].In [1], Chen and Nong developed two efficient, fully discrete schemes for solving Equation (1) based on the Galerkin finite method in space and the convolution quadrature in time.Error estimates and numerical experiments show that their schemes are efficient when dealing with different data regularity in Equation (1).It seems that the schemes mentioned above still have some limitations, especially in the high-dimensional and non-smooth solution problems.This motivated us to develop an efficient numerical scheme to cope with such situations.
For the high-dimensional time-fractional model, there are two ways to improve the computational efficiency of the numerical scheme.One involves the spatial direction, such as with the ADI technique, fast discrete Sine transforms (DSTs), and so on (see [6,10]).The other one involves the temporal direction, such as with the fast convolution quadrature and sum-of-exponentials (SOE) techniques (see [11][12][13][14]).In this paper, we focus on the DST and SOE techniques.
The DST technique has attracted the interest of many scholars in recent years, mainly because it avoids the direct inversion for the coefficient matrix of the linear discrete systems and, at the same time, ensures the convergence accuracy of the difference operators [15,16].In [15], Wang, et al. proposed an efficient fourth-order compact finite difference scheme for solving the Poisson equation based on the fast discrete Sine transform and greatly reduced the computational cost by avoiding matrix inversion for the discretized system.There has been some work in recent years on the application of the DST technique in numerically solving high-dimensional time-fractional equations (see [10,17,18]).
The SOE technique is an efficient way to reduce the computational complexity caused by the non-locality of fractional derivatives [13,19].In [13], Yan et al. presented the fast L2-1 σ (denoted by FL2-1 σ ) formula by employing the SOE technique with the kernel function in the Caputo derivative.Subsequently, some numerical studies of time-fractional models based on this approach have emerged (see [10,19] and the corresponding references therein).In [19], Liang et al. proposed a fast difference scheme for solving the one-dimensional time-fractional telegraph equation based on the FL2-1 σ formula.In [10], Li et al. derived a fast compact difference scheme for solving subdiffusion equations by applying the FL2-1 σ in time and the compact difference operator in space.It is worth mentioning that the compact difference operator is implemented by the DST via the fast Fourier transform (FFT) algorithm.Therefore, their scheme is computationally efficient in both the time direction and spatial dimensions.
Motivated by the work in [10], we aim to develop a fast compact difference scheme for solving Equation (1).To this end, we derive a compact difference scheme by applying the L2-1 σ formula in time and compact difference operator in space.Then, the SOE and DST techniques are both adopted to improve the computational performance of the derived scheme.The contributions of our paper are as follows.First, we construct a compact difference scheme with fast computation in both the temporal and spatial directions, which effectively handles the high computational cost caused by the high spatial dimension and the non-locality of the Caputo derivative (see the numerical scheme (10)).Second, we present the rigorous stability and error estimate for the uniform mesh-based fast compact difference scheme (see Theorems 1 and 2).Third, extensive numerical experiments are designed to verify the accuracy and efficiency of the fast compact difference scheme (see the Section 5, Tables 1 and 2 for the accuracy; Table 5 and Figure 1 for efficiency).Besides, the fast compact difference scheme based on graded meshes in time is adopted to handle non-smooth solution problems (see Tables 3 and 4 in Section 5).
This paper is organized as follows.In Section 2, we present the fully discrete difference scheme based on the compact difference operator for the Laplacian in space and the L2-1 σ formula for the Caputo derivative in time.In order to improve the computational performance of the derived scheme, we further adopt the DST and SOE techniques to obtain the fast compact difference scheme in Section 3. In Section 4, the stability and error estimate of the fast compact difference scheme based on the uniform meshes are provided rigorously.Numerical examples and the conclusions of this paper are given in Sections 5 and 6, respectively.

The Compact Difference Scheme
In this section, we first approximate the Laplacian ∆ in Equation ( 1) by applying the compact difference operator and then find the low-order system for the space's semidiscrete scheme with the aid of the reduced-order method so that the fully discrete scheme of the equation can be obtained by using the second-order L2-1 σ formula.
To this end, we introduce some useful notations in the following.Suppose the rectangle domain ) is a positive integer.The fully discrete grids on Ω are given by We will use M without a subscript to denote M 1 (or M 2 ) when M 1 = M 2 , unless otherwise specified.The inner and boundary grids are denoted by Ω h = Ω h ∩ Ω and ∂Ω h = Ω h ∩ ∂Ω, respectively.We denote the space of the grid function as For the grid function v h = v(x h ) with the index vector h = (i 1 , i 2 ) at the kth position Here, we have When employing the compact difference operator ∆ h for the discretization of the Laplacian ∆ in (1), one has from which we obtain the following space semi-discrete scheme: where u h (t) ≈ u(x h , t) and the truncation error R x (t) = O(h 4 ).
Next, we introduce the L2-1 σ formula for the approximation of the Caputo derivative C D β 0,t with β ∈ (0, 1).Let the temporal step size be τ = T/N t .Here, N t is a positive integer.Denote t n = nτ, where n ≥ 0 and t n+ 1 2 = (t n + t n+1 )/2.For a given function g(t), denote and For more details about the difference operator δ t , one may refer to Lemma 2.1 of [20].
We state the L2-1 σ formula as follows [21]. where )) and for n = 0, w Here , By setting φ h (t) = ∂ t u h (t) for ( 2), we find the following lower-order system: Therefore, by replacing C D α−1 0,t with the difference operator of the L2-1 σ Formula (3) in the space semi-discrete scheme (2), one has where the truncation error R n xt = O(τ 2 + h 4 ) and r 0 = r n = O(τ 2 ).By dropping the small terms R n xt , r 0 , and r n , we have the following compact difference scheme for solving Equation (1).We find the numerical solution where

The Derivation of the Fast Compact Difference Scheme
In this part, we use the DST and SOE techniques to improve the computational performance of the compact difference scheme (6).
By utilizing the discrete sine transform, one can derive that where [15] for more details about the derivation).
Therefore, the scheme in Equation ( 6) is equivalent to where ν is the index set defined by We can obtain the numerical solution to Equation ( 1) by the following steps: (a) Compute φ 0 ν , u 0 ν , and f n ν from φ 0 h , u 0 h , and f n h by means of DST, where n ≥ 1; (b) Compute φ n ν by solving the system (7), and then obtain u n ν , where n ≥ 1; (c) Obtain the numerical solutions φ n h and u n h from φ n ν and u n ν , respectively, by means of the inverse of DST, where n ≥ 1.It is worth mentioning that if we solve numerical scheme (6) directly, the inverse of the coefficient matrix needs to be solved at each time layer, and this leads to a computational cost of O(M 2 1 M 2 2 ).If the DST technique based on FFT, that is the scheme ( 7) is used, then the computational cost will be reduced to Next, we consider the SOE technique to further reduce the computational complexity of the scheme (7) caused by the non-locality of the Caputo derivative.
Without loss of generality, we adopt the graded meshes t n = T( n N t ) γ in time, where n = 0, • • • , N t and γ ≥ 1 is a proper chosen temporal mesh grading parameter.Denote the . The fast L2-1 σ formula based on the graded meshes is defined as where 0 = (στ n+1 ) 1−β Γ(2−β) and H k (t j ) (j ≥ 0, k ≥ 1) is obtained by the following recurrence: with the initial value H k (t 0 ) = 0 [10].Here, the coefficients The weights k (k ≥ 1) and the points s k (k ≥ 1) in Formula (8) are chosen to satisfy where is the absolute tolerance error and the number of exponentials N exp satisfies We remark that the letter F in the difference operator F D β τ n (see the left-hand side of the Formula (8)) represents the word "fast", which refers to the fact that the corresponding formula uses the SOE technique to improve the computational performance.It is shown in [13] that the number of exponentials N exp = O(log N t ) when T 1 and N exp = O log 2 N t when T ≈ 1.The total computational cost for the fast L2-1 σ Formula (8) is O(N t N exp ) with storage O(N exp ).
By replacing the classical L2-1 σ formula in the numerical scheme (7) with the above fast L2-1 σ formula based on graded meshes, we have the following fast compact difference scheme: where We will simply denote the fast compact difference scheme (10) as the FL2-1 σ (γ) scheme in what follows if no ambiguity arises.
One can observe that, compared with the numerical scheme ( 6), the FL2-1 σ (γ) scheme (10) reduces the overall computational cost from O( ), which greatly improves the computational efficiency of the original scheme (6).
When γ = 1, the graded meshes recover the uniform one.Thus, in this case, we shall use τ instead of τ n in (10) unless otherwise specified.The corresponding scheme based on uniform meshes is called the FL2-1 σ (1) scheme.When γ > 1, the grid points in the time direction are concentrated near the origin, which is beneficial for the numerical solution to capture the weak singularity solution, thus improving the accuracy loss of the scheme caused by the non-smooth solution.This will be verified by numerical examples in Section 5.

Stability and Error Estimates
We studied the stability and error estimates for the FL2-1 σ (1) scheme (10) in this section.For any given grid function v ∈ V h , the discrete L 2 norm is given by v = (v, v) h with the discrete inner product (u, v) h = (h 1 h 2 ) ∑ x h ∈Ω h u h v h .The discrete H 1 semi-norm and H 1 norm are denoted as 1 , respectively.Here, ∇ h = (δ 1 , δ 2 ).In light of the embedding theorem, we can conclude that for any given v ∈ V h , the discrete H 1 semi-norm |v| 1 and H 1 norm v 1 are equivalent.
Using the recursive relationship of the historical term H k (t j ) in Equation ( 9), one may reformulate the expression of the difference operator F D β τ based on the uniform meshes as follows (or see [10,13] more details on the derivation): 0 when n = 0, and for n ≥ 1, we have Here,
We need the following lemma: Lemma 1.The weights W (n+1) k in Equation (12) have the following properties: Proof.In light of Lemma 3.3 in [19] and the definitions of in Equation ( 12), one may readily obtain the desired results.
Let the symbol c be a positive constant which is independent of the temporal and spatial step sizes.The stability of the scheme ( 10) is described as follows: Theorem 1.The FL2-1 σ (γ) compact difference scheme (10) where ε n+σ h and ε h are the perturbations of h , respectively.
Proof.In light of the equivalence of schemes ( 6) and ( 7), one only needs to consider the following scheme instead of (10): We first consider the case n ≥ 1 for (13).By taking the discrete inner product on both sides of (13) with ϕ n+σ h , we have We can use Lemma 1 and the estimate for the discrete operator ∆ h presented in (3.14) of [10] such that From this, we derive that Note that δ t u n h = φ n+σ h + ε n+σ h , so it follows from Lemma 3.5 of [20] and the Cauchy-Schwarz inequality that By using the inequality ab ≤ a 2 + 1 4 b 2 with > 0, we obtain Therefore, the inequality (15) has the following estimate: In other words, we have We sum up n from n = 1 to n = m − 1 for the above inequality and obtain we have where Lemma 1 and the boundness of W (m) m−k presented in Lemma 4.1 of [13] are used.Here, The term E m has the estimate E m ≥ 1 σ ∇ h u m h 2 (see Lemma 3.5 of [20]).Thus, one has the further estimate for the inequality ( 16): by choosing suitable 1 and 2 values, the inequality ( 17) yields Next, we consider the bound of φ 1 2 on the right-hand side of (18).To this end, we operate the discrete inner product on both sides of the numerical scheme (6) with φ 1 2 h for n = 0 and find 0 (φ 1 h − φ 0 h ), and using the equality δ t u h , we obtain Following the idea in the proof of Theorem 4.1 of [6] for the case n = 0, we further have where we choose 1 = σ − 1 2 > 0.
For the first term on the right-hand side of Equation ( 19), we have where we let 2 = σ−0.5 2σ−1 > 0. The second term on the right-hand side of Equation ( 19) has the following estimate: By suitably choosing the parameters 3 and 4 ( e.g., let 3 = 1 σ ( 1 2τ − 1) > 0 and 4 = 1−σ 4 > 0), and noting that W (1) This together with (18) leads to By applying the Gronwall inequality to the above estimate and the equivalence of the H 1 norm and H 1 semi-norm, one has the desired result.All of this ends the proof.
Finally, we present the convergence result for the FL2-1 σ (1) scheme (10).Let Therefore, in light of (2), (6), and (13), one may find the following error equations: . By applying the stability result of Theorem 1, we obtain the error bound in the H 1 norm as follows: The above estimate leads to the following error estimate of the FL2-1 σ (1) scheme (10).The corresponding proof is thus omitted: Theorem 2. Let u(x h , t n ) and u n h be the solutions to Equation (1) and numerical scheme (10), respectively.Suppose the solution u(x, t) belongs to the function space C 6,3  x,t (Ω × [0, T]).Then, for γ = 1, we have Remark 1.In practice, the absolute tolerance error in the scheme (10) is always set to a very small number so as not to affect the temporal or spatial convergence accuracy.In this sense, the convergence result in Theorem 2 can be understood as O(τ 2 + h 4 ).
Remark 2. Although the nonuniform mesh-based FL2-1 σ (γ) scheme (i.e., γ > 1) can effectively handle non-smooth solution problems, its stability and error estimates are more difficult.In [10], Li et al. successfully gave a rigorous proof of the stability and convergence of the FL2-1 σ scheme based on the nonuniform meshes, but their derived scheme mainly deals with the subdiffusion problem.It seems that the application of their ideas here is not obvious, and further investigation is needed.

Numerical Examples
Two numerical examples will be presented to verify the convergence theoretical results and efficiency of the FL2-1 σ scheme (10).In the numerical examples, the computational domain is Ω = (0, 1) 2 , and the parameters κ and µ are both one.Using the equivalence of the H 1 norm and H 1 semi-norm, we calculated the H 1 norm error at t = t n by e(n, h) = |u(x h , t n ) − u n h | 1 .The corresponding convergence orders in time and in space were obtained by log(e(n, h)/e(2n, h)) and log(e(n, h)/e(n, h/2)), respectively.We always set the absolute tolerance error = 10 −10 in (10) so that it did not contaminate the temporal or spatial convergence orders.
We considered two cases, η = 3.5 and η = 1.5, to verify the accuracy of the FL2-1 σ (γ) scheme (10).For these two cases, we applied the uniform mesh-based and nonuniform mesh-based FL2-1 σ (γ) scheme (10), respectively, to obtain the numerical results.See Tables 1 and 2 for the first case and Tables 3 and 4 for the second case.We remark that the H 1 norm errors here and below were obtained at the final time T.   When η = 3.5, the solution of the equation satisfied the smoothness requirement of Theorem 2. From the numerical results in Tables 1 and 2, it can be observed that the accuracy of the scheme (10) was O(τ 2 + h 4 ) for different fractional orders α, which is consistent with the theoretical result of the convergence accuracy.
When η = 1.5, the second-order partial derivative of u with respect to t was ∂ 2 t u(x, y, t) = η(η − 1)t η−2 sin(πx) sin(πy), which was unbounded at t = 0 when η < 2. Therefore, the solution was not sufficiently smooth in this case.In order to capture the weak singularity solution at the initial point, we let the temporal mesh grading parameter be γ > 1.From the numerical results in Tables 3 and 4, one can see that by selecting a different γ value (i.e., γ = 1.5 and γ = 2), the theoretical convergence of the scheme (10) was restored, which shows that the scheme (10) based on graded meshes indeed captured the weak singularity solution and could significantly improve the convergence accuracy.
Example 2 (Computational efficiency).In this example, we studied the computational performance for the FL2-1 σ scheme (10).For the sake of simplicity, the case η = 3.5 in the above example was used.To better reveal that the scheme (10) had improved computational efficiency, we compared it with the two schemes (6) and (7).The numerical results are demonstrated in Table 5 and Figure 1.From Table 5, we can observe that when the spatial step size h was fixed, the accuracy of all three schemes increased as the number of temporal nodes increased.However, their computational cost also increased.With the fixed same number of spatial and temporal nodes, the accuracies of these threes scheme were almost the same, but the scheme (6) based on the direct solver took the most time, followed by the scheme (7), while the FL2-1 σ (γ) scheme (10) took the least time, and this time consumption gap increased with the increase in the number of temporal nodes.
For the case where α = 1.5, when the time step size τ was fixed to 1/4, and the spatial grid became denser, similar phenomena to those in Table 5 could be observed in the left subplot of Figure 1.
In particular, we compared the computational efficiency of the two numerical schemes (7) and (10) by plotting the CPU time versus N t in a log-log style.In order to numerically quantify the complexity of these two schemes, we expressed the CPU time in the form of N r t using least squares fitting and reported the corresponding order r in the legend of the right subplot of Figure 1.One can see that the order r of the FL2-1 σ (γ) scheme (10) was about one, while that of the DST-based scheme (7) was near two.This indicates that the proposed scheme (10) based on the DST and SOE techniques is more competitive in solving high-dimensional time-fractional problems.For other cases, such as α = 1.1 and 1.9, similar numerical results were also obtained, and they are not presented here for the sake of brevity.

Conclusions
In this paper, we proposed a fast compact difference scheme with O(τ 2 + h 4 ) convergence accuracy to solve the two-dimensional time-fractional Cattaneo Equation (1).This scheme is based on the DST and SOE techniques and has the ability of fast calculation in both time and space.The stability and error estimate of the numerical scheme (10) based on uniform meshes are presented rigorously.Numerical examples show that this scheme is efficient and suitable for numerically solving high-dimensional time-fractional problems.
We remark that the results obtained in this paper can easily be generalized to other high-dimensional problems.Aside from that, one may notice that although numerical experiments have shown that the scheme (10) based on graded meshes can improve the accuracy of solving high-dimensional non-smooth solution problems, the analysis of its corresponding numerical theory is not an easy job.Such an issue deserves further study and will be one of our upcoming research directions.
We note that the scheme presented in this paper is only for linear problems with constant coefficients, while in practical problems, the equations are often accompanied by variable coefficients or nonlinear terms.In the following work, we will try to extend the methods of this paper to solve variable coefficients and nonlinear problems.

Figure 1 .
Figure 1.Comparison of computational costs with fixed α = 1.5.(The dashed lines below the numerical curves are the corresponding fitted curves, which are shifted down here for the convenience of observation.)

Table 1 .
The H 1 -norm errors in time for the smooth case in Example 1 with h = 1/64.

Table 2 .
The H 1 -norm errors in space for the smooth case in Example 1 with τ = 1/4000.

Table 3 .
The H 1 -norm errors in time for the non-smooth case in Example 1 with h = 1/64 and α = 1.5.