Dissipativity and contractivity of the second-order averaged L1 method for fractional Volterra functional differential equations

1.   Introduction

Fractional Volterra functional differential equations (F-VFDEs), including fractional ordinary differential equations (F-ODEs), fractional delay differential equations (F-DDEs), fractional integro-differential equations (F-IDEs), fractional delay integro-differential equations (F-DIDEs) and other types, which appear in practice, are widely used to simulate some scientific problems in many fields of biology and finance ^[30,34]. Recently, F-VFDEs have received considerable attention because they can more accurately provide mathematical models of real-life problems with memory and hereditary characteristics than integer-order Volterra functional differential equations (VFDEs) due to the non-locality of the fractional derivative. Further, some studies (such as ^{[1,4,5,6,7,11,28]}) have been devoted to the existence and uniqueness of the solution for F-VFDEs and its special cases, which provides a theoretical foundation for its numerical computation and analysis.

Let $\mathbb{R}^{m}$ be an $m-$dimensional Euclidian space with the inner product $\langle\cdot, \cdot\rangle$ and the corresponding norm $\|\cdot\|$, and $C_{\mathbb{R}^{m}}(I)$ be a Banach space consisting of all continuous mapping $y:I\rightarrow \mathbb{R}^{m}$. In this paper, we study the long time behavior of averaged L1 method for the initial value problem of F-VFDEs:

where $T > 0$ and $0\leq\tau\leq +\infty$ are real constants, $\varphi \in C_{\mathbb{R}^{m}}[-\tau, 0]$ is a given initial function and $D_t^{\alpha}$ denotes the Caputo time fractional derivative of order $\alpha\in(0, 1)$, which is defined by

where the Riemann-Liouville fractional integral operator $\mathcal{I}_{t}^{\alpha}$ is given by

Throughout this paper, we assume that $f(t, y(t), y(\cdot))$ is independent of the values of the function $y(\xi)$ with $t < \xi\leq T$, i.e., $f(t, y(t), y(\cdot))$ as a Volterra functional, and that problem (1.1) has a unique true solution $y(t)$ with the required regularity. On this basis, we shall consider two types of problems. Firstly, we suppose that the continuous function $f:[0, +\infty)\times\mathbb{R}^{m}\times C_{\mathbb{R}^{m}}[-\tau, +\infty)\rightarrow \mathbb{R}^m$ satisfies the dissipative structural condition:

where $\gamma\geq 0, \; \beta < 0$ and $\lambda\geq 0$ are real constants, and the functions $\mu_{1}(t)$ and $\mu_{2}(t)$ are assumed to satisfy

For convenience, we shall always use the symbol $\mathcal{D}(\gamma, \beta, \lambda, \mu_{1}, \mu_{2})$ to denote the problem class consisting of all the problems (1.1) satisfying the condition $(1.2)$.

Secondly, we assume that the continuous function $f:[0, T]\times\mathbb{R}^{m}\times C_{\mathbb{R}^{m}}[-\tau, T]\rightarrow \mathbb{R}^m$ satisfies the one-sided Lipschitz condition

and the Lipschitz condition

where $p < 0$ and $q\geq 0$ are real constants, and the functions $\mu_{1}(t)$ and $\mu_{2}(t)$ also satisfy the condition (1.3). Similarly, we shall always use the symbol $\mathcal{D}(p, q, \mu_{1}, \mu_{2})$ to denote the problem class consisting of all the problems (1.1) satisfying the conditions (1.4) and (1.5). Note that the constant $T$ may also be $+\infty$, but in this case, the intervals $[0, T]$ and $[-\tau, T]$ should be replaced by $[0, +\infty)$ and $[-\tau, +\infty)$, respectively.

