Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Part of 6 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 $n$th 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.

[latex display=”true”] 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}) [/latex]

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.

[latex display=”true”] \large{f[x_i] = f(x_i)} [/latex]

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}$.

[latex display=”true”] \large{f[x_i, x_{i+1}] = \frac{f[x_{i+1}] – f[x_i]}{x_{i+1 – x_i}}} [/latex]

The second divided difference follows:

[latex display=”true”] \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}} [/latex]

This iteration continues until the $n$th divided difference:

[latex display=”true”] \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}} [/latex]

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

[latex display=”true”] 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}) [/latex]

###### 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.

[latex display=”true”] 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 [/latex]

Then, $f[x_1, x_2]$ is computed.

[latex display=”true”] 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 [/latex]

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

[latex display=”true”] 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 [/latex]

$f[x_2, x_3]$ is then found.

[latex display=”true”] 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 [/latex]

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

[latex display=”true”] 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 [/latex]

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

[latex display=”true”] 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 [/latex]

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

[latex display=”true”] 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) [/latex]

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)


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)
## $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.

The post Divided Differences Method of Polynomial Interpolation appeared first on Aaron Schlegel.