Convergence analysis and parity conservation of a new form of a quadratic explicit spline with applications to integral equations

Introduction In Mathematics, Physics, and Engineering, among other disciplines, there is a great need to adjust discrete sets of data and to approximate functions. In general, it is desired to know the values at intermediate points, which can be solved through interpolation polynomials. In practice, high-order polynomials can introduce significant errors due to various factors. Generally, these polynomials present fluctuations that are not present in the function to interpolate. For this reason, in this work, the piecewise interpolation will be considered, which is particularly useful when the data to adjust have a smooth behavior alternated with strong changes. The focus will be on quadratic piecewise interpolation S of continuous real functions. In the literature, methods and algorithms that depend on a determined criterion in the calculation of one of the coefficients of S are presented. In [1], a good development of interpolation methods is exposed, particularly spline methods. In [2], some optimization algorithms can be consulted. There is a very large and current body of literature on quadratic splines and collocation methods. In [3], quadratic spline interpolation is used where the coefficients of the polynomials are determined in matrix form, similarly in [4]. Typically, all of these methods require the resolution of algebraic systems or recursive equations. In this article, a variational alternative that minimizes the fluctuations of the interpolation polynomial S is presented. Its coefficients are determined explicitly through simple arithmetic computations.


Introduction
In Mathematics, Physics, and Engineering, among other disciplines, there is a great need to adjust discrete sets of data and to approximate functions. In general, it is desired to know the values at intermediate points, which can be solved through interpolation polynomials. In practice, high-order polynomials can introduce significant errors due to various factors. Generally, these polynomials present fluctuations that are not present in the function to interpolate. For this reason, in this work, the piecewise interpolation will be considered, which is particularly useful when the data to adjust have a smooth behavior alternated with strong changes. The focus will be on quadratic piecewise interpolation S of continuous real functions. In the literature, methods and algorithms that depend on a determined criterion in the calculation of one of the coefficients of S are presented. In [1], a good development of interpolation methods is exposed, particularly spline methods. In [2], some optimization algorithms can be consulted.
There is a very large and current body of literature on quadratic splines and collocation methods. In [3], quadratic spline interpolation is used where the coefficients of the polynomials are determined in matrix form, similarly in [4]. Typically, all of these methods require the resolution of algebraic systems or recursive equations. In this article, a variational alternative that minimizes the fluctuations of the interpolation polynomial S is presented. Its coefficients are determined explicitly through simple arithmetic computations. There is great interest of the Fredholm-Volterra integral equation in engineering, physics, and other disciplines, which drives a further need to find solutions. As is well known, explicit solutions are very difficult to obtain, so approximation methods should be in order here. There are many works which have analyzed and proposed numerical methods, as well as other based on the methods of splines. A recent interesting work is [5] in which they apply the collocation method using Chebyshev's polynomials to numerically evaluate the problem of Fredholm-Volterra; in [6] this kind of integral, equations are solved with B-Splines; in [7], fifth-order splines are used; in [8], the quadratic spline with end condition is utilized; and in [9], variable-order fractional functional differential equations are solved with Legendre collocation. In the present article, the spline obtained is used to solve Fredholm-Volterra integral equations and fractional differential equations.
An advantage of the spline presented in this work with respect to the cited spline methods is that an explicit law is determined for the computation of the coefficients, which avoids the resolution of algebraic systems that, although linear, are not exempt from method error, plus the possibility of not being well conditioned which requires a rescaling of the system. As for the application to integral equations, in particular the differential equations with fractional derivative, the need to solve an algebraic system lies in the lack of knowledge of the solution, but even so, the calculation of the coefficients of the algebraic system is simple and explicit.
The work is organized as follows: in "The quadratic spline" section, the interpolator is presented; in the "Convergence analysis" section, the convergence is studied; in the "Parity conservation" section, it is shown that S maintains the parity of the function to be interpolated; in the "Fredholm integral equation" section, the results are applied for Fredholm linear integral equations, similarly for Volterra linear integral equations in the "Volterra integral equation" section; in the "Examples" section, the numerical results are shown; and finally, in the "Conclusions" section, the conclusions are presented.
Let S(x) be the function in [ x 0 , x n ] defined through the polynomials P k (x) so that: Since S is an interpolator, P k (x) must verify: The continuity of S (x) imposes that: Lagrange polynomials p k (x) over I k are built as: which satisfy p k (x k−1 ) = y k−1 , p k (x k ) = y k , k = 1, . . . , n. Then: In [3,4], the coefficients a k are obtained by linear algebraic systems. Here, an explicit formula will be found. From (2) and through simple algebraic operations: from where a k = a k (a 1 ), linear in a 1 . The coefficient a 1 is determined from an additional condition.
From (5): Thus: In this way, an explicit expression is available in terms of elementary functions, in order to calculate the coefficients of the polynomial S. In addition, if we define r 1 = 0, it is easy to show that: Next, in order to minimize the fluctuations, a 1 is determined such that the sum of the quadratic errors between P k (x) and p k (x) in each I k is minimum. Let E be defined as E = n k=1 E k , where: Taking into account (3) and (4): from where: Then, recalling (6), results ∂ ∂a 1 a k = (−1) k+1 from where ∂ 2 E(a 1 ) ∂a 2 1 = 1 15 nh 5 > 0. Therefore, the solution a 1 that minimizes E is such that ∂E(a 1 ) ∂a 1 = 0. This leads to: From (7): Finally, combining (8) and (9): where Taking into account the definition of j : where c k,j is given by: ..,n . Then, the spline S can be written in the following matrix form: The importance of this equation lies in the fact that the matrix C is only once calculated because it only depends on x 0 , h, and n. It does not depend on the function to interpolate. For each data collection Y , the only new calculation is the scalar product C · Y .

