Divided Differences Method of Polynomial Interpolation

Part 2 of 9 in the series Numerical Analysis

The divided differences method is a numerical procedure for interpolating a polynomial given a set of points. Unlike Neville’s method, which is used to approximate the value of an interpolating polynomial at a given point, the divided differences method constructs the interpolating polynomial in Newton form.

Consider a table of values of x and corresponding values of a function f(x) at those points:

x y
x_0 f(x_0)
x_1 f(x_1)
\vdots \vdots
x_n f(x_n)

Also assume that P_n(x) is the nth Lagrangian polynomial that corresponds with the function f at these points. The polynomial P_n(x) can be expressed using the divided differences of the function f with respect to the x-values.

P_n(x) = a_0 + a_1(x - x_0) + a_2(x - x_0)(x - x_1) + \cdots + a_n(x - x_0) \cdots (x - x_{n-1})

Therefore the constants a_0, a_1, \cdots, a_n must be found to construct the polynomial. To find these constants, the divided differences are recursively generated until n iterations have been completed. We start with the zeroth divided difference of the function f with respect to x_i, which is the value of f at that point. Bracket notation is introduced to distinguish the divided differences.

\large{f[x_i] = f(x_i)}

Now begins the recursive generation of the divided differences. The first divided difference is then the function f with respect to the values x_i and x_{i+1}.

\large{f[x_i, x_{i+1}] = \frac{f[x_{i+1}] - f[x_i]}{x_{i+1 - x_i}}}

The second divided difference follows:

\large{f[x_i, x_{i+1}, x_{i+2}] = \frac{f[x_{i+1},x_{i+2}] - f[x_i, x_{i+1}]}{x_{i+2} - x_i}}

This iteration continues until the nth divided difference:

\large{f[x_0, x_1, \cdots, x_n] = \frac{f[x_1, x_2, \cdots, x_n] - f[x_0, x_1, \cdots, x_n]}{x_n - x_0}}

Thus the interpolating polynomial resulting from the divided differences method takes the form:

P_n(x) = f[x_0] + f[x_0, x_1](x - x_0) + f[x_0, x_1, x_2](x - x_0)(x - x_1) + \cdots + f[x_0, x_1, x_2, \cdots, x_n](x - x_0)(x - x_1) \cdots (x - x_{n-1})
Divided Differences Example

Consider the following set of x points and the corresponding values of a function f at those points. This data set is the same from the previous post on Neville’s method. We will perform the divided differences method to find the interpolating polynomial of these points and approximate the polynomial at f(8.4).

x f(x)
8.1 16.9446
8.3 17.56492
8.6 18.50515
8.7 18.82091

First, f[x_0, x_1] is found.

f[x_0, x_1] = \frac{f(x_1) - f(x_0)}{x_1 - x_0} = \frac{17.56492 - 16.9446}{8.3 - 8.1} = 3.1016

Then, f[x_1, x_2] is computed.

f[x_1, x_2] = \frac{f(x_2) - f(x_1)}{x_2 - x_1} = \frac{18.50515 - 17.56492}{8.6 - 8.3} = 3.1341

With f[x_0, x_1] and f[x_1, x_2], we find f[x_0, x_1, x_2].

f[x_0, x_1, x_2] = \frac{f[x_1, x_2] - f[x_0, x_1]}{x_2 - x_0} = \frac{3.1341 - 3.1016}{8.6-8.1} = 0.065

f[x_2, x_3] is then found.

f[x_2, x_3] = \frac{f(x_3) - f(x_2)}{x_3 - x_2} = \frac{18.82091 - 18.50515}{8.7 - 8.6} = 3.1576

We can then compute f[x_1, x_2, x_3]:

f[x_1, x_2, x_3] = \frac{f[x_2, x_3] - f[x_1, x_2]}{x_3 - x_1} = \frac{3.1576 - 3.1341}{8.7 - 8.3} = 0.05875

Lastly, we calculate f[x_1, x_2, x_3, x_4]:

