## Lagrangian Polynomial Interpolation with R

- 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

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')

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)

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

## Lagrangian Polynomial Interpolation with R – Mubashir Qasim

July 14, 2017 at 12:18 am[…] article was first published on R – Aaron Schlegel, and kindly contributed to R-bloggers) Part 2 of 4 in the series Numerical […]