Convergence analysis
Here, we will prove that the convergence order obtained for the spline presented here is O(h) when the interpolated function has a bounded second derivative.

be the spline presented in this work that interpolates y(x). Then, S(x) converges uniformly to y(x) when |h| −→ 0 with a convergence order of O(h).
Demonstration: Let: Recalling (1), results: From (10): Noting that | j n + s j − 1| ≤ 1 ∀j = 1, . . . , n − 1 and as the second derivative of y(x) is bounded, namely | j | ≤ M for h sufficiently small, results: In this way, back to (13): Using the result of [1]: It results:

Parity conservation
Next, it will be shown that the developed spline preserves the parity of the interpolated function.
Theorem 2 Coefficients a k of the spline S that interpolates y = y(x), −a < x < a, in n + 1 nodes, verify: 1) If y is even, then a k = a n−k+1 , k = 1, . . . , n.

Demonstration:
We will do the demonstration for y even. Since y is even, then Case 1: n odd. It will be proved that a k = a n−k+1 , k = 1, . . . , n−1 2 by induction on k.
Let us see that a k+1 = a n−k . From (5): a k+1 = k − a k = n−k − a n−k+1 = a n−k .
Therefore, a k = a n−k+1 , k = 1, . . . , n 2 if n is even. The proof for y odd and n even is similar to case 1 and the proof for y odd and n odd is similar to case 2.
Observation: In two of the cases (y even, n odd and y odd, n even), the property is verified independently of the election of a 1 , while in the other two cases (y even, n even and y odd, n odd), the election of a 1 is fundamental. It is easy to verify that for these cases, the property is not verified for every a 1 .

Demonstration:
Again, we will demonstrate for y even. Recalling (4) is easy to verify that p k (x) = p n−k+1 (−x), k = 1, . . . , n. Then, taking into account that a k = a n−k+1 , k = 1, . . . , n and that x k = −x n−k , k = 0, . . . , n, the result is effortlessly obtained. The demonstration for the case y odd is similar.

Fredholm integral equation
There exist numerous works to determine the numerical solution of Fredholm linear integral equations of the first and second kind. In [6], it is evaluated numerically using a mixed form between splines and Lagrange interpolation. This kind of problem is then solved using Taylor's expansion in [10]. In [11], the least squares method is utilized. The quadrature methods are multiple, see for example [12]. An interesting work is [5] where Chebyshev's polynomials are used. Finally, in [7], splines of order five are used. Here, the quadratic spline determined by our method is applied to numerically solve the linear problem. Consider the following equation: where , and λ is a real parameter. When f (x) ≡ 0, the equation is homogeneous and λ becomes an eigenvalue (or characteristic root) associated to the eigenfunction y(x) in (14).

