Pointwise-in-time $\alpha$-robust error estimate of the ADI difference scheme for three-dimensional fractional subdiffusion equations with variable coefficients

1.   Introduction

This paper focuses on the 3D subdiffusion equations with variable coefficients

in which: $a_{1}$, $a_{2}$, $a_{3}$, and $b$ are positive constants; $f$, $\varphi$, and $g$ are smooth functions; and $\Lambda = (0, {\Upsilon}_{x})\times(0, {\Upsilon}_{y})\times(0, {\Upsilon}_{z})$ is a rectangular field bounded by $\partial\Lambda$. The Caputo fractional derivative can be defined via

The fractional subdiffusion equation with variable coefficients is derived from the classical partial differential equation. Many scholars have studied the existence and uniqueness of the solution of this equation ^[1,2]. The high-dimensional subdiffusion problem can accurately describe the mechanical and physical processes with non-locality and spatial global correlation, it has become an important tool to describe a variety of complex mechanical and physical behaviors ^[3,4,5], and many numerical schemes have been proposed in ^{[6,7,8,9,10,11,12,13,14,15]}. Among them, the ADI method is favored by many scholars. Tian and Ge ^[16] offered the fourth-order ADI difference method for 2D unsteady convection diffusion problems. Wang et al. ^[17,18] considered the 2D fractional integro-differential equation by using an ADI compact difference method. Zhang et al. ^[19] provided the 2D semilinear multidelay parabolic equations by using Crank-Nicolson ADI compact difference schemes. Zhang et al. ^[20] considered the 2D space-fractional nonlinear Ginzburg-Landau equation via linearized ADI difference schemes. Wang et al. ^[21,22] presented an analysis of the convergence of the $\alpha$-robust $H^1$-norm for the ADI differece scheme for the 2D time fractional diffusion equation (TFDE). Chen et al. ^[23,24] considered the ADI difference scheme for a class of two-dimensional multi-term TFDEs and gave the point error estimation. Qiao et al. ^[25,26,27] presented the ADI difference scheme for 2D TFDE with a weakly singular kernel. Zhang et al. ^[28] presented an efficient ADI difference scheme for the nonlocal evolution problem in 3D space. Zhou et al. ^[29,30] solved the 3D TFDE via the ADI method.

Along with the previous research, many scholars considered the subdiffusion equation with variable coefficients. Ngondiep ^[31] discussed a class of convection-diffusion-reaction equations with variable coefficients and gave a two-level fourth-order scheme. Wei et al. ^[32] gave the anisotropic finite element method (FEM) for 2D time fractional variable coefficient diffusion equations on graded meshes. Zeng and Tan ^[33] proposed a two-grid FEM for nonlinear TFDEs with variable coefficients. Ma et al. ^[34] presented the L1 robust analysis of the fourth-order block-centered FDM for 2D TFDEs with variable coefficient. As we all know, the study of the variable coefficient subdiffusion equation is almost based on one or two dimensions. It is not easy to solve three-dimensional fractional subdiffusion equations with variable coefficients due to the existence of the weak singularity near the initial time. In this paper, we put forward an $\alpha$-robust L1 ADI difference scheme on graded mesh to overcome these difficulties. The L1 scheme on graded meshes is used to discretize the Caputo fractional derivative, and the difference method is used to approximate the spatial derivative. Using the modified Grönwall inequality, we prove the stability and $\alpha$-robust convergence of the proposed method.

In this article, the main achievements are as follows:

● We construct an ADI difference method for the 3D subdiffusion equation with variable coefficients. Most previous work has focused on 1D or 2D problems.

● We introduce a new norm to prove the stability and convergence of the proposed method.

● The $\alpha$-robust convergence is proved strictly, which is able to avoid the blow-up phenomenon in the case of $\alpha\rightarrow 1^{-}$.

The structure of this paper is as follows. In Section 2, the ADI scheme is developed for the 3D variable coefficient subdiffusion problems. In Section 3, we demonstrate the stability and convergence of the method. In Section 4, the three-dimensional numerical example is given to verify the effectiveness and accuracy of the proposed method. We conclude the paper in Section 5.

