A second order numerical method for a Volterra integro-differential equation with a weakly singular kernel

1.   Introduction

This article aims to study a second-order numerical method for the following Volterra integro-differential equation (VIDE) with a weakly singular kernel

where $\alpha\in(0, 1)$. $a(x), \; b(x)$, and $f(x)$ are smooth functions, and $\mu$ is an initial data. We assume that there exist two positive constants $\beta_1, \; \beta_2$ such that

To simplify the presentation, we assume that $\beta_1, \; \beta_2\leq\beta^*$, where $\beta^*$ is a positive constant. Based on the above assumptions, the problem (1.1) has a unique solution, $u(x)$, which satisfies the following lemma:

Lemma 1.1. [1, Theorem 4.1] If $f$ can be written $f(x) = f_1(x)+x^\beta f_2(x),$ where $\beta > 0$ and $\beta\neq1, 2, \ldots, N.$ Then there exist positive constants $C, d$ such that

where $\delta = \min\left(2-\alpha, 1+\beta\right)$.

It is well known that Volterra integro-differential equations widely exist in biology, finance, population growth models and other fields (see, e.g., [2,3]). In recent years, there has been tremendous interest in developing finite difference [4,5,6], finite element [7,8], and spectral methods [9,10,11] for first-order and second-order Volterra integro-differential equations. Among the existing numerical methods, most of the researchers focus their attention on high-accuracy finite element methods and spectral methods. Therefore, it is very necessary to study a class of high-order finite difference methods for VIDE(s).

As far as we know, linear multi-step methods with a uniform time grid are widely used in discretizing the first-order time derivative of partial differential equations (see [12,13], for example). It should be pointed out that the variable step-size linear multi-step methods allow us to take different step-sizes for different scales, i.e., small step-sizes for the domain with solution rapidly varying and large for the domain with solution slowly changing. Therefore, the variable step-size linear multi-step methods demonstrate the prominent advantages of high accuracy compared to the constant step-size linear multi-step methods. Recently, Liao and Zhang [14] developed the variable two-step backward differentiation formula (BDF2) to discretize the time derivative of diffusion equations and gave a new theoretical framework by using the positive semi-definiteness of BDF2 convolution kernels and a class of orthogonal convolution kernels for the first time. Furthermore, Liao et al. [15] derived the stability and convergence analysis of the second-order backward difference formula (BDF2) with variable steps for the molecular beam epitaxial model without slope selection. Wang et al. [16] gave stability and error estimates for time discretizations of linear and semi-linear parabolic equations by the two-step backward differentiation formula (BDF2) method with variable step sizes. In [17], the authors proposed linearly implicit backward differentiation formulas with variable step sizes to solve the two-dimensional incompressible Navier Stokes equations.

Inspired by the above references [14,15,16,17], the main purpose of this paper is to develop a second-order finite difference scheme for VIDE (1.1). This paper is organized as follows: In Section 2, we proposes a finite difference scheme on a graded mesh by using the variable step-size BDF2 to discretize the first derivative term and the linear interpolation technique to approximate the integral term and prove the stability of our proposed numerical method. In Section 3, we will show how to improve the point-wise accuracy to second order by selecting appropriate graded mesh parameters. Finally, the theoretical results are verified by numerical experiments in Section 4.

2.   Discretization scheme

Let $\bar{\Omega}^{N}: = \left\{0 = x_0 < x_1 < \cdots < x_N = L\right\}$ be a graded mesh (see [18]), where the grid points are given by $x_i = L\left(\frac{i}{N}\right)^\gamma, \; i = 0, 1, \cdots, N$ and $\gamma\in[1, \infty)$ is a given real number. For $i = 1, 2, \ldots, N$, let $h_i = x_i-x_{i-1}$ be the $i-$th mesh step and $h = \max\limits_{1\leq i\leq N}h_i$ be the corresponding maximum step size. Furthermore, we denote the $i-$th step-size ratios by $r_{i} = \frac{h_{i}}{h_{i-1}}, \; i = 2, \ldots, N.$ Obviously, $r_{\max} = \max\limits_{2\leq i\leq N-1} r_{i} = 2^{\gamma}-1$. Throughout this paper, for any continuous function $g(x)$, let $g_i = g(x_i)$, and $C$ represents a positive constant independent of the mesh parameter $N$.

Let $P_{i, k}g$ be the Lagrange interpolating polynomial of a function $g$ over points $x_{i}, x_{i-1}, \ldots, x_{i-k}$. Then the BDF2 formula can be given by

where $\nabla g_i: = g_i-g_{i-1}$. In addition, denote $D_2g_1: = \nabla g_1/h_{1}$. Then, on the above graded mesh $\bar{\Omega}^{N}$, we construct the following discretization scheme to approximate problem (1.1):

where $u_{i}^{N}$ is the approximation solution of $u(x)$ at point $x = x_{i}$ and

In the following numerical analysis, the method of the discrete orthogonal convolution (DOC) kernels $\theta_{i-j}^i$ will play an important role. Firstly, we rewrite the BDF2 formula (2.1) into the following formal

where $b_{i-k}^{i}$ are defined by $b_0^1: = 1/h_1$, and

The DOC kernels $\theta_{i-j}^i$ have the following the property (see [14])

where $\delta_{ik}$ is the Kronecker delta symbol. Furthermore, we have

We now introduce the following lemma for the DOC kernels $\theta_{i-j}^i$.

