Stability of polynomial interpolation
Contents
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.
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:
While the dependence on \(f\) is messy here, the error indicator \(\Phi(x)\) can be studied as a function of the nodes only.
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")