Fast finite difference/Legendre spectral collocation approximations for a tempered time-fractional diffusion equation

1.   Introduction

Fractional differential equations have attracted extensive attention for their ability to accurately describe complex physical phenomena in both micro and macro perspectives ^[1,2,3]. The differential equations with fractional derivative in time variable can better describe genetic and memory characteristics of some physical processes. For different physical situations, different memory function arise from various fractional derivative operators, such as the power memory kernel, exponential memory kernel, and Mittag-Leffler memory kernel functions ^[4].

In this work, we consider the following initial boundary value problem:

where $\Delta = \frac{\partial^2}{\partial x^{2}_{1}}+\frac{\partial^2}{\partial x^{2}_{2}}$ denotes the Laplace operator, $\Omega = [-1, 1]^{2}$ denotes the computational domain with boundary $\partial \Omega$, $\textbf{x} = ({x_{1}, x_{2}})\in \Omega$ stands for the space variables, and $\partial_{0, t}^{\alpha, \lambda(t)}$ denotes the generalized Caputo time-fractional derivative ^[5]

In fractional derivative (1.2), we assume that the kernel function satisfies $\lambda(t)\in C^{2}[0, T]$, $\lambda(t) > 0$ and $\lambda^{\prime}(t)\le 0$ for all $t\in [0, T]$. Particularly, if kernel takes $\lambda(t) = e^{-b t}(b > 0)$, fractional derivative (1.2) gives

provides the tempered fractional derivative and integral

Obviously, from a mathematical perspective, the Caputo-type fractional derivative (1.2) can be interpreted as a kind of generalized fractional derivative ^[4,6,7]. From Figure 1, we can see the influence of parameters on the solution. More application (in macroscopic and microscopic scales) of the model given in problem (1.1) with tempered fractional derivative, see references ^[3,8,9]. There are many research works on the discretization of fractional derivative (1.2). By splitting into two parts, Alikhanov first discussed the $L1$ formula of fractional derivative (1.2) in the work ^[5], and proposed a new difference scheme for the time-fractional diffusion equation with the fractional derivatives (1.2). Recently, in view of relation (1.4), Chen et al. ^[10] established two different $L1$ discretizations for tempered Caputo fractional derivative on the graded meshes, which developed in reference ^[11].

The nonlocality of the fractional derivative will make storage and computation very expensive if we adopt a lower-order method. For example, if we use the $L1$ scheme to calculate (1.1), the required computational and storage amounts are $O(NN_ {T})$ and $O(NN_{T}^{2})$, respectively, where $N$ is the number of spatial grid points and $N_{T}$ stands for the number of grid points in time. When the time step is small enough, it means that the numerical calculation needs expensive storage costs and computational costs, so it is particularly important to develop fast or high-order numerical schemes. The sum-of-exponential (SOE) technique is one of the effective techniques to approximate the fractional derivatives with the power kernel function ^[12,13]. In this topic, Jiang et al. ^[13] used the SOE method to approximate the kernel function $t^{-(\alpha+1)}(0 < \alpha < 1)$ and solved fractional diffusion equations on unbounded domains. Yan et al. ^[14] combined the $L2$-$1_{\sigma}$ scheme with SOE approximation successfully solved an initial boundary value problem of the time fractional diffusion equation. Combining the SOE approximation, Xu et al. ^[15] proposed a fast $L2$ discretization of Caputo fractional derivatives, and they get the error estimates by rearranging the coefficients of the fast $L2$ difference approximation. With the help of the SOE technique, Gu et al. ^[16] proposed a fast difference algorithm to solve the generalized time-space fractional diffusion equation with fractional derivative (1.3). It seems that the authors of ^[16] fail to give any convergence analysis of the fast $L1$ difference scheme.

