Lagrangian Polynomial Interpolation with R

Part 4 of 9 in the series Numerical Analysis

Polynomial interpolation is the method of determining a polynomial that fits a set of given points. There are several approaches to polynomial interpolation, of which one of the most well known is the Lagrangian method. This post will introduce the Lagrangian method of interpolating polynomials and how to perform the procedure in R.

Lagrangian Polynomial Interpolation

The Lagrangian method of polynomial interpolation uses Lagrangian polynomials to fit a polynomial to a given set of data points. The Lagrange interpolating polynomial is given by the following theorem:

For a set of data points (x_0, y_0), (x_1, y_1), \cdots, (x_n, y_n) with no duplicate x and there exists a function f which evaluates to these points, then there is a unique polynomial P(x) with degree \leq n also exists. The polynomial is given by:

P(x) = f(x_o)L_{n,0}(x) + \cdots + f(x_n)L_{n,n}(x) = \sum^n_{k=0} f(x_k) L_{n,k}(x)

Where each k in k = 0, 1, \cdots, n is:

L_{n,k} = \frac{(x - x_0)(x - x_1) \cdots (x - x_{k-1})(x - x_{k+1}) \cdots (x - x_n)}{(x_k - x_0)(x_k - x_1) \cdots (x_k - x_{k-1})(x_k - x_{k+1}) \cdots (x_k - x_n)} = \prod_{\stackrel{i \neq k}{i=0}}^{n} \frac{(x - x_i)}{(x_k - x_i)}
Lagrangian Polynomial Interpolation Example

Consider the following set of data points:

x y
0 7
2 11
3 28
4 63

The graph of the points looks like:

x <- c(0, 2, 3, 4)
y <- c(7, 11, 28, 63)

dat <- data.frame(cbind(x, y))

ggplot(dat, aes(x=x, y=y)) + 
  geom_point(size=5, col='blue')

plotted x-y points

We wish to find a polynomial that passes through these points. Note that much of the intermediate algebra will be skipped for the sake of brevity.

Start with L_1(x):

\large{L_1(x) = \frac{(x-2)(x-3)(x-4)}{(0-2)(0-3)(0-4)} = -\frac{1}{24}(x-2)(x-3)(x-4)}
L_2(x) is then:
\large{L_2(x) = \frac{(x-0)(x-3)(x-4)}{(2-0)(2-3)(2-4)} = \frac{1}{4}x(x-3)(x-4)}

Then finally, L_3(x) and L_4(x):

\large{L_3(x) = \frac{(x-0)(x-2)(x-4)}{(3-0)(3-2)(3-4)} = -\frac{1}{3}x(x-2)(x-4)}
\large{L_4(x) = \frac{(x-0)(x-2)(x-3)}{(4-0)(4-2)(4-3)} = \frac{1}{8}x(x-2)(x-3)}

These polynomials are then multiplied by they corresponding y value and added together, resulting in the following polynomial:

7\left(-\frac{1}{24}(x-2)(x-3)(x-4)\right) + 11\left(\frac{1}{4}x(x-3)(x-4)\right) - 28\left(-\frac{1}{3}x(x-2)(x-4)\right) + 63\left(\frac{1}{8}x(x-2)(x-3)\right)

Simplifying this equation yields:

\Large{x^3 - 2x + 7}
Polynomial Interpolation in R

We can implement a function in R to perform Lagrangian interpolation. The rSymPy package, which provides an R interface for the SymPy CAS (computer algebra system), will be used as we are interested in the symbolic form of the resulting polynomial.

library(rSymPy)

Our implementation of the Lagrangian polynomial interpolation method, lagrange.poly(), proceeds as follows:

lagrange.poly <- function(x, y) {

  l <- list() # List to store Lagrangian polynomials L_{1,2,3,4}
  k <- 1
  
  for (i in x) {
    # Set the numerator and denominator of the Lagrangian polynomials to 1 and build them up
    num <- 1
    denom <- 1

    # Remove the current x value from the iterated list
    p <- x[! x %in% i]

    # For the remaining points, construct the Lagrangian polynomial by successively 
    # appending each x value
    for (j in p) {
      num <- paste(num, "*", "(", 'x', " - ", as.character(j), ")", sep = ", collapse = ")
      denom <- paste(denom, "*", "(", as.character(i)," - ", as.character(j), ")", sep = ", collapse = ")
    }
  
    # Set each Lagrangian polynomial in rSymPy to simplify later.
    l[k] <- paste("(", num, ")", "/", "(", denom, ")", sep = ", collapse = ")
    k <- k + 1
  }
  
  # Similar to before, we construct the final Lagrangian polynomial by successively building 
  # up the equation by iterating through the polynomials L_{1,2,3,4} and the y values 
  # corresponding to the x values.
  eq <- 0
  
  for (i in 1:length(y)) {
    eq <- paste(eq, '+', as.character(y[i]), "*", l[[i]], sep = ", collapse = ")
  }
  
  # Define x variable for rSymPy to simplify
  x <- Var('x')
  
  # Simplify the result with rSymPy and return the polynomial
  return(sympy(paste("simplify(", eq, ")")))
}

With the x and y values used previously, perform Lagrangian polynomial interpolation to see if our function returns the same polynomial we found earlier.

lagrange.poly(x, y)
## [1] "7 - 2*x + x**3"

The polynomial x^3 - 2x + 7 resulting from our lagrange.poly() function matches what we found previously.

Lagrangian polynomial interpolation can also be performed in R using the poly.calc() function from the polynom package.

library(polynom)

Use the poly.calc() function to interpolate the data points using the Lagrangian method.

poly.calc(x, y)
## 7 - 2*x + x^3

We can plot this polynomial against the original data points to confirm visually that the polynomial does indeed pass through the points.

f <- function(x) {
  return(x^3 - 2 * x + 7)
}

ggplot(dat, aes(x=x, y=y)) + 
  geom_point(size=5, col='blue') + 
  stat_function(fun = f, size=1.25, alpha=0.4)

plotted points with polynomial resulting from polynomial interpolation method

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.

Series Navigation<< Neville’s Method of Polynomial InterpolationThe Newton-Raphson Root-Finding Algorithm in R >>
1 Comment

Post a Comment