When $\alpha = 1$, the F-VFDEs (1.1) are reduced to the classical integer-order VFDEs. Many papers have been obtained the dissipativity, contractivity and asymptotic stability results of the true and numerical solution to VFDEs, see ^{[13,14,15,29,31,32,33]}, and (related works for the special cases of VFDEs can be found in the references therein). Dissipative dynamical systems are widely used in physics and engineering with the property of possessing a bounded absorbing set that all trajectories enter in a finite time and, thereafter, remain inside ^[18,20,22]. So it's natural to ask whether the true and numerical solution to F-VFDEs (1.1) can also possess dissipativity, contractivity, and asymptotic stability properties similar to those of the VFDEs. There is no doubt that the answer is positive, but the literature about the long time behavior of solutions to F-VFDEs is limited.

Time delays occur in many interconnected real systems due to transportation of energy, materials and information. If the delay effects and memory characteristics of the real systems are taken into consideration at the same time, it could lead to various F-DDEs. In ^[10], Katja Krol discussed the asymptotic properties of $d$-dimensional linear F-DDEs and obtained the necessary and sufficient conditions for asymptotic stability of equations of this type using the inverse Laplace transform method, and proved polynomial decay of stable solutions. Kaslik and Sivasundaram ^[9] presented several analytical and numerical methods for the asymptotic stability and bounded-input, bounded-output stability analysis of linear F-DDEs. Additionally, the asymptotic stability of nonlinear discrete fractional pantograph equations with nonlocal initial conditions has been investigated in ^[2]. In 2015, Wang and Xiao ^[23] studied the dissipativity and contractivity of the F-ODEs, the particular cases of the F-VFDEs (1.1) without the functional term. ^[24] obtained the dissipativity and stability of the F-VFDEs (1.1) by using the fractional generalization of the Halanay-type inequality under almost the same assumptions as the classic integer-order VFDEs. Moreover, Wang et al. ^[25] established the contractivity and dissipativity of time fractional neutral functional differential equations and proved that the solutions have a polynomial decay rate. From the perspective of structure-preserving algorithms, it is worthwhile to investigate whether or not numerical methods retain the qualitative behavior of the underlying system. Motivated by this, Wang and Zou ^[27] not only discussed the asymptotic behavior with algebraic decay rate of the exact solution of the F-VFDEs (1.1), but also proved that the two schemes based on the Grünwald-Letnikov formula and L1 method are dissipative and contractive, and can preserve the exact algebraic decay rate as the continuous equations. It should be pointed out that the polynomial/algebraic decay rate of the solutions is a significant feature for fractional differential equations, and is also the essential difference between fractional and integer differential equations, mainly because of the nonlocal nature of fractional derivatives. Later, Wang et al. ^[26] improved the existing algebraic dissipativity and contractivity rates of the solutions to the scalar F-ODEs, and further established the contractivity and dissipativity with an algebraic decay rate of numerical solutions to fractional backward differential formulas under some assumptions on the weight coefficients. Li and Zhang ^[12] applied the L1 method with the linear interpolation procedure to solve the nonlinear fractional pantograph equations, and proved that the proposed numerical scheme can inherit the long time behavior of the underlying problems without any step size restrictions.

When considering the applicability of numerical methods for F-VFDEs, it is highly desirable to have numerical methods which, not only inherit properties of the underlying system, but also possess higher accuracy. However, we regretfully find that the existing research ^[12,26,27] devoted to the numerical dissipativity and contractivity analysis for the F-VFDEs and its special cases is still restricted within the limits of numerical methods that are of less than second-order accuracy. The reasons for this may be as follows: On the one hand, the nonlocal character of the fractional derivatives and the complex structure of the nonlinear right-hand side function in the problem (1.1) make the numerical analysis so much more complicated that some new analytical techniques or refined analyses may need to be developed. On the other hand, it can be seen from the literature ^[26,27] that the discrete version of the fractional generalization of the traditional Leibniz rule, which plays a central role in establishing numerical dissipativity and contractivity, holds only under the assumption that the weight coefficients have good signs and relationships as shown in Section 2, while the weight coefficients determined by the higher-order numerical methods generally do not satisfy these assumptions. Fortunately, the discrete convolution kernels obtained by using a novel second-order formula to approximate the Caputo fractional derivative have nice sign properties but not uniform monotonicity in ^[8]. Additionally, the uniform monotonicity for the coefficients of an averaged L1 scheme has been proved by Shen et al. ^[19] under some assumptions, which provide a possibility to break the aforementioned restriction. Inspired by this, we attempt to investigate the dissipativity and contractivity of the averaged L1 method for F-VFDEs, which is the aim of this paper.