f[x_1, x_2, x_3, x_4] = \frac{f[x_1, x_2, x_3] - f[x_0, x_1, x_2]}{x_3 - x_0} = \frac{0.05875 - 0.065}{8.7-8.1} = -0.010147

The interpolating polynomial, P_3(x), is then constructed.

P_3(x) = 16.9446 + 3.1016(x - 8.1) + 0.065(x - 8.1)(x - 8.3) - 0.010417(x - 8.1)(x - 8.3)(x - 8.6)

We can see the interpolated polynomial passes through the points with the following:

p3x <- function(x) {
f <- 16.9446 + 3.1016 * (x - 8.1) + 0.065 * (x - 8.1) * (x - 8.3) - 0.010417 * (x - 8.1) * (x - 8.3) * (x - 8.6)
return(f)
}
x <- c(8.1, 8.3, 8.6, 8.7)
fx <- c(16.9446, 17.56492, 18.50515, 18.82091)
dat <- data.frame(cbind(x, fx))
ggplot(dat, aes(x=x, y=fx)) + 
geom_point(size=5, col='blue') + 
stat_function(fun = p3x, size=1.25, alpha=0.4)

Divided differences interpolated function

Using the function above, we can also see the interpolated polynomial resulting from the divided differences method returns the same approximated value of the function f, f(x) as Neville’s method.

p3x(8.4)
## [1] 17.87709
Divided Differences in R

The following is an implementation of the divided differences method of polynomial interpolation in R. The function utilizes the rSymPy library to build the interpolating polynomial and approximate the value of the function f for a given value of x.

divided.differences <- function(x, y, x0) {
require(rSymPy)
n <- length(x)
q <- matrix(data = 0, n, n)
q[,1] <- y
f <- as.character(round(q[1,1], 5))
fi <- ''
for (i in 2:n) {
for (j in i:n) {
q[j,i] <- (q[j,i-1] - q[j-1,i-1]) / (x[j] - x[j-i+1])
}
fi <- paste(fi, '*(x - ', x[i-1], ')', sep = '', collapse = '')
f <- paste(f, ' + ', round(q[i,i], 5), fi, sep = '', collapse = '')
}
x <- Var('x')
sympy(paste('e = ', f, collapse = '', sep = ''))
approx <- sympy(paste('e.subs(x, ', as.character(x0), ')', sep = '', collapse = ''))
return(list('Approximation from Interpolation'=as.numeric(approx), 
'Interpolated Function'=f, 
'Divided Differences Table'=q))
}

Let’s use the function to interpolate a function given the values of x and the corresponding values of a function f, f(x) and approximate the function when x = 8.4.

divided.differences(x, fx, 8.4)
## Loading required package: rSymPy
## Loading required package: rJython
## Loading required package: rJava
## Loading required package: rjson
## $`Approximation from Interpolation`
## [1] 17.87709
## 
## $`Interpolated Function`
## [1] "16.9446 + 3.1016*(x - 8.1) + 0.065*(x - 8.1)*(x - 8.3) + -0.01042*(x - 8.1)*(x - 8.3)*(x - 8.6)"
## 
## $`Divided Differences Table`
##          [,1]   [,2]    [,3]        [,4]
## [1,] 16.94460 0.0000 0.00000  0.00000000
## [2,] 17.56492 3.1016 0.00000  0.00000000
## [3,] 18.50515 3.1341 0.06500  0.00000000
## [4,] 18.82091 3.1576 0.05875 -0.01041667

The function returns the approximated value, the interpolated polynomial use to approximate the function and the divided differences table which contains the intermediate calculations as we saw earlier. I wasn’t able to find a function in R or a package that performs the divided differences method, so, unfortunately, I have nothing for which to compare our function. However, the results agree with our manual calculations so that should be satisfactory for this introductory purpose.

References

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

Series Navigation<< Numerical Differentiation with Finite Differences in RNeville’s Method of Polynomial Interpolation >>
1 Comment

Post a Comment