**3. GMCk1**,**k2M 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.

#### *3.1. 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].

Consider the computation of

$$\mathbf{M}^{a,b}\_{\omega,n} := \int\_a^b \mathbf{K}(t\_n, \mathbf{s}) \mathbf{e}^{\mathbf{i}\omega \cdot \mathbf{y}(t\_n, \mathbf{s})} \boldsymbol{\Phi}(\mathbf{s}) d\mathbf{s}, \ n = 1, 2, \cdots, N. \tag{20}$$

When the phase has no stationary points, that is, *g* (*tn*,*s*) = 0 for any *<sup>s</sup>* <sup>∈</sup> [*a*, *<sup>b</sup>*], let {*ck*}*<sup>v</sup> <sup>k</sup>*=<sup>0</sup> be the equispaced nodes on the interval [*a*, *<sup>b</sup>*], that is, *ck* <sup>=</sup> *<sup>a</sup>* <sup>+</sup> *k v* (*<sup>b</sup>* <sup>−</sup> *<sup>a</sup>*) for *<sup>k</sup>* = 0, ··· , *<sup>v</sup>*. In addition, let {*mk*}*<sup>v</sup> <sup>k</sup>*=<sup>0</sup> denote a set of positive integers associated with nodes {*ck*}*<sup>v</sup> <sup>k</sup>*=0, which helps represent Hermite interplant later. Furthermore, define the function

$$\sigma\_k(s) = \begin{cases} \frac{K(t\_{n\prime}s)\phi(s)}{\mathcal{S}'(t\_{n\prime}s)}, & k=1, \\\frac{\sigma\_{k-1}^{\prime}(s)}{\mathcal{S}'(t\_{n\prime}s)}, & k \ge 2. \end{cases}$$

Then we can find a polynomial *<sup>p</sup>*(*s*) = *N*ˆ ∑ *<sup>q</sup>*=<sup>0</sup> *aqs <sup>q</sup>* with *<sup>N</sup>*<sup>ˆ</sup> = *v* ∑ *<sup>k</sup>*=<sup>0</sup> *mk* − 1 satisfying

$$\begin{cases} \begin{aligned} p(\mathcal{g}(c\_0)) &= \sigma\_1(c\_0) \\ \dots \\ p^{(m\_0-1)}(\mathcal{g}(c\_0)) &= \sigma\_{\mathfrak{m}\_0-1}(c\_0) \\ p(\mathcal{g}(c\_1)) &= \sigma\_1(c\_1) \\ \dots \\ p^{(m\_1-1)}(\mathcal{g}(c\_1)) &= \sigma\_{\mathfrak{m}\_1-1}(c\_1) \\ \dots \\ p(\mathcal{g}(c\_\mathbb{v})) &= \sigma\_1(c\_\mathbb{v}) \\ \dots \\ p^{(m\_\mathbb{v}-1)}(\mathcal{g}(c\_\mathbb{v})) &= \sigma\_{\mathfrak{m}\_\mathbb{v}-1}(c\_\mathbb{v}) \end{aligned} \end{cases}$$

With the coefficients *aq* by solving the above linear system, we can approximate M*a*,*<sup>b</sup> <sup>ω</sup>*,*<sup>n</sup>* by

$$\int\_{\mathcal{S}(t\_n,a)}^{\mathcal{S}(t\_n,b)} p(s) \mathbf{e}^{\mathrm{i}\omega s} ds = \sum\_{q=0}^{\hat{N}} a\_q \int\_{\mathcal{S}(t\_n,a)}^{\mathcal{S}(t\_n,b)} s^q \mathbf{e}^{\mathrm{i}\omega s} ds,$$

where *<sup>g</sup>*(*tn*,*b*) *<sup>g</sup>*(*tn*,*a*) *<sup>s</sup> q* ei*<sup>ω</sup><sup>s</sup> ds* can be calculated by incomplete Gamma function.

In the case of *g* (*tn*,*s*) = 0 for some *<sup>s</sup>* <sup>∈</sup> [*a*, *<sup>b</sup>*], suppose *<sup>s</sup>* = *<sup>a</sup>* without loss of generality. Then we insert the grid points

$$a, a + \frac{2^0}{\omega}, a + \frac{2}{\omega}, a + \frac{2^2}{\omega}, \dots, a + \frac{2^{\text{nr}}}{\omega}, b, 1$$

where *<sup>m</sup>* is the maximum integer less than log2 *<sup>ω</sup>*(*<sup>b</sup>* <sup>−</sup> *<sup>a</sup>*). 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.

#### *3.2. 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.

**Definition 7.** *Given any bivariate function <sup>g</sup>*(*t*,*s*) *defined on* [0, *<sup>T</sup>*] <sup>×</sup> [0, *<sup>T</sup>*], *suppose that <sup>g</sup>*(*t*,*s*) *has several stationary points <sup>ξ</sup>*1, ··· , *<sup>ξ</sup>nt over* [0, *<sup>T</sup>*] *for any fixed t*, *and*