The rest of this paper is organized as follows: In Section 2, we introduce the averaged L1 approximation of the Caputo fractional derivative and further construct a numerical method for F-VFDEs combining the averaged L1 scheme and appropriate piecewise interpolation. The dissipativity and contractivity with an algebraic decay rate of the numerical method are obtained in Section 3. In Section 4, some numerical experiments are carried out to verify our theoretical results. Finally, some conclusions are drawn in Section 5.

2.   Numerical method

Averaged L1 scheme, named as an averaged variant of the classic L1 scheme, is a fractional generalization of the Crank-Nicolson scheme, and has been proven to have second-order accuracy in ^{[8,17,19,21,35]}. In this paper, we use the notation in ^[19] to record this averaged scheme as the $\overline{L1}$ scheme, which is also called $L1^{+}$ formula in ^[8]. In this section, we will construct a numerical method based on the $\overline{L1}$ scheme to solve the initial value problem of F-VFDEs (1.1).

2.1.   $\overline{L1}$ approximation to the Caputo fractional derivative

Let $N$ be a positive integer, and consider the uniform time mesh on interval $[0, T]$ with the time stepsize $h = \frac{T}{N}$ and the mesh points $t_{n} = nh, \; n = 0, 1, \cdots, N$. For any function $y(t)\in C_{\mathbb{R}^{m}}[0, T]$, denote the piecewise linear interpolation function of $y(t)$ on each subinterval $[t_{k-1}, t_{k}]\; (1\leq k\leq N)$ as $(\Pi_{1}y)(t)$, i.e.,

thus

Taking the integral average for the Caputo fractional derivative $D_t^{\alpha}y(t)$ over each time interval $[t_{n-1}, t_{n}]$ and then approximating $y'(t)$ by $(\Pi_{1}y)'(t)$, leads to

where

and

Then, we can obtain the $\overline{L1}$ approximation for the Caputo fractional derivative $D_t^{\alpha}y(t)$ ^[19]:

From the formula (2.1), an easy computation gives rise to the precise expression of the coefficients $a_{n-k}^{(n)}$:

The following lemma gives some properties for the coefficients $a_{n-k}^{(n)}$ and $b_{n-k}^{(n)}$.

Lemma 2.1. Assume that $\alpha > 2-(\ln3/\ln2)\cong0.415$. Then, the coefficients $a_{n-k}^{(n)}$ and $b_{n-k}^{(n)}$ are defined, respectively, by formulas (2.1) and (2.2) to satisfy

${\rm(i)}$ $a_{0}^{(n)} > a_{1}^{(n)} > \cdots > a_{n-1}^{(n)} > 0$ for $n\geq 1$;

${{\rm(ii)}}$ $b_{0}^{(n)} = a_{0}^{(n)} = \frac{1}{\Gamma(3-\alpha)h^{\alpha}} > 0, \; b_{k}^{(n)} < 0, \; k = 1, 2, \cdots, n$;

${{\rm(iii)}}$ $\sum\limits_{k = 0}^{n}b_{k}^{(n)} = 0$.

Proof. Since the properties (ii) and (iii) can be acquired directly by the formula (2.2) and the property (i), we need only prove the property (i). For $n = 1$, it follows from formula (2.1) that $a_{0}^{(1)} = \frac{1}{\Gamma(3-\alpha)h^{\alpha}} > 0$. For $n\geq2, \; 1\leq k\leq n-1$, by the integral mean value theorem, there is $\xi_{k}\in(t_{k-1}, t_{k})$, such that

