Extrapolation Method for Non-Linear Weakly Singular Volterra Integral Equation with Time Delay †

: This paper proposes an extrapolation method to solve a class of non-linear weakly singular kernel Volterra integral equations with vanishing delay. After the existence and uniqueness of the solution to the original equation are proved, we combine an improved trapezoidal quadrature formula with an interpolation technique to obtain an approximate equation, and then we enhance the error accuracy of the approximate solution using the Richardson extrapolation, on the basis of the asymptotic error expansion. Simultaneously, a posteriori error estimate for the method is derived. Some illustrative examples demonstrating the efﬁciency of the method are given.


Introduction
Delay functional equations are often encountered in biological processes, such as the growth of the population and the spread of an epidemic with immigration into the population [1,2], and a time delay can cause the population to fluctuate. In general, some complicated dynamics systems are also modeled by delay integral equations since the delay argument could cause a stable equilibrium to become unstable. The motivation of our work is twofold: one of the reasons is based on the first-kind delay Volterra integral equation (VIE) of the form [3] t qt k(t, s)y(s)ds = f (t), t ∈ I := [0, T], which was discussed and transformed into the second-kind equivalent form k(t, t)y(t) − qk(t, qt)y(qt) + t qt ∂k(t, s) ∂t y(s)ds = f (t), if k(t, t) = 0 for t ∈ I, the normal form was given by y(t) = f (t) + y(qt) + There has been some research [4][5][6] to the following form y(t) = f (t) + Another source of motivation comes from the weakly singular delay VIE [7][8][9] y(t) = f (t) + qt 0 K(t, s) (qt − s) λ G(s, y(s))ds, t ∈ [0, 1], where λ ∈ (0, 1), K(t, s) is smooth and G(s, y(s)) is a smooth non-linear function. However, there has not yet been investigated for the case where two integral terms are presented, the first integral term is the weakly singular Volterra integral and the second integral terms not only has weak singularity in the left endpoint but also its upper limit is a delay function, which is challenging to calculate. It is the aim of this paper to fill this gap.
Then, Equation (1) possesses a unique solution (see Theorem 1). In this paper, we consider the case where the solution is smooth. Some numerical investigations of delay VIE have been conducted, such as discontinuous Galerkin methods [23], collocation methods [24][25][26], the iterative numerical method [27], and the least squares approximation method [28]. In [29], an h p version of the pseudospectral method was analyzed, based on the variational form of a non-linear VIE with vanishing variable delays. The algorithm increased the accuracy by refining the mesh and/or increasing the degree of the polynomial. Mokhtary et al. [7] used a well-conditioned Jacobi spectral Galerkin method for a VIE with weakly singular kernels and proportional delay by solving sparse upper triangular non-linear algebraic systems. In [8], the Chebyshev spectral-collocation method was investigated for the numerical solution of a class of weakly singular VIEs with proportional delay. An error analysis showed that the approximation method could obtain spectral accuracy. Zhang et al. [9] used some variable transformations to change the weakly singular VIE with pantograph delays into new equations defined on [−1, 1], and then combined it with the Jacobi orthogonal polynomial.
The extrapolation method has been used extensively [30,31]. We apply the extrapolation method for the solution of the non-linear weakly singular kernel VIE with proportional delay. We prove the existence of the solution to the original equation using an iterative method, while uniqueness is demonstrated by the Gronwall integral inequality. We obtain the approximate equation by using the quadrature method based on the improved trapezoidal quadrature formula, combining the floor technique and the interpolation technique. Then, we solve the approximate equation through an iterative method. The existence of the approximate solution is validated by analyzing the convergence of the iterative sequence, while uniqueness is shown using a discrete Gronwall inequality. In addition, we provide an analysis of the convergence of the approximate solution and obtain the asymptotic expansion of the error. Based on the error asymptotic expansion, the Richardson extrapolation method is applied to enhance the numerical accuracy of the approximate solution. Furthermore, we obtain the posterior error estimate of the method. Numerical experiments effectively support the theoretical analysis, and all the calculations can be easily implemented.
This paper is organized as follows: In Section 2, the existence and uniqueness of the solution for (1) are proven. The numerical algorithm is introduced in Section 3. In Section 4, we prove the existence and uniqueness of the approximate solution. In Section 5, we provide the convergence analysis of the approximate solution. In Section 6, we obtain the asymptotic expansion of error, the corresponding extrapolation technique is used for achieving high precision, and a posterior error estimate is derived. Numerical examples are described in Section 7. Finally, we outline the conclusions of the paper in Section 8.

