Revisiting the Ceschino Interpolation Method

well-known


Introduction
The Ceschino polynomial expansion method is a generalization of the Taylor polynomial expansion method where higher derivatives of a function are predicted in addition to the value of the function itself.This technique was first introduced by (Ceschino, 1956), but was largly forgotten afterward.An unsuccessfull atempt was tried in 1975 to apply the Ceschino coupling relations to the solution of an elliptic space-dependent differential equation, but the resulting spatial discretization was found to be less accurate than competing finite-element approaches, as presented by (Pageau, 1975).No further published work was reported after the Pageau thesis.
Here, we propose to apply the Ceschino coupling relations to the basic interpolation problem, as an alternative to existing univariate interpolation schemes, such as the cubic spline approach.The interpolation problem consists to evaluate a functional I{ f (x); ξ} of a continuous function (or dependent variable) f (x) at a specific point ξ in the case where function f (x) is only known at tabulated abscissa (or independent variables) {x m+1/2 ; m = 0, M}.We also introduce the concept of interpolation factors (a.k. a. , terp factors) that are useful for interpolating large databases with respect to a small number of independent variables, as presented by (MacFarlane, 1984).The Ceschino polynomial expansion method is the core component of the multiparameter reactor database system used in the reactor physics code DRAGON for performing cross section interpolation (Hébert, 2009).We will show that Ceschino polynomial expansion theory is an attractive choice for computing such interpolation factors and propose sample Matlab scripts for performing this task.

Ceschino polynomial expansion theory
The polynomial expansion theory is first applied over the one-dimensional domain depicted in Fig. 1.A continuous function f (x) is defined over this domain and is known at specific abscissa points x m+1/2 .A (J + 1)-th order Taylor series expansion of f (x) around x = x m−1/2 is written where the mesh width is equal to and where A Ceschino expansion is nothing but the Taylor's expansion for the derivatives where the binomial coefficients are defined as Our interpolation strategy is based on two-and three-point coupling relations obtained directly from the Ceschino polynomial expansion (4).Two points relations are used at the extremities of the domain and three-point relations are used inside.Cubic Hermite polynomials will also be introduced to perform the interpolation operation.

Two-points Ceschino coupling relations
Our relations are coupling the first N derivatives of f (x), with N = 1 leading to a cubic interpolation strategy.We set J = 2N in Eq. ( 4), leading to a truncation error of order 2N + 1 if k = 0. We next perform a linear combination of the first N components M (k) m+1/2 , introducing coefficients θ k .The linear combination permits to maintain the order of the truncation error to 2N + 1.We write After permutation of the two summations with the corresponding indices j and k in the right-hand-side, we get visiting We choose coefficients θ j in such a way that and we define coefficients θk as We have obtained our (2N + 1)-th order two-points Ceschino coupling relations as where the O(Δx m ) 2N+1 error term is not given.We need to determine a set of 2(N + 1) coefficients θ k and θk .Equations ( 8) and ( 9) permit to determine 2N + 1 of them, leaving θ 0 to be fixed.However, all values of θ 0 leads to valid solutions, making this choice arbitrary.We have chosen θ 0 = 1/(Δx m ) 2 in order to simplify the resulting mathematical formalism.
In the specific case of cubic Ceschino interpolation, we set N = 1, so that Eqs. ( 8) and ( 9) reduce to

Three-points Ceschino coupling relations
The three-points Ceschino coupling relations span two consecutive regions along the X axis, as depicted in Fig. 1.We set J = 3N in Eq. ( 4), leading to a truncation error of order 3N + 1 if k = 0.The Ceschino expansion are written where the mesh widths are equal to We next perform a linear combination of the first m+3/2 , introducing coefficients βk and β k .The linear combination permits to maintain the order of the truncation error to 3N + 1.We write where the truncation error is a linear combination of O(Δx m ) 3N+1 and O(Δx m+1 ) 3N+1 .
After permutation of the two summations with the corresponding indices j and k in the right-hand-side, we get We choose coefficients βj and β j in such a way that and we define coefficients βk as We have obtained our (3N + 1)-th order three-points Ceschino coupling relations as We need to determine a set of 3(N + 1) coefficients βk , βk and β k .Equations ( 18) and ( 19) permit to determine 3N + 1 of them, leaving β0 and β 0 to be fixed.A first set of coefficients can be obtained by setting β0 = −1/(Δx m ) 2 and β 0 = 1/(Δx m+1 ) 2 .A second independent set can be obtained by setting β 0 = 1/(Δx m ) 3 and β 0 = 1/(Δx m+1 ) 3 .Any other consistent set is a linear combination of these two.In the specific case of cubic Ceschino interpolation, we set N = 1, so that Eqs. ( 17) and ( 18) reduce to visiting so that our two independent sets of coefficients are (22)

