Analysis of Generalized Multistep Collocation Solutions for Oscillatory Volterra Integral Equations

Abstract: In this work, we introduce a class of generalized multistep collocation methods for solving oscillatory Volterra integral equations, and study two kinds of convergence analysis. The error estimate with respect to the stepsize is given based on the interpolation remainder, and the nonclassical convergence analysis with respect to oscillation is developed by investigating the asymptotic property of highly oscillatory integrals. Besides, the linear stability is analyzed with the help of generalized Schur polynomials. Several numerical tests are given to show that the numerical results coincide with our theoretical estimates.


Introduction
In many practical problems, such as epidemic diffusion, population dynamics and reaction processes, one may usually come across a class of Volterra integral equations (VIEs) (see [1] and references therein). Noting that most VIEs cannot be solved in closed forms, many researchers have made contributions to the numerical approaches to VIEs.
Particularly, the study of numerical solutions to VIEs with highly oscillatory Fourier or Bessel kernels has attracted much attention during the past decade. In [2], Xiang and Brunner first investigated Filon collocation approximations to highly oscillatory VIEs by employing the asymptotic property of oscillatory integrals. They found that errors of Filon collocation solutions decayed fast as the frequency increased. The third author presented an optimal convergence order for the direct Filon collocation solution to the first kind of oscillatory VIE arising in acoustic scattering in [3]. The convergence behavior of such kinds of numerical approaches was able to be revealed with the help of the detailed study of the remainder for the error function. Besides, it is noted that numerical analysis with respect to the frequency, which is usually done by solving error equations and extending van der Corput lemma (see [4] p. 333), is able to detect the ability of the numerical method to solve highly oscillatory VIEs. With these techniques in mind, several authors made great contributions to numerical solutions to highly oscillatory VIEs. For example, Galerkin and collocation solutions for VIEs with highly oscillatory trigonometric kernels were investigated in [5,6], highly oscillatory VIEs with weakly singular kernels were studied in [7], the Hermite-type Filon collocation method was presented in [8], and Clenshaw-Curtis-Filon qudrature for Cauchy singular integral equations was investigated in [9].
In this work, we consider the numerical computation of the following second-kind oscillatory VIE: where K(t, s), g(t, s) and f (t) are sufficiently smooth, u(t) is unknown, and ω denotes the oscillation parameter. When ω = 0, Equation (1) reduces to the classical VIEs. In the case of ω 1, the kernel in Equation (1) is highly oscillatory, and special quadrature rules should be employed in practical computation.
In the remaining part, we are restricted to the following problems. In the forthcoming section, we first develop a class of generalized multistep collocation methods (GMC k 1 ,k 2 M) for Equation (1) with non-oscillatory kernels, that is, ω = 0. Then, classical convergence analysis and linear stability analysis are implemented. In the third section, we study the numerical solution to VIE (1) when the kernel changes rapidly, that is, ω 1, and present the frequency-explicit convergence analysis. Some concluding remarks are given in Section 4.

GMC k 1 ,k 2 M in the Case of ω = 0
Frequently-used approaches for VIEs include collocation methods [10], the spectral collocation method [11,12], the spectral Galerkin method [13,14], the Nyström method [15,16], and so on. Among these numerical formulae, the collocation-based approach is one of the most important tools. In general, the collocation solution is obtained by making the polynomial or piecewise polynomial satisfy the collocation equation. For one-step collocation methods, one can find detailed analysis in [10]. To increase the convergence rate without adding collocation points, Conte and Paternoster studied multistep collocation solutions with the help of employing approximations to numerical solutions in computed steps in [17]. However, multistep methods usually tend to be unstable. Fazeli et al. further investigated the stability of multistep collocation methods in [18], and found some super implicit collocation solutions with wide stability regions. On the other hand, inspired by the study of boundary value methodology for solving ODE (see [19]), several authors made contributions to boundary value solutions to Volterra functional equations [20][21][22]. Based on interpolation outside the current subinterval and approximated end values, the third author and Xiang devised CBVM for second-kind VIEs in [22]. Furthermore, the third author extended said kind of methodology to VIEs with weakly singular kernels by employing the fractional polynomial interplant in [23], and the block CBVM for the first-kind VIE was investigated in [24]. In this section we first investigate the construction of GMC k 1 ,k 2 M with the help of local polynomial interpolation. Then, the convergence and linear stability analysis of GMC k 1 ,k 2 M are considered.