Suppose also that x, s ∈ R, f is continuous in I =[ a, b], and the kernel K(x, s) is continuous in the region
Consider a partition of the interval I with nodes x j = a + jh, h = b−a n , j = 0, 1, . . . , n. Using the quadratic spline S developed in "The quadratic spline" section, y(x) is interpolated in {y j } n j=0 , whose values must be determined. Evaluating (14) in nodes x j : Since a k is linear in y j , j = 0, . . . , n, this leads to a linear system. Taking into account (4) and (11): where δ i,k is the Kronecker delta, c k,i is that of (11), m j, From here: In the non-homogeneous case, (15) represents a linear algebraic system whose solution determines {y j } n j=0 and therefore S(x). Remark that in this case, if λ does not coincide with the eigenvalues of the associated homogeneous equation, the integral equation has a solution, but if it matches an eigenvalue, it does not always have a solution. Therefore, the non-homogeneous system of equations does not always have a solution.
When f (x) ≡ 0 in [ a, b], the system is homogeneous and the problem of eigenvalues λ and eigenfunctions y(x) is solved in the usual way.

Volterra integral equation of the second kind
Writing the Volterra integral equation of the second kind as and λ is a real parameter.
For the numerical resolution of (16), partite I using nodes x j = a + jh, h = (b − a)/n, j = 0, 1, . . . , n. The quadratic spline S interpolates y(x) in {y j } n j=0 , whose values must be determined. Then, evaluating (16) in the n + 1 nodes x j follows: Taking into account that y(x) is approximated by S(x): Similarly to the "Fredholm integral equation" section, these equations represent a nonhomogeneous system of order n. From its resolution, {y j } n j=0 is determined.

Volterra integral equation of the first kind
Writing the Volterra integral equation of the first kind as: In the particular case in which K(x, s), and K(x, x) does not vanish in 0 ≤ x ≤ b, it is possible to derivate (17) with respect to x, obtaining:  which is an equation of the second kind. A method is proposed to solve (17) which does not require the continuity conditions of f (x) and K(x, s) nor the restriction on the zeros of K(x, x). In the same way as in the previous cases, and F(x k ) = 0, k = 0, 1, . . . , n. If y(x) is the solution, F(x) ≡ 0. Taking into account that y is approximated by S, (18) writes as: To determine {y(x k )} n k=0 , y(x) is approximated by the spline S(x) given by (4). Evaluating (18) in x k , k = 1, . . . , n: and given the fact that P j (s) depends linearly on {y(x k )} n k=0 , a system of n nonhomogeneous linear equations is obtained which determines {y(x k )} n k=1 parametrized in y 0 = y(x 0 ), which is determined by the following ansatz: y(x 0 ) is the value that minimizes It is easy to see that it is sufficient that ∂ ∂y 0 G = 0. Derivability of G is guaranteed by the linearity of P j in y 0 .

