An efficient numerical method based on QSC for multi-term variable-order time fractional mobile-immobile diffusion equation with Neumann boundary condition

1.   Introduction

In recent years, the fractional partial differential equations (FPDEs) have been widely used to simulate various phenomena of anomalous diffusion, such as Brownian motion in fractal theory ^[1], and solute transport in porous media ^[2]. Samko and Ross ^[3] proposed an important kind of fractional operators, which are dependent on time and space. Recent studies showed that variable-order fractional derivatives can describe some complicated phenomena more precisely ^[4,5], such as the viscoelasticity oscillators ^[6] and the motion of particles in special circumstances ^[7]. The multi-term variable-order time fractional mobile-immobile diffusion (TF-MID) equation offers a more precise and realistic depiction of solute diffusion in heterogeneous porous media compared to its single-term counterpart, as referenced in ^[8,9].

Different kinds of numerical methods have been proposed for solving PDEs, such as the finite difference method ^[10,11,12] and finite element method ^[13,14,15]. The quadratic spline collocation (QSC) method stands out as an efficient discretization approach. The basis functions of the quadratic spline space are very smooth, and the smooth conditions at interpolation points are helpful to reduce the number of unknowns, which means the resulting algebraic equation has small scale. As a result, QSC has been widely adopted for solving both integer-order problems ^[16,17] and fractional-order differential equations ^[18,19].

Many efficient numerical schemes have been proposed to solve time fractional PDEs, including the classical $L1$ scheme ^[20,21] and $L2$-$1_{\sigma}$ scheme ^[22]. She et al. ^[23] introduced a transformed $L1$ method to solve the multi-term time-fractional initial-boundary value problem. Yuan et al. ^[24] introduced linearized transformed $L1$ Galerkin finite element method for nonlinear time fractional Schr$\ddot{{\rm{o}}}$dinger equations. Chen et al. ^[25] studied a Grunwald-Letnikov scheme on uniform mesh for reaction sub-diffusion equations with weakly singular solutions and presented a sharp error estimate. Zhou et al. ^[26] proposed the nonuniform Alikhanov schemes on the nonuniform meshes for nonlinear time fractional parabolic equations. Ren et al. ^[27] established sharp $H1$ norm error estimates for the $L1$ formula and a fractional Crank-Nicolson scheme on nonuniform meshes for reaction sub-diffusion problems. However, for variable-order models, since the fractional derivative has a variable-dependent kernel, it is more difficult to construct suitable numerical approximations. Zheng and Wang ^[28] presented an $L1$ scheme for variable-order time-fractional diffusion equations. Du et al. ^[29] developed an $L2$-$1_{\sigma}$ formula for the variable-order time-fractional wave equation. Zhang et al. ^[30] considered an implicit numerical method to solve space-time variable-order fractional advection-diffusion equations. Liu et al. ^[31] proposed the regularity of the solution to a variable-order time-fractional diffusion equation, and constructed a fully discrete numerical scheme. Building on this foundation, we ^[32] expand the application of the $L1$ scheme by developing a first-order numerical scheme for the variable-order TF-MID with variable coefficients. Ji et al. ^[33] proposed the $L1^+$ scheme of constant-order FPDEs, which can be derived by performing the classical $L1$ scheme to the integral form of the FPDEs. The $L1^+$ scheme has almost the same computational cost as the classical $L1$ scheme, while improving the convergence order $(2-\alpha)$ of the classical $L1$ scheme to second order, if the solution is sufficiently well-defined. In this paper, we will investigate the $L1^+$ formula for multi-term variable-order time fractional derivatives.

Furthermore, the fast implementation of numerical methods for time-fractional differential equations is another focus of researchers. Jiang et al. ^[34] introduced the sum-of-exponentials (SOE) technique, to greatly reduce the computational cost for the evaluation of the constant-order time fractional derivatives. Subsequently, Zhang et al. ^[35] employed the exponential sum approximation (ESA) method, adeptly handling the singular kernels of variable-order Caputo fractional derivatives. Building on this foundation, they developed a second-order fast evaluation method in ^[36].

In this paper, we combine the temporal $L1^+$ formula with the spatial QSC method to construct a QSC-$L1^+$ scheme, to solve two-dimensional multi-term variable-order TF-MID equations. The Neumann boundary conditions are incorporated into the framework of the energy method, which leads to a novel technique for the unconditionally stability and convergence of the QSC-$L1^+$ scheme. In addition, we conduct a comprehensive improvement and optimization of the ESA technique to fully meet the $L1^+$ formula. Such a fast evaluation demonstrates exceptional performance by significantly reducing the computational cost and memory requirements, particularly when combined with carefully selected parameters, which further highlights its notable advantages.

