Relationship between Interpolation and Differential Equations: A Class of Collocation Methods Relationship between Interpolation and Differential Equations: A Class of Collocation Methods

In this chapter, the connection between general linear interpolation and initial, bound- ary and multipoint value problems is explained. First, a result of a theoretical nature is given, which highlights the relationship between the interpolation problem and the Fredholm integral equation for high-order differential problems. After observing that the given problem is equivalent to a Fredholm integral equation, this relation is used in order to determine a general procedure for the numerical solution of high-order differential problems by means of appropriate collocation methods based on the integration of the Fredholm integral equation. The classical analysis of the class of the obtained methods is carried out. Some particular cases are illustrated. Numerical examples are given in order to illustrate the efficiency of the method.


Introduction
The relationship between interpolation and differential equations theories has already been considered. In Ref. ([1], p. 72), Davis observed that the Peano kernel in the interpolation problem is the Green's function of the differential problem φ″ðxÞ ¼ f ðxÞ where φðxÞ ¼ yðxÞ−P 1 ½yðxÞ, being P 1 ½yðxÞ the unique interpolatory polynomial for Eq. (1).
He observed that "these remarks indicate the close relationship between Peano kernels and Green's functions, and hence between interpolation theory and the theory of linear differential equations. Unfortunately, we shall not be able to pursue this relationship" [1].
Later, Agarwal ([2], p. 2), Agarwal and Wong ( [3], pp. 21, 151, 186) considered some separate boundary value problems and the related Fredholm integral equation, using only polynomial interpolation, without taking into account the related Peano kernel. They used Fredholm integral equation in order to obtain existence and uniqueness results for the solution of the considered boundary value problems.
Linear interpolation has an important role also in the numerical solution of differential problems. For example, finite difference methods (see, for instance, [4][5][6] and references therein) approximate the solution yðxÞ of a boundary value problem by a sequence of overlapping polynomials which interpolate yðxÞ in a set of grid points. This is obtained by replacing the differential equation with finite difference equations on a mesh of points that covers the range of integration. The resultant algebraic system of equations is often solved with iterative processes, such as relaxation methods.
Many authors (see [7][8][9][10] and references therein) used linear interpolation with spline functions for the numerical solution of boundary value problems.
Here, we take into account a more general nonlinear initial/boundary/multipoint value problems for high-order differential equations y ðrÞ ðxÞ ¼ f x, yðxÞ , x ∈ I ¼ ½a, b, r ≥1 L i ½yðxÞ ¼ w i , i ¼ 0;…, r−1; x ∈ I 8 < : (2) where yðxÞ ¼ ðyðxÞ, y ′ ðxÞ;…, y ðqÞ ðxÞÞ, 0 ≤ q < r, y ∈ C r ðIÞ, and L i are r linearly independent functionals on C r ðIÞ. Moreover, we suppose that the function f : ½a, b · R qþ1 ! R is continuous at least in the interior of the domain of interest, and it satisfies a uniform Lipschitz condition in y, which means that there exists a nonnegative constant Λ, such that, whenever ðx, y 0 , y 1 , …, y q Þ and ðx, y 0 , y 1 , …, y q Þ are in the domain of f , the following inequality holds (2) is a multipoint value problem (MVP).
In this chapter, -we assume that the conditions for the existence and uniqueness of solution of problem (2) in a certain appropriate domain of ½a, b · R qþ1 are satisfied and that the solution yðxÞ is differentiable with continuity up to what is necessary; -we get the Fredholm integral equation related to problem (2), by polynomial interpolation and the Peano kernel of the linear interpolation problem L i ½yðxÞ ¼ w i , i ¼ 0;…, r−1. In this way, we point out the close relationship between Green's function and Peano kernel; -then, we construct a class of spectral collocation (pseudospectral) methods which are derived by a linear interpolation process.
The reason for which we prefer collocation methods is their superior accuracy for problems whose solutions are sufficiently smooth functions. Recently, Boyd ([11], p. 8) observed that "When many decimal places of accuracy are needed, the contest between pseudospectral algorithms and finite difference and finite element methods is not an even battle but a rout: pseudospectral methods win hands-down."

The Fredholm integral equation for problem (2)
We consider the general differential problem (2), and we prove that it is equivalent to a Fredholm integral equation.
with L i , i ¼ 0;…, r−1, linearly independent functionals on C r ðIÞ, has the unique solution Proof. Since the L i , i ¼ 0;…, r−1 are linearly independent, the result follows from the general linear interpolation theory.
Theorem 1 With the above notations and under the mentioned hypothesis, problem (2) is equivalent to the Fredholm integral equation Proof. The result follows from the uniqueness of the Peano kernel and from Propositions 1 and 2.
From Theorem 1, general results on the existence and uniqueness of solution of problem (2) by standard techniques [2,3] can be obtained. In the following, we will not linger over them, but we will outline the close relationship between interpolation and differential equations. Particularly, we will use linear interpolation in order to determine a class of collocation methods for the numerical solution of problem (2).