As a higher-order numerical method, the spectral method has been widely used for solving many kinds of differential equations ^[17,18] including fractional differential equations ^[15]. Lin and Xu ^[19] solved the time-fractional diffusion equation based on the finite difference method in the time direction and the Legendre spectral method in space. Later, Li and Xu ^[20] improved their previous results and proposed a time-space spectral method for the time-fractional diffusion equation. Compared with one-dimensional problems, solving high-dimensional problems is more difficult due to the increase in storage and computational cost. Zeng et al. ^[21] proposed the Legendre Galerkin spectral method for two-dimensional Riesz space fractional nonlinear reaction-diffusion equations and obtained the optimal spatial error estimation. Guo et al. ^[22] used the Legendre spectral method to solve two-dimensional fractional nonlinear reaction-diffusion wave equations. Recently, combining the fast $L2$ discretization of Caputo fractional derivatives, Xu et al. ^[15] proposed a fast difference/spectral method for a time fractional equation in one dimension. In their numerical scheme, the convergence order in the temporal direction is $(3-\alpha)$ and the convergence order in the spatial direction is $O(N^{1-m})$, where $m$ is the regularity of the exact solution in the spatial direction. This work suggests that the fast difference approximation with high-order discretization schemes is an effective way to solve the time fractional differential equations. Following this line, in this paper, we try to develop a fast finite difference/Legendre spectral collocation method for problem (1.1) with smooth solution and discuss the error estimates for the present numerical schemes. First, the generalized fractional derivative is discretized by $L1$ discretization, and the time direction convergence order is analyzed. Secondly, by using SOE approximation of exponential kernel function, a fast $L1$ scheme with exponential memory kernel function is obtained. Therefore, it can be seen that the calculation and storage costs of the fast scheme are significantly reduced. With the similar technique given in ^[15], we will discuss the error estimates when the fast L1 method is used to solve problem (1.1).

The rest of this paper is organized as follows: In Section 2, we present time discretizations for problem (1.1) and discuss the stability and convergence of the present semi-discrete scheme. Combing the Legendre spectral collocation method for spatial variables, in Section 3, we discuss the finite difference/spectral collocation method for problem (1.1). The error estimates of the fully discrete scheme are also established in this section. For the generalized fractional derivative with exponential kernel, we develop the sum-of-exponentials technique developed in ^[13] to accelerate the calculation process. The error estimates of the related numerical schemes are discussed in Section 4. In Section 5, we present two numerical examples to verify the efficiency of our numerical algorithm and theoretical results. Finally, some conclusions are provided in Section 6.

2.   Time semi-discrete scheme

First, we discretize the fractional derivative with a generalized memory kernel, assuming the solution $u(t)\in C^{2}[0, T]$. The computational domain $[0, T]$ is uniformly divided into $N_{T}$ intervals with the uniform grids $t_{k} = k\tau$, $k = 0, 1, \ldots, N_{T}$. For $\lambda(t)\in C^{2}[0, T]$, the time fractional derivative with a generalized memory kernel is discretized as ^[5]

where

and the error

give ^[5]

where $\tilde{c}$ depends on $\alpha$.

We rewrite (2.1) as follows:

where $\alpha_{0} = \tau^{\alpha}\Gamma(2-\alpha), c_{l} = (\lambda_{l+1/2}a_{l}+(\lambda_{l}-\lambda_{l+1})b_{l})(l\ge 0)$, and the truncation error gives $r_{\tau}^{k+1} = R_{1}^{k}+R_{2}^{k}$. For the sake of simplification, we define the fractional differential operator $L_{t}^{\alpha, \lambda(t)}u(\textbf{x}, t_{k+1})$ by

Using $L_{t}^{\alpha, \lambda(t)}u(\textbf{x}, t_{k+1})$ to approximate $\partial_{0, t}^{\alpha, \lambda(t)}u(\textbf{x}, t)$ in problem (1.1) at $t_{k+1}$, and denoting $u^{k+1}$ is an approximation of $u(\textbf{x}, t_{k+1})$, $f^{k+1} = f(\textbf{x}, t_{k+1})$, we obtain the time semi-discrete scheme of problem (1.1), i.e.,

for $k = 0$.

For $1\leq k\leq N_T-1$,

After a simple calculation, the truncation error of the semi-discretized schemes (2.6) and (2.7) give $r^{k+1} = \alpha_{0}r_{\tau}^{k+1}$ and satisfies

Considering the variational formulation of (1.1), we define some functional spaces endowed with standard norms and inner products that will be used hereafter ^[17,18,26]

where $\beta = (\beta_{1}, \beta_{2})$, $|\beta| = \beta_{1}+\beta_{2}, D^{\beta} = \frac{\partial^{|\beta|}}{\partial x_{1}^{\beta_{1}}\partial x_{2}^{\beta_{2}}}$, $L^{2}(\Omega)$ is the space of measurable functions whose square is Lebesgue integrable in $\Omega$, which equip with the inner products

and the norm$H^{1}$ is defined by

The weak formulation of the (2.3) and (2.4) reads: Find $u^{k+1}\in H_{0}^{1}(\Omega)$, such that

for $k = 0$, and

for all $v\in H_{0}^{1}(\Omega)$ and $1\leq k\leq N_T-1$.

2.1.   Stability analysis of the semi-discrete scheme

Lemma 1. ^[5] If $\lambda(t) > 0$, $\lambda^{\prime}(t)\le 0,$ and $\lambda(t)\in C^{2}[0, T]$, the weighted coefficients $a_{l}$, $b_{l}$ and $c_{l}\; (l = 0, 1, \ldots)$ are given by (2.1), satisfy

