Neville algorithm

Constructing the interpolating polynomial is somewhat tedious. If our goal is merely to get the interpolated value, and we don't care to know the coefficients of the polynomial, we may use the Neville algorithm. This algorithm starts from the requested interpolation point $x$ and generates a table of the form

$x_0$	$P_0(x)$
		$P_{01}(x)$
$x_1$	$P_1(x)$		$P_{012}(x)$
		$P_{12}(x)$		$P_{0123}(x)$
$x_2$	$P_2(x)$		$P_{123}(x)$		$P_{01234}(x)$
		$P_{23}(x)$		$P_{1234}(x)$
$x_3$	$P_3(x)$		$P_{234}(x)$
		$P_{34}(x)$
$x_4$	$P_4(x)$

where $P_{i,\ldots{},k}(x)$ is the polynomial interpolating the points $x_i,x_{i+1}\ldots{},x_j$, evaluated at the point $x$. That is to say, the second column comes from interpolating only one point each, meaning that it reproduces trivially the original values: $P_i(x) = y_i$. The third column comes from interpolating pairs of adjacent points. Notice that $P_{0,1}$ is the polynomial we constructed above. The fourth column comes from interpolating three contiguous points. We illustrated $P_{123}$ above. The last column has only one entry and comes from interpolating all the points.

The identity

$$P_{i,\ldots{},\ell}(x) = \frac{(x - x_j)P_{i,\ldots{},j-1,j... ... -(x - x_k)P_{i,\ldots{},k-1,k+1,\ldots{}\ell}(x)}{x_k - x_j}$$ (1)

is used to calculate the values in any column from pairs of values in one column to the left. For example the rule gives

$$P_{123}(x) = \frac{(x - x_3)P_{12}(x) - (x - x_1)P_{23}(x)}{x_1 - x_3}$$

The student should verify that $P_{123}(x_i) = y_i$ for $i = 1,2,3$, as it should.

Here we apply the Neville scheme to the example of the exponential integral, interpolated at $x = 0.15$.

0.1	-1.6228
		-1.22230
0.2	-0.8218		-1.18706
		-1.08135		-1.17642
0.3	-0.3027		-1.12320		-1.17186
		-0.91395		-1.13992
0.4	0.1048		-1.02289
		-0.76870
0.5	0.4542

(We carry an extra decimal digit in the intermediate calculations to reduce round-off errors.)

Compare the two tables. As promised, the second column is a copy of the $y_i$ values. The third column gives the result of doing a linear interpolation between the two neighboring points and evaluating the resulting function at $x = 0.15$. Notice that only the first value, $-1.2223$, is really an interpolation. The others in the third column are actually linear extrapolations back to 0.15 from interpolations on the other intervals. The fourth column gives the result of a quadratic interpolation of the central point on that row and its two neighbors. For $x = 0.15$ only the first value in that column is really an interpolation. The others are extrapolations.

There can be two interpolating values in a column. For example, if we had constructed a similar table for interpolation at $x = 0.25$, the first two values in the fourth column would be interpolations, since the polynomials $P_{012}(x)$ and $P_{123}(x)$ both cover the $12$ interval.

Which value gives the best answer? Polynomial extrapolations can be risky, so clearly we should use only the interpolating values. For the example above, the interpolating values can be read from the first number in each column, starting with the third column: $-1.2223$, $-1.1871$, $-1.1764$, $-1.1719$. Since we expect the estimate to improve as we increase the degree of the polynomial, at least as long as the degree isn't too high, our best result should be $-1.1719$, coming from the quartic polynomial interpolating all the points.

It is not always necessary to carry out the Neville iteration to the bitter end. To avoid disasters with high degree polynomials on evenly spaced points near the edges of a table, it is best to do interpolations with roughly the same number of tabulated values on either side of the interpolation point, which may make stopping advisable, before exhausting all the points in a table.

Finally, we need an estimate of the interpolation error in our best result without peeking at the exact value. (If we know it, why interpolate?) A crude estimate of the error is given by the magnitude of the change between the best value and the next best. In our case it is $\vert-1.1719 + 1.1764\vert$, or about $0.005$. How far off are we, really? The actual value to six figures is $-1.16409$, so the error in the last value is really $0.008$, nearly twice our estimated error. Well, knowing what we do about the logarithmic singularity in this function at the origin and the hazards of high-degree polynomial interpolations, we might have expected trouble. The moral of the story: make your error estimates with a good measure of humility.