## Neville’s Method of Polynomial Interpolation

- Simpson’s Rule for Approximating Definite Integrals in R
- Numerical Differentiation with Finite Differences in R
- Divided Differences Method of Polynomial Interpolation
- Neville’s Method of Polynomial Interpolation
- Lagrangian Polynomial Interpolation with R
- The Newton-Raphson Root-Finding Algorithm in R
- The Secant Method Root-Finding Algorithm in R
- The Bisection Method of Root-Finding with R
- The Trapezoidal Rule of Numerical Integration in R

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.8770925Therefore 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

## Neville’s Method of Polynomial Interpolation – Mubashir Qasim

July 19, 2017 at 4:19 pm[…] article was first published on R – Aaron Schlegel, and kindly contributed to R-bloggers) Part 1 of 5 in the series Numerical […]