2.   Derivation of ADI scheme

Set $h_{x} = \frac{{\Upsilon}_{x}}{\ell_{1}}$, $h_{y} = \frac{{\Upsilon}_{y}}{\ell_{2}}$, $h_{z} = \frac{{\Upsilon}_{z}}{\ell_{3}}$ for the positive integers $\ell_{1}$, $\ell_{2}$, and $\ell_{3}$. Let $N\geq1$ be a positive integer number and set $t_{n} = T(n/N)^{\gamma}$ for $n = 0, 1, \cdots, N$, $\gamma\geq 1$. Denote time step $\iota_{n} = t_{n}-t_{n-1}$ for $1\leq n\leq N$.

Set $\overline{\Lambda}_{h} = \{(x_{i}, y_{r}, z_{k})|0\leq i\leq \ell_{1}, 0\leq r\leq \ell_{2}, 0\leq k\leq \ell_{3}\},$ and $\Lambda_{h} = \overline{\Lambda}_{h}\cap\Lambda,$ $\partial\Lambda_{h} = \Lambda_{h}\cap\partial\Lambda$. Let $u^{n}_{irk}$ be the approximation solution of (1.1)-(1.3).

For $n = 1, 2, \cdots, N$, by the L1 scheme ^[35,36], denote

where

Thus, the problem (1.1)-(1.3) can be approximated by

where

Set $S_{n} = d^{-1}_{n, 1} = \iota^{\alpha}_{n}\Gamma(2-\alpha)$, and to establish the ADI scheme, add a small term

onto the left side of Equation (2.2), which is order $O(N^{-2\alpha})$. Then, one has

Multiplying both sides of (2.6) by $S_{n}$, we get

Then, rewriting Eqs. (2.6)-(2.8) in an ADI form:

in which

Set $u^{n-\frac{1}{3}}_{irk} = (I-\frac{a_{3}S_{n}}{1+bS_{n}}\delta _{z}^{2})u_{irk}^{n}$ and $u^{n-\frac{2}{3}}_{irk} = (I-\frac{a_{2}S_{n}}{1+bS_{n}}\delta _{y}^{2})u_{irk}^{n}$. To get the solution $\{u^{n}_{irk}\}$ of problem (2.6)-(2.8), we just need to solve three of the one-dimensional equations.

First, for fixed $r\in \{1, 2, \cdots, \ell_{2}-1\}$ and $k\in \{1, 2, \cdots, \ell_{3}-1\}$, we can solve $\{u^{n-\frac{2}{3}}_{irk}\}$ in the $x$-direction as follows:

Then, for fixed $i\in \{1, 2, \cdots, \ell_{1}-1\}$ and $k\in \{1, 2, \cdots, \ell_{3}-1\}$, we can solve $\{u^{n-\frac{1}{3}}_{irk}\}$ in the $y$-direction as follows:

At last, for $i\in \{1, 2, \cdots, \ell_{1}-1\}$ and $r\in \{1, 2, \cdots, \ell_{2}-1\}$, we utilize $\{u^{n-\frac{2}{3}}_{irk}, u^{n-\frac{1}{3}}_{irk}\}$ in the $z$-direction to solve

where the solution $\{u^{n}_{ijk}\}$ is attained.

3.   Analysis of stability and convergence

For grid function $u = \{u^{n}_{irk}|0\leq i\leq \ell_{1}, 0\leq r\leq \ell_{2}, 0\leq k\leq \ell_{3}, 0\leq n\leq N\}$, we define

where $\delta_{x}u_{i-\frac{1}{2}, r, k} = \frac{1}{h_{x}}(u_{irk}-u_{i-1, r, k})$. The norms $\Vert \delta_{y}u^{n}\Vert$, $\Vert \delta_{z}u^{n}\Vert$, $\Vert \delta_{y}\delta_{z}u^{n}\Vert$, $\Vert \delta_{x}\delta_{z}u^{n}\Vert$, $\Vert \delta_{x}\delta_{z}\delta^{2}_{y}u^{n}\Vert$ and $\Vert \delta_{y}\delta_{z}\delta^{2}_{x}u^{n}\Vert$ can be defined similarly.

