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.

Polynomial interpolation methods

The simplest form of interpolation polynomial which you can find at the beginning of all numerical analysis textbooks is Lagrange's interpolation polynomial:

This interpolation polynomial is simple, but rarely used. In practice, Neville's algorithm is used for polynomial interpolation. This algorithm represents the N-th degree polynomial P_0,1,..,N-1,N , which passes through N+1 points (from 0-th to N-th) as the sum of two (N-1)-th degree polynomials:

This formula is used recursively until we get polynomials P_i  which are calculated by using formula P_i  = y_i . This algorithm calculates the N-th degree interpolation polynomial in O(N ²) operations. The advantages of this algorithm are its simplicity and the fact that the algorithm could be easily expanded to calculate derivate values as well as values of interpolation polynomial. The algorithm's disadvantage is that it has a relatively low performance.

It should be noted that Neville's algorithm isn't the latest thing in polynomial interpolation. From the new methods in this area we can point out the barycentric interpolation method:

Generally, this formula gives a rational function passing through (x_i , y_i ) independently of the weighting coefficients we have chosen. Thus, these coefficients affect only the way in which the interpolant approximates the function between the points. Two facts are important for us: first, when we choose

the rational function r(x) becomes the N-th degree polynomial. Secondly, for some important special cases (equidistant or Chebyshev nodes) there is an O(N) algorithm calculating all weighting coefficients. Thus, the value of the interpolation polynomial is calculated N times faster than by using Neville's algorithm.

Note #1
It should be noted that the straightforward use of the barycentric formula could cause bad consequences. If x_j  = 0, and x is close (but not equal) to it, it is possible to get an overflow because of dividing by a value which is too small. Such cases occur very seldom but, in spite of that, much more often that they should do according to probability theory. But don't worry, the algorithm from this module works correctly, switching to an overflow-stable type of formula when necessary.

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.

Neville's algorithm

The algorithm is implemented in two subroutines: NevilleInterpolation and NevilleDifferentiation. The first one calculates the value of the interpolation polynomial, the second - the polynomial value and the values of the first, second and the third derivatives. These subroutines perform the interpolation on an arbitrary set of nodes. Their calculation complexity is O(N ²).

Barycentric interpolation on special grids

To interpolate a function on an equidistant or Chebyshev grid, we can use the barycentric algorithm. To interpolate the function on the equidistant grid, use the EquidistantPolynomialInterpolation subroutine. To interpolate the function on the Chebyshev nodes of the first kind (in this case , i = 0,... , N), use the Chebyshev1Interpolation subroutine. To interpolate the function on the Chebyshev nodes of the second kind (in this case , i = 0,... , N), use the Chebyshev2Interpolation subroutine.

Barycentric interpolation on arbitrary grids

This problem is solved by using the BarycentricInterpolation subroutine. The subroutine arguments are a set of nodes, their function values and the weighting coefficients which should be calculated by using the following formula:

Links

1. 'Barycentric Lagrange Interpolation'. Jean-Paul Berrut, Lloyd N. Trefethen.

Report a bug

C#

C# 1.0 source.
polinterpolation.csharp.zip - Polynomial interpolation

C++

C++ source.
polinterpolation.cpp.zip - Polynomial interpolation
ablas.zip - optimized basic linear algebra subroutines with SSE2 support (for C++ sources only)

Delphi

Delphi source.
Can be compiled under FPC (in Delphi compatibility mode).
polinterpolation.delphi.zip - Polynomial interpolation

Visual Basic 6

Visual Basic 6 source.
polinterpolation.vb6.zip - Polynomial interpolation

Zonnon ^beta

Zonnon source.
Zonnon is an experimental language developed at ETH Zurich.
See www.zonnon.ethz.ch for more information.
polinterpolation.zonnon.zip - Polynomial interpolation