Piecewise cubic interpolation

Code to interpolate function values and corresponding derivatives.

Description

Piecewise cubic interpolation and approximated calculation of first and second derivative at the interpolation point.

• MATLAB release: MATLAB 5.2 (R10)

[John D'Errico]

This is a piecewise cubic Lagrange code, sort of.

What don't I like? Mathematically this is a poor choice. The author finds the 4 closest points to any point xi. Note that these 4 points need not even bracket the data point. So you may actually be extrapolating a given polynomial. This is not a good idea.

Next think about what happens as you move the point to be interpolated. You will find that the interpolant may actually be discontinuous. This discontinuity may occur in the interior of an interval if your data points are unequally spaced. Even in the best case, with equally spaced data, the interpolant will not be differentiable. This suggests that using this code to compute the first and even second derivatives of your curve is silly in the extreme.

Below is an example of what happens with this code. Note the discontinuity in the middle of an interval.

x = [0 .1 .2 .3 1 2 3];
y = sin(x);
xi = linspace(0,3,200);
yi = pcinterp(x,y,xi);
plot(x,y,'o',xi,yi,'r-')

The problem is, the basic algorithm is suspect. Its not as if this can be done better to eliminate the problem I've shown. A spline model or piecewise cubic Hermite is a far better choice for your interpolation problems. See interp1, spline, pchip, and the splines toolbox.

The help for pcinterp is lacking. In effect, there is no help. There is an example of use. But there are no explicit explanations of the input arguments, what size they may be, what type they are, etc.

There is no error checking here. Good code would check that the size of x and y are the same. Good code would look for replicate points in x, a serious error as this would cause a divide by zero.

The first line of the help is split over two lines. This should have been ONE line, so to enable lookfor to work properly.

The code itself has some problems. Odds are the author never thought to try his code with long vectors for both x and xi. If you do so, perhaps with 10000 points for each, you would very quickly see memory overflows when it creates 10000x10000 matrices. Why? Because the scheme used to find the closest points in x to any xi creates large arrays which are then sorted.

A few minor issues are pointed out by mlint. The author preallocates arrays for the output, but then reallocates them. In this case, there never was any need to preallocate the arrays, so it is a waste of time. Note that mlint points this fact out. Use it.

There is also a minor issue of interface I might suggest. The interpolated predictions are returned as a column vector, ALWAYS. A friendly code would return its results with the exact same shape as was the input vector xi.

What DO I like about this code? The author has gone to great pains to explain the methods used in the code. A great deal of nicely written documentation is found inside the file. I like that the author has left much white space in the code. It was very easy to read.

My rating here reflects the multiple problems in the code and the help. I hope the author does not take it badly, because the care taken in some areas of the code suggests some promise for the future in his work. But this particular code has serious flaws, too many for me to recommend it.

[V P]

John says:
This is a piecewise cubic Lagrange code, sort of.

With exception of xi in first or last interval in original x-data, this is true only if from 4 selected x-points 2 are to the left, and 2 to the right of xi.

Since for nonuniformly spaced data 4 x points closest to xi not always "doublebracket" xi, the program should be first corrected before it can be considered as a standard interpolation tool.

If we interpolate cubically (or in SPLIN mode, i.e. close to cubic spline) remaining within the class of uniformly spaced data, then INTERPOL will do the job faster.

[John D'Errico]

The author has improved this code, but it can still be better.

The help is now useful. There is error checking. The author has fixed the memory problems for large data sets, using the repmat code for small one, but a loop for large ones.

A fix has been introduced to look for replicate data, but it uses unique. So only the first point of a set of replicates is used. A better choice is to average the replicates, although this takes more work.

As Vassilli points out, this is not a traditional piecewise Lagrange code, but since I said it was only sort of such a code, I'll leave my statement as it is. When the data points are equally spaced in x, then this will indeed be a traditional piecewise cubic Lagrange. In fact, one can argue that this truly is a cubic Lagrange over a subset of the points. Its just that for unequally spaced data, that subset changes rather arbitrarily as you move along.

The example is extensive, as is the explanation of how cubic Lagrange works.

What remains that I'm not enthralled with is the use of the 4 nearest points for interpolation. This causes the result to be discontinuous. The basic algorithm is still suspect. At least the author makes a comment to that effect, so a user will know that fact. Even then, this is no better than a cubic Lagrange. Given these problems, I can't offer better than a 3 rating here.

A person who needs a serious tool for interpolation should still pass this over, in favor of the many better tools available for the job. Look for spline, pchip, interp1, or the splines toolbox. On the file exchange, you might even look at interpol.

[li wenlong]

4