The structure of the paper is as follows. In Section 2, we introduce the multi-term variable-order TF-MID equations and propose the QSC-$L1^+$ scheme. Subsequently, we meticulously construct an energy method aimed at rigorously demonstrating the unconditional stability and convergence of the QSC-$L1^+$ scheme in Section 3. In Section 4, we harness the ESA technique to achieve fast computations along the temporal direction. Section 5 presents detailed results from numerical experiments, which not only substantiate our theoretical findings but also highlight the efficiency and practicality of the proposed scheme. In Section 6, we provide a concise summary of the entire work. It is worth noting that, we use $C_{i}$ in this paper to denote constants which are independent of mesh sizes.

2.   Equation and discretization

In this section, we consider a kind of two-dimensional multi-term variable-order TF-MID equation as follows, and the equation can be used to understand and simulate diffusion behavior of solutes in heterogeneous porous media ^[9],

defined in a rectangular spatial domain $\Omega = (x_L, x_R)\times(y_L, y_R)$, and the temporal interval is denoted by $[0, T]$. At the initial time, the solution satisfies

and at the boundary $\partial \Omega$, the solution satisfies the Neumann condition

Here, the positive constant $\kappa$ signifies the diffusion coefficient and $f$ denotes a predefined smooth source function. Additionally, $\mathcal{L}$ stands as a spatial elliptic operator, which stands for:

$\alpha^{m}(t)$ for $m = 1, 2, \cdots, P$ are the variable time fractional orders which satisfy

For the solute transport problem in porous media, a small portion of solute transport follows the Fickian diffusion process, and it can be expressed as the term $u_{t}(x, y, t)$, while the majority of solute particle transport follows anomalous diffusion, which can be described by the variable-order Caputo fractional operator ${}_{0}^{C}D_{t}^{\alpha(t)}u(x, y, t)$. In order to distinguish anomalous diffusion with different speeds, we use multi-fractional orders ${}_{0}^{C}D_{t}^{\alpha^{m}(t)}u(x, y, t)$ in the modeling. The Caputo fractional operator ${}_{0}^{C}D_{t}^{\alpha^{m}(t)} u(x, y, t)$ is rigorously defined as

where the weight function $\omega_{\beta}(t)$ is specified by:

To characterize the initial singularity of the solution, one can introduce the weighted Banach space $C_{\mu}^{m}((0, T]; \mathcal{X})$, which incorporates with the norm $\|\cdot \|_{\mathcal{X}}$, where $m \geq 2$ and $0 \leq \mu < 1$, as defined in ^[37]. The space is defined as:

with the norm given by:

The eigenfunctions $\left\{\varphi_{i}\right\}_{i = 1}^{\infty}$ of the Sturm-Liouville problem, which satisfy the equations

constitute an orthogonal basis in the $L^{2}(\Omega)$ space. The eigenvalues $\left\{\lambda_{i}\right\}_{i = 1}^{\infty}$ are strictly positive and non-decreasing, approaching $\infty$ as $i$ increases. By harnessing the theory of sectorial operators, we define the fractional Sobolev space

equipped with the norm $\|v\|_{\breve{H}^{\gamma}}: = \left(\|v\|_{L^{2}}^{2}+|v|_{\breve{H}^{\gamma}}^{2}\right)^{1/2}$. Moreover, $\breve{H}^{\gamma}(\Omega)$ is a subset of the fractional Sobolev space ${H}^{\gamma}(\Omega)$, distinguished by the characteristics outlined in ^[38]:

for all $(x, y) \in \partial \Omega$, where $s < \gamma / 2$.

Lemma 2.1 (^[8,31]). If condition (2.4) holds, suppose that $u_{0} \in \breve{H}^{\gamma+6}$ for $\gamma > 1/2$, $f \in H^{1}\left([0, T]; \breve{H}^{s+4}\right) \cap H^{2}\left([0, T]; \breve{H}^{s+2}\right) \cap H^{3}\left([0, T]; \breve{H}^{s}\right)$ for $0 \leq s\leq \gamma$ and $\alpha \in C^{2}[0, T]$. If $\alpha^{m}(0) > 0$, for $m = 1, 2, \cdots, P$, we have $u \in C^{3}\left((0, T]; \breve{H}^{\gamma}(0, L)\right) \cap C_{1-\alpha^{m}(0)}^{3}\left((0, T]; \breve{H}^{\gamma}(0, L)\right)$ and

Next, we will conduct a comprehensive exploration of the $L1^{+}$ scheme applied to the temporal domain, accompanied by the analysis of the QSC discretization in the spatial domain.

2.1.   Temporal $L1^+$ scheme

Given a positive integer $N$, we divide the time interval $[0, T]$ uniformly as $0 = t_0 < t_1 < ... < t_N = T$. Therefore, we can conclude that $t_n = n\tau$, where $\tau$ is defined as $\tau = \frac{T}{N}$. For a given continuous function $v(t)$, we refer to its piecewise linear interpolation as $\Pi v(t)$. To quantify the discrepancy between the original function and its interpolation, we introduce the interpolation error $\theta_{e} v\left(t\right) = v\left(t\right)-\Pi v\left(t\right)$, which can be expressed as

