Neville’s Method of Polynomial Interpolation

Part 3 of 9 in the series Numerical Analysis

Neville’s method evaluates a polynomial that passes through a given set of x and y points for a particular x value using the Newton polynomial form. Neville’s method is similar to a now-defunct procedure named Aitken’s algorithm and is based on the divided differences recursion relation (“Neville’s Algorithm”, n.d).

It was stated before in a previous post on Lagrangian polynomial interpolation that there exists a Lagrange polynomial that passes through points y_1, y_2, \cdots, y_k where each is a distinct integer and 0 \leq y_i \leq n at corresponding x values x_0, x_1, x_2, \cdots, x_n. The k points y_1, y_2, \cdots, y_k are denoted P_{y_1, y_2, \cdots, y_k}(x).

Neville’s Method

Neville’s method can be stated as follows:

Let a function f be defined at points x_0, x_1, \cdots, x_k where x_j and x_i are two distinct members. For each k, there exists a Lagrange polynomial P that interpolates the function f at the k + 1 points x_0, x_1, \cdots, x_k. The kth Lagrange polynomial is defined as:

\large{P(x) = \frac{(x - x_j) P_{0,1,\cdots,j-1,j+1,\cdots,k}(x) - (x - x_i) P_{0,1,\cdots,i-1,i+1,\cdots,k}(x)}{(x_i - x_j)}}

The P_{0,1,\cdots,j-1,j+1,\cdots,k} and P_{0,1,\cdots,i-1,i+1,\cdots,k} are often denoted \hat{Q} and Q, respectively, for ease of notation.

\large{P(x) = \frac{(x - x_j) \hat{Q}(x) - (x - x_i) Q(x)}{(x_i - x_j)}}

The interpolating polynomials can thus be generated recursively, which we will see in the following example:

Neville’s Method Example

Consider the following table of x and corresponding y values.

x y
8.1 16.9446
8.3 17.56492
8.6 18.50515
8.7 18.82091

Suppose we are interested in interpolating a polynomial that passes through these points to approximate the resulting y value from an x value of 8.4.

We can construct the interpolating polynomial approximations using the function above:

Q_{1,1} = \frac{(8.4 - x_0)Q_{1,0} - (8.4 - x_1) Q_{0,0}}{x_1 - x_0} = \frac{(8.4 - 8.1)(17.56492) - (8.4 - 8.3)(16.9446)}{8.3 - 8.1} = 17.87508
Q_{2,1} = \frac{(8.4 - x_1)Q_{2,0} - (8.4 - x_2)Q_{1,0}}{(x_2 - x_1)} = \frac{(8.4 - 8.3)(18.50515) - (8.4 - 8.6)(17.56492)}{(8.6 - 8.3)} = 17.87833
Q_{3,1} = \frac{(8.4 - x_2)Q_{3,0} - (8.4 - x_3)Q_{2,0}}{(x_3 - x_2)} = \frac{(8.4 - 8.6)(18.82091) - (8.4 - 8.7)(18.50515)}{(8.7 - 8.6)} = 17.87363

The approximated values Q_{1,1} = 17.87508, Q_{2,1} = 17.87833, Q_{3,1} = 17.87363 are then used in the next iteration.

Q_{2,2} = \frac{(8.4 - x_0)Q_{2,1} - (8.4 - x_2)Q_{1,1}}{(x_2 - x_0)} = \frac{(8.4 - 8.1)(17.87833) - (8.4 - 8.6)(17.87508)}{(8.6 - 8.1)} = 17.87703
Q_{3,2} = \frac{(8.4 - x_1)Q_{3,1} - (8.4 - x_3)Q_{2,1}}{(x_3 - x_1)} = \frac{(8.4 - 8.3)(17.87363) - (8.4 - 8.7)(17.87833)}{(8.7 - 8.3)} = 17.877155

Then the final iteration yields the approximated y value for the given x value.

Q_{3,3} = \frac{(8.4 - x_0)Q_{3,2} - (8.4 - x_3)Q_{2,2}}{(x_3 - x_0)} = \frac{(8.4 - 8.1)(17.877155) - (8.4 - 8.7)(17.87703)}{(8.7 - 8.1)} = 17.8770925

Therefore 17.8770925 is the approximated value at the point 8.4.

Neville’s Method in R

The following function is an implementation of Neville’s method for interpolating and evaluating a polynomial.

poly.neville <- function(x, y, x0) {
  
  n <- length(x)
  q <- matrix(data = 0, n, n)
  q[,1] <- y
  
  for (i in 2:n) {
    for (j in i:n) {
      q[j,i] <- ((x0 - x[j-i+1]) * q[j,i-1] - (x0 - x[j]) * q[j-1,i-1]) / (x[j] - x[j-i+1])
    }
  }
  
  res <- list('Approximated value'=q[n,n], 'Neville iterations table'=q)
  return(res)
}

Let’s test this function to see if it reports the same result as what we found earlier.

x <- c(8.1, 8.3, 8.6, 8.7)
y <- c(16.9446, 17.56492, 18.50515, 18.82091)

poly.neville(x, y, 8.4)
## $`Approximated value`
## [1] 17.87709
## 
## $`Neville iterations table`
##          [,1]     [,2]     [,3]     [,4]
## [1,] 16.94460  0.00000  0.00000  0.00000
## [2,] 17.56492 17.87508  0.00000  0.00000
## [3,] 18.50515 17.87833 17.87703  0.00000
## [4,] 18.82091 17.87363 17.87716 17.87709

The approximated value is reported as 17.87709, the same value we calculated previously (minus a few decimal places). The function also outputs the iteration table that stores the intermediate results.

The pracma package contains the neville() function which also performs Neville’s method of polynomial interpolation and evaluation.

library(pracma)
neville(x, y, 8.4)
## [1] 17.87709

The neville() function reports the same approximated value that we found with our manual calculations and function.

References

Burden, R. L., & Faires, J. D. (2011). Numerical analysis (9th ed.). Boston, MA: Brooks/Cole, Cengage Learning.

Cheney, E. W., & Kincaid, D. (2013). Numerical mathematics and computing (6th ed.). Boston, MA: Brooks/Cole, Cengage Learning.

Neville’s algorithm. (2016, January 2). In Wikipedia, The Free Encyclopedia. From https://en.wikipedia.org/w/index.php?title=Neville%27s_algorithm&oldid=697870140

Series Navigation<< Divided Differences Method of Polynomial InterpolationLagrangian Polynomial Interpolation with R >>
1 Comment

Post a Comment