It is apparent, from the above equality, that $a_{n-k}^{(n)}$ is a positive function with respect to $\xi_{k}$. Consequently, differentiating both sides of the above equation with respect to $\xi_{k}$, we have

which implies that $a_{n-k}^{(n)}$ is a monotonically increasing positive function of $\xi_{k}$. It follows that $a_{1}^{(n)} > a_{2}^{(n)} > \cdots > a_{n-1}^{(n)} > 0$. Then, it remains to show that $a_{0}^{(n)} > a_{1}^{(n)}$ for $n\geq2$. By simple computation, we can derive from formula (2.4) that

Hence, the inequality $a_{0}^{(n)} > a_{1}^{(n)}$ is valid if, and only if, $\alpha > 2-(\ln3/\ln2)\cong0.415$. Thus we conclude that if $\alpha > 2-(\ln3/\ln2)\cong0.415$, then the property (i) holds. This completes the proof of Lemma 2.1.

It should be pointed out that $b_{k}^{(n)}$ is meaningless for $k > n$ by formula (2.2). Therefore, if there is no particular statement below, we always assume that $b_{k}^{(n)} = 0$ for $k > n$. From Lemma 2.1 and the references ^[26,27], we immediately obtain the following conclusion.

Lemma 2.2. Let $\left\{b_{k}^{(n)}\right\}_{k = 0}^{\infty}$ be weights obtained by the $\overline{L1}$ approximation (2.3) for the Caputo fractional derivative. If $\alpha\in(0.415, 1)$, then the following inequality holds:

2.2.   $\overline{L1}$ method for F-VFDEs

To discretize the initial value problem of F-VFDEs (1.1), taking the integral average over each subinterval $[t_{n-1}, t_{n}]$, we have

Then, the $\overline{L1}$ formula (2.3) approximates the Caputo derivative in Eq (2.6), the right rectangle rule deals with the integral, and an appropriate piecewise interpolation operator $\Pi^{h}$ treats the functional term. Thus we can propose the following $\overline{L1}$ method for the initial value problem (1.1) in F-VFDEs:

Here, the interpolation function $y^{h}(t)$ is an approximation to the true solution $y(t)$ of the problem (1.1), $y_{0}: = \varphi(0)$, and $y_{n}\in\mathbb{R}^{m}$ is an approximation to the value $y(t_{n})$ of $y(t)$ at the time point $t_{n}$. For simplicity, we always assume that the interpolation operator $\Pi^{h}$ satisfies the condition ^[29]:

where $-\tau\leq\bar{t}\leq t_{n}, y_{k}\in\mathbb{R}^{m}, k = 0, 1, \cdots, n$. The constant $c_{\pi}\geq 1$ is of moderate size and independent of $\bar{t}, \; n, \; y_{k}$ and $\varphi$, and the function $\eta(t)$ is defined by

where $\mathbb{Z_{+}}$ denotes a set which consists of all nonnegative integers, and $\bar{p} > 0$ be a positive integer depending only on the procedure of the interpolation.

From the condition (2.8), we can easily derive the canonical condition ^[14]:

Here and later, $\{z_{n}\}$ is a numerical solution sequence produced by applying the $\overline{L1}$ method (2.7) to any given perturbed problem

where $\phi(t)\in C_{\mathbb{R}^{m}}[-\tau, 0]$ is a given initial function.

3.   Dissipativity and contractivity of $\overline{L1}$ method

This section will focus on the dissipative and contractive analysis of $\overline{L1}$ method for the initial value problem (1.1). Before proceeding further, let us introduce two important lemmas.

Lemma 3.1 (^[16,26,27]). Consider the discrete Volterra equation

where the kernel $\left\{d_{k}\right\}_{k = 0}^{\infty}\in l^{1}$, i.e., $\sum\limits_{k = 0}^{\infty}|d_{k}| < \infty$. Then,

if, and only if,

Furthermore, let

