1. Introduction
Let
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
on
and
,
in
. Here,
.The parameters
and
are the positive constants related to the relaxation time and diffusion properties, respectively. The source term
f and the initial conditions
and
are given suitable functions. The Laplacian
is defined by
. The Caputo derivative
with
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
, where
. 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,
Table 1 and
Table 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
Table 3 and
Table 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
Section 5 and
Section 6, respectively.
2. 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 semi-discrete 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 . Let the spatial step size be and the grid points be , where is a positive integer. The fully discrete grids on are given by . We will use M without a subscript to denote (or ) when , unless otherwise specified. The inner and boundary grids are denoted by and , respectively. We denote the space of the grid function as .
For the grid function
with the index vector
at the
kth position
, the compact difference operator is given by
, where
Therefore, for , we have with and .
When employing the compact difference operator
for the discretization of the Laplacian
in (
1), one has
from which we obtain the following space semi-discrete scheme:
where
and the truncation error
.
Next, we introduce the L2-1
formula for the approximation of the Caputo derivative
with
. Let the temporal step size be
. Here,
is a positive integer. Denote
, where
and
. For a given function
, denote
and
, where
. For more details about the difference operator
, one may refer to Lemma 2.1 of [
20].
We state the L2-
formula as follows [
21]. If
, then
where
and for
, where
, we have
Here,
, and
By setting
for (
2), we find the following lower-order system:
Therefore, by replacing
with the difference operator
of the L2-1
Formula (
3) in the space semi-discrete scheme (
2), one has
where the truncation error
and
. By dropping the small terms
,
, and
, we have the following compact difference scheme for solving Equation (
1). We find the numerical solution
of
for
such that
where
,
, and
. From (
6), we obtain the numerical solution
of
to (
1):
and for
:
3. 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
and
(see [
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 , , and from , , and by means of DST, where ;
- (b)
Compute
by solving the system (
7), and then obtain
, where
;
- (c)
Obtain the numerical solutions and from and , respectively, by means of the inverse of DST, where .
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
. 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
in time, where
and
is a proper chosen temporal mesh grading parameter. Denote the time step size
. Set
and
. The fast L2-1
formula based on the graded meshes is defined as
where
and
is obtained by the following recurrence:
with the initial value
[
10]. Here, the coefficients
and
. The weights
and the points
in Formula (
8) are chosen to satisfy
where
is the absolute tolerance error and the number of exponentials
satisfies
We remark that the letter
F in the difference operator
(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
when
and
when
. The total computational cost for the fast L2-1
Formula (
8) is
with storage
.
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
and
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
to
, which greatly improves the computational efficiency of the original scheme (
6).
When
, the graded meshes recover the uniform one. Thus, in this case, we shall use
instead of
in (
10) unless otherwise specified. The corresponding scheme based on uniform meshes is called the FL2-1
(1) scheme. When
, 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.
4. 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 , the discrete norm is given by with the discrete inner product . The discrete semi-norm and norm are denoted as and , respectively. Here, . In light of the embedding theorem, we can conclude that for any given , the discrete semi-norm and norm are equivalent.
Using the recursive relationship of the historical term
in Equation (
9), one may reformulate the expression of the difference operator
based on the uniform meshes as follows (or see [
10,
13] more details on the derivation):
where
are defined as
when
, and for
, we have
Here,
,
, and
We need the following lemma:
Lemma 1. The weights in Equation (12) have the following properties:and 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) with is stable in the sense thatwhere and are the perturbations of and , 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
for (
13). By taking the discrete inner product on both sides of (
13) with
, we have
We can use Lemma 1 and the estimate for the discrete operator
presented in (3.14) of [
10] such that
From this, we derive that
Note that
, so it follows from Lemma 3.5 of [
20] and the Cauchy–Schwarz inequality that
where
.
By using the inequality
with
, we obtain
Therefore, the inequality (
15) has the following estimate:
In other words, we have
We sum up
n from
to
for the above inequality and obtain
Since
we have
where Lemma 1 and the boundness of
presented in Lemma 4.1 of [
13] are used. Here,
and
.
The term
has the estimate
(see Lemma 3.5 of [
20]). Thus, one has the further estimate for the inequality (
16):
Since
by choosing suitable
and
values, the inequality (
17) yields
Next, we consider the bound of
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
for
and find
Observing that
, and using the equality
, we obtain
Following the idea in the proof of Theorem 4.1 of [
6] for the case
, we further have
where we choose
.
For the first term on the right-hand side of Equation (
19), we have
where we let
.
The second term on the right-hand side of Equation (
19) has the following estimate:
By suitably choosing the parameters
and
( e.g., let
and
), and noting that
, we can derive that
This together with (
18) leads to
By applying the Gronwall inequality to the above estimate and the equivalence of the
norm and
semi-norm, one has the desired result. All of this ends the proof. □
Finally, we present the convergence result for the FL2-1
scheme (
10). Let
and
. From Lemma 2.3 of [
13], we have
Therefore, in light of (
2), (
6), and (
13), one may find the following error equations:
where
,
, and
. By applying the stability result of Theorem 1, we obtain the error bound in the
norm as follows:
The above estimate leads to the following error estimate of the FL2-1
scheme (
10). The corresponding proof is thus omitted:
Theorem 2. Let and be the solutions to Equation (1) and numerical scheme (10), respectively. Suppose the solution belongs to the function space . Then, for , 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 . Remark 2. Although the nonuniform mesh-based FL2-1 scheme (i.e., ) 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. 6. Conclusions
In this paper, we proposed a fast compact difference scheme with
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.