Ch4.1: Polynomial Interpolation

Polynomials

A common application is to approximate functions
Recall Taylor polynomials from Calculus II:

\[ f(x) \cong f(a) + f'(a)(x-a) + \frac{f''(a)}{2}(x-a)^2 + \cdots + \frac{f^{(n)}(a)}{n!}(x-a)^n \]

title

Taylor Polynomials

\[ f(x) \cong f(a) + f'(a)(x-a) + \frac{f''(a)}{2}(x-a)^2 + \cdots + \frac{f^{(n)}(a)}{n!}(x-a)^n \]

title

Linear Interpolation

Fit equation of line to two data points \( (x_1,y_1), ~(x_2,y_2) \)
For \( y = mx + b \), need \( m \) and \( b \)

linterp <- function (x1 , y1 , x2 , y2) {
m <- (y2 - y1) / (x2 - x1)
b <- y2 - m * x2
## Convert into a form suitable for horner
return (c(b, m))  }

Horner's Method (Ch1.3.2)

Nested polynomial evaluation

\[ \begin{align*} f(x) & = 1 - 2x + 3x^2 +4x^3 +5x^4 \\ & = 1 + x*(-2 + x*(3 + x*(4 + x*5)) \end{align*} \]

horner <- function(x, coefs) {
  y <- rep(0, length(x))
  for(i in length(coefs):1)
    y <- coefs[i] + x*y
  return (y)  }

Example: Linear Interpolation

u <- c(3, 5)
plot(u,type="l")

plot of chunk unnamed-chunk-3

(p<-linterp(1,3,2,5))

[1] 1 2

horner(10,p)

[1] 21

Higher Order Polynomial Interpolation

Fit degree \( n-1 \) polynomial \( p(x) \),

\[ p(x) = \beta_n x^n + \beta_{n-1} x^{n-1} + \cdots + \beta_1 x + \beta_0 \]

to \( n \) data points

\[ \{(x_1,y_1), (x_2,y_2), \ldots, (x_n, y_n) \} \]

Higher Order Polynomial Interpolation

Six data points, degree 5 polynomial

title

Higher Order Polynomial Interpolation

Our polynomial is

\[ p(x) = \beta_n x^n + \beta_{n-1} x^{n-1} + \cdots + \beta_1 x + \beta_0 \]

Then \( p(x_i) = y_i \) for each \( i \), which can be expressed as a matrix-vector equation:

title

Higher Order Polynomial Interpolation

polyinterp <- function(x,y) {
  if(length (x) != length (y))
   stop("Length of x & y vectors must be same")
  n <- length (x) - 1
  vandermonde <- rep(1,length (x))
  for(i in 1:n) {
    xi <- x^i
    vandermonde <- cbind(vandermonde,xi) }
  beta <- solve(vandermonde,y)

  names(beta) <- NULL
  return(beta) }

Higher Order Polynomial Interpolation

xdata <- c(0, 0.3, 0.8, 1.1, 1.6, 2.3)
ydata <- c(0.6 ,0.67, 1.01, 1.35, 1.47, 1.25)
(p <- polyinterp(xdata, ydata))

[1]  0.6000000  0.5910395 -2.5174715  5.4280961 -3.5892186  0.7303129

title

Higher Order Polynomial Interpolation

xdata <- c(0, 0.3, 0.8, 1.1, 1.6, 2.3)
ydata <- c(0.6 ,0.67, 1.01, 1.35, 1.47, 1.25)
(p <- polyinterp(xdata, ydata))

[1]  0.6000000  0.5910395 -2.5174715  5.4280961 -3.5892186  0.7303129

poly <- function(x) {p[1] + p[2]*x + p[3]*x^2 + p[4]*x^3 + p[5]*x^4 + p[6]*x^5}
poly(xdata)

[1] 0.60 0.67 1.01 1.35 1.47 1.25

Higher Order Polynomial Interpolation

polyinterpplot <- function(x,y) {
  n = length(x); d <- x[n] - x[1]
  p <- polyinterp(x, y) 
  N <- 100; t <- rep(0, N); g <- rep(0, N) 
  t[1] <- x[1]; g[1] <- y[1]; 
  for(i in 1:N) { 
     t[i+1] <- t[i] + d/N
     g[i+1] <- poly(t[i+1]) }
  ymin = min(g); ymax = max(g)
  plot(t,g,                
        main = "Interpolating Polynomial",
        xlab = "x", ylab = "y",     
        type = "l",col="blue", 
        xlim = c(x[1],x[n]),   
        ylim = c(ymin,ymax)   ) 
  lines(x,y,type = "p",col="red")}

Higher Order Polynomial Interpolation

xdata <- c(0, 0.3, 0.8, 1.1, 1.6, 2.3)
ydata <- c(0.6 ,0.67, 1.01, 1.35, 1.47, 1.25)

title

(p <- polyinterpplot(xdata, ydata))

plot of chunk unnamed-chunk-11

NULL

Polynomials

Recall Taylor polynomials from Calculus II:

\[ \begin{align*} f(x) & \cong \sum_{k=0}^{n} \frac{f^{(k)}(a)}{k!} (x-a)^k \\ & = f(a) + f'(a)(x-a) + \frac{f''(a)}{2}(x-a)^2 + \cdots + \frac{f^{(n)}(a)}{n!}(x-a)^n \end{align*} \]

Error Formula:

\[ E(x) = f(x) - p(x) = \frac{(x-a)^{n+1}}{(n+1)!}f^{(n+1)}(\xi(x)) \]

where \( \xi(x) \) between \( a \) and \( x \).

Higher Order Polynomial Interpolation

Fit degree \( n-1 \) polynomial \( p(x) \),

\[ p(x) = \beta_n x^n + \beta_{n-1} x^{n-1} + \cdots + \beta_1 x + \beta_0 \]

to \( n \) data points

\[ \{(x_1,y_1), (x_2,y_2), \ldots, (x_n, y_n) \} \]

Error Formula:

\[ E(x) = f(x) - p(x) = \frac{(x-x_1)(x-x_2)\cdots (x-x_n)}{n!}f^{(n)}(\xi(x)) \]

where \( x_1 \leq \xi(x) \leq x_n \).