Menu
Home Explore People Places Arts History Plants & Animals Science Life & Culture Technology
On this page
Neville's algorithm
Used for polynomial interpolation

In mathematics, Neville's algorithm is an algorithm used for polynomial interpolation that was derived by the mathematician Eric Harold Neville in 1934. Given n + 1 points, there is a unique polynomial of degree ≤ n which goes through the given points. Neville's algorithm evaluates this polynomial.

Neville's algorithm is based on the Newton form of the interpolating polynomial and the recursion relation for the divided differences. It is similar to Aitken's algorithm (named after Alexander Aitken), which is nowadays not used.

We don't have any images related to Neville's algorithm yet.
We don't have any YouTube videos related to Neville's algorithm yet.
We don't have any PDF documents related to Neville's algorithm yet.
We don't have any Books related to Neville's algorithm yet.
We don't have any archived web articles related to Neville's algorithm yet.

The algorithm

Given a set of n+1 data points (xi, yi) where no two xi are the same, the interpolating polynomial is the polynomial p of degree at most n with the property

p(xi) = yi for all i = 0,...,n

This polynomial exists and it is unique. Neville's algorithm evaluates the polynomial at some point x.

Let pi,j denote the polynomial of degree j − i which goes through the points (xk, yk) for k = i, i + 1, ..., j. The pi,j satisfy the recurrence relation

p i , i ( x ) = y i , {\displaystyle p_{i,i}(x)=y_{i},\,} 0 ≤ i ≤ n , {\displaystyle 0\leq i\leq n,\,}
p i , j ( x ) = ( x − x i ) p i + 1 , j ( x ) − ( x − x j ) p i , j − 1 ( x ) x j − x i , {\displaystyle p_{i,j}(x)={\frac {(x-x_{i})p_{i+1,j}(x)-(x-x_{j})p_{i,j-1}(x)}{x_{j}-x_{i}}},\,} 0 ≤ i < j ≤ n . {\displaystyle 0\leq i<j\leq n.\,}

This recurrence can calculate p0,n(x), which is the value being sought. This is Neville's algorithm.

For instance, for n = 4, one can use the recurrence to fill the triangular tableau below from the left to the right.

p 0 , 0 ( x ) = y 0 {\displaystyle p_{0,0}(x)=y_{0}\,}
p 0 , 1 ( x ) {\displaystyle p_{0,1}(x)\,}
p 1 , 1 ( x ) = y 1 {\displaystyle p_{1,1}(x)=y_{1}\,} p 0 , 2 ( x ) {\displaystyle p_{0,2}(x)\,}
p 1 , 2 ( x ) {\displaystyle p_{1,2}(x)\,} p 0 , 3 ( x ) {\displaystyle p_{0,3}(x)\,}
p 2 , 2 ( x ) = y 2 {\displaystyle p_{2,2}(x)=y_{2}\,} p 1 , 3 ( x ) {\displaystyle p_{1,3}(x)\,} p 0 , 4 ( x ) {\displaystyle p_{0,4}(x)\,}
p 2 , 3 ( x ) {\displaystyle p_{2,3}(x)\,} p 1 , 4 ( x ) {\displaystyle p_{1,4}(x)\,}
p 3 , 3 ( x ) = y 3 {\displaystyle p_{3,3}(x)=y_{3}\,} p 2 , 4 ( x ) {\displaystyle p_{2,4}(x)\,}
p 3 , 4 ( x ) {\displaystyle p_{3,4}(x)\,}
p 4 , 4 ( x ) = y 4 {\displaystyle p_{4,4}(x)=y_{4}\,}

This process yields p0,4(x), the value of the polynomial going through the n + 1 data points (xi, yi) at the point x.

This algorithm needs O(n2) floating point operations to interpolate a single point, and O(n3) floating point operations to interpolate a polynomial of degree n.

The derivative of the polynomial can be obtained in the same manner, i.e:

p i , i ′ ( x ) = 0 , {\displaystyle p'_{i,i}(x)=0,\,} 0 ≤ i ≤ n , {\displaystyle 0\leq i\leq n,\,}
p i , j ′ ( x ) = ( x − x i ) p i + 1 , j ′ ( x ) + p i + 1 , j ( x ) − ( x − x j ) p i , j − 1 ′ ( x ) − p i , j − 1 ( x ) x j − x i , {\displaystyle p'_{i,j}(x)={\frac {(x-x_{i})p'_{i+1,j}(x)+p_{i+1,j}(x)-(x-x_{j})p'_{i,j-1}(x)-p_{i,j-1}(x)}{x_{j}-x_{i}}},\,} 0 ≤ i < j ≤ n . {\displaystyle 0\leq i<j\leq n.\,}

Alternate notation easier for computer implementation

In the above formulae, if we take the degree of the successive interpolating polynomials d = j − i and change the notation to pd,i,

p 0 , i ( x ) = y i , {\displaystyle p_{0,i}(x)=y_{i},} d = 0 {\displaystyle d=0}
p d , i ( x ) = ( x − x i ) p d − 1 , i + 1 ( x ) − ( x − x i + d ) p d − 1 , i ( x ) x i + d − x i , {\displaystyle p_{d,i}(x)={\frac {(x-x_{i})p_{d-1,i+1}(x)-(x-x_{i+d})p_{d-1,i}(x)}{x_{i+d}-x_{i}}},} 1 ≤ d ≤ n , 0 ≤ i ≤ n − d {\displaystyle 1\leq d\leq n,0\leq i\leq n-d}

The final value pn,0 (in this notation) is the required interpolated value.

Since the number of computed items i.e. the range of i decreases with each successive d, a linear array can be used for memory efficiency with pi being overwritten and d being ignored. (For example: [1])

The derivative (using the product rule) can be computed likewise as:

p 0 , i ′ ( x ) = 0 , {\displaystyle p'_{0,i}(x)=0,} d = 0 {\displaystyle d=0}
p d , i ′ ( x ) = ( x − x i ) p d − 1 , i + 1 ′ ( x ) + p d − 1 , i + 1 ( x ) − ( x − x i + d ) p d − 1 , i ′ ( x ) − p d − 1 , i ( x ) x i + d − x i , {\displaystyle p'_{d,i}(x)={\frac {(x-x_{i})p'_{d-1,i+1}(x)+p_{d-1,i+1}(x)-(x-x_{i+d})p'_{d-1,i}(x)-p_{d-1,i}(x)}{x_{i+d}-x_{i}}},} 1 ≤ d ≤ n , 0 ≤ i ≤ n − d {\displaystyle 1\leq d\leq n,0\leq i\leq n-d}

As before, p′n,0 (in this notation) is the derivative.

As this depends on the successive computed values of p also for each d, it may be computed within the same loop. If linear arrays for p and p′ are used for efficiency, the p′ values should be computed before the p values are overwritten.

Application to numerical differentiation

Lyness and Moler showed in 1966 that using undetermined coefficients for the polynomials in Neville's algorithm, one can compute the Maclaurin expansion of the final interpolating polynomial, which yields numerical approximations for the derivatives of the function at the origin. While "this process requires more arithmetic operations than is required in finite difference methods", "the choice of points for function evaluation is not restricted in any way". They also show that their method can be applied directly to the solution of linear systems of the Vandermonde type.