Stability and pointwise-in-time convergence analysis of a finite difference scheme for a 2D nonlinear multi-term subdiffusion equation

1.   Introduction

The following two-dimensional nonlinear multi-term time-fractional subdiffusion equations with initial and boundary conditions are considered in this paper:

where $\Omega = (0, L_1)\times (0, L_2)$ with the boundary $\partial \Omega$, $f \; \in C^1(\mathbb{R})$, and $g$ is a continuous function in $\Omega$. The operator $D_{t}^{\boldsymbol{\alpha} }$ in (1.1) is defined by

where $D_{t}^{{\alpha }_{l}}$ denotes the Caputo derivative with respect to $t$, i.e.,

In order to describe the physical process better, some scholars use multi-term time fractional differential equations for modeling, such as the behavior of viscoelastic fluids^[1,2], the dispersion of pollutants^[3], and magnetic resonance imaging^[4]. The finite difference method^[5,6,7,8], the finite volume method ^[9,10], the fractional predictor-corrector method^[11], the finite element method^[12], the collocation method ^[13,14], and the spectral method^[15] have been developed for solving the multi-term time fractional equations. The work ^[16] proposed a fast linearized finite difference method for the nonlinear multi-term time-fractional wave equation using SOE approximation. Recently, a lot of work about the time fractional nonlinear subdiffusion equation has been published. For example, the paper ^[17] considered a Newton linearized Galerkin finite element method to solve the problem with non-smooth solutions in the time direction and provided an optimal error estimate by using the discrete fractional Gronwall-type inequality on the graded meshes. The paper ^[18] provided an efficient systematic framework for solving nonlinear fractional partial differential equations on unbounded domains and obtained an error estimate. Jiang et al. ^[19] proposed an efficient ADI scheme for the nonlinear subdiffusion equation with a weakly singular solution and obtained the pointwise-in-time error estimate. Li et al. ^[20] proposed a new tool, the refined discrete fractional-type Gr$\ddot{\text{o}}$nwall inequality, which is used to derive a sharp pointwise-in-time error estimate of the L1 scheme for the problem. However, as far as the authors know, at present there is no work dedicated to the pointwise-in-time error estimate of the L1 scheme for the nonlinear multi-term time-fractional subdiffusion equation with weakly singular solutions. The present work is designed to fill this gap. Following ^[21,22], in the remainder of our paper, we make the following assumption:

Assumption 1.1. For all $(x, y, t)\in \bar{\Omega}\times (0, T]$, we assume the solution satisfies

where $C$ is a positive constant. $C$ in this paper represents a constant independent of time and space step, and $C$ in different positions represents different values.

This paper consists of the following sections. In Section 2, we construct a fully discrete scheme for the problem (1.1)–(1.3). In Section 3, some lemmas are introduced that will be used in the subsequent analysis. In Section 4, convergence and stability of the fully discrete scheme proposed in Section 2 are given, and we obtain the sharp pointwise-in-time error estimate. In Section 5, the theoretical analysis is verified by four numerical examples.

2.   Fully discrete scheme

Let $M_1$, $M_2$, and $N$ be three positive integers. Divide space and time uniformly into $M_1 \times M_2$ and $N$ parts, respectively. Let $\{t_n|t_n = n\tau, \ 0\le n\le N\}$ be a uniform partition of $\left[0, T\right]$ with the time step $\tau = T/N$. Let $h_1 = L_1/M_1$ and $h_2 = L_2/M_2$ be the spatial steps. So spatial grid points consist of $x_i = ih_1$ and $y_j = jh_2$, where $i = 0, 1, \cdots, M_1$ and $j = 0, 1, \cdots, M_2$. The spatial grid is represented by $\Omega' = \{(x_i, y_j)|i = 0, 1, \cdots, M_1, j = 0, 1, \cdots, M_2\}$. Let $\Omega_h = \Omega'\cap\Omega$, and $\partial\Omega_h = \Omega'\cap\partial\Omega$. Defining $u_{i, j}^{n} = u(x_i, y_j, t_n)$, the previous equations (1.1)–(1.3) at the grid point $\left(x_i, y_j, t_n\right)$ can be transformed into