Existence and Uniqueness of Solution of the Original Equation
In this section, we discuss the existence and uniqueness of the solution of the original equation. There are two cases, 0 ≤ t ≤ T ≤ 1 and 1 < t ≤ T, that we will discuss in the following.
• Case I. For 0 ≤ s ≤ t ≤ T ≤ 1, by means of mathematical induction, when n = 1, Suppose that the following expression is established when n = k, Let n = k + 1; then, that is, the recurrence relation is established when n = k + 1, then the inequality (4) is also established. Next, we prove that the sequence y n (t) is a Cauchy sequence, i! is convergent, so the Cauchy sequence {y n } n∈N is convergent uniformly to y(t). Thus, y(t) is the solution to Equation (1), the existence is proved. • Case II. For 1 < s ≤ t ≤ T, the process is similar. Let γ = max{λ, µ}, when n = 1, Suppose that the following expression is established when n = k, Let n = k + 1. Then, we have i.e., the recurrence relation is established when n = k + 1, such that the inequality (6) is also established. For the sequence y n (t), Since the term i! is convergent, so the Cauchy sequence {y n } n∈N is convergent uniformly to y(t). Thus, y(t) is the solution to Equation (1), the existence is proved. Now, we prove that the solution to Equation (1) is unique. Let y(t) and v(t) be two distinct solutions to Equation (1), and denote the difference between them by w(t) = |y(t) − v(t)|. We obtain Let g(s) = Ls λ + Ls µ , then g(s) is a non-negative integrable function, according to Lemma 1. We obtain w(t) = 0, i.e., y(t) = v(t), the solution to Equation (1) is unique.

The Numerical Algorithm
In this section, we first provide some essential lemmas which are useful for the derivation of the approximate equation. Next, the discrete form of Equation (1) is obtained by combining an improved trapezoidal quadrature formula and linear interpolation. Finally, we solve the approximate equation using an iterative method. The process does not have to compute the integrals; hence, the method can be implemented easily.

Some Lemmas
Proof. The Taylor expansion of function u(x) at the point z is Similarly, the Taylor expansion of function u(y) at point z is combining (8) with (9), the proof is completed.
N , and t k = a + kh for k = 0, · · · , N, as for the integral b a G(t)dt. Then, the error of the modified trapezoidal integration rule has an asymptotic expansion where −1 < λ < 0, ζ is the Riemann-Zeta function and B 2j represents the Bernoulli numbers.

The Approximation Process
In this subsection, we describe the numerical method used to find the approximate solution to Equation (1). Let y(t) have continuous partial derivatives up to 3 on I, f (t), k 1 (t, s; y(s)), k 2 (t, s; y(s)) are four times continuously differentiable on I, D × R, D θ × R, respectively. Let y(t i ), y i denote the exact solution and approximate solution when t = t i , respectively. We divide (1). Then, where [qi] denotes the maximum integer less than qi. According to Lemma 3, we have For I 2 and I 3 , there are two cases.
• Case I. If [qi] = 0, then • Case II. If [qi] ≥ 1, we obtain y(qt i ) can be represented by linear interpolation of the adjacent points y(t [qi] ) and Then, the approximate expression of y(qt i ) is Then, (15) can be written as The approximation equations are as follows • Case II. When [qi] ≥ 1, where

