Ch4.1 Part 2: Polynomial Interpolation

• A common application is to approximate functions.
• Recall Taylor polynomials from previous section of notes.
• Taylor polynomials use derivative information to “shape-fit” an approximating polynomial.
• Focused on the point \( (a,f(a)) \) on graph; see next slide.

\[ \small{ \begin{aligned} 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 \\ & = \sum_{k=0}^{n} \frac{f^{(k)}(a)}{k!} (x-a)^k \\ & = P_n(x) \end{aligned} } \]

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

• Another application is polynomial interpolation to data.
• In this case we are interested in fitting a polynomial through multiple data points, rather than just one.

• For \( n \) data points, a degree \( n-1 \) polynomial \( p(x) \) can be found that passes through all of the data points.
  
• If we view the data points as being generated by a function \( f(x) \) on \( [x_1,x_n] \), then \( p(x) \) approximates \( f \) on this interval, agreeing with \( f \) at the nodes \( x_k \).
• This polynomial can be used to estimate \( f(x) \) at values of \( x \) in \( [x_1,x_n] \), but not outside of \( [x_1,x_n] \) (extrapolation).

What do you call a fake noodle?

An impasta.

Why can't you hear a Pterodactyl go to the bathroom?

Because the pee is silent.

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

• When computing with an interpolating polynomial, nested polynomial evaluation (Horner's method) is used.

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

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

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:

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

The polyinterp program computes the coefficients of the interpolating polynomial using the Vandermonde matrix.

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

• Now that we have the coefficients, we use nested evaluation to compute values of the polynomial.
  
• For example, as a check, we can plug in the \( x \)-data points to see if we obtain the \( y \)-data points.

horner(xdata,p)

[1] 0.60 0.67 1.01 1.35 1.47 1.25

ydata

[1] 0.60 0.67 1.01 1.35 1.47 1.25

The output of the polyinterp program is the vector of coefficients on the interpolating polynomial \( p(x) \).

(p <- polyinterp(xdata, ydata))

[1]  0.6000000  0.5910395 -2.5174715  5.4280961 -3.5892186  0.7303129

p[1]; p[2]

[1] 0.6

[1] 0.5910395

The naive method of polynomial evaluation uses the coefficients on powers of \( x \) to directly evaluate polynomial:

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

naivepoly(xdata)

[1] 0.60 0.67 1.01 1.35 1.47 1.25

ydata

[1] 0.60 0.67 1.01 1.35 1.47 1.25

The generate a plot of the data with interpolating polynomial, the program code on next slide does the following:

• Uses polyinterp to compute coefficients.
• Generates \( N = 100 \) new nodes on \( [x_1,x_n] \), with step size

\[ h = \frac{x_n-x_1}{N} \]

• Evaluates the interpolating polynomial at these 100 nodes using Horner's method (nested evaluation).
• Plots the original data together with the nested polynomial evaluation values on the 100 nodes.

HornerPlot <- function(x,y) {
p <- polyinterp(x, y)
n = length(x)
h <- (x[n] - x[1])/500
t <- seq(x[1],x[n],h)
f <- horner(t,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)))
points(x,y,type = "p",col="blue",lwd=3)
legend("topleft",
  legend = c("Data", "Polynomial"),
  col=c("blue","red"),  
  lty=c(1,1)  ) }      

\[ \small{ \begin{aligned} 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 \\ R(x) & = \frac{(x-a)^{n+1}}{(n+1)!}f^{(n+1)}(c(x)), \end{aligned} } \]

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

• Mathematically, \( p(x) \) fits the \( n \) data points perfectly.
• \( p(x) \) is typically used to evaluate \( x \)-values between \( [x_1,x_n] \), but not outside of this interval (extrapolation).
  
• In this way \( p(x) \) represents the function \( f(x) \) that generates the data, where \( f(x) \) typically doesn't have a formula.

Error formula is shown below, where \( c(x) \) between \( x_1 \) and \( x_n \).

• Note that the error formula is zero at the nodes \( x_k \).
• As with Taylor polynomials, calculus can be used to determine an upper bound on absolute error.

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