where the notations in (2.1) will be given afterwards.

Applying an L1 formula for the multi-term time Caputo derivatives, i.e., set

where $d_{k}^{{{\alpha }_{l}}} = \frac{{{\tau }^{-{{\alpha }_{l}}}}[{{k}^{1-{{\alpha }_{l}}}}-{{(k-1)}^{1-{{\alpha }_{l}}}}]}{\Gamma (2-{{\alpha }_{l}})}$. For simplicity, defining ${{d}_{k}} = \sum\limits_{l = 1}^{J}{{{b}_{l}}d_{k}^{{{\alpha }_{l}}}}$, it yields that

A standard second-order approximation is used to discretize $\Delta u_{i, j}^n$:

To approximate the nonlinear term $f(u_{i, j}^{n})$, we utilize the backward formula, expressed as

In (2.1), $({{r}_{l}})_{i, j}^{n}$ and $R_{i, j}^{n}$ are truncation errors, i.e.,

Leaving out the truncation errors in (2.1) and replacing $u_{i, j}^n$ by $U_{i, j}^n$, the fully discrete scheme is obtained:

3.   Preliminary lemmas

Lemma 3.1. For $(x_i, y_j)\in\Omega_h$, $1\le n\le N$, it holds that

Proof. Due to the continuity of $f'$ and the boundedness of $u$, $f'(\zeta)$ is bounded, where $\zeta$ is between $u_{i, j}^{n}$ and $u_{i, j}^{n-1}$. We obtain

Then we estimate the truncation error in two cases (^[19], page 6). Note the condition (1.5).

Case A: $n\ge 2$

where we have used $t_{n-1}\simeq t_{n}$ for $n \geq 2$ in the penultimate inequality.

Case B: $n = 1$

The proof is completed. □

Lemma 3.2. For $(x_i, y_j)\in\Omega_h$, $1\le n\le N$, we have $\left|R_{i, j}^n\right|\le C(h_1^2+h_2^2)$.

Proof. According to the Taylor expansion, we can see that

where ${{x}_{i-1}} < \xi_1 < {{x}_{i+1}}$, and ${{y}_{j-1}} < \xi_2 < {{y}_{j+1}}$. □

Lemma 3.3. [23, Lemma 6] $\left| \sum\limits_{l = 1}^{J}{{{b}_{l}}({{r}_{l}})_{i, j}^{n}} \right|\le C{{n}^{-\min \{2-{{\alpha }_{1}}, {{\alpha }_{1}}+1\}}}$, for $(x_i, y_j)\in\Omega_h$, $1\le n\le N$.

Define the positive stability multipliers ${\theta }_{0} = \mu, \ {{\theta }_{n}} = \mu \sum\limits_{k = 1}^{n}{(d_k-d_{k+1}){{\theta }_{n-k}}}$, where $\mu = d_{1}^{-1}$.

Lemma 3.4. (^[24], Lemma 5.1) ${\theta }_{n}$ is a monotonically decreasing sequence with respect to $n$. For $n = 0, 1, \cdots, N$,

Notation. $\theta _{n}$ here has the following relationship with $\sigma_n^{\boldsymbol{\alpha}}$ in Literature ^[24]: $\theta _{n} = \sigma_n^{\boldsymbol{\alpha}}$.

Lemma 3.5. (^[24], Lemma 5.2, Corollary 5.1) For $n = 1, 2, \cdots, N,$

where $\beta\ge 0$ and

If there exists $\underline{\alpha}$, which satisfies $0 < \underline{\alpha} \leqslant \alpha_1 < 1$, we have

where constant $\hat{C}$ depends on $\underline{\alpha}$ and $T$.

Lemma 3.6. (^[25], Lemma 5.1) If $\{y_n\}$ is a non-negative sequence, and it satisfies $y_n\le a_1t_n^{-\eta_1}+a_2t_n^{-\eta_2}+b\tau\sum\limits_{j = 1}^{n-1}{t_{n-j}^{\alpha-1}y_j}, \ 1\le n\le N$, we get $y_n\le C(a_1t_n^{-\eta_1}+a_2t_n^{-\eta_2}), \ 1\le n\le N$, where $0\le \eta_1, \eta_2 < 1$, $a_1, a_2, b > 0$, $0 < \alpha < 1$ and $N$ is a positive integer.