A class of Birkhoff-Lagrange collocation methods
Integral Eq. (8) allows to determine a very wide class of numerical methods for Eq. (2), which we call methods of collocation for integration.
Let fx i g m i¼1 be m distinct points in ½a, b and denoted by l i ðtÞ, i ¼ 1;…, m, the fundamental Lagrange polynomials on the nodes x i , that is ðt−x k Þ: Theorem 2 If the solution yðxÞ of Eq. (8) is in C rþm ðIÞ, then where p r, i, m ðxÞ ¼ and the remainder term T r, m ðy, xÞ is given by: being ξ x a suitable point of the smallest interval containing x and all x i , i ¼ 1;…, m.
Theorem 2 suggests to consider the implicitly defined polynomial y r, m ðxÞ ¼ P r−1 ½y r, m ðxÞ þ For polynomial (15), the following theorem holds.
that is, y r, m ðxÞ is a collocation polynomial for Eq. (2) at nodes x j , j ¼ 1;…, m.

A-priori estimation of error
In what follows, we consider the norm jf ðkÞ ðtÞj, ∀f ∈ C ðqÞ ðIÞ: Moreover, we define where R m ðy, tÞ is defined as in (14).

Algorithms and implementation
To calculate the approximate solution of problem (2) by polynomial y r, m ðxÞ at x ∈ I, we need the values y ðsÞ r, m ðx k Þ, k ¼ 1;…, m, s ¼ 0;…, q. In order to get these values, we propose the following algorithm: q and consider the following system System (26) can be written in the form From Eq. (27), we get or, putting GðYÞ ¼ AFðYÞ þ C, For the existence and uniqueness of solution of system (34), we can prove, with standard technique, the following theorem.
Theorem 5 If T ¼ Λ∥A∥ ∞ < 1, system (34) has a unique solution which can be calculated by an iterative method with a fixed ðY m Þ 0 ∈ R mðqþ1Þ and GðY m Þ ¼ AFðY m Þ þ C: Remark 3 If f is linear, then system (27) is a linear system which can be solved by a more suitable method. (27) can be considered as a discrete method for the numerical solution of (2). (15) can generate the polynomial sequence

Remark 5 Method
which is equivalent to the discretization of Picard method for differential equations.

Numerical computation of the entries of matrix
i, j ¼ 1;…m. To this aim, it suffices to compute For the computation of the integral (41), we use the recursive algorithm introduced in Ref. [13]: Remark 6 An alternative approach for the exact computation of integrals (39) and (40) is to use a quadrature formula with a suitable degree of precision.

Outline of the method
Summarizing the proposed method consists of the following steps: 1. determine the interpolation polynomial P r−1 ½yðxÞ satisfying the boundary conditions and compute the Peano remainder; 2. approximate y ðrÞ ðxÞ by Lagrange interpolation on a set of fixed nodal point; 3. compute the elements of matrix A (28) and solve system (26); 4. obtain polynomial (15).

Initial value problems
In the case of initial value problems, in Refs. [13,17,25] If fx i g m i¼1 are the zeros of Chebyshev polynomials of first and second kind, the explicit expression for polynomials p r, i, m ðxÞ can be obtained [13,17,25] for some values of r.
Particularly, for r ¼ 1 and r ¼ 2, in the case of zeros of Chebyshev polynomials of first kind, we get where T k−1 ðxÞ and T kþ1 ðxÞ are the Chebyshev polynomials of the first kind and degree k−1 and k þ 1, respectively, and In the case of zeros of Chebyshev polynomials of second kind and In Refs. [13,25], the authors presented the corresponding implicit Runge-Kutta methods too.
In Ref. [26], Coleman and Booth used also a polynomial interpolant of degree n for y″, but they started from an identity different to Eq. (8) and derived a collocation method for which the nodes {x i } m i¼1 are the zeros of Chebyshev polynomials of second kind.
where v k ðxÞ are the complementary Lidstone polynomials of degree k [27]. The kernel is In Ref. [19], the proposed method applied to problem (2) with conditions (64) and (67), respectively, has been examined in detail.

Multipoint boundary value problems
Let us now consider [15] the following conditions in I ¼ ½−1; 1 In this case P r−1 ½yðxÞ ¼ with p r, k ðxÞ ¼ and l k ðtÞ are the fundamental Lagrange polynomials on the points x j , j ¼ 1;…, r−s. P r−1 ðxÞ is the unique polynomial of degree ≤ r−1 which satisfies the Birkhoff interpolation problem Hence, the solution of problem (2), with multipoint conditions (71), is with P r−1 ½yðxÞ given in Eq. (72) and Observe that Eq. (74) is a special type of Birkhoff interpolation problem with incidence matrix E ¼ ðe ij Þ defined by e 1j ¼ e is ¼ 1;j ¼ 0;⋯, s−1; i ¼ 2;…, r−s þ 1; e ij ¼ 0 otherwise and r ≥ 1.
In Ref. [23], P r−1 ½yðxÞ is presented in a little different form: where E s ðx, l k ðxÞÞ ¼ Let us now consider the following conditions [12,20] yð−1Þ ¼ ω 0 , The solution to the Birkhoff interpolation problem with q r, i ðxÞ ¼ and Hence, the solution of problem (2) is with P r−1 ½yðxÞ given in Eq. (80) and

Numerical examples
In this section, we present some numerical results obtained by applying method (15), which we call CGN method, to find numerical approximations y r, m ðxÞ to the solution of some test problems. In order to solve the nonlinear system (19), we use the so-called modified Newton method [29] (the same Jacobian matrix is used for more than one iteration) and we use algorithm (44) for the computation of the entries of the matrix, when polynomials p r, i, m ðxÞ are not explicitly known. Since the true solutions of the analyzed problems are known, we consider the error function e m ðxÞ ¼ jyðxÞ−y r, m ðxÞj.
The maximum values of e m ðxÞ over the interval ½a, b have also been calculated by using Matlab, particularly the built-in solvers • ode15s, a variable-step, variable-order multistep solver based on the numerical differentiation formulas of orders 1-5; • ode45, a single-step solver, based on an explicit Runge-Kutta (4,5) formula, the Dormand-Prince pair for initial value problems, and the finite difference codes; • bvp4c (with an optional mesh of 200 points) that implements the three-stage Lobatto IIIa formula; • bvp5c that implements the four-stage Lobatto IIIa formula.
for boundary value problems.
Moreover, the powerful tool Chebfun [30] has been used.
Example 1 Consider the following linear ninth-order BVP [28] y From Eq. (7), we get K x 9 ðx, tÞ ¼ Now we calculate the values of the integrals (39) by using Eq. (45), and we solve system (26). Thus, we obtain the approximate solution (15) to problem (85). Table 1 shows the numerical results. The absolute errors are compared with those obtained in Ref. [28], where a modified decomposition method is applied for the solution of problem (85). The second and third columns of Table 1 show the error, respectively, in the method in Ref. [28] and in the CGN method, using in both cases polynomials of degree 12. The last column contains the error in the approximation by a polynomial of degree 14 using CGN method. As collocation points, equidistant nodes in ½0; 1 are chosen. Analogous results are obtained by using Chebyshev nodes of first and second kind, and Legendre-Gauss-Lobatto points.
The maximum absolute errors calculated by using Matlab are displayed in Table 4.
Example 3 Consider now the following nonlinear problem [31] y ð4Þ ðxÞ ¼ sin x þ sin 2 x− y″ðxÞ 2 x ∈ ½0; 1 This kind of problems models several nonlinear phenomena such as traveling waves in suspension bridges [32] or the bending of an elastic beam [33].
Suspension bridges are generally susceptible to visible oscillations, due to the forces acting on the bridge (including the force due to the cables which are considered as a spring with a one-sided restoring, the gravitation force and the external force due to the wind or other external sources). f represents the forcing term, while y represents the vertical displacement when the bridge is bending.
In the case of elastic beam, f represents the force exerted on the beam by the supports. x measures the position along the beam (x ¼ 0 is the left-hand endpoint of the beam), y and y ′ indicate, respectively, the height and the slope of the beam at x. y ″ measures the curvature of the graph of y, and, in physical terms, it measures the bending moment of the beam at x, that is, the torque that the load places on the beam at x.
The considered boundary conditions state that the beam has both endpoints simply supported. Moreover, the derivative of the deflection function is not zero at those points, and it indicates that the beam at the wall is not horizontal. Table 5 shows the comparison between the NMD method presented in Ref. [31] and the CGN method with m ¼ 5 and m ¼ 9, respectively. The approximating polynomial of NMD method has degree 11, while the polynomial considered in CGN method for m ¼ 5 has degree 8.