Discretization of VIE
Let the interval [0, T] be divided uniformly, that is, For the first k 1 subintervals, that is, for any t ∈ [t 0 , t k 1 ], the collocation polynomial is represented by For In the last subinterval [t N−k 2 , t N ], we rewrite u h (t) as Finally, the collocation equation follows: A direct calculation leads to Denoting we have for k = −k 1 , · · · , k 2 + 1, Now we are able to rewrite Equation (6) in the closed form: where I denotes the identity matrix, A = , · · · , f (t N )] T . By employing proper numerical integration approaches such as Clenshaw-Curtis quadrature and applying iterative solvers to Equation (9), we are able to obtain the collocation solution at the grid.

Convergence Analysis with Respect to Stepsize
Now we turn to studying the convergence behavior of the piecewise collocation polynomial computed by Equation (9). Firstly, we revisit some helpful results from approximation theory. Lemma 1 ([10] p. 43). Consider the following assumption.
• Defining abscissa a ≤ ξ 1 < ... < ξ m ≤ b, we obtain the error between f (x) and the Lagrange interpolation polynomial of degree m − 1 with respect to the given points {ξ j }.
where L j (x) denotes Lagrange basis.
Then we can represent the error function ε m ( f ; x) as follows.
Here the kernel function κ d (x, t) can be obtained by Existing studies show that we cannot compute collocation boundary value solutions by recurrences. All numerical values should be computed simultaneously through solving linear systems. Note that the element of hA(1 : N, 2 : N + 1) is bounded by whereK denotes the maximum of the kernel function K(t, s), and the above inequality is derived from the Lesbegue constant of the polynomial interpolant (see [25]). We obtain hA(1 : N, 2 : N + 1) < 1 whenever h < (K2 k 1 +k 2 +4 ) −1 , which enables us to compute det(I − hA(1 : N, 2 : N + 1)) = 0 by Gaussian elimination, as is done in [22]. Therefore, the well-posedess of the solution computed by GMC k 1 ,k 2 M is guaranteed. It is noted that when we encounter stiff problems, the maximumK may be large, which implies we have to apply a particularly small stepsize h and restricts the application of the collocation method. However, due to the compactness of Volterra integral operator, the spectrum of hA(1 : N, 2 : N + 1) will be found in the neighborhood of 0 with a tolerance stepsize, and the multistep collocation method is feasible in practical uses. In Figure 1, we show the discretized spectrum of hA(1 : N, 2 : N + 1) by considering the kernel function K(t, s) = 50e iω(t−s) with the maximum 50. It can be seen that eigenvalues are bounded by the unit circle when the stepsize decreases to 1/64, which guarantees the solvability of the linear system. Furthermore, we arrive at the following theorem by employing Lemmas 1 and 2.
Theorem 3. Suppose that K(t, s), f (t) in VIE (1) are sufficiently smooth, that is, K ∈ C k 1 +k 2 +2 (D) and f ∈ C k 1 +k 2 +2 (I). Furthermore, let u h (t) denote the collocation polynomial computed by GMC k 1 ,k 2 M with a stepsize h. Then the collocation error e h (t) = u(t) − u h (t) in the collocation grid is bounded by where the constant C is independent of the stepsize h but depends on T.
Proof. Note that by we obtain the collocation error function e h (t) satisfying Let Furthermore, by letting we obtain v n := n−k 1 ,k 1 k 1 ,n,0 + h k 1 +k 2 +3 RES n−k 1 ,k 1 k 1 ,n,0 , n = 1, · · · , k 1 , Suppose that MOM b,d a,c,i and RES b,d a,c,i are bounded by the constant B. It is easily noted from Equations (8) and (13) Hence, we have According to Lemma 2, we have or equivalently, for sufficiently small stepsize h.
In this example, we test the performance of GMC k 1 ,k 2 M. We mainly focus on two terms of data, the maximum of error functions (INAE), and the convergence order. Computed results are shown in Tables 1-3. It can be seen from these tables that as the quantity of nodes increases, absolute errors decay fast, and as k 1 and k 2 get bigger, the convergence order enlarges. Besides, numerical results illustrate that GMC k 1 ,k 2 M achieves the expected order of the estimate given in Theorem 3.

