Mittag-Leffler stability of complete monotonicity-preserving schemes for time-dependent coefficients sub-diffusion equations

A typical application of time fractional equations is to describe various anomalous diffusion models, namely sub-diffusion and super-diffusion phenomena. For the standard diffusion process, which is usually described by the classical heat equation, where the mean square displacement (MSD) of particles is a linear function of time, i.e., $\langle x^{2}(t)\rangle\sim Ct$. ¹¹1As usual, our notation $\sim$ here represents the asymptotic limit $\lim_{t\to\infty}\frac{\langle x^{2}(t)\rangle}{t}=C$, where $C$ represents a general positive constant, which may take different values in different places, but is always independent of $t$ or $n$. We also use the standard notation $O$ to represent a certain asymptotic order, for example $y(t)=O(t^{-\alpha})$ means $\lim_{t\to\infty}\left|\frac{y(t)}{t^{-\alpha}}\right|\leq C$. However, many actual model data indicate that the MSD of anomalous diffusion processes may have a form $\langle x^{2}(t)\rangle\sim Ct^{\alpha}$, where $\alpha\in(0,2)$. Based on random walk theory, we can rigorously derive fractional diffusion equations, where Caputo derivatives are typically used in the time direction instead of classical first-order derivatives [3, 8]. Specifically, when $\alpha\in(0,1)$, it is referred to as the sub-diffusion process.

For the time fractional order model, its solutions exhibit asymptotic behavior completely different from the standard diffusion equation over a long time interval. In fact, solutions to time fractional equations typically have algebraic decay rates, while standard first-order equations have exponential decay rates. This can be seen from the simple example below. Consider the fractional ODE: $D_{t}^{\alpha}y(t)=\lambda y(t)$ with $y(0)=y_{0}>0$, where the Caputo fractional derivative of order $\alpha$ is defined by $D_{t}^{\alpha}y(t):=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{t}(t-s)^{-\alpha}y^{\prime}(s)ds$ [8]. The solution is given by $y(t)=y_{0}E_{\alpha,1}(\lambda t^{\alpha})$, where $E_{\alpha,\beta}(t)$ is the Mittag-Leffler funciton. When $\lambda<0$, we have $y(t)\sim Ct^{-\alpha}$ by the asymptotic expansion formula of Mittag-Leffler funciton. This long-term algebraic decay rate is commonly referred to as Mittag-Leffler stability. Another surprising result is presented in [19], which shows that for the fractional ODE: $D_{t}^{\alpha}y(t)=\lambda y^{\gamma}(t)$ with $y_{0}>0,\gamma>0$ and $\lambda<0$, the solution has the asymptotic decay rate $y(t)\sim Ct^{-\frac{\alpha}{\gamma}}$.

Once we have the above two results and combine appropriate energy methods, we can obtain the long-term decay rate estimation of the standard $L^{s}(\Omega)$-norm $(s>1)$ for the solution of linear or some nonlinear anomalous diffusion model. See for example [4, 10, 19, 2]. Currently, most of these analyses are focused on constant coefficient equations. However, the anomalous diffusion problems can be both constant coefficients and variable coefficients. Compared with constant coefficient equations, the variable coefficient equations are appropriate for modeling of anomalous diffusive processes in turbulent media. For instance, modeling sub-diffusive transport inside a fractured porous media [6]. Diffusion representation of tracer particles in heterogeneous biological ambient fluids [1]. In particular, the application of the time dependence coefficient is a new perspective of coupling the spatial and temporal variables, which avoids the divergence caused by the long tail of the temporal probability distribution at a fixed time [5]. For more information on the application of variable coefficients to more complex physical and biological systems, see also [3, 5, 7].

The aim of this paper is to study both theoretically and numerically the long time Mittag-Leffler stability of the solutions for a class of nonlinear anomalous diffusion models with time-dependent coefficients

where $0<\alpha<1$, $\beta>-\alpha$, $\nu>0$, and $\Omega\subset\mathbb{R}^{d}~{}(d\geq 1)$ is a bounded domain with smooth boundary $\partial\Omega$. And $\mathcal{N}(u)$ can be linear or nonlinear operator, some typical examples including the Laplace operator $\mathcal{N}(u)=-\Delta u$; the $p$-Laplace operator $\mathcal{N}(u)=-\Delta_{p}u$ with $\Delta_{p}u=\mathrm{div}(|\nabla u|^{p-2}\nabla u)$ and the mean curvature operator $\mathcal{N}(u)=-\mathrm{div}\Big{(}\frac{\nabla u}{\sqrt{1+|\nabla u|^{2}}}\Big{)}$. Following [4, 16],we assume that $\mathcal{N}(u)$ satisfies the following structural assumption

