Lidstone–Euler interpolation and related high even order boundary value problem

We consider the Lidstone–Euler interpolation problem and the associated Lidstone–Euler boundary value problem, in both theoretical and computational aspects. After a theorem of existence and uniqueness of the solution to the Lidstone–Euler boundary value problem, we present a numerical method for solving it. This method uses the extrapolated Bernstein polynomials and produces an approximating convergent polynomial sequence. Particularly, we consider the fourth-order case, arising in various physical models. Finally, we present some numerical examples and we compare the proposed method with a modified decomposition method for a tenth-order problem. The numerical results confirm the theoretical and computational ones.

One of these is the sequence {ε k } k , with ε k satisfying the conditions Polynomials ε k are connected to the classic Euler polynomials E k by the following relation [14] ε k (x) = 2 2k+1 E 2k+1 x + 1 2 .
Also, they are fundamental polynomials for the interpolation problem In fact, the polynomial Q r , given by [16] Q r (x) = r −1 is the unique polynomial that solves the interpolation problem (6). Therefore we call {ε k } k odd Lidstone-Euler polynomial sequence and (6) Lidstone-Euler interpolation conditions, or boundary conditions. The sequence {ε k } k is not known in the literature, but the polynomials have already been used [14,17,36,38,44]. In fact these polynomials are fundamental polynomials in the modified Abel expansion [36,38,44]. We can consider [16] the Lidstone-Euler boundary value problem (LEbvp in the following): ⎧ ⎨ ⎩ (−1) r y (2r ) (x) = f x, y, y , . . . , y (q) (7a) with 0 ≤ q ≤ 2r − 1 fixed, α i , β i , i = 0, . . . , r − 1, finite real constants. Problem (7a)-(7b), in the general case, does not appear in the literature. In [9] the following problem is considered: and in [29] the more general problem is studied. Sufficient conditions for the existence of at least one solution for problems (8) and (9) are established under appropriate assumptions on f . Some special cases of the LEbvp are important for their applications in numerous branches of applied sciences [31,32,45,46]. This paper is intended to be a contribution to the study of Lidstone-Euler interpolation and of LEbvp.
The paper is organized as follows: in Sect. 2 we consider the Lidstone-Euler interpolation problem and we get some known and new relations that will be used in the sequel. In Sect. 3 we give some theoretical results on Lidstone-Euler boundary value problems. We provide sufficient conditions for the existence and uniqueness of the solution. Section 4 is devoted to presenting a method for the numerical solution. In Sect. 5 we propose an algorithm for effective calculation of the numerical solution of the considered boundary value problem. In Sect. 6 we consider the special case r = 2, and we provide illustrative numerical examples. In order to compare the proposed method with other existing ones, we apply the new method and the Modified Decom-position method presented in [43] to the numerical solution of a LEbvp of order 10. The last Section contains some conclusions and notations on future works.

Lidstone-Euler polynomials and related interpolation
The so-called Lidstone-Euler polynomials have been introduced in [14,16,17] in a broader theoretical and applicative framework. They, however, at least in the even case, are known as Abel polynomials since they are fundamental in the modified Abel series [36,38,44]. In this Section we will recall some known results on odd Lidstone-Euler polynomials sequence {ε k } k , and introduce new properties in order to solve the LEbvp.

Proposition 1 The Lidstone-Euler polynomials ε r can be expressed at any
where Proof The proof follows by induction. From the theory of differential equations the solution ε 1 of the boundary value problem can be written as Next, suppose that (10) is true for r ≥ 1. Then the solution ε r +1 of the boundary value problem can be written as Thus, in view of (12), Hence, from (10) we get Proposition 2 For the function K r the following relations hold Proof The (i) and (ii) follow from Proposition 1 and the boundary conditions ε k (0) = 0 = ε k (1) = 0. The (iii) follows by induction. From (11)-(12) By differentiating, we get By differentiating again we obtain and this proves (iii) for s = 1. Now, suppose that (iii) is true for s. Then Relation (iv) easily follows from (iii).

Proposition 3 For the function K r the following inequalities hold
Proof For r = 1 a direct computation gives Next, for r > 1, relation (12) yields In order to prove the (16), from (14) we get From (15) we obtain The polynomial sequence {ε k } k is useful in the interpolation problem [16]. In fact it is fundamental in the solution of Birkoff interpolation problem given by The following two theorems have been proved in [16]: is the unique solution of the interpolation problem (17).
is the unique polynomial of degree 2r − 1 such that The polynomial P r [ f ] is called the Lidstone-Euler interpolating polynomial of the function f .

Proposition 4 For the derivatives of the Lidstone-Euler interpolating polynomial
there exist constants C 2i and C 2i+1 such that and Proof The proof follows after easy calculations from Theorem 2, relations (4) and from the inequality for the n-th From the theory of differential equation it is sufficient to prove that R (2r )

