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 \(k\)th 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

The post Neville’s Method of Polynomial Interpolation appeared first on Aaron Schlegel.