Interpolation with cubic Hermite polynomials
Knowledge of M (0) m+1/2 and the capability to easily obtain M (1) m+1/2 on each tabulated point x m+1/2 makes possible the interpolation of function f (x) at each values of the independent variable x with a cubic Hermite polynomial in x.Such polynomial guarantee that the interpolated value and first derivative of the dependent variable remains continuous in x over the complete domain.As pointed out by (Rozon et al., 1981), this continuity property of the first derivative is often required in numerical applications such as those based on perturbation theory.The first operation consists to solve a tridiagonal linear matrix system for obtaining the unknown vector M (1) = col{M (1) m+1/2 ; m = 0, M} over a M-region domain, considering the known values M (0) m+1/2 of f (x) at tabulation points x m+1/2 .The linear matrix system is made with the first independent set of coefficients from Eq. ( 21) for linking the unknowns inside the domain.We have selected the first set in order to obtain a symmetric C matrix with minimum powers of Δx m as coefficients.The first and last line coefficients are obtained from Eq. ( 12).Using coefficients from Eq. ( 12) with those from Eq. ( 22) leads to a singular C matrix.This last observation gives an additional clue for selecting three-point coefficients from Eq. ( 21).The linear system is written where the symmetric tridiagonal matrix is written and where the source term is written The solution of the linear system in Eq. ( 23) can be performed without pivoting, as matrix C is diagonally dominant.
We next introduce the cubic Hermite polynomials defined over a reference region −1/2 ≤ u ≤ 1/2.They are so that a function f (u) defined over this domain can be expressed as where −1/2 ≤ u ≤ 1/2.The above relation can be generalized to the interpolation of function f (x) at ξ over region m.We first perform the change of variable where visiting

Introduction of interpolation factors
Interpolation factors are useful to interpolate a large number of dependent variables at a unique value ξ of the independent variable.The interpolation factors are function only of the tabulated abscissas {x m+1/2 ; m = 0, M} and on the interpolation abscissa x.Using interpolation factors {t m+1/2 (ξ) ; m = 0, M}, an interpolated dependent variable I{ f (x); ξ} of f (ξ) is obtained from where Interpolation factors can be obtained if the interpolation operation is distributive, that is, if it can be distributed to the sum of two functions f (x) and g(h) according to The simplest form of interpolation factors are those corresponding to linear Lagrange interpolation.In this case, the interpolated value of f (x), with x m−1/2 ≤ ξ ≤ x m+1/2 , is given by Eq. ( 30) with Similar interpolation factors exist for cubic Ceschino interpolation and can be obtained with the following procedure.The source term defined in Eq. ( 25) can be written in matrix form as where The interpolated value of f (ξ), with x m−1/2 ≤ ξ ≤ x m+1/2 , is therefore given by the relation where and The vector T(ξ) = {t m+1/2 (ξ) ; m = 0, M} of interpolation factors is obtained after transposition of Eq. ( 36), leading to

Matlab scripts and numerical examples
Two Matlab scripts are proposed in Appendices A and B as prototypes of the cubic Ceschino interpolation method.The first script, alterp() is used to obtain the terp factors corresponding to an interpolation (if lderiv=false) or to a derivation (if lderiv=true).
The second script, alteri() is used to obtain the terp factors corresponding to the definite integration of f (x).The following Matlab session is an example of interpolation similar to the spline Matlab tutorial.Execution of the above script leads to Fig. 2. Similarly, the first derivative of f (x) = sin(x) can be computed by setting lderiv = true, as described in the following Matlab session.Execution of the above script leads to Fig. 3.We observe that the order of the numerical derivation approximation is less than the order of the interpolation, as expected.The higher derivation errors are observed at extremities of the domain, where two-point Ceschino coupling relation are used.

Conclusion
We have presented a straightforward numerical technique based on Ceschino polynomial expansion.Three applications of this approach permit to perform interpolation, derivation and definite integration of tabulated data.Equation ( 36) is efficient to interpolate few dependent variables over a large number of points ξ.Equation (39) introduces the concept of interpolation factors and is efficient to interpolate a large number of dependent variables over a few number of points ξ.Matlab scripts are provided as basic implementation of the Ceschino interpolating method.