Why is the pchip/cubic interpolation error the same as linear interpolation?

I was trying to estimate the error bounds for a numerical simulation, and found that pchip/cubic interpolation were not behaving as I expected. The error on an interpolating polynomial should go like ~(dx)^(n+1) where dx is the grid spacing and n is the order of the polynomial. So, for linear interpolation the error should go like dx^2, and for a cubic spline it should go like dx^4. I find that both of these thing are true using Matlab's linear and spline interpolation function. But, for some reason the 'cubic/pchip' interpolation error also scales like dx^2---as if its linear. Is this expected?!?

Here's the code I'm using, which actually does the analysis by changing the wavenumber of a sine function on a fixed grid (rather than fixing the function and changing the grid).

L = 128;
N = 128;
dx = L/N;
x = dx*(0:N-1)';
x_fine = linspace(0,L-dx,100*N)';

mgrid = 0:.1:7;
error = zeros(length(mgrid),1);
error_theory = zeros(length(mgrid),1);
method = 'spline'; n=3;
% method = 'pchip'; n=1;
% method = 'linear'; n=1;
for i=1:length(mgrid)
  m = mgrid(i);
  k = 2*pi*(2^m)/L;
  f_interp = interp1(x,sin(k*x),x_fine,method);
  f_exact = sin(k*x_fine);
  error(i) = max(abs(f_interp-f_exact));
  error_theory(i) = ((k*dx)^(n+1))/(4*(n+1));
end

figure
plot(mgrid,error), ylog, hold on
plot(mgrid,error_theory)
legend('actual', 'theory')
title(sprintf('Error with %s interpolation',method))

Answer 1

[John D'Errico]

Because you need to understand how pchip works. As well, you chose a poor metric to evaluate the error, at least for this specific problem.

You used the maximum error to measure the prediction error, on an oscillatory function. Why is this a problem when using pchip as an interpolant? Pchip was designed with very specific goals in mind. One big factor was to avoid introducing ringing behavior in an interpolant. Think of it as Gibbs phenomena, if you recognize the term. Splines tend to be terribly bad in that respect.

x = 0:10;
y = double(x > 5);
xfine = 0:.1:10;
yfinespl = interp1(x,y,xfine,'spline');
yfinepchip = interp1(x,y,xfine,'pchip');
plot(x,y,'ro',xfine,[yfinespl;yfinepchip]','-')

As you can see, the blue curve is a cubic spline, but green is pchip.

The function we are trying to interpolate is a simple step function. A spline interpolates that like crap.

So, why is that a problem for you? PCHIP is designed to NEVER add an extremum in your function that was not there already. PCHIP is a great tool for functions with sharp transitions, where monotonicity is of primal importance, yet still be able to interpolate moderately smoothly.

But what happens when you use pchip to interpolate an oscillatory function? Remember, this is the antithetical problem to the design of PCHIP.

x = linspace(0,pi,6);
y = sin(x);
xfine = linspace(0,pi,1000);
yfinepchip = interp1(x,y,xfine,'pchip');
yfinespl = interp1(x,y,xfine,'spline');
plot(x,y,'ro',xfine,[yfinespl;yfinepchip]','-')

Hmm. As you can see, the spline interpolant fit very well. But PCHIP "saw" the flat spot on top, and refused to introduce an extremum BETWEEN the points at the peak.

So, for oscillatory functions, IF you use a maximum error measure, AND you use PCHIP, it will indeed behave much like a linear interpolant, by a maximum error metric. Remember that the maximum error for a pchip interpolant will probably always be across a peak on an oscillatory function, since pchip will act much like a linear interpolant there.

Had you instead tried to interpolate a function like exp(x), you would find that PCHIP has better performance than linear, although probably not as good as a spline, because PCHIP is not a C2 (twice continuously differentiable) interpolant, as is a cubic spline.

The point is, rather than just throwing numbers at a function and wondering why it does as it does, you need to understand the tools you are using. When you see something strange, the first thing you need to do is plot things. PLOT EVERYTHING. In this case, plotting the various interpolants might have helped you to understand what is happening.

When used to solve a problem from the class it was designed to solve, PCHIP is a great thing. But not every tool is perfect for every problem. If that was true, we would only be able to buy one kind of screwdriver, and we would only ever need one such tool.