where $\left\{s_{k}\right\}_{k = 0}^{\infty}$ are coefficients. If $\left\{d_{k}\right\}_{k = 0}^{\infty}\in l^{1}$ and the condition (3.2) holds, then we can get $\left\{s_{k}\right\}_{k = 0}^{\infty}\in l^{1}$ and the estimate $\|x\|_{l^{\infty}}\leq \|s\|_{l^{1}}\|c\|_{l^{\infty}}$.

Lemma 3.2 (^[27]). Consider the linear Volterra convolution equation

where $\left\{G_{k}\right\}_{k = 0}^{\infty}$ satisfies the spectral condition $\rho = \sum\limits_{k = 0}^{\infty}|G_{k}| < 1$.

${\rm(i)}$ If $\lim\limits_{n\rightarrow \infty}g_{n} = g_{\infty}$, then $\lim\limits_{n\rightarrow \infty}x_{n} = (1-\rho)^{-1}g_{\infty}$ ^[3].

${{\rm(ii)}}$ Let $C$ be a positive constant and $\alpha\in(0, 1)$. If $g_{n}\rightarrow \frac{C}{n^{\alpha}}$ as $n\rightarrow \infty$, then $x_{n}\rightarrow \frac{C(1-\rho)^{-1}}{n^{\alpha}}$ as $n\rightarrow \infty$ ^[26].

Based on Lemma 2.2 and the above two lemmas, we now give the main results of this section.

Theorem 3.1. Suppose $\{y_{n}\}$ be an approximation sequence produced by using the method (2.7) to solve the problem (1.1) $\in \mathcal{D}(\gamma, \beta, \lambda, \mu_{1}, \mu_{2})$ with $\beta+\lambda c_{\pi}^{2} < 0$, $\left\{b_{k}^{(n)}\right\}_{k = 0}^{\infty}$ are weights determined by the $\overline{L1}$ method, and $y^{h} = \Pi^{h}$ is the piecewise $\sigma-$degree $(\sigma\geq 1)$ Lagrangian interpolation operator satisfying the condition (2.8). Assume that $\alpha\in(0.415, 1)$, and $\eta(\hat{t}_{n})\geq 1$ with $\hat{t}_{n} = t_{n}-\mu_{2}(t_{n})$. Then, for any given $\varepsilon > 0$, there exists a positive integer $n_{0}$, such that

It follows that the method (2.7) is dissipative with an absorbing set

Proof. From the method (2.7) and the condition (1.2), we can get

In view of $\eta(\hat{t}_{n})\geq 1$, it follows further from the condition (2.8) and the inequality (3.6) that

Since $\alpha\in(0.415, 1)$, we can know from Lemma 2.2 that the inequality (2.5) holds. Hence, combining the inequalities (2.5) and (3.6), we have

Because of the fact that $\beta < 0$, the above inequality can be further deduced as

For the sake of simplicity, let

then the inequality (3.8) can be rewritten as the equivalent form

Now we consider the following two cases successively:

For the case of $\max\limits_{1\leq k\leq n}\|y_{k}\|^{2} = \|y_{n}\|^{2}$, the inequality (3.9) can be rewritten as

Set $u_{0}: = \|y_{0}\|^2, \; d_{0}: = R, \; u_{n}: = \|y_{n}\|^2, \; c_{n} = P$ and $d_{n}: = Q_{n}$ for $n\geq 1$, then the above inequality is equivalent to

To get a bound of $\|y_{n}\|$, we define a sequence $\{x_{n}\}$ by

Due to the fact that $\beta+\lambda c_{\pi}^{2} < 0$, we can derive from simple calculation that