Corollary 1 For the kernel K r the following explicit representation holds
Proof This can be proved from the comparison between Theorem 3 and Theorem 5.2 in [16], after observing the uniqueness of Peano's kernel.
The result follows from the mean value Theorem for integrals, since f ∈ C 2r [0, 1].

The Lidstone-Euler boundary value problem
As we said, the Lidstone-Euler boundary value problem consists of a general nonlinear differential equation of order 2r as in (7a) with boundary conditions as in (7b).
If f ≡ 0, the problem (7a)-(7b) has the unique solution y = P r , where P r is defined in (18).
So far we are not aware of any proof of the existence and uniqueness of problem (7a)-(7b). For similar problems there are existence theorems in the literature. For example, in [9] the existence of positive solutions is studied for the problem (8) and also for the same differential equation as in (8) but with different boundary conditions. In [29] the existence of solutions of the two-point boundary value problem (9) is analyzed (for other examples see [8,48] and references therein).
The following theorem provides sufficient conditions for the existence and uniqueness of the solution of general LEbvp (7a)-(7b).  (20) and (21) respectively; Then the boundary value problem (7a)-(7b) has a unique solution on D.
Proof From the results in [20] it follows that the boundary value problem (7a)-(7b) is equivalent to the Fredholm type integral equation with y = y, y , . . . , y (q) . If we introduce in C q [0, 1] the finite norm it becomes a Banach space. Let's define B[0, 1] as the set We will show that the operator T : (23) and hence, from hypothesis (i), Propositions 2 and 3, we obtain Similarly, Relations (24) and (25) Moreover, the inequalities (24)- (25) imply that the sets with θ < 1. Hence the uniqueness of the solution follows.

The numerical solution
High order boundary value problems arise in the mathematical modelling of viscoelastic and inelastic flows, deformation of beams, plates deflection theory and in many other applications in engeneering, physics and science.
Only a limited number of this type of problems can be solved analytically. For this reason many researchers have proposed numerical solutions (see, for instance, [1,13,[19][20][21]33,34,43]). Different approaches have been considered, among which collocation, also with spline functions and Lagrange interpolation, Galerkin weighted residual, decomposition, variational techniques.
Here we present a general method for higher order LEbvp, inspired to the idea in [20]. Next we consider in more detail the important case r = 2, that arises in many applications. We consider also a tenth-order LEbvp and compare the proposed method with a Modified Decomposition method [43].
The proposed method is based on extrapolated Bernstein polynomials [10].

The extrapolated Bernstein method
The extrapolated Bernstein method is based on the asymptotic expansion for Bernstein polynomials [10] and the classical Richardson extrapolation process [39].
are the Bernstein basis polynomials.

S i [y] does not dependent on h, and E h [y] vanishes as h → 0.
Proof As usual, the solution y of the LEbvp can be written as where P r [y] is the polynomial interpolating the boundary conditions and K r is the related Peano kernel. For the function y (2r ) Theorem 6 holds. Hence from (28) we get where B n y (2r ) is as in (27), S i y (2r ) (t) , i = 0, . . . , m are coefficients independent on h, and E h y (2r ) → 0 as h → 0. After easy calculations we get the relations (29), (30), (31).

Corollary 2 For the solution y of the LEbvp we get
uniformly in x ∈ [0, 1].
The expansion of Theorem 7 suggests the following extrapolation procedure.

Theorem 8
Let y ∈ C 2(r +m) [0, 1], with m a given integer. Let {n k } k be an increasing sequence of positive integers and h k = n −1 k . We define a sequence of polynomials of degree n i+k as follows For a fixed i Moreover, the following representations of the error and of T Proof See [10,40].
The proposed numerical method consists in approximating the solution y of the LEbvp by the extrapolated polynomial T (0) m−1 [y], being n m the last element of the considered numerical sequence {n i } i . We observe that our method is a continuos method since it provides an approximating convergent polynomial sequence. Of course, for any z ∈ [0, 1], y(z) is approximated by T In [15] the five sequences above have been compared.
In the next Section we outline an algorithm for computing of the approximating extrapolated polynomial.

An algorithm for computing the extrapolated polynomial
In order to calculate the approximate solution we observe that from (30) φ n [y] is an implicitly defined polynomial. Moreover These values can be obtained by solving the following algebraic non-linear system with k, s as in (35), and y i = y i , y i , . . . , y (q) i . The system (36) can be written in the form Let L = q i=0 L i , being L i the Lipschitz constant for f as in Theorem 5. Then for the existence and uniqueness of solution of (37) the following theorem can be proved with standard techniques [20].

) has a unique solution which can be calculated by an iterative method
Moreover, if Y is the exact solution of the system (37), Hence the result follows by applying the contraction mapping theorem.

