Divided Differences Method of Polynomial Interpolation

July 27, 2017
By

[This article was first published on R – Aaron Schlegel, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
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 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.

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

To leave a comment for the author, please follow the link and comment on their blog: R – Aaron Schlegel.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Search R-bloggers

Sponsors

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)