Based on Lemma 2.1, we have $\left\|\frac{\partial^{2} v}{\partial t^{2}}\right\|_{\mathcal{X}} \leq C_{2}t^{-\alpha^{m}(0)}$. Therefore, we can verify that

We represent $v^{n}$ as the numerical approximation of $v(t)$ at the specific time instant $t = t_{n}$, and define $\mathfrak{T} = \big\{v^{n}, \ n = 0, 1, \cdots, N\big\}$ as a finite-dimensional function space that spans over the temporal grid. For convenience, we further define

Then, averaging the integral of ${}_{0}^{C}D_{t}^{\alpha^{m}(t)} v(t)$ over $[t_{n-1}, t_{n}]$, and approximating $\alpha^{m}(t)$ at the midpoint of this subinterval, we have

where $\tilde{\alpha}_{n}^{m} : = \alpha_{n-\frac{1}{2}}^{m}$, for $n = 1, 2, \cdots, \; N$ and $m = 1, 2, \cdots, P$. Based on the trapezoidal formula and Lemma 2.1, we can verify that

where $C_{4}$ is the upper bound for $v(t)$, and $C_{5}$, $C_{6}$ are positive constants related to $C_{4}$. According to the boundedness of $\Gamma(x)$ and applying Taylor's expansion for two numerators above, we have

where $C_{i}$ for $i = 7, 8, \cdots, 12$ are the upper bounds of some functions. Since the integration of linear terms are equal to zero, we can obtain $\left|r_{1, n, m}\right| = \mathcal{O}\big(\tau^{2}\big)$. Furthermore, the first term on the righthand side of (2.6) can be approximated by the $L1$ formula as

According to the truncation error analysis in ^[33,39], we have $r_{2, n, m} = \mathcal{O}\Big(\tau^{2}t_{n}^{-\tilde{\alpha}_{n}^m-\alpha^{m}(0)}\Big)$. Substituting Eq (2.7) into (2.6), we have the approximation for $n = 1, 2, \cdots, \; N$ that

where

The discretization (2.8) for the variable-order Caputo fractional derivative ${}_{0}^{C}D_{t}^{\alpha^{m}(t)} v(t)$, with the coefficients (2.9), is called the $L1^+$ formula.

2.2.   Spatial QSC method

For given positive integer $M_x$, we divide spatial domain $[x_L, x_R]$ uniformly as

with mesh size $\Delta x = \frac{x_R-x_L}{M_x}$. Then, we define the quadratic spline space with the variable $x$ as

where ${\mathbf P}_{2}(\cdot)$ represents the set of piecewise quadratic polynomials with the variable $x$. Similarly, we can define the quadratic spline space $\mathcal{V}_{y}$ for the variable $y$. Furthermore, we denote by $\triangle : = \triangle_x \times \triangle_y$ the mesh partition of the spatial domain, and denote by $\mathcal{V}: = \mathcal{V}_{x} \otimes \mathcal{V}_{y}$ the space of piecewise biquadratic polynomials.

Next, we consider the basis functions of the space $\mathcal{V}$. We let

and define the quadratic B-splines

as the basis function of $\mathcal{V}_{x}$. Similarly, we get $\big\{\phi_{j}(y), \ j = 0, \cdots, M_{y}+1 \big\}$ as the basis functions for $\mathcal{V}_{y}$. Thus, the quadratic spline solution $u_{h}^{n} \in \mathcal{V}$ of Eq (2.1) is represented as

where $c_{i, j}^{n}$ are the coefficients to be solved.

In order to solve these unknowns, we choose the centers of the rectangular mesh $\triangle$ as the collocation points, denoted by

where $\big\{\xi_{i}^{x} = \frac{1}{2}\left(x_{i-1}+x_{i}\right), \ i = 1, 2, \cdots, M_{x}\big\}$ and $\big\{\xi_{j}^{y} = \frac{1}{2}(y_{j-1}+y_{j}), \ j = 1, 2, \cdots, M_{y}\big\}$. We denote $\xi_{0}^{x} = x_{L}$, $\xi_{M_{x}+1}^{x} = x_{R}$ and $\xi_{0}^{y} = y_{L}$, $\xi_{M_{y}+1}^{y} = y_{R}$, and denote $\partial{\xi} = \big\{ \big(x_{L}, \xi _{j}^{y}\big), \ j = 0, 1, \cdots, M_{y}+1 \big\} \cup \big\{\big(x_{R}, \xi _{j}^{y}\big), \ j = 0, 1, \cdots, M_{y}+1 \big\} \cup \big\{\big(\xi _{i}^{x}, y_{L}\big), \ i = 0, 1, \cdots, M_{x}+1 \big\} \cup \big\{\big(\xi _{i}^{x}, y_{R}\big), \ i = 0, 1, \cdots, M_{x}+1 \big\}$ as the boundary collocation points. Consequently, we select $\bar{\xi} = \xi \cup \partial\xi$ as the set encompassing both interior and boundary collocation points. For simplicity, we further define index sets: $\Lambda = \big\{ (i, j), \ (\xi _{i}^{x}, \xi _{j}^{y}) \in \xi \big\}$, $\partial\Lambda = \big\{ (i, j), \ (\xi _{i}^{x}, \xi _{j}^{y}) \in \partial \xi \big\}$, and $\bar{\Lambda} = \Lambda \cup \partial\Lambda$.