Theorem 1. For all $\tau > 0$, the semi-discretized schemes (2.6) and (2.7) are unconditionally stable, i.e.,

Proof. For $k = 0$, substituting $v = u^{1}$ in (2.6) and using $\|u^{1}\|\leq \|u^{1}\|_{1}$, we have

From (2.9), we have $\|u^{1}\|_{1}\le \|u^{0}\|+\tilde{c} \alpha_{0}c^{-1}_{0}\max\limits_{0\le j\le 1}\|f^{j}\|$. For $k > 1$, we prove (2.8) by mathematical induction. Assume that there holds

Now prove (2.8) using inequality (2.10). Substituting $v = u^{k+1}$ into (2.7), we obtain

which implies

Using Lemma 1 and the inequality $\|u^{k+1}\|\le \|u^{k+1}\|_{1}$, we have

Note that $\|u^{0}\|+\tilde{c}\alpha_{0}c^{-1}_{0}\max\limits_{0\le j\le k-j}\|f^{j}\|\le \|u^{0}\|+\tilde{c}\alpha_{0}c^{-1}_{0}\max\limits_{0\leq j\leq k}\|f^{j}\|(j = 0, 1, \ldots, k-1)$. Therefore

Finally, we obtain (2.8). □

2.2.   Convergence analysis of the semi-discrete scheme

Here we provide some lemmas to analyze coefficients $c_{l}$ and $a_{l}$, which will help to study the convergence analysis of the solution. In the following sections, $\tilde{c}$ is used to represent a constant, which may be different under different situations. Here, we employ the technique used in reference ^[19] to do the error estimates of schemes (2.6) and (2.7).

Theorem 2. Let $u(\textbf{x}, t_{k})$ be the exact solution of problem (1.1), $\{u^{k}\}_{k = 0}^{N_{T}}$ be the solution of schemes (2.6) and (2.7), there holds

where $\tilde{c}$ is independent of $u, \tau$ and $T$.

Proof. Firstly, we prove that the following inequality

Let $e^{j} = u(\textbf{x}, t_{j})-u^{j}$, for $j = 1$, we obtain

If we assume the initial value is exact (i.e., $e^{0} = 0$), let $v = e^{1}$ in (2.14), using (2.5), we obtain

Suppose (2.13) holds for all $j = 1, 2, \ldots, k$, we need to prove it also holds for $j = k+1$. In view of

Let $v = e^{k+1}$ in (2.15), we have

We rewrite (2.16) as follows:

Using Lemma 1, we have $c_{k-j}^{-1} < c_{k}^{-1}$, it follows that

which implies (2.13). Lemma 1 suggests that

Hence, we have

Next we consider the case $\alpha\rightarrow 1$, we firstly prove that

Inequality (2.17) obviously holds for $j = 1$. Suppose (2.17) holds for $j = 1, 2, \ldots, k$. Now we will prove that (2.17) still holds when $j = k+1$. In fact, with the similar method given in Theorem 1, we have

Simple calculations, we can check that

With the help of (2.18), Eq (2.19) provides that

Since

we have

Substituting (2.21) into (2.20), we obtain

We finally obtain our conclusion. □

3.   The fully discrete scheme

Let $P_{N}(I_{x})$ and $P_{N}(I_{y})$ be the spaces of all polynomials of degree less than or equal to $N$ defined on domains $I_{x} = (-1, 1)$ and $I_{y} = (-1, 1)$, respectively. Let $V_{N}(\Omega) = P_{N}(I_{x})\bigotimes P_{N}(I_{y})$, $\{\xi_{j}\}_{j = 0}^{N}$, $\{\xi_{i}\}_{i = 0}^{N}$ are Legendre-Gauss-Labatoo (LGL) points, and $\{\omega_{j}\}_{j = 0}^{N}$ and $\{\omega_{i}\}_{i = 0}^{N}$ are the weights, which satisfy the following quadrature formula ^[25,26]

Define the inner product and norm as follows ^[18,25,26]

We will use the $H_{0}^{1}$-orthogonal projection operator $\pi_{N}$ such that ^[24]

for all $u \in H_{0}^{1}(\Omega)$. Here $V_{N}^{0}(\Omega)$ is defined as

For this projection, there have estimates ^[25,26]

Lemma 2. ^[25] For all $\varphi \in H^{m}(\Omega)(\ m\ge 1)$, $v_{N}\in P_{N}(\Omega)$, there holds

Finite difference/Legendre collocation approximation of problem (1.1) gives: Find $u_{N}^{k}\in V_{N}^{0}(\Omega)$,

