9.3. Stability of polynomial interpolation

With barycentric interpolation available in the form of Function 9.2.3, we can explore polynomial interpolation using a numerically stable algorithm. Any remaining sensitivity to error is due to the conditioning of the interpolation process itself.

Demo 9.3.1

We choose a function over the interval \([0,1]\).

f = x -> sin(exp(2*x));

Here is a graph of \(f\) and its polynomial interpolant using seven equally spaced nodes.

plot(f,0,1,label="function",legend=:bottomleft)
t = (0:6)/6
y = f.(t)
scatter!(t,y,label="nodes")

p = FNC.polyinterp(t,y)
plot!(p,0,1,label="interpolant",title="Equispaced interpolant, n=6")

This looks pretty good. We want to track the behavior of the error as \(n\) increases. We will estimate the error in the continuous interpolant by sampling it at a large number of points and taking the max-norm.

n = 5:5:60;  err = zeros(size(n))
x = range(0,1,length=2001)      # for measuring error
for (i,n) in enumerate(n) 
  t = (0:n)/n                   # equally spaced nodes
  y = f.(t)                     # interpolation data
  p = FNC.polyinterp(t,y)
  err[i] = norm( (@. f(x)-p(x)), Inf )
end
plot(n,err,m=:o,title="Interpolation error for equispaced nodes",
    xaxis=(L"n"),yaxis=(:log10,"max error"),)

The error initially decreases as one would expect but then begins to grow. Both phases occur at rates that are exponential in \(n\), i.e., \(O(K^n\)) for a constant \(K\), appearing linear on a semi-log plot.

Runge phenomenon

The disappointing loss of convergence in Demo 9.3.1 is a sign of ill conditioning due to the use of equally spaced nodes. We will examine this effect using the error formula (9.1.5) as a guide:

\[ f(x) - p(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!} \Phi(x), \qquad \Phi(x) = \prod_{i=0}^n (x-t_i). \]

While the dependence on \(f\) is messy here, the error indicator \(\Phi(x)\) can be studied as a function of the nodes only.

Demo 9.3.2

We plot \(|\Phi(x)|\) over the interval \([-1,1]\) with equispaced nodes for different values of \(n\).

plot(xaxis=(L"x"),yaxis=(:log10,L"|\Phi(x)|",[1e-25,1]),legend=:bottomleft)

x = range(-1,1,length=2001)
for n in 10:10:50
    t = range(-1,1,length=n+1)
    Φ = [ prod(xₖ.-t) for xₖ in x ]
    scatter!(x,abs.(Φ),m=(1,stroke(0)),label="n=$n")
end
title!("Error indicator for equispaced nodes")