where $s\in(1,\infty)$, $\gamma\in(0,+\infty)$ and $C_{N}>0$. One can verify that the above example of $\mathcal{N}(u)$ indeed satisfies the elliptical structure assumption [4, 14]. The non-zero source terms $f(t,x)$ and satisfies the decay estimate

Let $u$ be the solution of (1). Under the above assumptions, we can establish the long-term estimation of $\|u(t,\cdot)\|_{L^{s}(\Omega)}$, i.e. Mittag-Leffler stability, by combining the energy method with the long-term estimation of time fractional ODEs. For the homogeneous linear constant coefficient equation with $\beta=0,f=0$ and $\mathcal{N}(u)=-\Delta u$, we can apply standard eigenvalue and eigenfunction pairs methods to obtain that $\|u(t,\cdot)\|_{L^{s}(\Omega)}\sim Ct^{-\alpha}$. For the constant coefficient nonlinear equation with $\beta=0$ and $\mathcal{N}(u)$ satisfying (2), it is shown that $\|u(t,\cdot)\|_{L^{s}(\Omega)}\sim Ct^{-\alpha/\gamma}$; see [16, 4, 2]. And for the homogeneous varying coefficient nonlinear equation with $f=0$ and $\mathcal{N}(u)$ satisfying (2), the asymptotic estimate $\|u(t,\cdot)\|_{L^{s}(\Omega)}\sim Ct^{-(\alpha+\beta)/\gamma}$ for weak solution is obtained in [14], from which can see from that the time-dependent factor $t^{\beta}$ has a direct impact on the long-term decay rate of the solution.

At present, most numerical methods for time fractional equations focus on the reduction in numerical convergence rate caused by weak singularity of the solution near the initial moment. Various effective techniques have been developed to restore high-order convergence of numerical methods; see [13, 11, 15, 9]. From the perspective of the structure-preserving algorithm, we naturally hope that the numerical solution can preserve the long-term optimal decay rate of the sub-diffusion equation, that is, establish numerical Mittag-Leffler stability. At present, there are relatively few research results on qualitative analysis of numerical solutions to anomalous diffusion equations, especially long-term behavior analysis. For the linear time fractional evolutional equation, the numerical ML stability, i.e. $\|u_{n}\|_{L^{s}(\Omega)}\sim Ct_{n}^{-\alpha}$ is established in [18] through singularity analysis of the generating function of the numerical solution, where $u_{n}$ is the numerical approximation of $u$ at $t_{n}$. For the constant coefficient nonlinear equation with $\beta=0$ and $\mathcal{N}(u)$ satisfying (2), the numerical ML stability as $\|u_{n}\|_{L^{s}(\Omega)}\sim Ct_{n}^{-\alpha/\gamma}$ was recently established in [17], using the discrete comparison principle and the method of constructing fine discrete sup and sub solutions. The key to success here lies in a class of time fractional derivative structure-preserving discrete schemes, namely the $\mathcal{CM}$-preserving schemes developed in [12]. The discrete coefficients of this type of numerical method have many excellent properties.

The objective of this article is twofold.

• 1.

  We derive the long-term decay estimate for the solution of the general non-homogeneous variable coefficient nonlinear equation (1) with $\beta\neq 0,f\neq 0$ and $\mathcal{N}(u)$ satisfying (2) and prove that $\|u(t,\cdot)\|_{L^{s}(\Omega)}\sim Ct^{-(\alpha+\beta)/\gamma}$. Our results extend the work of [14, 2].

• 2.

  We establish the numerical ML stability for (1) based on the $\mathcal{CM}$-preserving schemes and show that $\|u_{n}\|_{L^{s}(\Omega)}\sim Ct_{n}^{-(\alpha+\beta)/\gamma}$. Firstly, we prove that the $\mathcal{CM}$-preserving schemes also satisfies the discrete comparison principle for variable coefficient equations. Then, applying this result, we construct discrete sup and sub solutions with variable coefficients to obtain the expected optimal numerical decay rate estimate. Our work in this article extend the results in [17] from homogeneous constant coefficient equations to non-homogeneous, variable coefficient equations.

The rest of this article is organized as follows. In Section 2, we derive the long time decay rate of continuous solutions for (1) and provide detailed proofs of the estimates. In Section 3, we recall $\mathcal{CM}$-preserving schemes on uniform mesh $t_{n}=n\tau$ with time step size $\tau>0$ and their main properties. In Section 4, we establish the discrete numerical ML stability for (1). We first derive the discrete comparison principle of $\mathcal{CM}$-preserving schemes for a class of time-dependent fractional ODEs. Then, by constructing discrete sub and sup solution we show the numerical ML stability for fractional non-homogeneous, variable coefficient ODEs. Then, by establishing a variable coefficient discrete energy inequality, we obtain the estimate of the long-term decay rate of the numerical solution of (1). Finally, we present several typical numerical examples to illustrate and verify the optical numerical decay rate in Section 5.