for all $v_{N}\in V_{N}^{0}(\Omega)$. The discrete $H^{1}$-norm is defined by

To do the error estimate, we split the error as $u(t_{k})-u_{N}^{k} = (u(t_{k})-\pi_{N}u(t_{k}))-(u_{N}^{k}-\pi_{N}u(t_{k})): = \epsilon_{N}^{k}-e_{N}^{k}$. Without loss of generality, we only consider the homogeneous case of problem (1.1), i.e., $f = 0$.

Theorem 3. Let $u(\textbf{x}, t)$ be the exact solution of problem (1.1), and $u_{N}^{k}$ be the solution of fully discrete scheme (3.4) with the initial condition $u_{N}^{0} = \pi_{N}u^{0}$. Suppose $\partial_{t}^{2}u \in L^{\infty}((0, T]; H^{m}(\Omega))(m\ge 1)$, then

where $\|\partial_{t}^{\alpha, \lambda(t)}u\|_{L^{\infty}(H^{m})} = \mathrm{sup}_{t\in (0, T)}\|\partial_{t}^{\alpha, \lambda(t)}u(\textbf{x}, t)\|_{m},$ $\|u\|_{L^{\infty}(H^{m})} = \mathrm{sup}_{t\in (0, T)}\|u(\textbf{x}, t)\|_{m},$ $\|\partial_{t}^{2}u\|_{L^{\infty}(H^{m})} = \mathrm{sup}_{t\in (0, T)}\|\partial_{t}^{2}u(\textbf{x}, t)\|_{m},$ and $\|\partial_{t}^{2}u\|_{L^{\infty}(L^{2})} = \mathrm{sup}_{t\in (0, T)}\|\partial_{t}^{2}u(\textbf{x}, t)\|.$

Proof. Subtracting $e_{N}^{k} = u_{N}^{k}-\pi_{N}u(t_{k})$ on both sides of Eq (3.4), we obtain

Rewriting the above formula, we obtain

where

and

For term $\left(\varepsilon_{1}^{k+1}, v_{N}\right)_{N}$, we have

where $I_{d}$ is the identity operator. Using (3.3), we have

By using of (3.2), we obtain

Note that

the facts $(\nabla \pi_{N}u(t_{k}), \nabla v_{N})_{N} = (\nabla \pi_{N}u(t_{k}), \nabla v_{N})$ and $\left(\partial_{t}^{\alpha, \lambda(t)}u(t_{k}), v_{N}\right) = -\left(\nabla u(t_{k}), \nabla v_{N}\right)$, we obtain

Furthermore, we obtain

Combining the estimates of $\left(\varepsilon_{1}^{k}, v_{N}\right)_{N}$ and $\left(\varepsilon_{2}^{k}, v_{N}\right)_{N}$, we have

Moreover, with the same technique used in Theorem 2, using inequality (3.6), we can prove that

Hence, applying triangle inequality, we obtain

For $\alpha\rightarrow 1$, similar to the case discussed in Theorem 2, we obtain

Applying triangle inequality once again, we have

which implies the desired result. □

4.   Fast time discretization

4.1.   Construction of the fast difference operators

Applying the same technique given in reference ^[13], we use the sum-of-exponentials (SOE) approximation to approximate the power function $t^{-\alpha}(0 < \alpha < 1)$ in fractional derivative (2.1). And for the approach error, we have the following results:

Lemma 3. ^[13] For a given absolute error $\varepsilon$, the kernel functions $\frac{1}{t^{\alpha}}(\alpha \in (0, 1))$, there exist positive real numbers $\{\eta_{i}^{\alpha}\}_{i = 1}^{N_{exp}}$ and $\{w_{i}^{\alpha}\}_{i = 1}^{N_{exp}}$ such that

where $N_{exp} = O\left(\log\frac{1}{\varepsilon}\left(\log\log\frac{1}{\varepsilon}+\log\frac{T}{\tau}\right)+\log\frac{1}{\tau}\left(\log\log\frac{1}{\varepsilon}+\log\frac{1}{\tau}\right)\right).$

For the details of the proof of Lemma 3, see ^[13]. Following the same idea given in references ^[13,16], we split the Caputo fractional derivative (1.3) into a local term $L(t_{k})$ and a history term $H(t_{k})$ in discretized fractional derivative (2.1) with kernel $\lambda(t) = e^{-bt}(b\ge 0)$, i.e.,

For the local term, we discrete it as the follows

The history term, using SOE approximation, we obtain

where $\widetilde{\eta}_{i}^{\alpha} = \eta_{i}^{\alpha}+b$. For simplicity, we rewrite $H(t_{k})$ as