This implies that $\left\{d_{k}\right\}_{k = 0}^{\infty}\in l^{1}$ and the condition (3.2) in Lemma 3.1 holds for $\zeta = 1$. Let $\sum\limits_{k = 0}^{\infty}s_{k} = \left(1-\sum\limits_{k = 0}^{\infty}d_{k}\right)^{-1}$, then it follows from $\rho_{1}\in(0, 1)$ that $\sum\limits_{k = 0}^{\infty}s_{k} = (1-\rho_{1})^{-1} = \sum\limits_{j = 0}^{\infty}(\rho_{1})^{j}$. Thus, it is easy to acquire that $s_{k}\geq 0$ for $k\geq 0$ and $\|s\|_{l^{1}} = \sum\limits_{k = 0}^{\infty}s_{k} = (1-\rho_{1})^{-1}$. Therefore, using Lemma 3.1 can yield the estimate

So, if we can prove the inequalities $0\leq u_{n}\leq x_{n}$ holds for all $n\geq 0$, then the bound of $\|y_{n}\|$ follows immediately. Next, we give the proof of the above statement by mathematical induction.

Since $u_{0} = \|y_{0}\|^2$ is given, we can choose $x_{0}$, such that $0\leq u_{0}\leq x_{0}$. Let $c_{0} = (1-d_{0})x_{0}$, then

which indicates the inequality (3.10) is valid for $n = 0$. When $n = 1$, subtracting equality (3.11) from inequality (3.10) yields

Hence, it can be shown that

By the definition of $d_{0}$ and $\beta+\lambda c_{\pi}^{2} < 0$, we get $0 < d_{0} < 1$. Then, the inequality $u_{1}\leq x_{1}$ can be deduced immediately. Similarly, for $n\geq 2$, taking the difference between inequality (3.10) and equality (3.11) leads to

i.e.,

Thus, we obtain $u_{n}\leq x_{n}$. In conclusion, we have proved that $0\leq u_{n}\leq x_{n}$ holds for all $n\geq 0$. By this proven result and the inequality (3.12), we have the estimate

As a result, we acquire

For the case of $\max\limits_{1\leq k\leq n}\|y_{k}\|^{2} = \max\limits_{1\leq k\leq n-1}\|y_{k}\|^{2}$, we can easily observe from the inequality (3.9) that

For simplicity, define

Then, the inequality (3.14) can be rewritten as

Set $v_{0} = \|y_{0}\|^2, \; g_{n} = P, \; v_{n} = \|y_{n}\|^2$ and $G_{n} = Q_{n+1}+RI_{k}$ for $n\geq1$. Then, the inequality (3.15) is equivalent to

Define a sequence by

Note that $\beta+\lambda c_{\pi}^{2} < 0$. By simple computation, we can attain

and

Therefore, it follows from the result (i) in Lemma 3.2 that

Proceeding as in the proof of the case of $\max\limits_{1\leq k\leq n}\|y_{k}\|^{2} = \|y_{n}\|^{2}$, we can show that $0\leq v_{n}\leq x_{n}$ holds for all $n\geq 0$. Consequently, we can get the estimate

which implies that

According to the definition of limit, a combination of inequalities (3.13) and (3.16) shows that the open ball $B = B\left(0, \sqrt{-\frac{\gamma}{\beta+\lambda c_{\pi}^{2}}+\varepsilon}\right)$ is an absorbing set for any $\varepsilon > 0$, and the method (2.7) is dissipative. Therefore, we complete the proof of Theorem 3.1.

Theorem 3.2. Assume that the piecewise $\sigma-$degree $(\sigma\geq 1)$ Lagrangian interpolation operator $\Pi^{h}$ satisfies the canonical condition (2.9), and $\eta(\hat{t}_{n})\geq 1$ with $\hat{t}_{n} = t_{n}-\mu_{2}(t_{n})$. Then, the numerical solutions $\{y_{n}\}$ and $\{z_{n}\}$, attained by using the method (2.7) to solve respectively the problems (1.1) and (2.10), which belong to the class $\mathcal{D}(p, q, \mu_{1}, \mu_{2})$ with $p+qc_{\pi} < 0$, satisfy the contractivity estimate

provided that $\alpha\in(0.415, 1)$, where $C_{\alpha}$ is a positive constant independent of n.