Vandermonde approach is convenient way of theoretically describing how to find interpolating polynomial, but computing with it is ill-conditioned; see weblink.

\[ p(x) = c_0 + c_1 x + \cdots + c_{N-1} x^{N-1} \]

• Small deviation in entries causes large deviation in solution.

• These systems are ill-conditioned becayse second row is very close to being a multiple of first row.
  
• This means that the two lines are nearly parallel.
  
• If two lines are truly parallel, then they will never intersect.
• If two lines are nearly parallel, then they intersect but it may take a while relative to one of the axes.

• If two lines are nearly parallel, then they intersect but it may take a while relative to one of the axes.
• Thus a slight change in a nearly parallel system sends the intersection point a disproportionately larger distance away.
• Row operations in this case also creates inaccuracies due to cancellation errors and loss of significant digits.
  

• Vandermonde matrices should be avoided numerically, since they are ill-conditioned.
• We will instead use Lagrange Interpolating Polynomials.
• These can be constructed with pencil and paper; click here for YouTube video.

Suppose we are given \( n \) data points

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

Define the Lagrange helping function \( L_k \) as follows:

\[ \small{ L_k = \frac{(x-x_0)\cdots (x-x_{k-1})(x-x_{k+1})\cdots (x-x_n)}{(x_k-x_0)\cdots (x_k-x_{k-1})(x_k-x_{k+1})\cdots (x_k-x_n)} } \]

Suppose we are given \( n \) data points

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

With the Lagrange helping functions

\[ \small{ L_k = \frac{(x-x_0)\cdots (x-x_{k-1})(x-x_{k+1})\cdots (x-x_n)}{(x_k-x_0)\cdots (x_k-x_{k-1})(x_k-x_{k+1})\cdots (x_k-x_n)}, } \]

the Lagrange interpolating polynomial is defined as

\[ \small{ P(x) = y_1 L_1(x) + y_2 L_2(x) + \ldots + y_n L_n(x) } \]

Suppose we are given the \( n = 3 \) data points below:

\[ \small{ (1,5), (2,6), (3,7) } \]

Lagrange helping functions:

\[ \small{ L_1 = \frac{(x-2)(x-3)}{(1-2)(1-3)}, \,\, L_2 = \frac{(x-1)(x-3)}{(2-1)(2-3)}, \,\, L_3 = \frac{(x-1)(x-2)}{(3-1)(3-2)} } \]

Lagrange interpolating polynomial:

\[ \small{ \begin{aligned} P(x) & = 5 \frac{(x-2)(x-3)}{(1-2)(1-3)} + 6 \frac{(x-1)(x-3)}{(2-1)(2-3)} + 7 \frac{(x-1)(x-2)}{(3-1)(3-2)} \\ & = x + 4 \,\, \mathrm{(after \, some \, algebra)} \end{aligned} } \]

Suppose we are given the \( n = 3 \) data points below:

\[ \small{ (2,16), (5,12), (7,25) } \]

Lagrange helping functions:

\( L_1= \)

\( L_2 = \)

\( L_3 = \)

Lagrange interpolating polynomial (do not simplify):

\( P(x)= \)

We are interested in estimating missing value at \( x = 0.9 \).

\( (a) \) Use \( x_3, x_4 \) to construct a degree one interpolation.

We are interested in estimating missing value at \( x = 0.9 \).

\( (b) \) Use \( x_2, x_3, x_4 \) to construct a degree two interpolation.

We are interested in estimating missing value at \( x = 0.9 \).

\( (c) \) Use \( x_1, x_2, x_3, x_4 \) to construct degree three interpolation.

library(pracma)
LagrangePlot <- function(x,y) {

  n = length(x)
N <- 500 

t <- linspace(x[1], x[n], N)
p <- lagrangeInterp(x, y, t)

bottom = min(0,min(p)); top = max(0,max(p))

plot(t,p, 
  main = "Lagrange Interpolating Polynomial",
  xlab="x",ylab="y",type="l",lwd=2,col="red", 
  xlim = c(x[1],x[n]), ylim = c(bottom, top) )
points(x,y,type = "p",col="blue",lwd=3)}