with

Using linear interpolation operators to approximate $u^{\prime}(s)$, we have

where $\hat{a}_{i} = \frac{1-e^{-(\eta_{i}^{\alpha}+b)}\tau}{\tau(\eta_{i}^{\alpha}+b)e^{(\eta_{i}^{\alpha}+b)\tau}}$, $\hat{b}_{i} = -\hat{a}_{i}$. Hence, history term has recurrence relation

Collecting the above formulas, we obtain the fast finite difference operator of fractional derivative (1.3)

where $U_{h, i}^{\alpha}(t_{1}) = 0$.

4.2.   The fast semi-discrete problem

If we use the difference operator (4.3) to approximate fractional derivatives (1.3) at time level $t_{k}$, we get the semi-discrete scheme

The weak form in $L^{2}$-inner product gives: Find $u^{k}\in H_{0}^{1}(\Omega)(1\le k\le N_{T})$, such that

For the sake of simplification, (2.2) can be rewritten as

where $\rho_{k-j}^{k, \alpha} = \frac{c_{j-1}-c_{j}}{c_{0}}$ as $1\le j\le k-1$, $\rho_{k-j}^{k, \alpha} = \frac{c_{k-1}}{c_{0}}$ as $j = k$. By comparing the two operators $F_{t}^{\alpha, \lambda(t)}$ and $L_{t}^{\alpha, \lambda(t)}$, we have

where $\Pi_{1, j}u(s) = u(t_{j-1})\frac{t_{j}-s}{\tau}+u(t_{j})\frac{s-t_{j-1}}{\tau}$, $s \in [t_{j-1}, t_{j}]$,

and

Simple calculations, we can check that

For the truncation error of the fast difference operator $F_{t}^{\alpha, \lambda(t)}$ gives

Lemma 4. Suppose that $u(t)\in C^{2}[0, T]$, for any $0 < \alpha < 1$, there holds

Proof. By the triangle inequality

□

Furthermore, we can rewrite the fast difference operator as

where

The coefficients $\rho_{k-j}^{k, \alpha}$ have the following properties.

Lemma 5. For $k\geq 1$, the coefficients $\rho_{k-j}^{k, \alpha}$ and $c_{0}^{-1}$ satisfy

(1) $\rho_{k-j}^{k, \alpha} > 0, j = 1, 2\ldots, k-1$; (2) $\sum\limits_{j = 1}^{k}\rho_{k-j}^{k, \alpha} = 1, j = 1, 2\ldots, k-1$;

(3) $\rho_{0}^{k, \alpha} > \frac{2(1-\alpha)\lambda_{k-1/2}}{3k^{\alpha}}$;    (4) $0 < \rho_{k-1}^{k, \alpha} < 1$;   (5) $1 < c_{0}^{-1} < e^{\frac{b\tau}{2}}$.

Proof. See the Appendix A. □

To do the error estimates of the fast difference scheme, we follow the technique used in reference ^[24]. Introducing a parameter $\sigma$, then $u^{k}-\sum_{j = 1}^{k}\rho_{k-j}^{k, \alpha}u^{k-j}$ can be rewritten as following

where $\sigma = \frac{1}{2}\rho_{k-1}^{k, \alpha}$, $\sigma^{j}$ means power of $\sigma$, and

Similarly, the sum-of-exponentials part can be rewritten as

where $\bar{\zeta}_{k-j}^{k, \alpha} = \sum_{l = 0}^{j}\sigma^{j-l}{\zeta}_{k-l}^{k, \alpha}$, $j = 0, 1, \ldots, k.$ The above statements show that the fast difference operator can be rewritten as

Thus, the semi-discrete scheme (4.4) can be rewritten as

For the weighted coefficients $\bar{\rho}_{k-j}^{k, \alpha}$ and $\bar{\zeta}_{k-j}^{k, \alpha}$, there holds

Lemma 6. For $k\geq 2$, the coefficients $\bar{\rho}_{k-j}^{k, \alpha}$ and $\bar{\zeta}_{k-j}^{k, \alpha}$ satisfy

(1) $\bar{\rho}_{k-j}^{k, \alpha} > 0$, $j = 1, 2\ldots, k$; (2) $\frac{1}{\bar{\rho}_{0}^{k, \alpha}} < \frac{1}{\rho_{0}^{k, \alpha}} < \frac{3k^{\alpha}}{2(1-\alpha)\lambda_{k-1/2}}$;

