**2. GMCk1**,**k2M 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–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 *GMCk*1,*k*<sup>2</sup> *M* with the help of local polynomial interpolation. Then, the convergence and linear stability analysis of *GMCk*1,*k*<sup>2</sup> *M* are considered.

#### *2.1. Discretization of VIE*

Let the interval [0, *T*] be divided uniformly, that is,

$$X\_h = \{ t\_{\bar{j}} \colon t\_{\bar{j}} = jh, \ j = 0, 1, \dots, \text{''}, \text{''} = T/h \} \dots$$

Then define local basic functions

$$\phi\_j^{k\_1,k\_2}(s) = \prod\_{i=-k\_1, j \neq j}^{k\_2+1} \frac{s-i}{j-i}, j = -k\_1, \dots, k\_2+1,\tag{2}$$

For the first *<sup>k</sup>*<sup>1</sup> subintervals, that is, for any *<sup>t</sup>* <sup>∈</sup> [*t*0, *tk*<sup>1</sup> ], the collocation polynomial is represented by

$$u\_h(t\_{k\_1} + sh) = \sum\_{i=-k\_1}^{k\_2+1} y\_{k\_1+i} \phi\_i^{k\_1,k\_2}(s), \text{ s} \in (-k\_1, 0], \tag{3}$$

For *<sup>k</sup>*<sup>1</sup> <sup>≤</sup> *<sup>n</sup>* <sup>≤</sup> *<sup>N</sup>* <sup>−</sup> *<sup>k</sup>*<sup>2</sup> <sup>−</sup> 1, *uh*(*t*) over the interval [*tn*, *tn*+1] is rewritten as

$$u\_h(t\_h + sh) = \sum\_{i=-k\_1}^{k\_2+1} y\_{n+i} \phi\_i^{k\_1, k\_2}(s), \; s \in (0, 1], \tag{4}$$

In the last subinterval [*tN*−*k*<sup>2</sup> , *tN*], we rewrite *uh*(*t*) as

$$u\_h(t\_{N-k\_2-1} + sh) = \sum\_{i=-k\_1}^{k\_2+1} y\_{N-k\_2-1+i} \phi\_i^{k\_1,k\_2}(s), \; s \in (1, k\_2+1]. \tag{5}$$

Finally, the collocation equation follows:

$$
\mu\_h(t\_n) = f(t\_n) + \int\_0^{t\_n} K(t\_n, s) \mu\_h(s) ds, \ t\_n \in X\_h. \tag{6}
$$

A direct calculation leads to

$$y\_n - f(t\_n) = \begin{cases} h \int\_{-k\_1}^{n-k\_1} K(t\_n, t\_{k\_1} + sh) \left( \sum\_{i=-k\_1}^{k\_1 + 1} y\_{i, t\_i} \phi\_i^{k\_1, k\_2}(s) \right) ds, & n = 1, \dots, k\_1, \\\ h \int\_{-k\_1}^0 K(t\_n, t\_{k\_1} + sh) \left( \sum\_{i=-k\_1}^{k\_1 + 1} y\_{i, t\_i} \phi\_i^{k\_1, k\_2}(s) \right) ds \\\ + \begin{aligned} h \int\_{-k\_1 + 1}^0 \int\_0^1 K(t\_n, t\_{j-1} + sh) \left( \sum\_{i=-k\_1 - 1}^{k\_1 + 1} y\_{j-1, i} \phi\_i^{k\_1, k\_2}(s) \right) ds, & n = k\_1 + 1, \dots, N - k\_2, \\\ h \int\_{-k\_1}^0 K(t\_n, t\_{k\_1} + shh) \left( \sum\_{i=-k\_1 - 1}^{k\_1 + 1} y\_{i, t\_i} \phi\_i^{k\_1, k\_2}(s) \right) ds \\\ + \begin{aligned} h \sum\_{j=-k\_1 + 1}^{N - k\_1 - 1} \int\_0^1 K(t\_n, t\_{j-1} + sh) \left( \sum\_{i=-k\_1 - 1}^{k\_1 + 1} y\_{j-i, t\_i} \phi\_i^{k\_1, k\_2}(s) \right) ds \\\ + h \int\_{1}^{n - N + k\_2 + 1} K(t\_n, t\_{N - k\_2 - 1} + sh) \sum\_{i=-k\_1 - 1}^{k\_1 + 1} y\_{N - k\_2 - 1, i} \phi\_i^{k\_1, k\_2}(s) ds, & n = N - k\_2 + 1, \dots, N. \end{cases} \end{cases} \tag{7}$$