Taking the collocation points into expressions (2.11), we get

and the Neumann boundary condition

where operators $\boldsymbol{\theta}_{x}$, $\boldsymbol{\varsigma}_{x}$, and $\boldsymbol{\eta}_{x}$ are defined as

We further define

which leads to

In addition, the operators $\boldsymbol{\theta}_{y}$, $\boldsymbol{\eta}_{y}$, and $\boldsymbol{\vartheta}_{y}$ are defined along the $y$ direction similarly. We need to notice that $\boldsymbol{\varsigma}_{x}$ and $\boldsymbol{\varsigma}_{y}$ are discretizations for the Neumann boundary conditions, and their definitions are compatible with $\boldsymbol{\vartheta}_{x}$ and $\boldsymbol{\vartheta}_{y}$, respectively. In this paper, we keep both definitions without confusion. Next, we will delve into the full discretization scheme, which is built upon the spatial and temporal discretizations, to construct a robust numerical method for Eqs (2.1)–(2.3).

2.3.   QSC-L1$^+$ scheme

We consider Eq (2.1) on the time subinterval $[t_{n-1}, t_{n}]$, and take the integral average to obtain

Through calculations, we can validate the first term of Eq (2.15),

Regarding the first term at the righthand side of Eq (2.15),

where

and based on estimations (2.5), it satisfies

Similarly, for the last term of Eq (2.15), we have

where $r_{4, n}$ satisfies

with $C_{14}$ as an upper bound for $f_{tt}(x, y, t)$.

Utilizing Eqs (2.16) and (2.17), in conjunction with the discretization method outlined in (2.8), Eq (2.15) can be reformulated as

with the truncation errors

Subsequently, we will proceed to approximate the solution of Eq (2.1) with the QSC method. Specifically, we substitute $u_{h}^{n}(x, y)$ of the form given by (2.11) into Eq (2.18), neglect the truncation errors, and obtain

Taking the collocation points $(\xi_{i}^{x}, \xi_{j}^{y})$ into Eq (2.20), we get the QSC-$L1^+$ scheme,

For the initial condition (2.2), we have

and on the boundary collocation points, we have

Remark 2.1. If the equation includes additional convection terms, such as $u_x + u_y$, the corresponding QSC-$L1^+$ scheme (2.21) will also be appended an additional term $(\boldsymbol{\vartheta}_x \boldsymbol{\theta}_y + \boldsymbol{\theta}_x\boldsymbol{\vartheta}_y)c_{i, j}^{n-\frac{1}{2}}$, which leads to a sparse matrix-vector multiplication. Therefore, the QSC-$L1^+$ scheme can be applied for TF-MID equations with convection terms, with additional little computational cost.

Next, we will delve into a comprehensive analysis of the stability and convergence properties of the proposed numerical scheme.

3.   Numerical analysis of the QSC-$L1^+$ method

Before conducting the numerical analysis, it is imperative to establish some fundamental definitions concerning inner products and norms. Specifically, we define $\mathcal{M}_h = \big\{\boldsymbol w, \ \boldsymbol w = \{w_{i, j}, \ (i, j)\in {\Lambda}\} \big\}$ as the spatial grid function space and further introduce $\mathring{\mathcal{M}}_{h} = \{ \boldsymbol w\in\mathcal{M}_h, \ \boldsymbol{\varsigma}_{x}w_{i, j} = \boldsymbol{\varsigma}_{y} w_{i, j} = 0 \ \text{for}\ (i, j)\in \partial{\Lambda} \}$. For two functions $w, v \in \mathring{\mathcal{M}}_{h}$, we proceed to define the inner product,

which induce the associated discrete norms

We can reformulate scheme (2.21) as

To facilitate clarity and consistency in our notation, we uniformly define the coefficients for Eq (3.1) as follows:

With these new coefficients, the QSC-$L1^+$ method (2.21) can be rewritten in a more compact form as