(3) $\sum\limits_{j = 1}^{k}\bar{\rho}_{k-j}^{k, \alpha}\le 1-\frac{2\sigma(1-\alpha)\lambda_{k-1/2}}{3(1-\sigma)k^{\alpha}}$; (4) $|\bar{\zeta}_{k-j}^{k, \alpha}|\le 4\varepsilon\tau^{\alpha}$, $j = 1, 2\ldots, k$; (5) $0 < \sigma < \frac{1}{2}$.

Proof. See the Appendix B. □

Theorem 4. Assume that the permissible error $\varepsilon$ satisfies $\varepsilon \le \min\left\{\frac{1}{4}\sqrt{\frac{\Gamma(1-\alpha)\sigma\lambda(T)}{3\hat{c}(1-\sigma)T^{\alpha}N_{T}e^{b\tau/2}}}, \frac{\Gamma(1-\alpha)}{4\hat{c}}\right\}$, then the semi-discrete scheme (4.4) is stable, and its solution satisfies

Proof. For $k = 1$, taking $v = u^{1}$ in (4.4), there holds

Using Schwarz inequality, we arrive at $\|u^{1}\|^{2}_1\le \|u^{1}\|\|u^{0}\|+\|f^{1}\|\|u^{1}\|,$ which produces $\|u^{1}\|_1\le \|u^{0}\|+\|f^{1}\|.$ Using (5) of Lemma 5, we have

For $k\geq2$, using mathematical induction on the index $k$, we will prove the following result:

Firstly, substituting $v = 2\bar{u}^{k}$ into semi-discrete scheme (4.5), it yields

Using the identity $2(\nabla u^{k}, \nabla \bar{u}^{k}) = \|\nabla u^{k}\|^{2}+\|\nabla \bar{u}^{k}\|^{2}-\sigma^{2}\|\nabla u^{k-1}\|^{2}$ and (3) of Lemma 6, we have

Furthermore, using Poincaré inequality ^[26]

and (5) of Lemma 5, if $\varepsilon \le \frac{\Gamma(1-\alpha)}{4\hat{c}}$, we obtain

where we used the relation $\bar{\zeta}_{k}^{k, \alpha} = \zeta_{k}^{k, \alpha} = \tilde{b}_{0}\leq \tau^{\alpha} \varepsilon$ in above derivation. Setting $\theta_{k} = 1-\sum\limits_{j = 1}^{k}\bar{\rho}_{k-j}^{k, \alpha}$ and using Young inequality, we obtain

According to (3) and (4) of Lemma 6, we obtain $\theta_{k} = 1-\sum\limits_{j = 1}^{k}\bar{\rho}_{k-j}^{k, \alpha}\ge \frac{2\sigma (1-\alpha)\lambda_{k-1/2}}{3(1-\sigma)k^{\alpha}}$. Moreover, with the help of Poincaré inequality, we obtain

If $\varepsilon\le \frac{1}{4}\sqrt{\frac{\Gamma(1-\alpha)\sigma\lambda(T)}{3\hat{c}(1-\sigma)T^{\alpha}N_{T}e^{b\tau/2}}}$, there holds

Rewriting (4.8) as follows:

In view of relation $\sigma+\theta_{2}+\bar{\rho}_{0}^{2, \alpha} = 1$, we obtain (4.7) for $k = 2$. Assuming that (4.7) is already correct for all $k = 2, 3, \ldots, n-1$, we can check that

which implies

Furthermore, using (4.10) and Poincaré inequality (4.9), we obtain

Combining (4.10) and (4.11), we obtain the desired results. □

For the convergence of time semi-discrete scheme (4.4), we have

Theorem 5. Let $u(\textbf{x}, t)$ be the exact solution of problem (1.1), and $\{u^{k}\}_{k = 1}^{N_{T}}$ be the semi-discrete solution of scheme (4.4) with the initial $u^{0} = u(0)$. Under the assumption given in Theorem 4 and $\partial_{t}^{2}u(\textbf{x}, t)\in L^{\infty}((0, T]; L^{2}(\Omega))$, there holds

where $\tilde{c}_{\alpha, T}$ depends only on $\alpha$ and $T$.

Proof. Let $e^{k} = u(t_{k})-u^{k}$. Subtracting (4.4) from problem (1.1) at time level $t_{k}$, it follows that

With the similar technique in Theorem 4, we obtain

Applying Lemma 4, we obtain (4.12). □

4.3.   Fast difference/Legendre spectral collocation

We consider the spectral collocation method in space as follows: Find $u_{N}^{k}\in V_{N}^{0}(\Omega)$, such that

Theorem 6. Let $u(\textbf{x}, t)$ be the exact solution of (1.1), and $u_{N}^{k}$ be the solution of scheme (4.13) with the initial condition $u_{N}^{0} = \pi_{N}u^{0}$. Suppose $\partial_{t}^{2}u(\textbf{x}, t)\in L^{\infty}((0, T]; H^{m}(\Omega))(m\ge 1)$, for $k = 2, \ldots, N_{T}$, there holds

