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.
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 2) 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 2).
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 larger 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 significantly lower too.
You can use following functions to build barycentryc representation of the polynomial interpolant:
Functions mentioned above returnstructure as result. This structure can be used with following functions from subpackage:
Following functions are used to switch between barycentric representation and other representations:
Finally, following functions can be used for quick - O(N) - evaluation of interpolating polynomial on one of the special grids without creation of barycentric model:
ALGLIB Reference Manual includes following examples on polynomial interpolation:
ALGLIB package supports polynomial curve fitting, either unconstrained (chapter of ALGLIB User Guide.function) or constrained ( function). Arbitrary number of constraints on function value - f(xc)=yc - or its derivative - df(xc)/dx=yc - is supported. This functionality is described in more details in the corresponding
This article is intended for personal use only.
C++ source. MPFR/GMP is used.
GMP source is available from gmplib.org. MPFR source is available from www.mpfr.org.
Python version (CPython and IronPython are supported).
ALGLIB® - numerical analysis library, 1999-2013.