We notice that the orders of magnitude of the coefficients can be estimated as $b_{1}^{(n)} = \mathcal{O}(1)$, $b_{2}^{(n)} = \mathcal{O}(\tau^{1-\tilde{\alpha}_{n}^{m}})$. Considering the coefficient properties outlined in ^[33,39], the newly defined set of coefficients $\big\{b_{n-k+1}^{(n)}, k = 1, 2, \cdots, n\big\}$ fulfills the properties stated in the following lemma.

Lemma 3.1. For sufficiently small values of $\tau$, the coefficients $b_{k}^{(n)}$ are monotonically decreasing for $k = 1, 2, \cdots, n$, i.e.,

where $n$ is a fixed time index.

3.1.   Basic preparations

Before the numerical analysis, we need to provide some basic lemmas which are necessary in the stability and convergence analysis. Here, we first introduce some lemmas on the coefficients of the QSC-$L1^+$ method.

Lemma 3.2. We assume that the fractional order ${\alpha^m}(t)$ satisfies $({\alpha^m}(t))^{\prime} \leq 0$, for $0\leq t\leq T$ and $m = 1, 2, \cdots, P$, then the coefficients $b_{n-k}^{(n)}$ in the QSC-$L1^+$ scheme (3.2) fulfill the estimation

Lemma 3.3. The coefficients of the QSC-$L1^+$ scheme have a positive lower bound, i.e.,

where $n\geq 2$.

Lemma 3.4. A special summation of the coefficients of the QSC-$L1^+$ scheme is bounded, i.e.,

The proof of above lemmas can be referred to ^[40]. Furthermore, we also require the subsequent lemmas concerning the operators that have been previously defined.

Lemma 3.5 (^[22]). If the coefficients $b_{k}^{(n)}$ are monotonically decreasing, for $k = 1, 2, \cdots, n$, then we have the estimate,

where $\boldsymbol c^k = \big\{ c_{i, j}^k, \ (i, j)\in \bar{\Lambda} \big\}$ are the coefficients in the quadratic spline approximation $u_h^{k}$.

Lemma 3.6. For function $\boldsymbol w \in \mathring{\mathcal{M}}_{h}$, we have

Proof. We prove the first estimation for simplicity. According to reference ^[40], the operator $\boldsymbol{\theta}_{x}$ can be decomposed as $\boldsymbol{\theta}_{x} = \boldsymbol{\zeta}_{x}^2$, where $\boldsymbol{\zeta}_{x}$ is also a spatial operator. Similarly, $\boldsymbol{\theta}_{y} = \boldsymbol{\zeta}_{y}^2$. Since $\left\|\boldsymbol\zeta_{x} \boldsymbol w\right\|^{2} = \langle\boldsymbol\zeta_{x} \boldsymbol w, \boldsymbol\zeta_{x} \boldsymbol w\rangle = \langle \boldsymbol\theta_{x} \boldsymbol w, \boldsymbol w \rangle$, we can get from the definition of $\boldsymbol\theta_{x}$ in (2.12) that

According to the homogeneous Neumann boundary conditions, we get $w_{0, j} = w_{1, j}$ and $w_{M_{x}+1, j} = w_{M_{x}, j}$. Then, we have

We utilize the inequality $2ab \leq a^2+b^2$ in equality (3.4) to get

Then, using the inequality $2ab \geq -a^2-b^2$ in equality (3.4), we can get

The second estimation can be obtained in the similar way.

Lemma 3.7. For any $\boldsymbol v^{n} \in \mathring{\mathcal{M}}_{h}$, $n = 1, 2, \cdots, N$, we have

Proof. Recalling the notation $\boldsymbol v^{n-\frac{1}{2}} = \frac{1}{2}\left(\boldsymbol v^n+\boldsymbol v^{n-1} \right)$, we have

We first present the estimation for $P_{1}$,

Applying the homogeneous Neumann boundary conditions $\boldsymbol{\varsigma}_{x}\boldsymbol{\theta}_{y} v_{0, j}^{n-1} = 0$ and $\boldsymbol{\varsigma}_{x}\boldsymbol{\theta}_{y} v_{M_{x}+1, j}^{n-1} = 0$, for $j = 1, 2, \cdots, M_y$, we rearrange the terms in $P_{1}$ and obtain

Likewise, we have

Furthermore, we use the Cauchy-Schwarz inequality and $-2ab \leq a^2 + b^2$ to obtain

The terms $P_3$ and $P_4$ have similar results. Thus, we can complete the proof.

Lemma 3.8 (^[41]). We suppose $v^{n}$, $w^{n}$ $\in \mathfrak{T}$ satisfy the inequality $v^{n} \leq (1 + \tau C_{18})v^{n-1} + \tau w^{n-1}$, with $n = 1, 2, \ldots, N$, then we can obtain

where $C_{18}$ is a positive constant.

Given the lemmas presented above, we will focus on the stability of the QSC-$L1^+$ scheme (2.21)–(2.23), or its equivalent form (3.2).