Denoting

$$\text{MO}\mathbf{M}\_{a,c,i}^{b,d} = \int\_{a}^{b} K(t\_c, t\_d + sh)\phi\_i^{k\_1,k\_2}(s)ds,\tag{8}$$

we have for *<sup>k</sup>* <sup>=</sup> <sup>−</sup>*k*1, ··· , *<sup>k</sup>*<sup>2</sup> <sup>+</sup> 1,

$$A\_{k}^{initial} = (a\_{i,j}^{(k)}) = \begin{cases} \text{MOM}\_{-k\_1, n, j-k\_1-1}^{0, k\_1} & i \le N, j = k + k\_1 + 1, \\ 0, & \text{others}, \end{cases}$$

$$A\_{k}^{main} = (b\_{i,j}^{(k)}) = \begin{cases} \text{MOM}\_{0, n, k}^{1, j-1}, & k\_1 < i \le N, i + k \le j \le i + k + \min\{0, n - N + k\_2 + 1\}, \\ 0, & \text{others}, \end{cases}$$

$$A\_{k}^{end} = (c\_{i,j}^{(k)}) = \begin{cases} \text{MOM}\_{i, n, j-N + k\_2 + 1, N - k\_2 - 1}^{n-N + k\_2 + 1, N - k\_2 - 1}, & i > N - k\_2, j = N - k\_2 + k\_r \\ 0, & \text{others}. \end{cases}$$

Now we are able to rewrite Equation (6) in the closed form:

$$(\mathbf{I} - h\mathbf{A}(\mathbf{1}:N, \mathbf{2}:N+1))\mathbf{Y} = \mathbf{F} + h\mathbf{y}\_0\mathbf{A}(\mathbf{1}:N, \mathbf{1}),\tag{9}$$

where **<sup>I</sup>** denotes the identity matrix, **<sup>A</sup>** = *<sup>k</sup>*2+<sup>1</sup> ∑ *<sup>k</sup>*=−*k*<sup>1</sup> *Ainitial <sup>k</sup>* + *<sup>k</sup>*2+<sup>1</sup> ∑ *<sup>k</sup>*=−*k*<sup>1</sup> *Amain <sup>k</sup>* + *<sup>k</sup>*2+<sup>1</sup> ∑ *<sup>k</sup>*=−*k*<sup>1</sup> *Aend <sup>k</sup>* , **<sup>Y</sup>** <sup>=</sup>

[*y*1, *<sup>y</sup>*2, ··· , *yN*] *<sup>T</sup>*, and **<sup>F</sup>** = [ *<sup>f</sup>*(*t*1), *<sup>f</sup>*(*t*2), ··· , *<sup>f</sup>*(*tN*)]*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.

#### *2.2. 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 <sup>a</sup>* <sup>≤</sup> *<sup>ξ</sup>*<sup>1</sup> <sup>&</sup>lt; ... <sup>&</sup>lt; *<sup>ξ</sup><sup>m</sup>* <sup>≤</sup> *<sup>b</sup>*, *we obtain the error between <sup>f</sup>*(*x*) *and the Lagrange interpolation polynomial of degree m* − 1 *with respect to the given points* {*ξj*}.

$$\varepsilon\_{\mathfrak{m}}(f; \mathfrak{x}) = f(\mathfrak{x}) - \sum\_{j=1}^{m} L\_j(\mathfrak{x}) f(\mathfrak{x}\_j), \ \mathfrak{x} \in [a, b]\_{\mathfrak{m}}$$

*where Lj*(*x*) *denotes Lagrange basis.*

• *Letting* <sup>1</sup> <sup>≤</sup> *<sup>d</sup>* <sup>≤</sup> *<sup>m</sup>*, *we suppose f*(*x*) *belongs to the space Cd*[*a*, *<sup>b</sup>*].