Remark 1.
When numerical solutions of evolution equations are considered, Courant proposes that the combination of a consistent and stable numerical approach led to its convergence, which contributes to the foundation of classical numerical analysis theory of numerical studies on differential equations. On the other hand, the above convergence analysis is based on a fixed integration interval [0, T], which differs from the convergence analysis for evolution problems where we usually consider the case of T → ∞. In addition, it should be noted that the convergence result in Theorem 3 does not guarantee a feasible approximation in practical computation for long-time integration, especially when we are met with stiff problems. Therefore, we give linear stability analysis of the presented collocation method in the forthcoming subsection.

Linear Stability Analysis
For a long-time integration problem, round-off errors may dramatically affect the numerical solution. In this subsection, we analyze the collocation solution's linear stability originating from the study of numerical solutions of ordinary differential equations, where one usually considers the test equation y (t) = λy(t), Re(λ) < 0.
Particularly, Brugnano and Trigiante investigated multistep methods for solving differential problems with the above scalar equation in [19]. For the general linear multistep formula k ∑ j=0 α j y n+j − hλ k ∑ j=0 β j y n+j = 0, we can introduce two polynomials and define the associated characteristic polynomial π(z, q) = ρ(z) − qσ(z) with q = hλ. When π(z, q) is a Schur polynomial for fixed q, the method is absolutely stable at q. For the moment the definition of the region of absolute stability is Since both of discretization of ODE and VIE result in difference equations, we can investigate the generalized multistep collocation method with the help of stability studies of ODE. Consider the following test equation: We turn to study the linear stability of the collocation solution by investigating Equation (15). By applying GMC k 1 ,k 2 M we have Next, noting the difference between y j and y j−1 in Equation (15) leads to Then the characteristic polynomial is defined by Before investigating the linear stability region, we introduce some helpful definitions and theorems in the version of GMC k 1 ,k 2 M.

Definition 4 ([19]
). For any complex number q := hλ, if the collocation solution u h to Equation (15) computed by GMC k 1 ,k 2 M goes to 0 as T goes ∞ for fixed stepsize, then GMC k 1 ,k 2 M is said to be absolutely stable at q.

Definition 5 ([19]
). For any z ∈ S, if GMC k 1 ,k 2 M is absolutely stable at z, then the set S is said to be the linear stability region of GMC k 1 ,k 2 M. Particularly, if the left part of the complex plane is contained in S, then GMC k 1 ,k 2 M is said to be A−stable.

Theorem 6 ([19]). For any complex number q, if roots of Equation (18) satisfy
then GMC k 1 ,k 2 M is stable at q.
By a direct calculation, we find that roots of π k 1 ,k 2 (z, q) do not satisfy the condition given in Theorem 6 in the case of k 1 = k 2 . Hence, the region of stability cannot be shown. In Figures 2 and  3, we list the boundary locus corresponding to various multistep collocation methods with k 1 = k 2 , where the boundary Γ is defined by It can be seen that these trajectories are Jordan curves, which implies Γ is the boundary of corresponding absolute stability region. The stability region in Figure 2 is the part outside the boundary curves, while that in Figure 3 is the inside part. Therefore, we can conclude that GMC k 1 ,k 2 M has wide stability region in the case of k 2 > k 1 . In addition, the boundary trajectories of GMC k1,k2 M and GMC k2,k1 M are symmetric with respect to virtual axis.

GMC k 1 ,k 2 M in the Case of ω 1
When the oscillation parameter ω 1 in Equation (1), classical quadrature usually results in time-consuming algorithms. Hence, we first give an efficient numerical approach for moments in Equation (6) in this section. Then the frequency-explicit convergence analysis is presented.