3.2.   Stability analysis

Theorem 3.1. We denote by $\boldsymbol c^n = \{c_{i, j}^n, \ (i, j)\in {\Lambda}, 0\leq n\leq N\}$ the coefficients of the approximation (2.11) solved by the QSC-$L1^+$ method. If the fractional orders $\alpha^{m} (t)$ are monotonically decreasing functions, then the following estimate holds:

Proof. We multiply Eq (3.2) by $2\Delta x\Delta y\boldsymbol\theta_{x}\boldsymbol\theta_{y}c_{i, j}^{n}$. Taking summation for the indices $i$ from $1$ to $M_{x}$ and for $j$ from $1$ to $M_{y}$ yield

According to Lemmas 3.5 and 3.7, we achieve the following estimation

Then, based on Lemma 3.2, we obtain

Denote

where $1\leq n \leq N$. With inequality (3.7), the inequality (3.6) can be simplified as

We utilize Lemma 3.8 to derive that

According to the definition of $\Big\{ b_{n-k+1}^{(n)}, \ k = 1, 2, \cdots, n \Big\}$, we have by Lemma 3.3 that $b_{1}^{(n)} = 1+\tau a_{1}^{(n)} > 1+C_{16}\tau$, and $b_{n-k+1}^{(n)} > C_{16}\tau$, for $k = 1, 2, \cdots, n-1$. Then, $G^{n}$ has the lower bound

We combine estimates (3.8) and (3.9) to conclude that for $1\leq n\leq N$,

Next, we will analyze the terms in (3.10) individually. First, based on Lemma 3.6, we obtain the estimates

By the inequality $2ab\leq\ 2\varepsilon a^2 +(1/2\varepsilon) b^2$,

Taking the estimates (3.11) – (3.12) into (3.10), we can derive the conclusion of the theorem.

3.3.   Convergence analysis

Recall that the quadratic spline interpolation $\mathcal I w(x, y)$ of the function $w(x, y)$ satisfies

and satisfies

on the left and right boundaries, as well as

on the other two boundaries. For $w(x, y)\in C^4\big(\bar{\Omega}\big)$, we denote a norm $\| w\|_c = \max_{(i, j)\in{\Lambda}}|w\big(\xi _{i}^{x}, \xi _{j}^{y}\big)|$, and the interpolation error satisfies ^[17,42]

We denote $\boldsymbol u^n = \big\{ u^n(\xi _{i}^{x}, \xi _{j}^{y}), \ (i, j)\in \bar{\Lambda} \big\}$ as the true solution, and $\boldsymbol u_{h}^n = \big\{ u_{h}^n(\xi _{i}^{x}, \xi _{j}^{y}), \ (i, j)\in \bar{\Lambda} \big\}$ as the numerical solution by the the QSC-$L1^+$ scheme. Then, the convergence of the QSC-$L1^+$ scheme can be derived as follows.

Theorem 3.2. If the fractional orders $\alpha^{m}(t)$ satisfy $0 < \alpha_{*} \leq\alpha^{m}(t)\leq \alpha^{*} < 1$, then the numerical solution of the QSC-$L1^+$ scheme converges and satisfies

Proof. With Eq (2.18), for any fixed $n$, we can construct the following difference equation on the interpolation $\mathcal I u^{n}(x, y)$:

where

with the definition of $R^n$ in (2.19). We take the collocation points $(\xi_{i}^{x}, \xi_{j}^{y})$ for $(i, j)\in {\Lambda}$ into (3.15)–(3.16), and obtain

Based on (2.19) and properties (3.14), $g^{n-\frac{1}{2}}$ can be estimated by

Since $\mathcal I u^n(x, y)\in \mathcal{V}$, we assume $\mathcal I u^n(x, y)$ can be written as

where $d_{i, j}^{n}$ are degrees of freedom (DOFs) of $\mathcal I u^{n}(x, y)$. Then, Eq (3.17) can be rewritten as

where $g_{i, j}^{n-\frac{1}{2}} = g^{n-\frac{1}{2}}\big(\xi _{i}^{x}, \xi _{j}^{y}\big)$, and the boundary conditions for $(i, j)\in \partial{\Lambda}$, $1\leq n \leq N$,

Denote $\boldsymbol e^n = \boldsymbol d^n - \boldsymbol c^n$, and we substitute (3.2) and (2.23) from (3.19) and (3.20) to obtain

and for $(i, j)\in \partial{\Lambda}, \quad 1\leq n \leq N$ such that

Applying estimation (3.10) in the proof of Theorem 3.1, together with $\boldsymbol e^0 = 0$, we can get

Defining $E^{n} = \max _{1 \leq l \leq n}\left\|\boldsymbol{\theta}_{x} \boldsymbol{\theta}_{y} \boldsymbol e^{l}\right\|$, for $n = 1, 2, \cdots, N$, we have from (3.21) that