Suppose ${g^j}$ is an arbitrary mesh function; we define an integral operator $B_{t}^{\boldsymbol{\alpha} }$, which satisfies $\ B_{t}^{\boldsymbol{\alpha} }({{g}^{0}}) = 0, \ B_{t}^{\boldsymbol{\alpha} }({{g}^{n}}) = \sum\limits_{j = 1}^{n}{{{\theta }_{n-j}}{{g}^{j}}}$, for $n = 1, 2, \cdots, N$.

Lemma 3.7. (^[23], Lemma 2) For any mesh function $\{U^j\}_{j = 0}^N$, the following formula holds.

4.   Convergence and stability analysis

Theorem 4.1. If $U_{i, j}^n$ is the solution of (2.5)–(2.7) at point $(x_i, y_j, t_n)$, and suppose $\left|U_{i_0, j_0}^{n}\right| = {{\left\| U^{n} \right\|}_{\infty }}$, it holds that

Proof. From (2.5), it is easy to see that $U_{i, j}^{n}$ satisfies

Based on (2.7) and the fact that $d_{k-1}^{\alpha_l} > d_k^{\alpha_l}$, we have

Then, we obtain

which, according to the definition of $D_{\tau }^{\boldsymbol{\alpha} }$, is equivalent to the following expression:

By applying the integral operator to both sides of the above formula, we obtain

Due to Lemma 3.7 and the above definition of $B_{t}^{\boldsymbol{\alpha} }$, (4.1) holds for $n\ge 1$. Obviously, (4.1) also holds for $n = 0$. □

Theorem 4.2. Suppose that $u$ is the solution of (1.1)–(1.3) and satisfies Assumption 1.1, and that $U$ is the solution of (2.5)–(2.7). There exist positive constants $\tau_0$ and $h_0$. When $\tau < \tau_0$ and $h_1, h_2 < h_0$, it holds for $m > 0$ that

Proof. For $m = 0$, (4.3) holds obviously. Assuming that (4.3) holds for $m = 0, 1, 2\cdots n-1$($n\ge 1$), then for sufficiently small $\tau$, $h_1$, $h_2$, and $1\leq k \leq n$, we obtain

Let us discuss whether the inequality (4.3) holds when $m = n$. Let $e_{i, j}^m = u_{i, j}^m-U_{i, j}^m$ for $0\leq m\leq n$. Subtracting (2.5)–(2.7) from (2.1)–(2.3), we obtain

where $({{R}_{f}})_{i, j}^{m} = f(u_{i, j}^{m-1})-f(U_{i, j}^{m-1})$. Considering the continuity of $f'$, the boundedness of $u_{i, j}^{m-1}$ and $U_{i, j}^{m-1}$, we have

Let $\left|e_{{{i}_{0}}, {{j}_{0}}}^{m}\right| = {{\left\| {{e}^{m}} \right\|}_{\infty }}$. Using Lemmas 3.1–3.3 and inequality (4.4), similar to (4.2), it can be obtained that for $1\leq m\leq n$,

Applying the definition of $B_{t }^{\boldsymbol{\alpha}}$ again, and Lemmas 3.4 and 3.5, we can further obtain

where $\sigma_1 = \min \{2-{{\alpha }_{1}}, {{\alpha }_{1}}+1\}$. From Lemma 3.6, we have

In summary, the inequality (4.3) holds for $m = n$, which finishes the mathematical induction. The proof is complete. □

Remark 4.1. The global maximum error of the numerical solution is

When $t_n$ is away from 0, the local maximum error is

Theorem 4.3. The fully discrete scheme (2.5)–(2.7) is stable with respect to the initial value. If $\hat{U}_{i, j}^n$ satisfies the following equations,

and ${{\left\| g- \hat{g}\right\|}_{\infty }}$ is sufficiently small, then

where $\hat{e}_{i, j}^n = U_{i, j}^n-\hat{U}_{i, j}^n$.

Proof. Subtracting (4.7)–(4.9) from (2.5)–(2.7) and according to Theorem 4.1, we have