Set $\overline{U} = \{u_{irk}|u_{irk} = 0 \text{ if } (x_{i}, y_{r}, z_{k})\in\partial\Lambda_{h} \text{ and } (x_{i}, y_{r}, z_{k})\in \overline{\Lambda}_{h}\}$, for $\forall u$, $w\in \overline{U}$, we define

3.1.   Stability analysis

This section mainly discusses the stability analysis of ADI difference scheme (2.6)-(2.8).

Lemma 3.1. For grid functions $u, \omega\in \overline{U}$, one has

where the inequality is equal when $\omega = u$.

Proof. Using the Cauchy-Schwartz inequality, one has

The proof is completed. □

For $n = 1, 2, \cdots, N$ and $j = 1, 2, \cdots, n-1$, define a positive sequence

Lemma 3.2. ^[37,38] Assume that the sequences $\{m_{n}\}^{\infty}_{n = 1}$, $\{g_{n}\}^{\infty}_{n = 1}$ are nonnegative, and assume that the grid function $\{u^{n}:n = 0, 1, \cdots, N\}$ satisfies $u_{0}\geq 0$ and

Then

Lemma 3.3. ^[39] When $n = 1, 2, \cdots, N$, one has

Next, the stability theorem is presented.

Theorem 3.4. For $n = 1, 2, \cdots, N$, the solution ${u^{n}_{irk}}$ of ADI scheme (2.6)-(2.8) satisfies

Proof. Multiplying both sides of $-h_{x}h_{y}h_{z}(L_{1}\delta _{x}^{2}+L_{2}\delta _{y}^{2}+L_{3}\delta _{z}^{2})u_{irk}^{n}$, and summing over $i, r, k$ for $(x_{i}, y_{r}, z_{k})\in \Lambda_{h}$, we have

Then, using the Cauchy-Schwartz inequality, Young's inequality, and Lemma 3.1, one has

and

Thus, we get

which obtains

Finally, Lemma 3.2 and Lemma 3.3 are applied in (3.3), and the proof is completed. □

3.2.   $H^{1}$-norm convergence analysis

Lemma 3.5. ^[24] Assume that $|\partial^{p}_{t}u(x, y, z, t)|\leq C(1+t^{\alpha-1})$, $p = 0, 1, 2$. Then, for $(x_{i}, y_{r}, z_{k}, t_{n})\in \Lambda$, we have

Lemma 3.6. ^[40] For grid function $v\in\nu_{h}$, one has $\parallel v\parallel\leq\frac{Ł_{1}Ł_{2}}{\sqrt{6(Ł_{1}^{2}+Ł_{2}^{2})}}\parallel \nabla_{h}v\parallel.$

Define

and

According to the definition of the small term (2.5), one has

Using a Taylor expansion, one has

Applying Lemma 3.5 and Eq. (3.8), one has

We now show the $H^{1}$-norm convergence analysis.

Theorem 3.7. Suppose that $u(x, y, z, \cdot)\in C^{4, 4, 4}(\Lambda\cap\partial\Lambda)$, $|\partial^{p}_{t}u(x, y, z, t)|\leq C(1+t^{\alpha-1})$, $p = 0, 1, 2$. one has

where $C$ is an $\alpha$-robust constant that does not blow up as $\alpha\rightarrow1^{-}$.

Proof. Subtracting Eqs. (1.1)-(1.3) from Eqs. (2.6)-(2.8), the error equations are presented via

Multiplying both sides of (3.11) by $-h_{x}h_{y}h_{z}(L_{1}\delta _{x}^{2}+L_{2}\delta _{y}^{2}+L_{3}\delta _{z}^{2})e_{irk}^{n}$ and summing over $i, r, k$ for $(x_{i}, y_{r}, z_{k})\in \Lambda_{h}$, one has

Since

and

Then, we get

That is,

Thus

By the Poincaré inequality, Eq. (3.9), and Lemma 3.5, we have

From (3.6), (3.14), and (3.16), using Lemma 3.2 and Lemma 3.3, and noting that $\Vert e^{0}\Vert_{A} = 0$, one has