$$\begin{cases} \mathbf{g}'(t, \xi\_1) = \dots = \mathbf{g}^{(r\_1)}(t, \xi\_1^{\mathbf{r}}) = 0, \mathbf{g}^{(r\_1+1)}(t, \xi\_1^{\mathbf{r}}) \neq 0 \\ \mathbf{g}'(t, \xi\_2^{\mathbf{r}}) = \dots = \mathbf{g}^{(r\_2)}(t, \xi\_2^{\mathbf{r}}) = 0, \mathbf{g}^{(r\_2+1)}(t, \xi\_2^{\mathbf{r}}) \neq 0, \\ \dots \\ \mathbf{g}'(t, \xi\_{n\_l}^{\mathbf{r}}) = \dots = \mathbf{g}^{(r\_{ll})}(t, \xi\_{n\_l}^{\mathbf{r}}) = 0, \mathbf{g}^{(r\_{n\_l}+1)}(t, \xi\_{n\_l}^{\mathbf{r}}) \neq 0. \end{cases}$$

*Let <sup>ρ</sup>*(*t*) = max *<sup>i</sup>*=1,··· ,*nt* {*ri*} *and r* = sup *<sup>t</sup>*∈[0,*T*] {*ρ*(*t*)}. *Then g*(*t*,*s*) *is said to be in* <sup>A</sup>(*r*).

Secondly, we give a slight extension of the classical van der Corput Lemma (see [4] p. 333).

**Lemma 8.** *Suppose that <sup>g</sup>*(*t*,*s*) ∈ A(*r*). *Moreover, suppose <sup>φ</sup>*(*s*) <sup>∈</sup> *<sup>C</sup>*1(*a*, *<sup>b</sup>*) *and <sup>φ</sup>* (*s*) *is integrable. We can conclude that*

$$\left| \int\_{a}^{b} \phi(s) \mathbf{e}^{\mathrm{i}\omega \cdot \mathbf{y} \cdot (t\_{n}, s)} ds \right| \leq \mathsf{C} \omega^{-1/(r+1)} . \, n = 1, 2, \cdots \text{,} N.$$

*Here the constant C is independent of ω*.

Finally, we are able to develop the convergence behavior of collocation polynomials computed by *GMCk*1,*k*<sup>2</sup> *M* in the highly oscillatory case.

**Theorem 9.** *Assume both of <sup>g</sup>*(*t*,*s*) ∈ A(*r*) *and <sup>f</sup> are sufficiently smooth. Then the numerical solutions derived from GMCk*1,*k*<sup>2</sup> *M for VIE* (1) *satisfy*

$$\max\_{t \in I\_h} \{|\mu(t) - \mu\_h(t)|\} = O(\omega^{-1/(r+1)}), \omega \to \infty. \tag{21}$$

**Proof.** To begin with, we explore the boundedness of the solution *<sup>u</sup>*(*t*) to Equation (1) and its derivative. By applying Picard iteration, we can rewrite *<sup>u</sup>*(*t*) as

$$
\mu = f + \sum\_{j=1}^{\infty} (\mathbb{K}^j f). \tag{22}
$$

Here **K** denotes the integral operator

$$(\mathbf{K}\boldsymbol{\phi})(t) := \int\_0^t \boldsymbol{K}(t, \mathbf{s}) \mathbf{e}^{\mathbf{i}\omega \cdot \mathbf{g}(t, \mathbf{s})} \boldsymbol{\phi}(\mathbf{s}) d\mathbf{s}$$

*Mathematics* **2020**, *8*, 2004

According to Lemma 8, we get that *<sup>u</sup>*(*t*) is bounded as *<sup>ω</sup>* <sup>→</sup> <sup>∞</sup>. On the other hand, the derivative can be rewritten by a direct calculation