Numerical examples
First, we consider the case r = 2. In this case the LEbvp becomes with the boundary conditions We observe that (38) is a general non-linear fourth-order differential equation. It often arises from mathematical modeling of many physical phenomena in numerous branches of applied sciences. For example it models the equilibrium states of traveling waves in a suspension bridge [26]. Different boundary conditions represent types of bridge being modeled. Often fourth-order differential equations are called beam equations [32,47], for their relevance in beam theory. In this case the boundary value problem (38)-(39) describes the equilibrium state of the deformation of an elastic beam column with one of its end simply supported and the other end is clamped by sliding clamps. The deformation is caused by a load f ; y represents the bending moment stiffness, y represents the shear force and y (4) is the load density stiffness. Equation (38) can be equipped also with different boundary conditions, for example [31]: The method proposed in this paper provides the explicit expression of the approximating polynomials φ n [y] and of the extrapolated polynomials T with The elements of each block matrix A k , k = 0, . . . , 3, are Observe that the integrals that appear in (40) are easily calculable exactly since the integral functions in (40) are polynomials. Therefore also the elements of the matrix can be easily calculated.
So far we are not aware of any specific methods for the numerical solution of problem (1)-(2) even in the case r = 2 (see [7,22,23,25,[30][31][32]35,37,41,[45][46][47] and references therein). The reason for this could be the non-symmetry of the boundary conditions. Not even specific numerical examples we have found in the literature.
Pending further clarification of these aspects, in the following we report two examples to validate the theoretical results previously given. Since the analytical solutions of the considered examples are known, ∀x ∈ [0, 1] we compute the true absolute errors

Example 1 Consider the following problem
The theoretical solution is y ( The first approximating polynomials, obtained by using equidistant nodes, are − 0.00581x 7 − 0.000482x 8 − 0.0000148x 9 − 6.04527 · 10 −8 x 10 . Figure 1 shows the graphs of the error functions e n i with n i = 2 i , i = 2, . . . , 5 ( Figure 1a) and the graph of E 3 (Figure 1b). Figure 2a contains the plots of the error functions by using the first terms of the Bulirsh sequence for the degree of the approximating polynomials φ n i [y]. Figure 2b shows the graph of E 4 . The errors in x = 1 2 by using extrapolation for different sequences n i are displayed in Table 1.

Example 2 Consider now the following nonlinear problem
The theoretical solution for this problem is y (x) = sin(x). Figure 3a contains the plots of the error functions e n i with n i = 7 + i, i = 1, . . . , 5; Fig. 3b shows the graph of E 4 . Table 2 contains the errors in x = 1 2 by using extrapolation for different sequences n i .

Comparisons
Now we apply the Bernstein extrapolation method proposed previously for the numerical solution of a higher-order LEbvp. To show the reliability of our method, in the absence of existing specific methods for problem (7a)-(7b), we compare it with the modified decomposition method for boundary value problems presented in [43]. The method in [43] is a modified Adomian decomposition method [2,3], providing the solution of boundary value problems in the form of a rapidly convergent series with components that are recursively computed.

Example 3
Consider the following tenth-order nonlinear problem The theoretical solution for this problem is y (x) = e x . In [43] the same differential equation is considered, but with different boundary conditions.
In Table 3 we show the absolute errors in equidistant points in [0, 1], obtained by using the two 12−th degree polynomials T (0) 1 (x) and y W (x). The Table 3 shows that the numerical performance is better for the extrapolated Bernstein method. Both methods provide a succession of converging approximations. The computational cost is not comparable, as the modified decomposition method requires preliminary analytical work and, therefore, does not produce automatic codes, unlike our method.
We can have better results by considering approximating polynomials of higher degree. In Table 4 we show the absolute errors in equidistant points in [0, 1], obtained by using the extrapolation technique for several sequences n i .

Conclusions
In this paper we considered the Lidstone-Euler interpolation problem and we obtained new properties, for example the Cauchy representation of the error, a new representation of the Peano Kernel, bounds for the Kernel and for its derivatives. Moreover we considered the associated Lidstone-Euler boundary value problem, in both theoretical and computational aspects. We introduced a method for the numerical solution of the LEbvp. This method uses the extrapolated Bernstein polynomials. Hence we gave an approximating, convergent polynomial sequence for the numerical solution of the LEbvp. An algorithm for effective calculation is given too. Numerical examples support theoretical results and show that high accuracy in the approximation is achieved by using extrapolation, both in theoretical and computational aspects. A comparison with a modified decomposition method is given for a tenth-order problem.
For future developments we note that the existence and uniqueness theorem can be improved, especially in inequalities. Other numerical approaches are possible and a global comparison on performance can be made. Finally, it is necessary to study the complementary case that id odd degree equations and suitable boundary conditions can be analized.