Lemma 2.1. [14, Lemma 2.3] For $i\geq1$, the DOC kernels $\theta_{i-j}^i$ satisfy the following formula

Next, the stability result of the discretization scheme (2.2) is given as follows:

Lemma 2.2. Assume that $h\leq\left(1-\alpha\right)/(16\beta^*)$. Then the discrete solution $u_i^N$ of the discretization scheme (2.2) satisfies the following formula

Proof. Firstly, for the integral term of (2.2), we have

For $i\geq1$, multiplying both sides of the discretization scheme (2.2) by the DOC kernels $\theta_{i-j}^{i}$ and summing $j$ from 1 to $i$ yields,

Applying Eq (2.5) to Eq (2.8), one has

Multiplying both sides of the above equation (2.9) by $2u_i^N$ and summing the resulting equality from 1 to $n$, one has

Furthermore, we have

Set $\left|u_m^N\right|: = \max\limits_{0\leq i\leq n}\left|u_i^N\right|$, letting $n = m$ in the inequality (2.10), we get

Using Lemma 2.1 to Eq (2.11), one has

If $h\leq\left(1-\alpha\right)/(16\beta^*)$, we have

Based on the discrete Grönwall inequality [19, Lemma 3.2], we have

which completes the proof. □

3.   Error analysis

For $i = 0, 1, \cdots, N$, let $e_i = u_i-u_i^N$ be the absolution error between numerical solution $u_i^N$ and exact solution $u(x)$ at point $x = x_i$. Then the error equation can be written by

where $R_{1, i}, \; R_{2, i}$ are characterized by

Lemma 3.1. The truncation error $R_{1, i}$ and $R_{2, i}$ estimations can be given by the following inequality:

Proof. For the truncation error $R_{1, 1}$, it follows from Taylor's expansion formula and Lemma 1.1 that

Similarly, for $i\geq2$, based on [14], one has

Furthermore, based on Eq (3.5), yields,

and

For $R_{2, i}$, $1\leq i\leq N$, it is easy to obtain

where we have used the following inequality:

for any decreasing function $\phi > 0$ on $[a, b]$, see [20]. When graded mesh parameter $\gamma\leq2/\delta$, Eq (3.8), the following estimates can be given:

where the following inequality is used

Conversely, if $\gamma > 2/\delta$,

where $\xi_k\in\left(\frac{k-1}{N}, \frac{k}{N}\right)$. To sum up, this lemma is proved. □

Now, we derive the main result of this paper as follows:

Theorem 3.1. Let $u_n$ be the exact solution of problem (1.1) at the point $x = x_n$ and $u_n^N$ be the solution of problem (2.2) on the mesh $\bar{\Omega}^{N}$. Then, under the condition $h\leq\left(1-\alpha\right)/(16\beta^*)$, we have

Proof. For $n\geq1$, Lemma 2.2 and Eq (3.1) imply that

where

Based on Lemma 2.2 and Lemma 3.1, it is easy to obtain

Next, for $P_n$, the following estimate can be obtained by exchanging the summation order

where the fact $\sum\limits_{i = j}^n\theta_{i-j}^i\leq Ch_j$ is used. Then it is easy to get the following estimations:

and

Furthermore, by using $h_j\leq T\gamma N^{-1}\left(j/N\right)^{\gamma-1}$, yields,

If $\gamma < 2/\delta$, we have

If $\gamma = 2/\delta$, by using $\sum\limits_{j = 3}^n j^{-1}\leq\int_1^n\frac{1}{x}dx\leq \ln N$, yields,

If $\gamma > 2/\delta$, the following estimation will be obtained through the definition of the definite integral

According to Eqs (3.14)–(3.20), the estimation of $P_n$ in Eq (3.12) can be obtained. The desirable result can be followed by Eq (3.13).

□

4.   Numerical results and discussion

In order to verify our theoretical results, we consider the following test problem [1]

where $f = \left(2-\alpha\right)x^{1-\alpha}+\frac{\Gamma(1-\alpha)\Gamma(3-\alpha)}{\Gamma(4-2\alpha)}x^{3-2\alpha}$. The exact solution of this problem is $u(x) = x^{2-\alpha}e^{-x}$. Since the exact solution is given, the maximum absolute error and the convergence order are calculated as follows:

For different mesh parameters $\gamma$, $N$ and $\alpha$, Table 1 shows the maximum absolute errors and convergence orders. It is shown from these numerical experiments that they complement the theoretical results given in Theorem 3.1.

5.   Concluding remarks

Based on the variable step size BDF2, this paper proposes a second-order numerical method on the graded mesh to solve a Volterra integro-differential equation with a weak singular kernel and gives rigorous stability and convergence analysis for our presented numerical method. In the further work, by using the analysis of BDF3 given in [21], we shall study a third-order numerical method for the Volterra integro-differential equation with a weakly singular kernel.

Author contributions

Li-Bin Liu: Methodolog, writing review and editing; Limin Ye: software and writing original draft preparation; Xiaobing Bao: software, formal analysis and writing review and editing; Yong Zhang: numerical analysis and check the English writing. Further, all the authors have checked and approved the final version of the manuscript.

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.

Acknowledgement

This work was supported by the National Science Foundation of China (12361087), the Guangxi Science Technology Base and special talents (Guike AD20238065, Guike AD23023003), the projects of Excellent Young Talents Fund in Universities of Anhui Province (gxyq2021225) and Innovation Project of Guangxi Graduate Education (YCSW2024463).