Reliable determination of interpolating polynomials

This paper addresses a fundamental problem in mathematics and numerical analysis, that of determining a polynomial interpolant to specified data. The data is taken as consisting of a set of points (abscissae), at each of which is specified a function value. Additionally, at each point, any number of leading derivative values of the function may be given. Mathematically, this problem is solved. The classical Lagrangian interpolation formula applies in the derivative-free case, and the Newton form of the interpolating polynomial in general.Numerically, few reliable algorithms are available; most published algorithms concentrate on speed of computation. This paper describes an algorithm that delivers the required polynomial in Chebyshev form. It is based on the repeated use of the Newton representation, with a data ordering strategy and iterative refinement to improve accuracy, using a carefully devised merit function to measure success. The algorithm attempts to provide a polynomial that is stable in the sense of backward error analysis, i.e. that is exact for slightly perturbed data.Implementations of the algorithm have been in use since the early 1980s in the NAG Library and NPL's Data Approximation Subroutine Library (DASL). In addition to providing polynomial interpolants in their own right, these implementations are used as computational modules in the NAG and DASL routines for constrained least-squares polynomial data fitting.This paper constitutes the first detailed presentation of the algorithm.