Ch4.1: Polynomial Interpolation

• A common application is to approximate functions.
• Recall Taylor polynomials from Calculus II.
• Focused on the point \( (a,f(a)) \) on graph; see next slide.

\[ \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*} \]

• 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))  }

• Nested polynomial evaluation

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

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

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

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

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

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

horner(xdata,p)

[1] 0.60 0.67 1.01 1.35 1.47 1.25

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

HornerPlot <- function(x,y) {
  n = length(x)
  N <- 100; h <- (x[n] - x[1])/N
  t <- rep(0, N); f <- rep(0, N) 
  t[1] <- x[1]; f[1] <- y[1]; 
  p <- polyinterp(x, y)

  for(i in 1:N) { 
     t[i+1] <- t[i] + h 
     f[i+1] <- horner(t[i+1],p) }

  plot(t,f, 
    main = "Horner Interpolating Polynomial",
    xlab="x",ylab="y",type="l",lwd=2,col="red", 
    xlim = c(x[1],x[n]), ylim = c(min(f),max(f)) )
  lines(x,y,type = "p",col="blue",lwd=3)}

\[ \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 \).

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 \( \xi(x) \) between \( x_1 \) and \( x_n \).