Iterative Scheme
Now, the solution of the approximate equation can be solved by an iterative algorithm.
Step 2. Letỹ 0 i =ỹ i−1 , m := 0, then we compute y m+1 i (i ≤ N) as follows: • Case I. When [qi] = 0, • Case II. When [qi] ≥ 1, where Step 3. If |y m+1 i − y m i | ≤ , then letỹ i := y m+1 i and i := i + 1, and return to step 2. If otherwise, let m := m + 1, and return to step 2. Remark 1. In Section 3.2, we considered the regularity of k i (t, s; y(s))(i = 1, 2) only up to r = 2 in Lemma 3, since the desired accuracy has been obtained, and it is sufficient for the subsequent convergence analysis and extrapolation algorithm.

Existence and Uniqueness of the Solution to the Approximate Equation
In this section, we investigate the existence and uniqueness of the solution to the approximate equation. We first introduce the following discrete Gronwall inequality.
Proof. We discuss the existence of the approximate solution under two cases.
When h is sufficiently small, such that L 1 Therefore, the iterative algorithm is convergent and the limit is the solution to the approximation equation. The existence of approximation is proved when [qi] = 0. Now, we prove the uniqueness of approximation. Suppose y i and x i are both solutions to Equation (20). Denote the absolute differences as w i = |y i − x i |. We have where L = max{L 1 , L 2 }. When h is sufficiently small, such that L 1 According to Lemma 4 with A = 0, we have w i = 0, i.e., y i = x i , the solution of Equation (20)  (1) The first situation is [qi] + 1 = i, namely, when i ≤ 1 1−q , we have Let the step size h be small enough, such that L h 2 t λ i + (qt i ) µ (1 − β i ) ≤ 1 2 . Then, we can determine that |y m+1 The second situation is [qi] + 1 < i, namely, when i > 1 1−q , we obtain The above two situations show that the iterative algorithm is convergent and that the limit is the solution to Equation (21).
Next, we prove that the solution to Equation (21) is unique. Suppose y i and x i are both solutions to Equation (21). Denote the differences as w i = |y i − x i |, i = 1, · · · , N. Then, we have (1) The first situation is [qi] + 1 = i (i.e., when i ≤ 1 1−q ). Then, (24) entails By letting h be so small that L h 2 t . According to Lemma 4 with A = 0, we have w i = 0, and the solution of Equation (21) is unique.
(2) The second situation is [qi] + 1 < i (i.e., when i > 1 1−q ). Then, (24) can imply Letting h be so small that L h 2 t γ i ≤ 1 2 , then According to Lemma 4 with A = 0, we have w i = 0, i.e., y i = x i , the solution of Equation (21) is unique. Combining the above situations, the proof of Theorem 2 is completed.

Convergence Analysis
In this section, we will discuss errors caused by the process of obtaining discrete equations using a quadrature formula and interpolation technique and the errors caused by solving the discrete equation using iterative algorithms. According to the quadrature rule, Equation (12) can be expressed as From Lemmas 2 and 3, the remainders are In order to investigate the error between the exact solution and the approximate solution of Equation (1), we first give the following theorem. Proof. Subtracting (19) from (27), Letting h be so small, that h 2 t λ i L 1 ≤ 1 2 , it is easy to derive By Lemma 4, we have max The proof is complete.
Next, we evaluate the error arising from the iterative process.

Theorem 4.
Under the conditions of Theorem 2, y i is the solution of Equation (19) and y i is the approximate solution of Equation (1), and y i is defined by (21). The absolute error is denoted by e 2,i = |y i − y i |. Assume that h is sufficiently small, then, there exist two positive constants, C 1 and C 2 , which are independent of h = T N , such that Proof. Subtracting (21) from (19), we have e 2,0 = 0. We consider two cases.
(1) The first case is [qi] + 1 = i (i.e., when i ≤ 1 1−q ). Then, we have According to Lemma 4, we have e 2,i ≤ C 1 h . (2) The second case is [qi] + 1 ≤ i (i.e., when i > 1 1−q ). Then, we obtain Proof. By Theorems 3 and 4, the absolute error between y(t i ) and y i has the expression We obtain the conclusion of Theorem 5.

