Polynomial interpolation

Polynomial interpolation is the most known one-dimensional interpolation method. Its advantages lies in its simplicity of realization and the good quality of interpolants obtained from it. However, it has several disadvantages (some of them will be considered later) and is lately hard-pressed by alternative interpolation methods: splines and rational functions. But, in spite of that, polynomial interpolation is still one of the main tools of numerical analysis.

    Barycentric representation of polynomial
    Choosing nodes
    Polynomial interpolation
    Polynomial fitting
    Manual entries
    Links

Barycentric representation of polynomial

It is known that any rational function may be represented in barycentric form:

The same formula may be used to represent polynomial, which is a special case of a rational function. Such representation have following properties: a) it is more stable than the power basis, b) after barycentric representation is built it takes only O(N) time to interpolate (instead of O(N²) for Neville-type algorithm), c) in several important cases (equidistant or Chebyshev grids) barycentric model can be built in O(N) time instead of O(N²).

For these reasons polynomial interpolation in ALGLIB is implemented on top of ratint unit, which implements rational interpolation/fitting. polint unit is used only to build barycentric model, which are processed with ratint subroutines (evaluated, differentiated, etc.).

Choosing nodes

The interpolation accuracy depends on the interpolation nodes. The equidistant grid, despite its simplicity and convenience, should not be used for interpolation because of two related reasons.

To understand the first reason, let's consider the function graph. Here are the function

(black) and the interpolation polynomial constructed by using 11 points of the equidistant grid (red). We can see that the interpolation accuracy is good in the center. The closer to the edges the less accurate the interpolation becomes. Although the polynomial passes through all the points, between them it deviates wildly from the function. We could expect that the interpolation error will decrease when we take more points. But in practice, the more n we take, the wilder the deviation of the interpolation polynomial becomes from the function (near the edges). If the number of points goes to infinity, the interpolation error will tend towards infinity too.

It is proved that there is a class of functions which cannot be interpolated by the polynomial on an equidistant grid. These functions have poles on the complex plane in the neighborhood of the interpolation interval (in particular, the function mentioned above has poles in points x = +i è x = -i). It should be understood that this effect is not a bug in the algorithm and not a consequence of the floating point errors. It is a fundamental property of an interpolation polynomial.

The second reason not to use the equidistant grid is associated with the first one. When using the equidistant grid, floating point errors could be accumulated and cause low interpolation quality. The reason is even if the function to be interpolated is a "good function" (i.e. hasn't any poles in the neighborhood of the interpolation interval), floating point errors disfigure its graph. These small distortions can make this function look like a "bad function", which will increase an error dramatically.

This problem has two solutions. If we cannot back away from the equidistant grid, we can use cubic splines or rational functions. If we can choose the nodes, we can interpolate the function using Chebyshev nodes (the node arrangement is presented later). On this grid in the overwhelming majority of cases the interpolation error decreases when increasing the number of points (specifically, it is true for all smooth functions). The floating point errors are also less tends to accumulate.

Polynomial interpolation

Interpolation polynomial, i.e. polynomial which passes through all points of data set, can be built with several subroutines:

PolynomialBuild, which builds interpolation polynomial using arbitrary set of nodes in O(N²) time (but interpolation has only O(N) complexity). Its result is BarycentricInterpolant structure, which may be used with subroutines from ratint unit.
PolynomialBuildEqDist, which builds interpolation polynomial using equidistant nodes. Its complexity is O(N).
PolynomialBuildCheb1 and PolynomialBuildCheb2, which build interpolation polynomial using Chebyshev nodes. Their complexity is O(N).

Polynomial fitting

Polynomial fitting is supported by PolynomialFit (unconstrained fitting without individual weights) and PolynomialFitWC (fitting with individual weights and constraints) subroutines. Arbitrary number of constraints on function value - f(xc)=yc - or its derivative - df(xc)/dx=yc - is supported.

Note #1
Excessive constraints may be inconsistent, especially when you work with simple models. See Reference manual, PolynomialFitWC comments, for more information on this subject.

Both subroutines use Chebychev basis with automatic scaling to reduce condition number. Because such fitting task is linear, Linear Least Squares solver with O(N·M²) complexity is used (here N is a number of points, M is a basis size).