Fast Calculation of Moments
Numerical studies of highly oscillatory integrals (HOIs) have been intensively focused on in the past few decades. High-order algorithms, such as Filon-type quadrature [26], Levin quadrature [27], and the numerical steepest decent method [28], have been proposed. In this subsection, we consider a composite quadrature rule based on Xiang's modified Filon-type quadrature developed in [29].
When the phase has no stationary points, that is, g (t n , s) = 0 for any s ∈ [a, b], let {c k } v k=0 be the equispaced nodes on the interval [a, b], that is, In addition, let {m k } v k=0 denote a set of positive integers associated with nodes {c k } v k=0 , which helps represent Hermite interplant later. Furthermore, define the function Then we can find a polynomial p(s) =N ∑ q=0 a q s q withN = v ∑ k=0 m k − 1 satisfying With the coefficients a q by solving the above linear system, we can approximate M a,b ω,n by g(t n ,b) g(t n ,a) p(s)e iωs ds =N ∑ s q e iωs ds can be calculated by incomplete Gamma function.
In the case of g (t n , s) = 0 for some s ∈ [a, b], suppose s = a without loss of generality. Then we insert the grid points where m is the maximum integer less than log 2 ω(b − a). Integration with Xiang's Filon quadrature in each subintervals results in the composite Filon quadrature. It is noted that the integral over the first interval is non-oscillatory and we can employ classical quadrature such as Gauss or Clenshaw-Curtis instead to avoid the stationary problem.

Convergence Analysis with Respect to the Frequency
Collocation methods with high-order quadrature usually lead to a class of fascinating algorithms, which are able to provide high-precision collocation solutions in the case of high frequency. In this subsection, we consider the general oscillator and investigate the convergence analysis for multistep collocation solutions, where the convergence order is represented by the frequency parameter ω.
Firstly, let us restrict ourselves to considering the following set of functions.
Secondly, we give a slight extension of the classical van der Corput Lemma (see [4] p. 333).
Here the constant C is independent of ω.
Finally, we are able to develop the convergence behavior of collocation polynomials computed by GMC k 1 ,k 2 M in the highly oscillatory case. Theorem 9. Assume both of g(t, s) ∈ A(r) and f are sufficiently smooth. Then the numerical solutions derived from GMC k 1 ,k 2 M for VIE (1) satisfy Proof. To begin with, we explore the boundedness of the solution u(t) to Equation (1) and its derivative. By applying Picard iteration, we can rewrite u(t) as Here K denotes the integral operator (Kφ)(t) := t 0 K(t, s)e iωg(t,s) φ(s)ds According to Lemma 8, we get that u(t) is bounded as ω → ∞. On the other hand, the derivative can be rewritten by a direct calculation By letting ω → ∞, I j is bounded due to Lemma 8, and II j is bounded by noting that g (t, s) vanishes at s = 0.
When noting that the collocation error function e h (t) defined in the previous section satisfies we obtain where Since both of u(t) and u h (t) are bounded as ω → ∞, employing Lemma 8 implies Hence for fixed stepsize h, I − hA(1 : N, 2 : N + 1) is invertible for sufficiently large ω, and we can represent e h by e h = (I − hA(1 : N, 2 : N + 1)) −1 R.
By noting that maximum of R goes to 0 with a speed of O(ω −1/(r+1) ) as ω goes to ∞, we obtain the estimate (21).
In the following example, we test the convergence rate of GCM 1,2 M in the case of high frequency.

Example 2.
In this example, we solve the following VIE with GCM 1,2 M, The exact solution is u(t) = t 0 (−ce s )e −cs ds + 1 e ct , c = iω − 1.
In Figure 4, we plot the scaled infinite norm of absolute error according to the corresponding order by letting N = 32, and ω varies from 50 to 1000. The left part shows the infinite norm of the error and the right part shows the absolute error scaled by corresponding rates. It can be seen that the increase of the frequency parameter ω makes the absolute error get smaller. This indicates as the kernel becomes more highly oscillatory, computed approximation becomes more accurate. Considering the right part of Figure 4, we find that when the frequency parameter ω reaches 150, the curve turns to a horizontal straight line, which is in agreement with the estimate given in Theorem 9.

Final Remark
For VIEs with oscillatory and non-oscillatory kernels, we have investigated the generalized multistep collocation solution to VIE (1). Detailed convergence properties with respect to the stepsize and oscillation are presented. Noting that the new approach coupled with mild composite oscillatory quadrature rules is able to produce high-order approximation as the frequency goes to infinity, we could expect it is valuable to conduct further studies in related highly oscillatory problems, such as oscillatory Riemann-Hilbert problems, spectral calculation of oscillatory Fredholm operators, and so on.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

VIE
Volterra integral equation GMC k 1 ,k 2 M generalized multistep collocation method ODE ordinary differential equation CBVM collocation boundary value method HOI highly oscillatory integral