Taking the maximum on the left-hand side for $l = 1, 2, \cdots, \; n$, we have

which leads to

Based on the estimation (3.18), we see

Next, we give further discussions on the values of $t_{n}$.

(Ⅰ) If $t_{n} \leq 1$, we have $t_{n}^{1-\alpha^{m}(0)} \leq t_{n}^{1-{\alpha}(0)}$, and $t_{k}^{-\tilde{\alpha}_k^{m}-\alpha^{m}(0)} \leq t_{k}^{-\alpha^{*}-\alpha(0)}$ for $k\leq n$, where $\tilde{\alpha}_n = \max_{1 \leq m \leq P} \tilde{\alpha}^{m}_n$, and $\alpha(0) = \max_{1 \leq m \leq P} \alpha^{m}(0)$. Therefore, we have

(Ⅱ) If $t_{n} > 1$, we have $t_{n}^{1-\alpha^{m}(0)} \leq t_{n}^{1-\bar{\alpha}(0)}$ for $\bar{\alpha}(0) = \min_{1 \leq m \leq P} \alpha^{m}(0)$. Meanwhile, we define an integer $\bar{k} = \lfloor 1/\tau \rfloor$, which leads to $t_{k} \leq 1$ for $k \leq \bar{k}$, and we can refer to the case (Ⅰ) above for the summation. For $k \geq \bar{k} + 1$, we have $t_{k} > 1$ and $t_{k}^{-\tilde{\alpha}_k^{m}-\alpha^{m}(0)} \leq t_{k}^{-\alpha_{*}-\bar{\alpha}(0)}$, then we have

Thus, we can substitute (3.23) and (3.24) into (3.22) to obtain

Since $\big(\mathcal I u - u \big)^n\left(\xi _{i}^{x}, \xi _{j}^{y}\right) = 0$, we can get

With the estimate provided in (3.25), we can obtain the convergence result.

Remark 3.1. If the fractional orders $\alpha^{m}(t)$ satisfy $\alpha^{*}+ \alpha(0) < 1$, one can obtain

Remark 3.2. The convergence of the numerical solution always depends on the well-posedness of the solution. If the fractional orders do not satisfy the restriction (2.4), the well-posedness of the solution in Lemma 2.1 will not hold anymore. Under some weak regularity of the solution, such as

for a positive parameter $\sigma$, then it can be obtained similarly that

Remark 3.3. If the solution satisfies $u \in C^2[0, T]$, the QSC-$L1^+$ scheme can achieve second temporal convergence order, which fits well with the example confirmation in Ref. ^[33].

4.   Acceleration techniques

The Caputo time fractional operator is a nonlocal operator, and its discretization typically involves the numerical approximations at all the historical time points, which can be prohibitively expensive, especially for complicated systems and long-term simulations. To mitigate this computational burden, we utilize the ESA technique introduced in ^[35] to reduce the computational cost during the implementation of the $L1^+$ formula. The primary task is to efficiently approximate the function $t^{-\tilde{\alpha}_{n }^m}$ in the integration (2.7) over the subinterval $[\tau, T]$, while the implementation on the subinterval $[0, \tau]$ does not require acceleration.

It is reported in ^[35] that for a small tolerance $\varepsilon > 0$, one can define

For $t \in\left[t_{n-1}, t_{n}\right]$, $s \in\left[0, t_{n-2}\right]$, such that $t-s\geq \tau$, and the exponential function in the definition of the variable order Caputo time fractional derivative can be approximated by

where

Now for any $v \in \mathfrak{T}$ and $1\leq m \leq P$, the $L1^+$ operator $\bar{\boldsymbol\delta}_{t}^{\tilde{\alpha}_{n}^m} v(t_{n-\frac{1}{2}})$ with $n\geq2$ defined in (2.8) can be decomposed as

The first term in (4.1) represents the historical summation, which can be written as

where

We can see that $V^{(2, r)} = 0$, and $V^{(n, r)}$ can be expressed as obtained recursively by

This means that $V^{(n, r)}$ can be recursively obtained with the value $V^{(n-1, r)}$ at the previous time instant. Taking expressions (4.2) into (4.1), we can get the fast evaluation method of the $L1^+$ formula for the variable order fractional operator,

Using the fast evaluation (4.3), we can get an accelerated scheme, named the QSC-F$L1^+$ scheme. Actually, when $n = 1, 2$, we still employ $(2.21)$ to compute Eq (2.1). For $n \geq 3$, we take the following QSC-F$L1^+$ scheme,