*Then we can represent the error function <sup>ε</sup>m*(*<sup>f</sup>* ; *<sup>x</sup>*) *as follows.*

$$\kappa\_{\mathfrak{M}}(f; \mathfrak{x}) = \int\_{a}^{b} \kappa\_{d}(\mathfrak{x}, t) f^{(d)}(t) dt, \ \mathfrak{x} \in [a, b]. \tag{10}$$

*Here the kernel function <sup>κ</sup>d*(*x*, *<sup>t</sup>*) *can be obtained by*

$$\kappa\_d(\mathbf{x}, t) := \frac{1}{(d - 1)!} \left\{ (\mathbf{x} - t)\_+^{d - 1} - \sum\_{j = 1}^m L\_j(\mathbf{x}) (\boldsymbol{\xi}\_j - t)\_+^{d - 1} \right\},$$

*and*

$$(\mathbf{x} - t)^p\_+ := \begin{cases} \mathbf{0}, & \mathbf{x} < t, \\ (\mathbf{x} - t)^p, & \mathbf{x} \ge t. \end{cases}$$

**Lemma 2** ([10] p. 81)**.** *Suppose that there a sequence* {*ki*} *with ki* ≥ 0 *and another sequence* { *<sup>i</sup>*} *with* <sup>0</sup> ≤ *ρ*0. *Moreover,* {*ki*} *and* { *<sup>i</sup>*} *satisfy*

$$
\mathfrak{e}\_{\mathcal{U}} \le \rho\_0 + \sum\_{i=0}^{n-1} q\_i + \sum\_{i=0}^{n-1} k\_i \mathfrak{e}\_{i\prime} \ n \ge 1\_{\prime \prime}
$$

*with ρ*<sup>0</sup> ≥ 0, *qi* ≥ 0, *i* ≥ 0. *Then*

$$\epsilon\_n \le \left(\rho\_0 + \sum\_{i=0}^{n-1} q\_i\right) \mathbf{e}^{\sum\_{i=0}^{n-1} k\_i} \text{ } n \ge 1.$$

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 *<sup>h</sup>***A**(1 : *<sup>N</sup>*,2: *<sup>N</sup>* + <sup>1</sup>) is bounded by

$$\left\| h(k\_1 + k\_2 + 1)\bar{\mathcal{K}} \left\| \sum\_{i=-k\_1}^{k\_2+1} \phi\_i^{k\_1,k\_2}(t) \right\|\_{\infty} \leq h\bar{\mathcal{K}} 2^{k\_1 + k\_2 + 4}, \bar{\mathcal{K}}$$