Examples
Quadratic spline Consider f (x) = |x|, −1 ≤ x ≤ 1, interpolated with the Lagrange polynomial L n (x) with equidistant nodes. It turns out as is well known, . Let e n be: It can be seen that for n > 15 the Lagrange interpolation I presents sharp fluctuations. In Table 1, some errors are calculated.
However, with the S spline from this work, the undesired fluctuations are markedly reduced, as seen in Table 2. Remark that despite the fact f is not a function of class C 1 , a good interpolator is achieved.
Regarding the space convergence order, using D defined on (12), it is possible to see that In this example, for n = 20, D = 5 · 10 −2 , then 2 D 2 = 5 · 10 −3 which is effectively larger than e 20 = 6.6 · 10 −4 . The space convergence order is linear as shown in Fig. 1. The dots were numerically obtained, and the curve was obtained by a least squares method.
As a second example, consider g(x) = sin 2πx, −1 ≤ x ≤ 1. Interpolations are made with the Quadratic Spline end Condition routine (Q) from [8], with the natural cubic spline of Mathematica 9.0.1.0 software (M) and with the S spline from this work. Results are shown in Table 3.  Therefore, in this example, the spline presented here is a better approximation than the Quadratic Spline end Condition routine.
Regarding the space convergence order, for n = 50, D = 1.3 · 10 −4 then 2 D 2 = 3.38 · 10 −8 which is effectively bigger than e 50 (S) = 9 · 10 −9 . In this case, Fig. 2 shows that the space convergence order is cubic. The dots were numerically obtained, and the curve was obtained by the least squares method using a third degree polynomial.
Observe that in the first example a linear space convergence order was obtained while in the second it was cubic. The difference in behavior is due to the fact that f (x) = |x| has no continuous derivative at x = 0 while the derivatives of g(x) = sin 2πx are continuous at every point.

Fredholm equation
Consider two extra examples from [13]. The first of them being whose solution is y(x) = e x (2x − 2 3 ). The solution obtained using the quadratic spline gives the results shown in Table 4, where the error e n is defined in (19) and The second example corresponds to an equation of the first kind: whose solution is y(x) = x(1 − 2x), being the eigenvalue λ = −3, of multiplicity 2. The results are shown in Table 5. Finally, consider from [11]: whose solution is y(x) = cos x. With the S spline, for n = 5 and n = 10, e 5 (S) = 1.17·10 −3 and e 10 (S) = 1.79 · 10 −5 are obtained respectively. In [11], y is approximated by {1, x, x 2 } using the least squares method (L) where it is obtained that e 2 (L) = 1.49 · 10 −3 , which is very similar to our results.
In [6], the error is not specified and the solutions are only compared graphically. For the developed spline, the curves corresponding to the numerical and analytical solutions are completely overlapped.
Again from [13]: whose solution is y(x) ≡ 3. The results are shown in Table 8.

Fractional differential equation
Let α be a real positive number and denote by m = α the smaller integer bigger than α. Let us define the Caputo fractional derivative, D α a , of a n times differentiable function of real variable, y(x), as [14] D α a y( where y (m) (t) means the mth derivative of the function y(t).
Let y(x) :[ 0, X] −→ R be a differentiable real function of the real variable x and f (x) : [ 0, X] −→ R continuous. Let 1 < α ≤ 2. Then, let us consider the following fractional differential equation:  Converting the initial value problem for the differential equation into an equivalent Volterra integral equation [14]: Consider the fractional oscillator [15]: where 1 < α ≤ 2, whose exact solution is: where E α,β (z) is the so-called Mittag-Leffler function .
Taking into account the results of the "Volterra integral equation" section, the equation of Volterra associated to the equation (27) is solved with the S spline and compared with the exact solution. The parameters used are α = 3/2, X = 10, y(0) = 1, and y (0) = 0. Results are shown in Table 9.

Conclusions
Using variational calculus, a quadratic spline method that minimizes the spline fluctuations has been developed. Even though piecewise interpolation has several decades of study, since the 1940s when it was developed by the mathematician I. J. Schoenberg, the main advantage of this scheme is that the coefficients of the spline are explicitly determined through simple arithmetic calculations without needing recursive equations nor solving algebraic systems. The reason of using a quadratic spline is that in this case the explicit law of the coefficients of the interpolating segmental polynomial is simply obtained. For higher order splines, although they improve the interpolation, we were not able to find a simple and explicit expression of the coefficients. This variational method may be considered as an adaptive scheme, due to its simplicity in the determination of the coefficients. The spline has a linear space convergence order for functions with bounded second derivative and it preserves the parity of the function. It is also useful for solving Fredholm-Volterra linear integral equations and fractional differential equations given its simplicity. Some weaknesses of the method are that the nodes must be equally spaced and we could not extend this form to the general high-order splines and multiple dimensions. In a second stage, the results will be extended to fractional ordinary differential equations. Table 9 Errors e n and E T with the S spline for (27)