$$\begin{split} u'(t) &= f'(t) + \sum\_{j=1}^{\infty} \left( \mathcal{K}(t,t) f(t) \mathbf{e}^{i\omega\chi(t,t)} (\mathbf{K}^{j-1} f)(t) + \mathbf{i}\omega \int\_{0}^{t} \mathcal{K}(t,s) f(s) g'(t,s) \mathbf{e}^{i\omega\chi(t,s)} (\mathbf{K}^{j-1} f)(s) ds \right) \\ &= f'(t) + \sum\_{j=1}^{\infty} \left( \mathcal{Z}\_{j} + \mathcal{Z}\mathcal{Z}\_{j} \right) \,, \end{split}$$

where

$$\mathcal{L}\_j := \mathcal{K}(t, t) f(t) \mathrm{e}^{\mathrm{i}\omega \cdot \mathbf{g}(t, t)}(\mathbf{K}^{j-1} f)(t),\\\mathcal{L}\mathcal{T}\_j := \mathrm{i}\omega \int\_0^t \mathcal{K}(t, s) f(s) \mathbf{g}'(t, s) \mathrm{e}^{\mathrm{i}\omega \cdot \mathbf{g}(t, s)}(\mathbf{K}^{j-1} f)(s) ds.$$

By letting *ω* → ∞, I*<sup>j</sup>* is bounded due to Lemma 8, and II*<sup>j</sup>* is bounded by noting that *g* (*t*,*s*) vanishes at *<sup>s</sup>* = 0.

When noting that the collocation error function *eh*(*t*) defined in the previous section satisfies

$$c\_{\rm li}(t) = \int\_0^t K(t, s) \mathbf{e}^{\rm i} \boldsymbol{\omega}^{\rm j}(t, s) \, d\boldsymbol{s} \,\,\, t \in \mathcal{X}\_{\rm li} \tag{23}$$

we obtain

$$(\mathbf{I} - h\mathbf{A}(\mathbf{1}:N, \mathbf{2}:N+1))\mathbf{e}\_h = \mathbf{R}\_\prime \tag{24}$$

.

where

$$\mathbf{e}\_{h} = \begin{pmatrix} c\_{h}(t\_{1}) \\ c\_{h}(t\_{2}) \\ \dots \\ c\_{h}(t\_{N}) \end{pmatrix}, \mathbf{R} = \begin{pmatrix} h^{k\_{1}+k\_{2}+3} \operatorname{Res}\_{0,0}^{1,0} \\ h^{k\_{1}+k\_{2}+3} \operatorname{Res}\_{0,0,0}^{0,0} \\ \dots \\ h^{k\_{1}+k\_{2}+3} \operatorname{Res}\_{0,0,0}^{k\_{1},0} + h^{k\_{1}+k+2+3} \sum\_{l=k\_{1}}^{N-k\_{2}-1} \operatorname{Res}\_{0,N,N}^{1,l} + h^{k\_{1}+k\_{2}+3} \operatorname{Res}\_{0,n,N-k\_{2}-1}^{k\_{2}+1,N-k\_{1}-1} \end{pmatrix}$$

Since both of *<sup>u</sup>*(*t*) and *uh*(*t*) are bounded as *<sup>ω</sup>* <sup>→</sup> <sup>∞</sup>, employing Lemma <sup>8</sup> implies

$$|\mathsf{M}\mathsf{OM}\_{a,c,i}^{b,d}| \leq \mathsf{C}\omega^{-1/(r+1)}, |\mathsf{R}\mathsf{ES}\_{a,c,i}^{b,d}| \leq \mathsf{C}\omega^{-1/(r+1)}$$

Hence for fixed stepsize *<sup>h</sup>*, **<sup>I</sup>** <sup>−</sup> *<sup>h</sup>***A**(<sup>1</sup> : *<sup>N</sup>*, 2 : *<sup>N</sup>* + <sup>1</sup>) is invertible for sufficiently large *<sup>ω</sup>*, and we can represent **e***<sup>h</sup>* by

$$\mathbf{e}\_{\mathrm{lt}} = (\mathbf{I} - h\mathbf{A}(\mathbf{1}:\mathcal{N}, \mathbf{2}:\mathcal{N} + \mathbf{1}))^{-1}\mathbf{R}.$$

By noting that maximum of **<sup>R</sup>** goes to 0 with a speed of *<sup>O</sup>*(*ω*−1/(*r*+<sup>1</sup>)) as *<sup>ω</sup>* goes to <sup>∞</sup>, 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*,

$$u(t) + \int\_0^t e^{i\omega(t-s)} u(s)ds = \mathfrak{e}^t, t \in [0, 1]. \tag{25}$$

*The exact solution is u*(*t*) = *<sup>t</sup>* 0 (−*ce<sup>s</sup>* )*e* <sup>−</sup>*csds* + <sup>1</sup> *e ct*, *<sup>c</sup>* = *<sup>i</sup><sup>ω</sup>* <sup>−</sup> 1.

In Figure 4, we plot the scaled infinite norm of absolute error according to the corresponding order by letting *<sup>N</sup>* = 32, and *<sup>ω</sup>* 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.

**Figure 4.** *GMCk*1,*k*<sup>2</sup> *M* for the highly oscillatory problem.

#### **4. 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.

**Author Contributions:** H.C. and J.M. conceived and designed the experiments; H.C. and L.L. performed the experiments; H.C. and J.M. analyzed the data; J.M.contributed reagents/materials/analysis tools; H.C. and J.M. wrote the paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by National Natural Science Foundation of China (number 11901133) and the Science and Technology Foundation of Guizhou Province (number QKHJC[2020]1Y014).

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

#### **Abbreviations**

The following abbreviations are used in this manuscript:

