Collocation Methods for Volterra Integral and Integro-Differential Equations : A Review

We present a collection of recent results on the numerical approximation of Volterra integral equations and integro-differential equations by means of collocation type methods, which are able to provide better balances between accuracy and stability demanding. We consider both exact and discretized one-step and multistep collocation methods, and illustrate main convergence results, making some comparisons in terms of accuracy and efficiency. Some numerical experiments complete the paper.


Introduction
It is the purpose of this paper to illustrate recent results on collocation methods for Volterra integral equations (VIEs) and Volterra integro-differential equations (VIDEs), mainly due to the authors.Such equations model evolutionary problems with memory in many applications, such as dynamics of viscoelastic materials with memory, electrodynamics with memory, heat conduction in materials with memory [1][2][3][4][5][6].The numerical solution of these equations has a high computational cost due both to the nonlinearity of the advancing term and to the evaluation of the lag term, which contains the past history of the solution.Therefore, a crucial point is finding accurate and efficient numerical methods.
Collocation methods have several desirable properties.They provide an approximation over the entire integration interval to the solution of the equation, which reveals to be quite useful in a variable-stepsize implementation: indeed, it is easy to recover the missing past values when the stepsize is changed, by evaluating the collocation polynomial.Other good properties of collocation methods are their high order of convergence, strong stability properties and flexibility.As a matter of fact, if some information is known on the behavior of the exact solution, then it is possible to choose the collocation functions in order to better follow such behavior, so giving rise to mixed collocation methods, see for example [7] in the case of ordinary differential equations (ODEs), and [8] in the case of VIEs.It is also worthwhile mentioning that collocation also has an important theoretical relevance: in fact, many numerical methods are difficult to be analyzed as discrete schemes while, re-casted as collocation-based methods, their analysis is reasonably simplified and can be carried out in a very elegant way.There is, however, a remarkable drawback of one-step collocation methods: they suffer from order reduction phenomenon when applied to stiff problems [9][10][11], since the order of convergence is not uniform (for instance, in the case of s-stage collocation based Runge-Kutta methods on Gauss-Legendre collocation points, the order is p = 2s in the grid points, but it degenerates to p = s for stiff problems, since the order is s in the internal stages).Such a drawback is successfully solved by two-step collocation methods [12], having high uniform order on the overall integration interval.On the side of computational cost, collocation methods are usually more expensive than other classes of methods.In fact, a collocation method with m collocation parameters requires at each time-step the solution of a nonlinear system of dimension m.To face this drawback, multistep collocation methods can be adopted which increase the order of convergence at the same computational cost of one-step ones.When a collocation method is applied to an integral equation, several integrals must be computed, thus suitable quadrature rules are needed to complete the discretization, with the introduction of an additional error.Lastly, a reliable error estimation for collocation methods for integral equations is still missing: there have been some advances (compare [13] and references therein contained), however considerable work needs to be done.
One-step collocation methods first appeared in the literature and main results are collected in the monographs [2,3].Recently, we have proposed multistep collocation methods [13][14][15][16] and two step almost collocation methods [13,17,18], where the collocation polynomial depends on the approximate solution in a fixed number of previous time steps, with the aim of increasing the order of convergence of classical one-step collocation methods, without additional computational cost at each time step, and at the same time obtaining highly stable methods.This idea has been already proposed for the numerical solution of ODEs [19][20][21] (see also [11], Section V.3), and afterward modified in [12], by also using the inherent quadratic technique [22][23][24].We also underline that they have high uniform order, thus they do not suffer from the order reduction phenomenon, well-known in the ODEs context [9].Other approaches, based on multistep collocation, have been proposed in [25][26][27][28][29][30][31][32].
Here we briefly introduce one-step collocation methods and illustrate with more detail the construction and analysis of multistep collocation methods for VIEs and VIDEs, with the aim of giving a complete idea on the recent developments in this context.We give practical indications on how to choose the quadrature formulas in the discretized methods for an efficient implementation.In this review, we consider VIEs and VIDEs with smooth kernel and solution.We illustrate methods with a uniform mesh, however they could easily be applied to a non-uniform mesh (compare [2] for one-step collocation methods).
The paper is organized as follows.Sections 2 and 3 deal with one-step and multistep collocation methods for VIEs, respectively.Section 4 illustrates two-step almost collocation methods for VIEs.Sections 5 and 6 focus on one-step and multistep collocation methods for VIDEs, respectively.

One Step Collocation Methods for VIES
We consider VIEs of the form where k ∈ C(D × R), with D := {(t, τ) : 0 ≤ τ ≤ t ≤ T}, and g ∈ C(I).In the literature, many authors (see [2,3] and references therein contained) have analyzed one step collocation methods for VIEs.As it is well known, a collocation method is based on the idea of approximating the exact solution of a given integral equation with a suitable function belonging to a chosen finite dimensional space, usually a piecewise algebraic polynomial which satisfies the integral equation exactly on a certain subset of the integration interval (called the set of collocation points).
Let us discretize the interval I by introducing a uniform mesh The Equation ( 1) can be rewritten, by relating it to this mesh, as where represent respectively the lag term and the increment function.
Collocation methods provide an approximation P(t) to the solution y(t) of ( 1) on [0, T], such that its restriction to each interval (t n , t n+1 ] is a polynomial:

Exact One-Step Collocation Methods
Let us fix m collocation parameters 0 ≤ c 1 < ... < c m ≤ 1 and denote by t nj = t n + c j h the collocation points.The collocation polynomial, restricted to the interval [t n , t n+1 ], is of the form: where L j (s) is the j-th Lagrange fundamental polynomial with respect to the collocation parameters and Y nj := P n (t nj ).Exact collocation methods are obtained by imposing that the collocation polynomial (2) exactly satisfies the VIE (1) in the collocation points t n,i and by computing y n+1 = P n (t n+1 ): where i = 1, ..., m.Note that the first equation in (3) represents a system of m nonlinear equations in the m unknowns Y ni .We recall that generally P(t) is not continuous in the mesh points, as where Here, Π µ denotes the space of (real) polynomials of degree not exceeding µ.
The classical collocation methods have uniform order O(h m ) for any choice of the collocation parameters, and can achieve local superconvergence in the mesh points by opportunely choosing the collocation parameters, i.e., order 2m − 2 with m Lobatto points or m − 1 Gauss points with c m = 1 and order 2m − 1 with m Radau II points.The optimal superconvergence order O(h 2m ) in the mesh points can be achieved with Gauss nodes in the iterated collocation methods [2,3].

Discretized One-Step Collocation Methods
The collocation Equation (3) is not yet in a form amenable to numerical computation: another discretization step, based on quadrature formulas Fni F ni and Φni Φ ni for the approximation of (4) and ( 5) are needed in order to obtain the fully discretised collocation schemes, thus leading to Discretized collocation methods.
The discretized collocation polynomial is of the form where Ỹnj := Pn (t nj ).The m unknowns Ỹnj are determined by imposing that the collocation polynomial (7) satisfies exactly the integral equation at the collocation points and by using quadrature formulas of the form i = 1, ..., m, for approximating the lag term (4) and the increment function (5).The Formulas (8) and ( 9) are obtained by using quadrature formulas of the form where the quadrature nodes ξ l and d il satisfy 0 ≤ ξ 1 < ... < ξ µ 1 ≤ 1 and 0 ≤ d i1 < ... < d iµ 0 ≤ 1, µ 0 and µ 1 are positive integers and w il , b l are suitable weights.The numerical method is then of the form: where Φn (t ni ) and Fn (t ni ) are given by ( 8) and ( 9).Note that the first equation in (9) represents a system of m nonlinear equations in the m unknowns Ỹni .
Such methods preserve, under suitable hypothesis on the quadrature Formulas (8) and ( 9), the same order of the exact collocation methods [3].
A collocation method for VIEs is equivalent to an implicit Runge-Kutta method for VIEs (VRK method) if and only if c m = 1 (see Theorem 5.2.2 of [3]).As the lag-term computation is the most expensive part in the numerical solution of VIEs, fast collocation and Runge-Kutta methods have been constructed for convolution VIEs of Hammerstein type [33,34] in order to reduce the computational effort in the lag-term computation.The stability analysis of collocation methods for VIEs can be found in [3,35] and the related bibliography.

Multistep Collocation Methods for VIEs
Multistep collocation methods for VIEs have been introduced in [16] by adding interpolation conditions in r previous step points, with the aim of increasing the uniform order of convergence of one step collocation methods without increasing the computational cost.The multistep collocation polynomial, restricted to the interval [t n , t n+1 ], is of the form n = r, ..., N − 1, where again and ϕ k (s), ψ j (s) are the following polynomials of degree m + r − 1 The collocation parameters are assumed to satisfy c i = c j and c 1 = 0.

Exact Multistep Collocation
The exact multistep collocation methods are obtained by imposing that the collocation polynomial (12) exactly satisfies the VIE (1) at the collocation points t ni , and by computing where the lag-term F ni and increment-term Φ ni are given by ( 4) and ( 5) respectively.The r-step m-point exact collocation method ( 12)-( 15) has uniform convergence order of at least p = m + r, for any choice of distinct collocation abscissas 0 < c 1 < ... < c m ≤ 1, as stated in the following theorem proved in [16].
Theorem 1.Let ε(t) = y(t) − P(t) be the error of the exact collocation method ( 12)-( 15) and p = m + r.Suppose that i. the given functions describing the VIE (1 and ρ denotes the spectral radius. Then Moreover, a suitable choice of collocation parameters can ensure superconvergence in the mesh points, as pointed out in the following theorem [16].

Theorem 2. Let us suppose that
• the hypothesis of the Theorem 1 hold with p = 2m + r − 1.

•
the collocation parameters c 1 , ..., c m are the solution of the system with

Discretized Multistep Collocation
The discretized multistep collocation methods are obtained by using quadrature formulas of the form (8) to (9) for approximating the lag term and the increment function.The discretized multistep collocation polynomial, denoted by Pn (t), is then of the form n = 0, ..., N − 1,where the functions ϕ k (s) and ψ j (s) are given by ( 14), and Ỹnj := Pn (t nj ) are determined by the solution of the following nonlinear system The following theorem [16] shows that, as in the exact case, the r-step m-point discretized collocation method ( 19) and ( 20) has convergence order of at least p = m + r, for any choice of distinct collocation abscissas 0 < c 1 < ... < c m ≤ 1.
Theorem 3. Let ε(t) := y(t) − P(t) be the error of the discretized collocation method (19) and (20) and let p = m + r.Suppose that i.
An analogous result holds concerning the local superconvergence:

Two Step Almost Collocation Collocation Methods for VIEs
Within the class of multistep collocation methods, although methods with unbounded stability regions exist, no A-stable methods have been found [16].In order to determine A-stable methods, two step almost collocation (TSAC) methods have been introduced in [18] and further analyzed in [13,17].
As regards the global error, the method has uniform order of convergence p * = min{l + 1, q, p + 1}, where l and q are the order of the starting procedure (for the computation of the starting values y 1 and Y [1] i , i = 1, 2, ..., m) and the order of the quadrature Formulas ( 23) and ( 24) respectively (see Theorem 2.5 in [18]).Then we use as starting procedure a one step collocation method having uniform order of convergence l = p.Two-step collocation methods are obtained by solving the system of order conditions up to the maximum uniform attainable order p = 2m + 1, and, in this way, all the basis functions are determined as the unique solution of such system.However, as observed in [18], it is not convenient to impose all the order conditions because it is not possible to achieve high stability properties (e.g., A-stability) without getting rid of some of them.Therefore, almost collocation methods have been introduced by relaxing a specified number r of order conditions, i.e., by a priori opportunely fixing r basis functions, and determining the remaining ones as the unique solution of the system of order conditions up to p = 2m + 1 − r.Within the class of TSAC methods, A-stable methods have been constructed in [18] by fixing one (case r = 1) or both (case r = 2) of the polynomials ϕ 0 (s) and ϕ 1 (s) as where α j and β j , j = 0, 1, . . ., p − m, are free parameters, which have to be determined in order to obtain desired stability properties.
A error estimation of the local discretization error for TSAC methods has been derived in [13].

Diagonally Implicit TSAC Methods for VIEs
The computational cost associated to the solution of the nonlinear system ( 22) can be reduced by making the coefficient matrix have a structured shape, e.g., lower triangular or diagonal.This strategy, in the field of Runge-Kutta methods for ODEs, leads to the raise of the famous classes of Diagonally Implicit and Singly Diagonally Implicit Runge-Kutta methods (DIRK and SDIRK), see [10,11] and bibliography therein contained.Moreover, in the field of collocation-based methods for ODEs, an analogous strategy has been applied, obtaining TSAC methods having structured coefficient matrix [12].
In fact, a lower triangular matrix allows to solve the equations in m successive stages, with only a d-dimensional system to be solved at each stage.Moreover, if all the elements on the diagonal are equal, in solving the nonlinear systems by means of Newton-type iterations, one may hope to use repeatedly the stored LU factorization of the Jacobian.If the structure is diagonal, the problem reduces to the solution of m independent systems of dimension d, and can therefore be solved in a parallel environment.
Methods of this type have been derived in [17], where first of all it was assumed w j,m+1 = 0, j = 1, . . ., m, in such a way that ( 22) becomes a nonlinear system of dimension md only depending on the stage values Y [n+1] i , i = 1, . . ., m, and assumes the following form where By defining , the nonlinear system in ( 28) takes the form where • denotes the usual Hadamard product, I is the identity matrix of dimension d and e is the unit vector of dimension md.The tensor form (30) clearly shows as the matrices which determine the structure of the nonlinear system (28) are Ψ and W. In [17] a strategy was described to obtain lower triangular or diagonal structures for the matrices Ψ and W: in particular a quadrature formula of the form was proposed for the increment in addition to the quadrature formula for the approximation of the lag term We observe that in Formula (31), in case of triangular structure, wjl = 0, l = 1, . . ., j while, in case of diagonal structure, wj1 = 0 and w jl = 0, l = 1, . . ., j − 1.The determination of the weights in Formulas ( 31) and (33) has been described in [17].
Assuming that Ψ and W are lower triangular, we obtain the diagonally implicit TSAC methods (DITSAC) ), where B [n] i is given by ( 29), and F [n] j , Φ [n+1] j are approximations of ( 34) by means of the quadrature Formulas ( 31) and ( 33).

Numerical Results
We present some numerical results which confirm that, differently from one step collocation methods, the TSAC methods do not suffer form the order reduction in the integration of stiff systems, as we expect from the uniform order of convergence stated in Theorem 5.In order to illustrate this phenomenon, we show the results obtained on both a non stiff and a stiff equation: with exact solution y(t) ≡ 1; with λ = −10 4 and exact solution y(t) = sin(t).This is a stiff problem because it is equivalent to the Prothero-Robinson problem for ODEs.
We compare TSAC methods with superconvergent one step collocation methods of [2,3], where m denotes the number of collocation points and p denotes the order of the method: The method TSAC2 is the two-stage TSAC method described in Example 1.The accuracy is defined by the number of correct significant digits cd at the end point (the maximal absolute end point error is written as 10 −cd ).For each test we plot in Figure 1 the number of cd versus the number of mesh points N. We observe as for non stiff Problem (37) the effective order of the all methods is coherent with the theoretical order, while for stiff Problem (38) the one step methods show order reduction as the effective order reduces to p = 2.

One-Step Collocation Methods for VIDEs
We concentrate on VIDEs of type: where g(t, y) : For sake of completeness we report the theorem of existence and uniqueness of solution for (39) [3].Theorem 6.Let g(t, y) and k(t, s, y) be continuous functions and satisfy a uniform Lipschitz condition with respect to y.Then there exists a unique solution y ∈ C 1 ([0, T]) of the problem (39).
T} be a partition of the time interval [0, T] with constant stepsize h = t n+1 − t n , n = 0, . . ., N − 1.The integro-differential Equation (39) can be written as follows: where represent respectively the lag term and the increment function.

Exact One-Step Collocation Methods
Here we briefly expose the classical one-step collocation methods for VIDEs and their main properties [2,3].
The numerical approximation at the point t n+1 is then given by y n+1 = P n (t n+1 ).The final form of an exact collocation method is n = 0, . . ., N, where the lag term and the increment function can be written as The first equation in (41) requires, at each time step, the solution of an m-dimensional nonlinear system in the unknowns {U ni } m i=1 .For every choice of the collocation parameters c 1 , . . ., c m , the collocation polynomial P(t) is continuous on [0, T] and provides a uniform approximation of order O(h m ).Moreover, if c 1 , . . ., c m are suitably chosen, the order of convergence at the mesh points increases (local superconvergence): is 2m − 2 for the Lobatto points, 2m − 1 for the Radau points and 2m for the Gauss ones [2,3].

Discretized One-Step Collocation Methods
In the general case, the integrals appearing in (42) and ( 43) cannot be exactly evaluated, so a further approximation is needed in order to fully discretize the method.Let us suppose to approximate these integrals by quadrature formulae of the type: These formulae are then used to define the discretized collocation methods as where the collocation polynomial is now of the form The discretized collocation methods are a special class of the Runge-Kutta extended methods and preserve the order of convergence and superconvergence of exact collocation methods, if the quadrature Formulae (44) and ( 45) are sufficiently accurate [3].Some relevant stability results for one-step collocation methods are derived in [36,37].

Exact Multistep Collocation
Recently, in order to obtain an higher order of convergence at the same computational effort, multistep collocation methods have been introduced: the solution y(t) is approximated by a piecewise algebraic polynomial P(t): where again and the functions ϕ k (s), ψ j (s) are polynomials of degree m + r − 1 which are determined by imposing that the polynomial ( 48) satisfies ( 49) and the interpolation conditions: For any fixed set of collocation parameters c 1 , . . ., c m , conditions ( 49) and ( 50) lead to the following non linear system of (r + m) 2 equations, where the (r + m) 2 unknowns are the coefficients of the polynomials ϕ k (s) and ψ j (s): l, k = 0, ..., r − 1, i, j = 1, ..., m.Exact multistep collocation methods are obtained by imposing that the collocation polynomial (48) satisfies the VIDE at the collocation points t ni , and by computing y n+1 = P n (t n+1 ): n = r − 1, . . ., N, where now the lag term and the increment function can be written as We note that at each time step, the approximations y n−k , k = 0, ..., r − 1 are already known, so only the unknowns {U ni } m i=1 need to be computed, by solving the nonlinear system given by the first equation of (52).
Observe that we are able to give an approximate value P(t) of the solution y(t) at each point t of the integration interval, therefore we have a uniform approximation of the solution on [0, T].
The classical one-step collocation methods described in the previous section can be seen as a particular case of multistep methods with r = 1 and where L j (τ) is the j-th Lagrange fundamental polynomial with respect to the collocation parameters.We observe that, at each time step, both one-step and multistep collocation methods require the solution of a non linear system of dimension m for the stages U ni , i = 1, . . ., m.The multistep methods only need in addition the computation of the starting values y 1 , . . ., y r−1 .

Discretized Multistep Collocation
As in the case of one-step collocation methods, it is evident that the exact multistep collocation methods (52) are not directly applicable for the implementation, since approximations of the integrals F n (t ni , P(•)) and Φ n (t ni , P(•)) are needed.With the aim of fully discretizing the multistep collocation methods we consider the following quadrature formulas where the weights and nodes are suitably chosen, as it will be illustrated later.
The discretized multistep collocation method for the problem (39) approximates the solution y(t) with a piecewise polynomial P(t), with where the polynomials {ϕ k (s)} r−1 k=0 , {ψ j (s)} m j=1 are the same as in the exact collocation, and can be computed by solving the system (51).
We impose that at the collocation points P(t) satisfies the VIDE (39), where the integrals appearing in both the lag term (53) and the increment function (54) are approximated by the quadrature formulae defined in (55), and we set ỹn+1 = P(t n + h).Thus the discretized multistep method is where Fn (t ni , P(•)) and Φn (t ni , P(•)) are of the form i.e., they are obtained by applying the quadrature Formula (55) to the integrals appearing in ( 53) and (54).

Convergence Analysis
The multivalue nature of the multistep methods imposes to analyze first the zero-stability of the methods.When h → 0, second equation of (52) reduces to Therefore, the method (48) and ( 52) is said to be zero-stable, if all of the roots of the polynomial with g(t, y) such that y(t) = exp(−t); and on the nonlinear problem with g(t, y) such that y(t) = cos(t).We consider three methods • TS3: superconvergent discretized two-step collocation method, with r = 2 and m = 1, with order p = 3; • TS3b: two-step discretized collocation method, with r = 2 and m = 2, c 1 = 0.9, c 2 = 1, with uniform order 3; • TS5: superconvergent discretized two-step collocation method, with r = 2 and m = 2, with order p = 5.
Method TS3b has an unbounded stability region, while TS3 and TS5 have a bounded stability region.The exact expression of the methods and their stability region can be found in [14].To confirm the theoretical order of convergence, in Figure 2 the error (in logarithmic scale) produced by methods TS3, TS3b and TS5 when applied to problems ( 61) and (62), and the slopes corresponding to order 3 and 5.We see that the effective order is equal to the theoretical one.

Conclusions
We have illustrated multistep collocation methods for VIEs and VIDEs and gave an overview of their convergence and superconvergence properties.This idea may be exploited to obtain high order methods for solving other types of equations as well.For example, recently two-step collocation methods have been proposed for fractional differential equations [44], and further developments may be achieved for other fractional models, as time fractional differential equations [45].Further issues of this research will focus on oscillatory problems [46,47] and in particular on the application of multistep collocation methods to periodic integral equations [48,49].Moreover, it seems reasonable to consider the possibility of employing collocation spaces based on functions other than polynomials, as in [50][51][52] and similarly as in the case of oscillatory problems [53], and merge into the numerical scheme as many known qualitative properties of the continuous problem as possible, in a structure-preserving perspective [54].
The literature on the numerical treatment of VIEs is quite rich and goes beyond the results considered in this review.Here we would like to mention some other results, at least.In [55] the modified Newton-Kantorovich method combined with collocation were applied non linear and nonlinear VIE with piecewise smooth kernels.Such VIE were introduced in [56] and asymptotic approximations to parametric families of solutions were constructed and the existence of continuous solutions was proved.The review of the numerical methods of optimal accuracy (spline-collocation technique) for multidimensional weakly singular VIEs is given in [57].Some other interesting papers regard the distance between the approximate and exact solutions of various generalizations of the Volterra equations [58][59][60][61][62][63].Lastly, we underline that in the practical applications of VIE based models it is extremely important to have the numerical method to be stable with respect to the measurement errors both in the source function and in the kernel.It is well known that the 1st kind of VIEs enjoy self-regularization property when the mesh step serves as the regularisation parameter.In addition, the Lavrentiev type regularisation is a good option [64,65].These issues have been discussed in [66].

Figure 1 .
Figure 1.Number of correct significant digits with respect to the number of mesh points.(a) Problem (37); (b) Problem (38).