where $\hat{e}_{i_0, j_0}^n = {{\left\| \hat{e}^n \right\|}_{\infty }}$. (4.10) holds for $n = 0$ obviously. Suppose that (4.10) holds when $n = 0, 1, 2, \cdots, m-1(m\ge 1)$, we have ${{\left\| \hat{U}^r \right\|}_{\infty }}\le {{\left\| \hat{e}^r \right\|}_{\infty }}+{{\left\| U^r \right\|}_{\infty }}\le C{{\left\| \hat{e}^0 \right\|}_{\infty }}+{{\left\| U^r \right\|}_{\infty }}\le 1+{{\left\| U^r \right\|}_{\infty }}$ on the condition that ${{\left\| \hat{e}^0 \right\|}_{\infty }}$ is small, for $r\le m-1$. According to (4.3), $U_{i, j}^r$ is bounded. Then, combining the continuity of $f'$ and the boundedness of $\hat{U}_{i, j}^r$, it holds that

Due to Lemma 3.6, we have

Therefore, the mathematical induction ends, and the proof is complete. □

5.   Numerical experiments

Example 5.1. We first consider the following two-term time-fractional nonlinear subdiffusion equation with $b_1 = b_2 = 1$.

where $h(x, y, t)$ and $g(x, y)$ are up to the exact solution. We set the exact solution to $t^{\alpha_1}sin(\pi x)sin(\pi y)$, which satisfies Assumption 1.1. We consider the spatial domain $\Omega = (0, 1)\times (0, 1)$ and set $T = 1$.

Above all, we compute the convergence order in spatial direction. We set $N = 1000$ so that the influence of errors in temporal direction can be ignored compared with errors in spatial direction. The maximum errors at $t_n = 1$ and rates, when $\alpha_1 = 0.4$ and $\alpha_2 = 0.3$, are presented in Table 5.1. Numerical results show that the spatial accuracy is $O(h_1^2+h_2^2)$. For studying temporal convergence rates, we define global errors $E_G$ and local errors $E_L$ by

Global errors and rates in Table 5.2 show that the global temporal convergence order is $\alpha_1$. When $t_n$ is far away from 0, results are shown in Table 5.3, and we get the local temporal accuracy $O(\tau)$. All in all, numerical results are consistent with the theoretical analysis in Remark 1.

Example 5.2. Secondly, we consider a two-dimensional three-term time-fractional nonlinear subdiffusion equation with $b_1 = b_2 = b_3 = 1$.

where $\Omega = (0, 1)\times (0, 1)$ and $T = 1$. We calculate $h(x, y, t)$ and $g(x, y)$ based on the exact solution

Numerical results are shown in Tables 5.4 and 5.5, which verify the theoretical analysis as well.

Example 5.3. Thirdly, we investigate the scenario in which the nonlinear term is represented by $f(u) = sin(u)$, which satisfies Lipschitz condition. The corresponding equation is formulated as follows:

where $\Omega = (0, 1)\times (0, 1)$ and $T = 1$. Similarly, we derive $h(x, y, t)$ and $g(x, y)$ based on the exact solution

The corresponding numerical results are presented in Tables 5.6 and 5.7. The global convergence order is $\alpha_1$, and the local convergence order is 1 in the temporal direction.

Example 5.4. Finally, we consider a two-dimensional two-term time-fractional nonlinear subdiffusion equation with $b_1 = b_2 = 1$, whose exact solution is unknown.

where $\Omega = (0, 1)\times (0, 1)$ and $T = 1$.

The two-mesh method^[26] is applied to compute errors and convergence rates. We take $M_1 = M_2 = 60$. $E_L$ is redefined by

in which $W^n$ is the numerical solution of Example 5.4 with $\tau = T/2N$. The local errors are shown in Table 5.8. The local temporal convergence rate $O(\tau)$ is consistent with Remark 1.

Use of AI tools declaration

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

Acknowledgments

The research is supported in part by the Natural Science Foundation of Shandong Province under Grant ZR2023MA077, Fundamental Research Funds for the Central Universities (No. 202264006), and the National Natural Science Foundation of China under Grant 11801026.

Conflict of interest

The authors declare there are no conflicts of interest.