where *<sup>K</sup>*¯ denotes the maximum of the kernel function *<sup>K</sup>*(*t*,*s*), and the above inequality is derived from the Lesbegue constant of the polynomial interpolant (see [25]). We obtain *<sup>h</sup>***A**(<sup>1</sup> : *<sup>N</sup>*, 2 : *<sup>N</sup>* + <sup>1</sup>) <sup>&</sup>lt; <sup>1</sup> whenever *<sup>h</sup>* <sup>&</sup>lt; (*K*¯2*<sup>k</sup>*1+*k*2+4)<sup>−</sup>1, which enables us to compute det(**<sup>I</sup>** <sup>−</sup> *<sup>h</sup>***A**(<sup>1</sup> : *<sup>N</sup>*, 2 : *<sup>N</sup>* <sup>+</sup> <sup>1</sup>)) <sup>=</sup> 0 by Gaussian elimination, as is done in [22]. Therefore, the well-posedess of the solution computed by *GMCk*1,*k*<sup>2</sup> *<sup>M</sup>* is guaranteed. It is noted that when we encounter stiff problems, the maximum *<sup>K</sup>*¯ 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 *<sup>h</sup>***A**(<sup>1</sup> : *<sup>N</sup>*, 2 : *<sup>N</sup>* + <sup>1</sup>) 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 *<sup>h</sup>***A**(<sup>1</sup> : *<sup>N</sup>*, 2 : *<sup>N</sup>* + <sup>1</sup>) by considering the kernel function *<sup>K</sup>*(*t*,*s*) = 50ei*ω*(*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.

**Figure 1.** The spectrum of *<sup>h</sup>***A**(1 : *<sup>N</sup>*,2: *<sup>N</sup>* + <sup>1</sup>) for various stepsizes *<sup>h</sup>*.

Furthermore, we arrive at the following theorem by employing Lemmas 1 and 2.

**Theorem 3.** *Suppose that <sup>K</sup>*(*t*,*s*), *<sup>f</sup>*(*t*) *in VIE* (1) *are sufficiently smooth, that is, <sup>K</sup>* <sup>∈</sup> *<sup>C</sup>k*1+*k*2+2(*D*) *and <sup>f</sup>* <sup>∈</sup> *<sup>C</sup>k*1+*k*2+2(*I*). *Furthermore, let uh*(*t*) *denote the collocation polynomial computed by GMCk*1,*k*<sup>2</sup> *<sup>M</sup> with a stepsize h*. *Then the collocation error eh*(*t*) = *<sup>u</sup>*(*t*) <sup>−</sup> *uh*(*t*) *in the collocation grid is bounded by*

$$\max\_{t \in X\_h} |e(t)| \le Ch^{k\_1 + k\_2 + 2},\tag{11}$$

*where the constant C is independent of the stepsize h but depends on T.*

**Proof.** Note that by

$$u(t) = f(t) + \int\_0^t K(t, s)u(s)ds, \ t \in [0, T]\_{\prime}$$

and

$$\mu\_h(t) = f(t) + \int\_0^t K(t,s)\mu\_h(s)ds, \ t \in X\_{h\nu}$$

we obtain the collocation error function *eh*(*t*) satisfying

$$
\sigma\_h(t) = \int\_0^t K(t, s)\sigma\_h(s)ds, \; t \in X\_{\text{lt}}.\tag{12}
$$

Let

$$R\_{\mathfrak{n}}^{k\_1,k\_2}(v) := \int\_0^{k\_1+k\_2+1} \kappa\_{k\_1}^{k\_2}(v,z) u^{(k\_1+k\_2+2)}(t\_{\mathfrak{n}}+zh) dz\_{\mathfrak{n}}$$

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

where

$$\kappa\_{k\_1}^{k\_2}(\upsilon, z) := \frac{1}{(k\_1 + k\_2 + 1)!} \left( (\upsilon - z)\_+^{k\_1 + k\_2 + 1} - \sum\_{j = -k\_1}^{k\_2 + 1} \phi\_j^{k\_1 k\_2}(\upsilon)(j - z)\_+ \right).$$

For *<sup>t</sup>* <sup>∈</sup> [0, *tk*<sup>1</sup> ], we have

$$c\_h(t) = c\_h(t\_{k\_1} + sh) = \sum\_{i=-k\_1}^{k\_2+1} v\_{k\_1+i} \phi\_i^{k\_1,k\_2}(s) + h^{k\_1+k\_2+2} R\_0^{k\_1,k\_2}(s), \; s \in [-k\_1, 0],$$

where *vn* :<sup>=</sup> *eh*(*tn*). For *<sup>t</sup>* <sup>∈</sup> [*tj*, *tj*+1], *<sup>j</sup>* <sup>=</sup> *<sup>k</sup>*1, ··· , *<sup>N</sup>* <sup>−</sup> *<sup>k</sup>*<sup>2</sup> <sup>−</sup> 1, we have

$$\varepsilon\_h(t) = \varepsilon\_h(t\_j + sh) = \sum\_{i=-k\_1}^{k\_2+1} \upsilon\_{j+i} \phi\_i^{k\_1 k\_2}(s) + h^{k\_1+k\_2+2} R\_j^{k\_1, k\_2}(s), \; s \in [0, 1].$$

For *<sup>t</sup>* <sup>∈</sup> [*tN*−*k*<sup>2</sup> , *tN*], we have

$$c\_h c\_h(t) = c\_h(t\_{N-k\_2-1} + sh) = \sum\_{i=-k\_1}^{k\_2+1} v\_{N-k\_2+i} \phi\_i^{k\_1,k\_2}(s) + h^{k\_1+k\_2+2} R\_{N-k\_2}^{k\_1,k\_2}(s), \; s \in [1, k\_2+1].$$

Furthermore, by letting

$$\text{Res}\_{a,c,i}^{b,d} = \int\_{a}^{b} \mathcal{K}(t\_{c\prime}, t\_d + sh) \mathcal{R}\_i^{k\_1,k\_2}(s) ds,\tag{13}$$

we obtain

$$v\_{n} := \begin{cases} h \sum\_{\substack{i=-k\_{1} \\ i=-k\_{1} \\ \text{ ${i}$ }}}^{k\_{1}+1} \text{Mon}\_{k\_{1},n,0}^{\text{nR}} + h^{k\_{1}+k\_{2}+3} \text{REs}\_{k\_{1},n,0}^{n-k\_{1},k\_{1}}, & n=1,\cdots,k\_{1}, \\ h \sum\_{\substack{i=-k\_{1} \\ i=-k\_{1} \\ \text{ ${i}$ }}}^{k\_{1}+1} \text{Mon}\_{k\_{1},n,i}^{\text{M}\_{k\_{1},n,i}^{\text{M}}} + h^{k\_{1}+k\_{2}+3} \text{REs}\_{k\_{1},n,0}^{\text{M}\_{k\_{1},n,0}^{\text{M}}} \\ \quad + h \sum\_{\substack{i=k\_{1} \\ i=k\_{1} \\ \text{ ${i}$ }}}^{n-1} \sum\_{\substack{i=-k\_{1} \\ i=-k\_{1} \\ \text{ ${i}$ }}}^{k\_{1}+1} v\_{i+1} \text{Mon}\_{k\_{1},n,i}^{\text{M}\_{k\_{1},n,i}^{\text{M}}} + h^{k\_{1}+k\_{2}+3} \text{REs}\_{0,n,n-1}^{\text{M}\_{k\_{1},n,0}^{\text{M}}} & n=k\_{1}+1,\cdots,N-k\_{2}, \\ h \sum\_{\substack{i=-k\_{1} \\ i=-k\_{1} \\ \text{ ${i}$ }}}^{N\_{i}+1} \text{Mon}\_{0,n,i}^{\text{M}\_{k\_{1},n,i}^{\text{M}}} + h^{k\_{1}+k\_{2}+3} \text{REs}\_{0,n,0}^{\text{M}\_{k\_{1},n,i}^{\text{M}}} \\ \quad + h \sum\_{\substack{i=k\_{1} \\ i=k\_{1}$$

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

Suppose that MOM*b*,*<sup>d</sup> <sup>a</sup>*,*c*,*<sup>i</sup>* and RES*b*,*<sup>d</sup> <sup>a</sup>*,*c*,*<sup>i</sup>* are bounded by the constant *B*. It is easily noted from Equations (8) and (13) that *B* does not depend on the stepsize. A direct calculation leads to

$$|v\_{n}| \leq \begin{cases} h\mathcal{B} & \sum\_{i=-k\_{1}j\neq n-k\_{1}}^{k\_{2}+1} |v\_{i+k\_{1}}| + h|v\_{n}| + h^{k\_{1}+k\_{2}+3}B\_{\prime} & n=1, \cdots, k\_{1}, \\ h\mathcal{B} & \sum\_{i=-k\_{1}j\neq n-k\_{1}}^{k\_{2}+1} |v\_{i+k\_{1}}| + h^{k\_{1}+k\_{2}+3}B + h(k\_{1}+k\_{2}+2)|v\_{n}| \\ & \quad + h\mathcal{B} \sum\_{i=-k\_{1}j\neq n-l}^{n-1} \left(\sum\_{i=-k\_{1}j\neq n-l}^{k\_{2}+1} |v\_{i+l}| + h^{k\_{1}+k\_{2}+3}B\right), & n=k\_{1}+1, \cdots, N-k\_{2}, \\ & \quad \sum\_{i=-k\_{1}j\neq n-k\_{1}}^{k\_{2}+1} |v\_{i+k\_{1}}| + h^{k\_{1}+k\_{2}+3}B + h(k\_{1}+k\_{2}+2)|v\_{n}| \\ & \quad + h\mathcal{B} \sum\_{l=k\_{1}}^{N-k\_{1}-1} \left(\sum\_{i=-k\_{1}j\neq n-l}^{k\_{2}+1} |v\_{l+i}| + h^{k\_{1}+k\_{2}+3}B\right) \\ & \quad + h\mathcal{B} \sum\_{i=-k\_{1}j\neq n-N+k\_{2}+1}^{k\_{2}+1} |v\_{i+N-k\_{2}}| + h|v\_{n}| + h^{k\_{1}+k\_{2}+3}B, & n=N-k\_{2}+1, \cdots, N. \end{cases}$$

Hence, we have

$$(1 - h(k\_1 + k\_2 + 2))|v\_n| \le h^{k\_1 + k\_2 + 2}B + h(k\_1 + k\_2 + 2)B \sum\_{i=n+1}^{n+k\_2} |v\_i| + h(k\_1 + k\_2 + 2)B \sum\_{i=1}^{n-1} |v\_i|.$$

Since 1 <sup>−</sup> *<sup>h</sup>*(*k*<sup>1</sup> <sup>+</sup> *<sup>k</sup>*<sup>2</sup> <sup>+</sup> <sup>2</sup>) <sup>≈</sup> 1 for sufficiently small stepsize *<sup>h</sup>*, we obtain

$$|v\_{\mathfrak{n}}| \le h^{k\_1+k\_2+2}\bar{B} + hk\_2(k\_1+k\_2+2)\bar{B}||\mathbf{e}\_{\mathfrak{h}}||\_{\infty} + h(k\_1+k\_2+2)\bar{B}\sum\_{i=1}^{n-1}|v\_i|.$$

According to Lemma 2, we have

$$||\mathbf{e}\_h||\_{\infty} \le \mathbf{e}^{(k\_1+k\_2+2)\not B} \mathcal{B}h^{k\_1+k\_2+2} + h \mathbf{e}^{(k\_1+k\_2+2)\not B} (k\_1+k\_2+2)k\_2 \mathcal{B} ||\mathbf{e}\_h||\_{\infty}.$$

or equivalently,

$$||\mathbf{e}\_{\mathrm{h}}||\_{\infty} \le \frac{\mathbf{e}^{(k\_1+k\_2+2)\bar{B}}\mathcal{B}}{1 - h\mathbf{e}^{(k\_1+k\_2+2)\bar{B}}(k\_1+k\_2+2)k\_2\bar{B}}h^{k\_1+k\_2+2}$$

for sufficiently small stepsize *h*.

**Example 1.** *Let us solve VIE with GMCk*1,*k*<sup>2</sup> *M*

$$u(t) = \mathbf{e}^t + \int\_0^t 2\cos(t-s)u(s)ds, t \in [0,2] \tag{14}$$

*with the exact solution u*(*t*)=(<sup>1</sup> + *<sup>t</sup>*)2e*<sup>t</sup>* .

In this example, we test the performance of *GMCk*1,*k*<sup>2</sup> *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.


**Table 1.** Collocation error and convergence order of *GMCk*1,*k*<sup>2</sup> *M* for Example 1.

**Table 2.** Collocation error and convergence order of *GMCk*1,*k*<sup>2</sup> *M* for Example 1.


**Table 3.** Collocation error and convergence order of *GMCk*1,*k*<sup>2</sup> *M* for Example 1.


It can be seen from these tables that as the quantity of nodes increases, absolute errors decay fast, and as *k*<sup>1</sup> and *k*<sup>2</sup> get bigger, the convergence order enlarges. Besides, numerical results illustrate that *GMCk*1,*k*<sup>2</sup> *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, *<sup>T</sup>*], *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.*

#### *2.3. 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) = \lambda y(t), \text{Re}(\lambda) < 0.1$$

Particularly, Brugnano and Trigiante investigated multistep methods for solving differential problems with the above scalar equation in [19]. For the general linear multistep formula

$$\sum\_{j=0}^{k} \alpha\_j y\_{n+j} - h\lambda \sum\_{j=0}^{k} \beta\_j y\_{n+j} = 0\_{\prime \prime}$$

we can introduce two polynomials

$$\rho(z) = \sum\_{j=0}^{k} \alpha\_j z^j, \sigma(z) = \sum\_{j=0}^{k} \beta\_j z^j.$$

and define the associated characteristic polynomial *<sup>π</sup>*(*z*, *<sup>q</sup>*) = *<sup>ρ</sup>*(*z*) <sup>−</sup> *<sup>q</sup>σ*(*z*) with *<sup>q</sup>* = *<sup>h</sup>λ*. When *<sup>π</sup>*(*z*, *<sup>q</sup>*) 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

$$\mathbb{D} := \{ q \in \mathbb{C} : \pi(z, q) \text{ is a Schur polynomial} \}$$

If <sup>C</sup><sup>−</sup> <sup>⊆</sup> <sup>D</sup>, the method is said to be *<sup>A</sup>*−stable.

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:

$$u(t) = 1 + \lambda \int\_0^t u(s)ds, \ t \in [0, T], \ \text{Re}(\lambda) < 0. \tag{15}$$

We turn to study the linear stability of the collocation solution by investigating Equation (15). By applying *GMCk*1,*k*<sup>2</sup> *M* we have

$$y\_j = 1 + \lambda \int\_0^{jh} u\_h(s)ds, j = k\_1 + 1, \dots, N - k\_2. \tag{16}$$

Next, noting the difference between *yj* and *yj*−<sup>1</sup> in Equation (15) leads to

$$y\_j - y\_{j-1} = h\lambda \sum\_{i=-k\_1}^{k\_2+1} y\_{j-1+i} \int\_0^1 \phi\_i^k(s) ds, \; j = k\_1 + 1, \dots, N - k\_2. \tag{17}$$

Then the characteristic polynomial is defined by

$$
\pi^{k\_1,k\_2}(z,q) = z^{k\_1+1} - z^{k\_1} - q \sum\_{i=0}^{k\_1+k\_2+1} z^i \int\_0^1 \phi\_{i-k\_1}^{k\_1,k\_2}(s)ds = \rho(z) - q\sigma(z). \tag{18}
$$

Before investigating the linear stability region, we introduce some helpful definitions and theorems in the version of *GMCk*1,*k*<sup>2</sup> *M*.

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

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

**Theorem 6** ([19])**.** *For any complex number q*, *if roots of Equation* (18) *satisfy*

$$\left|z\_1^k\right| \le \dots \le \left|z\_{k\_1}^k\right| < 1 < \left|z\_{k\_1+1}^k\right| \le \dots \le \left|z\_{k\_1+k\_2+1}^k\right| \tag{19}$$

*then GMCk*1,*k*<sup>2</sup> *M is stable at q*.

By a direct calculation, we find that roots of *<sup>π</sup>k*1,*k*<sup>2</sup> (*z*, *<sup>q</sup>*) do not satisfy the condition given in Theorem <sup>6</sup> in the case of *<sup>k</sup>*<sup>1</sup> <sup>=</sup> *<sup>k</sup>*2. Hence, the region of stability cannot be shown. In Figures <sup>2</sup> and 3, we list the boundary locus corresponding to various multistep collocation methods with *<sup>k</sup>*<sup>1</sup> <sup>=</sup> *<sup>k</sup>*2, where the boundary Γ is defined by

$$\Gamma := \{ z \in \mathbb{C}, z = \frac{\rho(e^{i\theta})}{\sigma(e^{i\theta})}, 0 \le \theta < 2\pi \}.$$

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 *GMCk*1,*k*<sup>2</sup> *M* has wide stability region in the case of *k*<sup>2</sup> > *k*1. In addition, the boundary trajectories of *GMCk*1,*k*2*M* and *GMCk*2,*k*1*M* are symmetric with respect to virtual axis.

**Figure 2.** Linear stability region for *GMC*1,2*M*, *GMC*1,3*M*, *GMC*2,3*M*.

**Figure 3.** Linear stability region for *GMC*2,1*M*, *GMC*3,1*M*, *GMC*3,2*M*.