Compared with the computational cost $\mathcal{O}(M_{x}M_{y}N^2)$ for the QSC-$L1^+$ scheme (2.21), the cost for the QSC-F$L1^+$ scheme (4.4) is reduced to $\mathcal{O}(M_{x}M_{y}N\log^{2} N)$. Moreover, the memory requirement is also significantly decreased from $\mathcal{O}(M_{x}M_{y}N)$ to $\mathcal{O}(M_{x}M_{y}\log^{2} N)$.

5.   Numerical experiments

Example 5.1. We limit Eq (2.1) within the spatial domain $\Omega = (0, 1)\times (0, 1)$ and over the temporal interval $[0, T] = [0, 1]$, and consider various combinations of the subsequent variable fractional time orders,

We take $\kappa = 1$, $u^{0}(x, y) = \cos (\pi x) \cos (\pi y)$, and impose the homogeneous Neumann boundary conditions and a homogeneous source function.

Since the true solution is unknown, the numerical solution on a finer space-time mesh is chosen as the reference solution for comparison, then the error can be measured by

Now, we compute the convergence orders of the new proposed scheme. We first choose $M_{x} = M_{y} = 2^{12}$, which is small enough, and change the value of $N$ from $2^{9}$ to $2^{12}$. The observed errors near the initial time point and the corresponding orders of convergence in time are shown in Table 1. It can be seen that the temporal convergence orders fit the theoretical order $\mathcal{O}\big(\tau^{\min{\{3-\alpha^*-\alpha(0), \ 2\}}} \big)$ very well. The errors at the final point and corresponding orders are shown in Table 2. It seems that the initial weak singularity almost does not effect the convergence orders at the instances far away. Then, we choose $N = 2^{10}$ such that temporal mesh size is small enough, and change the value of $M_x$ and $M_y$ from $2^{6}$ to $2^{9}$. The observed errors and the spatial orders are shown in Table 3. We can see that the numerical results are consistent with the theoretical results.

Besides, in order to thoroughly investigate the effectiveness and efficiency of the ESA technique, we compare the observed errors and the running time for both the QSC-$L1^+$ scheme and the QSC-F$L1^+$ scheme in Table 4. We can clearly observe that the QSC-F$L1^+$ scheme requires significantly less running time compared to the QSC-$L1^+$ scheme, while achieving comparable levels of errors with the same discretization parameters. This finding exhibits the superiority of the ESA technique in terms of computational efficiency.

Example 5.2. We extend Eq (2.1) into three-dimensional spatial domain $\Omega = (0, 1) \times (0, 1) \times (0, 1)$, and still on the temporal interval $[0, T] = [0, 1]$, with the same variable fractional orders as Example 5.1. We choose the proper source function, such that the true solution is $u(x, y, z, t) = (1+ t^3)\cos \pi x \cos \pi y\cos \pi z$, which satisfies the homogeneous Neumann boundary conditions and the initial value $u^{0}(x, y) = \cos \pi x \cos \pi y \cos \pi z$. For such a problem, there is no singularity at the initial point, and the theoretical convergence order is $\mathcal{O}\big(\tau^{2} +\Delta{x}^2+\Delta{y}^2\big)$. To verify such a result, we first choose $M_{x} = M_{y} = M_{z} = 2^{6}$, and change the value of $N$ from $2^{2}$ to $2^{5}$. The observed errors at the final point and the corresponding temporal convergence orders are shown in Table 5. Then, we choose $N = 2^{8}$, and change the value of $M_x$, $M_y$, and $M_z$ from $2^{2}$ to $2^{5}$. The observed errors and the spatial convergence orders are shown in Table 6. We can see that the numerical results are consistent with the theoretical results.

6.   Concluding remark

In this paper, we have proposed a novel QSC-$L1^+$ scheme specifically designed for solving the multi-term variable-order TF-MID equation with Neumann boundary conditions, which is often used to model solute transport in porous media. The new scheme has been proved to be unconditionally stable and convergent with order $\mathcal{O}\big(\tau^{\min{\{3-\alpha^*-\alpha(0), \ 2\}}} +\Delta{x}^2+\Delta{y}^2\big)$, with some proper assumptions on $\alpha^{m}(t)$ and without any restrictions on the solution of the original model. The numerical experiments have demonstrated that the results obtained using the proposed QSC-$L1^+$ scheme align well with the theoretical analysis, even when the variable-order function $\alpha^{m}(t)$ is not monotonically decreasing or even not smooth. Furthermore, to enhance computational efficiency, we incorporate a fast evaluation method based on the ESA technique into the QSC-$L1^+$ scheme, resulting in the QSC-F$L1^+$ scheme. This refined approach significantly reduces the computation cost, making it more practical for large-scale and time-consuming simulations.

Use of AI tools declaration

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

Acknowledgments

The work is supported by Shandong Provincial Natural Science Foundation (No. ZR2021MA020), and the Fundamental Research Funds for the Central Universities (No. 22CX03016A).

Conflict of interest

The authors declare there are no conflicts of interest.