Proof. Let $e_{N}^{k} = u_{N}^{k}-\pi_{N}u(t_{k})$, using (4.13), we have

where

Using equality $\bar{F}_{t}^{\alpha, \lambda(t)}\bar{u}^{k} = F_{t}^{\alpha, \lambda(t)}u^{k}$, scheme (4.15) can be rewritten as

where $\bar{e}_{N}^{k} = e_{N}^{k}-\sigma e_{N}^{k-1}$. Moreover, using Lemma 2, we have

According to Lemma 4, for $l = 0, 1$, we have

Combing (3.2), Eq (4.17) produces

According to the definition of the $\|\cdot\|_{0, N}$ and projection operator $\pi_{N}$, we have

Applying (1.1), (4.16), and (4.18), we obtain

Furthermore, using Lemmas 2 and 4, we obtain

Combining the estimates of $|\varepsilon_{1}^{k}(v_{N})|$, $|\varepsilon_{2}^{k}(v_{N})|$ and the Poincaré inequality, we have

With the similar argument given in Theorem 4, we obtain

By the triangle inequality $\|u(t_{k})-u_{N}^{k}\|_{1, N}\leq \|u(t_{k})-\pi_{N}u(t_{k})\|_{1, N}+\|u_{N}^{k}-\pi_{N}u(t_{k})\|_{1, N}$, we finally obtain our conclusion. □

5.   Numerical experiments

In this section, we will test our theory results given in previous sections. The error in time direction at terminal time $T = 1$ is measured by the point-wise maximum norm. In the following numerical experiments, we used two examples with the exponential kernel function $\lambda(t) = e^{-t}$, which modify the numerical examples given in reference ^[5]. In our numerical experiments, all the algorithms are implemented by MATLAB R2022b, which were run on a 2.90 GHz PC with 32 GB of RAM and the Windows 11 operating system.