Finally, using Lemma 3.6 and the definition of $\Vert e^{n}\Vert_{A}$, the proof is completed. □

4.   Numerical experiment

Example 4.1. For problem (1.1)-(1.3), considering ${\Upsilon}_{x} = {\Upsilon}_{y} = {\Upsilon}_{z} = \pi$, $\ell_{1} = \ell_{2} = \ell_{3} = \ell$, $T = 1$, $\varphi(x, y, z, t) = 0$, $g(x, y, z) = 0$, and

The exact solution is $u(x, y, z, t) = t^{\alpha}\sin x\sin y\sin z$. Set $a_{1} = 1$, $a_{2} = 0.1$, $a_{3} = 0.5$, and $b = 1$. In Tables 1–3, we fix $\ell = 128$ to test the error and convergence in time. In Table 1, choosing different $\alpha$, we give errors, convergence orders, and CPU time in time for $\gamma = \frac{2-\alpha}{\alpha}$. In Table 2, set $a_{1} = 1$, $a_{2} = 0.1$, $a_{3} = 0.5$, and $b = 100$, and choosing different $\alpha$ we give errors, convergence orders and CPU time in time for $\gamma = \frac{2-\alpha}{\alpha}$. In Table 2, the convergence order goes down for $N = 32$ and $\alpha = 0.2$, which can be caused by the round-off error computing the coefficients of the L1 scheme (2.1), see Remark 2.1 in ^[41]. In Table 3, set $a_{1} = 1$, $a_{2} = 0.1$, $a_{3} = 0.5$, and $b = 0.001$, and choosing different $\alpha$ we give errors, convergence orders, and CPU time in time for $\gamma = \frac{2-\alpha}{\alpha}$. Due to the influence of multiple variable coefficient values such as $a_{1}$, $a_{2}$, $a_{3}$, $b$, the degree of freedom of the equation fitting is higher, and there may be some differences for CPU times. From Tables 1–3, it is obtained that the order of convergence in time is $(N^{-min\{2\alpha, 2-\alpha, \gamma\alpha\}})$, which is consistent with the theoretical analysis.

In Tables 4–6, we fix $N = 512$ to present errors and convergence in space for $\gamma = (2-\alpha)/\alpha$. In Table 4, set $a_{1} = 1$, $a_{2} = 0.1$, $a_{3} = 0.5$, and $b = 1$, and by choosing different $\alpha$ we give $H^{1}$-norm errors, convergence orders and CPU time in space. In Table 5, choosing different $\alpha$ we give $H^{1}$-norm errors, convergence orders and CPU time in space for $a_{1} = 1$, $a_{2} = 0.1$, $a_{3} = 0.5$, and $b = 100$. In Table 6, choosing different $\alpha$ we give $H^{1}$-norm errors, convergence orders, and CPU time in space for the case where $a_{1} = 1$, $a_{2} = 0.1$, $a_{3} = 0.5$, and $b = 0.001$. In Tables 4-6, we can find that the spatial convergence order is order 2, which is consistent with our theory. From Tables 1-6 we can notice that our scheme is computationally efficient and consistent with the analysis.

5.   Conclusions

In this paper, an ADI scheme is proposed to solve the 3D variable coefficient subdiffusion problem, and our theoretical analysis shows that the method is unconditionally stable, its temporal convergence order is order $\min\{2-\alpha, \gamma\alpha, 2\alpha\}$, and its spatial convergence order is order 2. We present numerical results when taking different coefficients, which show that our ADI scheme is consistent with the theoretical analysis and is very efficient in solving such problems. The ADI technique proposed can reduce the computational cost from spatial discretization. In the future, we will consider some fast and parallel numerical methods ^[42,43,44] to improve computational efficiency in time.

Acknowledgments

The work was supported by National Natural Science Foundation of China Mathematics Tianyuan Foundation (12226337, 12226340, 12126321, 12126307), Scientific Research Fund of Hunan Provincial Education Department (21B0550, 22C0323, 23C0193), Hunan Provincial Natural Science Foundation of China (2022JJ50083, 2023JJ50164).

Use of AI tools declaration

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

Conflict of interest

The authors declare there is no conflict of interest.