Consider the time-dependent coefficients fractional ODE, which highlight the influence of $t^{\beta}$:

where $0<\alpha<1$, $\beta>-\alpha$, $\gamma>0$ and $\nu>0$. The decay result concerning (4) reads as follows.

Note that when $\beta=0$, the fractional ODE (4) is reduced to the power nonlinear equation and its solution has decay rate $y(t)\sim Ct^{-\alpha/\gamma}$ [16]. Further, when $\beta=0$ and $\gamma=1$, it becomes fractional linear ODE and the solution has the decay rate $y(t)\sim t^{-\alpha}$. Our first result below is to extend the estimation of the long-term decay rate of homogeneous equation (4) to non-homogeneous cases with source terms.

The proof the Mittag-Leffler stability of fractional ODEs in [16, 14] is mainly based on the comparison principle and the construction of appropriate barrier functions.

It is noted that in (4), the nonlinear function $g(y):=\nu t^{\beta}y^{\gamma}$ is positive and nondecreasing, which allows the utilization of Lemma 2.2 to prove Lemma 2.1. Now we extend Lemma 2.1 to non-homogeneous cases with source terms.

Our main idea of proof is as follows: by constructing new fine upper and lower solutions, we can control the source function $f$ and obtain inequality of the corresponding homogeneous equation, which can be applied to Lemma 2.2.

As an immediate application of Theorem 2.1, we can derive the long-term decay rate for $\|u(x,t)\|_{L^{s}(\Omega)}$ for the solution of the time fractional PDE model (3) by combination of appropriate energy estimate.

From the perspective of structural preservation algorithms, it is natural to develop structural preservation algorithms for time fractional order equations, so that numerical solutions can accurately preserve the long-term optimal decay rate of the solution in Lemma 2.1 for nonlinear fractional ODEs and in Corollary 2.1 for nonlinear fractional PDEs. That is to develop numerical Mittag-Leffler stability for the model (1).

The Caputo fractional derivative $D_{t}^{\alpha}(y(t))$ has a convolutional structure, so we consider that its numerical approximation also has a corresponding discrete convolution structure, which takes the form of

where $\delta_{n}:=-\sum_{k=0}^{n-1}\omega_{k}$ and $y_{k}$ denotes the numerical approximation of $y(t_{k})$ on the uniform meshes $t_{n}=n\tau$ with step size $\tau>0$. The weight coefficients $\{\omega_{k}\}_{k=0}^{\infty}$ can be computed explicitly. The two most popular methods are Grünwald-Letnikov scheme and L1 method, both of which can be seen as fractional extensions of the classical backward Euler scheme.

By using Riemann-Liouville fractional integration, the initial value problem of the Caputo derivative can be transformed into the corresponding Volterra integral equation, which also has a convolutional structure and its kernel function is $k_{\alpha}(t):=t^{\alpha-1}/\Gamma(\alpha)$ for $t>0$. The standard kernel function $k_{\alpha}(t)$ is a very important class of completely monotonic ($\mathcal{CM}$) functions ²²2A function $g(t):(0,\infty)\to\mathbb{R}$ is said to be completely monotone if $(-1)^{n}g^{(n)}(t)\geq 0$ for all $n\geq 0$ and $t>0$. . Based on the basic idea of preserving structure algorithm, the $\mathcal{CM}$-preserving numerical methods for Riemann-Liouville fractional integration is developed in [12]. The new concept and framework of $\mathcal{CM}$-preserving numerical methods has many excellent properties in the discrete coefficients. This kind of numerical methods has recently been successfully applied to a class of nonlinear anomalous diffusion models [17]. It is proved in [12] that for $\mathcal{CM}$-preserving numerical methods the weights $\omega_{k}$ and $\delta_{n}$ have the convenient properties

These properties imply that there exist constants $0<C_{3}\leq C_{4}$, which are independent of $n$, sunch that

It is easy to verify the Grünwald-Letnikov scheme and L1 method in the above two examples are $\mathcal{CM}$-preserving and their coefficients meet the above properties. See more details in [12].

In the section, we study the numerical Mittag-Leffler stability for $\mathcal{CM}$-preserving numerical methods for nonlinear sub-diffusion models with variable coefficients. Specifically, we prove that for the $\mathcal{CM}$-preserving schemes, its numerical solution has a long-term optimal decay rate that is completely consistent with the corresponding continuous equation. Our proof mainly consists of the following three steps. Firstly, we establish the discrete comparison principle for variable coefficient equations. Then we apply this result to homogeneous equations with variable coefficients and establish the Mittag-Leffler stability of numerical solutions. Finally, for general non-homogeneous equations with source terms, we also establish the numerical Mittag-Leffler stability.