Example 1. First, we consider problem (1.1) with exact solution $u(x, y, t) = (1+(6-(6+6t+3t^{2}+t^{3})e^{-t})\sin(\pi x)\sin(\pi y).$ The source term and initial values are given by the exact solution.

In Table 1, we list the errors and convergence orders of the present numerical scheme (3.4). Here, we fix the degree of polynomial $N_{x} = N_{y} = 20$. It can be observed that the time direction convergence order is $O(\tau^{2-\alpha})$, which is consistent with our theoretical analysis. To test the convergence order in spatial direction, we set $\tau = 1/1000$ for $\alpha = 0.3, 0.5$, and $0.9$. The relationship between the errors and polynomial degree $N$ is shown in the semi-log graph in Figure 2. We find that the numerical solution has the exponential accuracy.

Example 2. To test the effectiveness of the fast difference/spectral collocation scheme (4.13), we consider problem (1.1) with exact solution $u(x, y, t) = (3-(2+2t+t^{2})e^{-t})\sin(\pi x)\sin(2\pi y),$ and the source term and initial values are calculated by the exact solution.

The numerical results of schemes (3.4) and (4.13) are listed in Tables 2 and 3. In these tables, we compared the errors, convergence orders, and CPU times of the two schemes in Tables 2 and 3 for $\alpha = 0.5$ and $\alpha = 0.7$, respectively. In the implementation of fast difference/spectral collocation, we take $\varepsilon = 1e-9$, and it can be seen that the time convergence order of fast solver is $O(\tau^{2-\alpha})$, which is consistent with our theoretical results. In this example, the CPU(s) are measured, and the total times of numerical errors for the time steps $\tau$ vary in $[1/40, 1/80, 1/160, 1/320, 1/640, 1/1280]$. We also observe that the fast $L1$ scheme/spectral collocation has lower complexity and achieves the same accuracy as the $L1$ scheme ones.

6.   Conclusions

In this paper, the two-dimensional fractional diffusion equation with generalized memory kernel is analyzed and approximated in time and space. We used the spectral collocation method in the spatial direction, and the $L1$ formula and the fast $L1$ formula are used in time to approximate the fractional derivative. Finally, we obtained estimates of the temporal and spatial errors of the present schemes, the time and space error estimates are $O(\tau^ {2-\alpha}+O(N^{1-m})$. By comparing the $L1$ scheme with the fast $L1$ scheme, we found that the fast $L1$ scheme can significantly reduce storage costs and computational costs. There were also many improvements in this article, such as the use of graded grids instead of uniform grids (e.g., Stynes et al. ^[11]) for time grids when the initial singularity of solution is considered. Our future work will focus on this issue with the help of the existing techniques developed in references ^[10,27,28].

Author contributions

Zunyuan Hu: Methodology, software, writing-original draft, writing-review & editing; Can Li: Conceptualization, methodology, writing-original draft, writing-review & editing; Shimin Guo: Methodology, writing-review & editing. All authors have read and approved the final version of the manuscript for publication.

Use of Generative-AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

In this research, Can Li was partially supported by the Science Basic Research Plan in Shaanxi Province of China under Grant No. 2023-JC-YB-045. Shimin Guo was partially supported by the National Natural Science Foundation of China under Grant No. 12271427.

Conflict of interest

It is declared that none of the authors have any competing interests in this manuscript.

Appendix A.   Proof of Lemma 5

(1) According to the definitions of $\rho_{k-j}^{k, \alpha}$ and $c_{0}$, we obtain

By simple calculation, and using Lemma 1, we have

(2) After calculating term by term, we obtain

(3) Using Lemma 1 once again, we have

(4) Since (1) provides that $\rho_{k-1}^{k, \alpha} > 0$, then $\rho_{k-1}^{k, \alpha} = (c_{0}-c_{1})c_{0}^{-1} = 1-\frac{c_{1}}{c_{0}} < 1.$

(5) Note that

where $Z = \frac{1}{2-\alpha}-\frac{1}{2}$, $0\le Z\le \frac{1}{2}$. We may rewrite $c_{0}$ as follows

By checking the partial derivatives of $g$ with respect to variable $Z$, we obtain $\frac{\partial g(\tau, Z)}{\partial Z} = 1-e^{-b\tau}\ge 0,$ which means $g(\tau, Z)$ is monotone increasing function with respect to variable $Z$, thus we get $c_{0}\le e^{-\frac{1}{2}b\tau}\left[1-\frac{1}{2}e^{-\frac{1}{2}b\tau}\right]+\frac{1}{2},$ and $c_{0}\geq e^{-\frac{1}{2}b\tau}$. Furthermore, we define $\tilde{g}(\tau) = e^{-\frac{1}{2}b\tau}\left[1-\frac{1}{2}e^{-\frac{1}{2}b\tau}\right]+\frac{1}{2}$ and take the derivative of $\tilde{g}(\tau)$, we obtain

Hence, $\tilde{g}(\tau)\le 1$, $0 < c_{0} = g(\tau, Z)\le \tilde{g}(\tau)\le 1$, it provides that $c_{0}^{-1} = \frac{1}{g(\tau, Z)}\ge \frac{1}{\tilde{g}(\tau)}\ge 1.$ Finally we get $1\leq c_{0}^{-1}\leq e^{\frac{1}{2}b\tau}$.

Appendix B.   Proof of Lemma 6

(1) From Lemma 5, for $j = 1, 2, ..., k$, we have $\rho_{k-j}^{k, \alpha} > 0$, thus $\sigma = \frac{1}{2}\rho_{k-1}^{k, \alpha} > 0.$ By using $\bar{\rho}_{k-j}^{k, \alpha} = \sum_{l = 1}^{j}\sigma^{j-l}\bar{\rho}_{k-l}^{k, \alpha}-\sigma^{j}$, we further have

(2) According to Lemma 5, we have $\frac{1}{\rho_{0}^{k, \alpha}} < \frac{3k^{\alpha}}{2(1-\alpha)\lambda_{k-1/2}}$. Moreover, note that $\rho_{0}^{k, \alpha} < \bar{\rho}_{0}^{k, \alpha}$, we have $\frac{1}{\bar{\rho}_{0}^{k, \alpha}} < \frac{1}{\rho_{0}^{k, \alpha}} < \frac{3k^{\alpha}}{2(1-\alpha)\lambda_{k-1/2}}$.

(3) Note that $\frac{1}{\bar{\rho}_{0}^{k, \alpha}} < \frac{3k^{\alpha}}{2(1-\alpha)\lambda_{k-1/2}}$, we may check that

(4) In view of $|\zeta_{k-j}^{k, \alpha}|\le 2\tau^{\alpha}\varepsilon$, we obtain $|\bar{\zeta}_{k-j}^{k, \alpha}| = |\sum_{l = 0}^{j}\sigma^{j-l}\zeta_{k-j}^{k, \alpha}|\le 2\tau^{\alpha}\varepsilon\sum_{l = 0}^{j}\sigma^{l}\le 4\tau^{\alpha}\varepsilon.$

(5) From Lemma 5, we have $\sigma = \frac{1}{2}\rho_{k-1}^{k, \alpha} < \frac{1}{2}$.