Proof. Let $e_{k} = y_{k}-z_{k}$. By the method (2.7) and conditions (1.4) and (1.5) can yield

Since $\alpha\in(0.415, 1)$, it follows from Lemma 2.2 that the inequality (2.5) holds. Therefore, using inequalities (2.5) and (3.18) together with the canonical condition (2.9), we can get

Thanks to the fact that $b_{0}^{(n)} > 0$ and $p < 0$, the above inequality can be rewritten as

Set $F_{n-k} = \left|b_{n-k}^{(n)}\right|\Big/\left(b_{0}^{(n)}-p\right)$ and $H = qc_{\pi}\Big/\left(b_{0}^{(n)}-p\right)$, such that the above inequality is equivalent to

Next, we consider the following two cases successively.

For the case of $\max\limits_{1\leq k\leq n}\|e_{k}\|^{2} = \|e_{n}\|^{2}$, the inequality (3.19) can be further deduced as

In view of $p+qc_{\pi} < 0$, it is easy to check that $H\in(0, 1)$. Thus, we can further obtain

Define a sequence $\{x_{n}\}$ satisfying

By some routine calculations, we have

According to the Taylor formula, it holds that

Hence, the result (ii) in Lemma 3.2 leads to

From the mathematical induction method and a similar proof process in Theorem 3.1, it can be proven that $0\leq \|e_{n}\|^{2}\leq \|x_{n}\|^{2}$ is valid for all $n\geq0$. As a result, we can get the estimate

For the case of $\max\limits_{1\leq k\leq n}\|e_{k}\|^{2} = \max\limits_{1\leq k\leq n-1}\|e_{k}\|^{2}$, the inequality (3.19) has the equivalent form

Then, some simple computations yield

Similarly, combining the formula (3.15) and the result (ii) in Lemma 3.2 yields the estimate

Based on the estimates (3.17) and (3.22), it can easily be concluded that the method (2.7) is contractive and there exists a positive constant $C_{\alpha}$ independent of $n$, such that the contractive estimate (3.17) holds. The proof of Theorem 3.2 is now completed.

4.   Numerical experiments

In this section, we will present some numerical experiments to verify our theoretical results in previous sections.

Example 4.1. Consider the nonlinear F-DIDE:

where $y_{1}(t), y_{2}(t)$ are real-valued scalar functions.

Let $u = (u_{1}, u_{2})^{T}, \psi(t) = (\psi_{1}(t), \psi_{2}(t))^{T}$, and

Then

Thus, we get

which means that the problem (4.1) belongs to the class $\mathcal{D}(\gamma, \beta, \lambda, \mu_{1}, \mu_{2})$ with $\gamma = 2, \; \beta = -3.5, \; \lambda = 0.6, \; \mu_{1}(t) = 0$ and $\mu_{2}(t) = 1$. Therefore, the F-DIDE (4.1) is dissipative according to Theorem 2.1 in ^[24].

Let $v = (v_{1}, v_{2})^{T}, \chi = (\chi_{1}(t), \chi_{2}(t))^{T}$. Then,

which gives

Further,

which leads to

Consequently, it follows from the inequalities (4.2) and (4.3) that the system (4.1) belongs to the class $\mathcal{D}(p, q, \mu_{1}, \mu_{2})$ with $p = -3.6, \; q = 2\sqrt{1.9}, \; \mu_{1}(t) = 0$ and $\mu_{2}(t) = 1$. Therefore, the F-DIDE (4.1) is contractive and asymptotically stable according to Theorem 2.2 in ^[24].

Now, we apply the $\overline{L1}$ method (2.7) with piecewise linear interpolation operator (i.e., $\sigma = 1$ and $c_{\pi} = 1$ ^[14]) to solve the system (4.1). It is easy to check that

and

We choose the step size $h = 0.01$, and successively take $\alpha = 0.5, \; 0.7, \; 0.9, \; 1$ to plot the numerical solutions of the system (4.1) with the following two different initial functions:

$\text{(I)}$ $y_{1}(t) = \sin(t), \; y_{2}(t) = \cos(t), \; t\in[-1, 0]$;

$\text{(II)}$ $y_{1}(t) = 5\cos(t), \; y_{2}(t) = 2\sin(t), \; t\in[-1, 0]$,

respectively. The numerical results are given in Figures 1–2. It can be seen from the first to third subfigures in Figures 1–2 that the numerical solutions can preserve the dissipativity of the problem (4.1), which confirms the result of Theorem 3.1. Comparing the four subfigures in Figures 1–2, we find that the solutions of the integer-order DIDE ($\alpha = 1$) decay exponentially into a given ball. However, the solutions of the F-DIDE no longer decay exponentially but with a polynomial rate into a bounded absorbing set because of the nonlocal nature of the fractional derivative.

Let

and

be the difference between two solutions $y(t)$ and $z(t)$ with the different initial functions $\varphi(t)$ and $\phi(t)$. In the numerical simulations, we take the initial functions (I) as $\varphi(t)$, and initial functions (II) as $\phi(t)$. To observe the contractivity behavior of numerical solutions more intuitively, Figure 3 draws the error curves of the numerical solutions of F-DIDE (4.1) with the different initial functions $(\mathrm{I})$ and $(\mathrm{II})$ for $\alpha = 0.5, \; 0.7, \; 0.9, \; 1$. We can observe from Figure 3 that the numerical solutions are contractive, and, the larger the order $\alpha$ is, the faster the contractive rate becomes. Furthermore, we also find that the contractivity of the numerical solutions for the integer-order DIDE has an exponential decay rate, while the contractivity for the F-DIDE is polynomial.

From Theorem 6 in ^[27], we have asymptotic contractive rate

To further measure the quantitative behavior of the contractivity rate corresponding to two different initial functions $\varphi(t)$ and $\phi(t)$, we define an index function $I_{\alpha}$ as in ^[26,27]:

Clearly, the index function

as $t\rightarrow \infty$, and it is independent of the initial value $\max\limits_{-1\leq\xi\leq0}\|\varphi(\xi)-\phi(\xi)\|^2C_{\alpha}$. Thus, we can take

in our numerical experiments. Table 1 shows the values of the index function $I_{\alpha}(t)$ at

for

with $h = 0.01$. We find that the contractivity rate is about $\|e(t)\|^2 = O(t^{-2\alpha})$, which is about twice as much as our theoretical prediction for numerical contractivity rate given in Theorem 3.2. For scalar F-ODE or essentially decoupled linear systems, Wang et al. ^[26] obtained the optimal contractivity rate $\|e(t)\|^2 = O(t^{-2\alpha})$ by directly estimating the decay rate of $\|y(t)-z(t)\|$ to avoid the square-root operation of the Mittag-Leffler function. Based on this and the results shown in Table 1, we believe that the optimal contractivity rate for F-VFDEs can be achieved theoretically and numerically with some new analytical techniques developed in the future.

5.   Conclusions

In this paper, we mainly investigate the long time behavior of the $\overline{L1}$ method for the initial value problem of F-VFDEs. With the help of the fact that the weight coefficients of $\overline{L1}$ scheme for Caputo fractional derivative have good signs and uniform monotonicity for $\alpha\in(0.415, 1)$, we prove that the $\overline{L1}$ method is dissipative and contractive, and can preserve the algebraic contractive rate. Finally, the numerical experiments are conducted to illustrate our theoretical results.

Acknowledgments

The authors are grateful to Professor Dongling Wang (Xiangtan University) for very helpful discussions about some technical issues, and to the two anonymous referees for valuable and useful comments and suggestions. This research is supported by the National Natural Science Foundation of China (Grant No. 12071403), the Research Foundation of Education Department of Hunan Province of China (Grant No. 21A0108), and the Postgraduate Scientific Research Innovation Project of Xiangtan University (Grant No. XDCX2021B108).

Conflict of interest

We declare that there are no conflicts of interest.