Extrapolation Method
In this section, we first describe the asymptotic error expansion and then present an extrapolation technique for achieving high precision. Finally, a posterior error estimate is derived.
Theorem 6. Let f (t), k 1 (t, s; y(s)), k 2 (t, s; y(s)) are four times continuously differentiable on I, D × R, D θ × R, respectively. Additionally, y(t) has continuous partial derivatives up to 3 on I and k i (t, s; y(s)) (i = 1, 2) satisfy Lipschitz conditions (2). There exist functionsŴ i (t)(i = 1, 2, 3) independent of h, such that we have the following asymptotic expansions: Proof. Assume that {Ŵ k (t), k = 1, 2, 3} satisfy the auxiliary delay equationŝ W k (t) = W k (t) + t 0 s λ k 1 (t, s; y(s))Ŵ k (s)ds + qt 0 s µ k 2 (t, s; y(s))Ŵ k (s)ds, andŴ k (t i ), i = 1, · · · , N satisfy the approximation equationŝ The analysis procedure is similar to the proof of Theorem 3. We obtain Then, we obtain According to Lemma 4, there exists a constant d such that The asymptotic expansion is From Theorem 6, we consider the Richardson extrapolation method to achieve higher accuracy.

Numerical Experiments
In this section, we illustrate the performance and accuracy of the quadrature method using the improved trapezoid formula. For ease of notation, we define where y h i is the approximate solution of Equation (1), y k,h i is the approximate solution of k-th extrapolation, E k,i is the absolute error between the exact solution and the approximate solution of k-th extrapolation when t = t i . The procedure was implemented in MATLAB. Example 1. Consider the following equation with T = 1, λ = − 1 2 , and q = 0.95. The exact solution is given by y(t) = t and f (t) is determined by the exact solution.
Applying the algorithm with N = 2 4 , 2 5 , 2 6 , 2 7 , 2 8 , the numerical results at t = 0.4 are presented in Table 1, the CPU time(s) are 0.34, 0.55, 0.98, 1.62, and 3.01 s, respectively. By comparing E h and E 1,i , we can observe that the accuracy was improved and the extrapolation algorithm was effective. In the third column, the rate values show that the convergence order was consistent with the theoretical analysis.
By applying the numerical method for N = 2 4 , 2 5 , 2 6 , 2 7 , and 2 8 , the obtained results at t = 0.4 are shown in Table 3. As λ was not equal to µ, we first applied the Richardson h 2+λ extrapolation, and then adopted the Richardson h 2+µ extrapolation. By comparing E h , E 1,i and E 2,i , these results verify the theoretical results, and we can see that the extrapolation improved the accuracy dramatically. When N = 8, 16, 32, 64, 128, the CPU time(s) are 1.43, 2.41, 3.99, 17.46, and 21.36 s, respectively. The exact solution and the approximation when N=8 are plotted in Figure 1.

Conclusions
In this paper, by using the improved trapezoidal quadrature formula and linear interpolation, we obtained the approximate equation for non-linear Volterra integral equations with vanishing delay and weak singular kernels. The approximate solutions were obtained by an iterative algorithm, which possessed a high accuracy order O(h 2+γ ). Additionally, we analyzed the existence and uniqueness of both the exact and approximate solutions. The significance of this work was that it demonstrated the efficiency and reliability of the Richardson extrapolation. The computational findings were compared with the exact solution: we found that our methods possess high accuracy and low computational complexity, and the results showed good agreement with the theoretical analysis. For future work, we can apply this method for solving two-dimensional delay integral equations.