Hi,

I'm trying to benchmark cubic spline interpolation function interp against MATLAB's spline function using the following pieces of code but there is a difference that I'm not able to explain.

/*helloworld*/
#include "nr3matlab.h"
#include "interp_1d.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {

VecDoub xx(prhs[0]);
VecDoub yy(prhs[1]);
VecDoub x(prhs[2]);
int n = x.size();
VecDoub y(n,plhs[0]);
int i;

Spline_interp * myfunc= new Spline_interp(xx,yy);

for (i=0;i<n;i++){
y[i] = myfunc->interp(x[i]);
}
}

%MATLAB code for benchmark
x = 0:0.1:10;
y = sin(x);
xx = 0:.025:10;
yspline = spline(x,y,xx);
yyhelloworld=helloworld(x,y,xx);
stem((yyspline-yyhelloworld))

interp as called this way retuns the natural cubic spline, right? What about spline.m?

Appreciate your comments, thanks

...there is a difference

What is the nature of the difference? Large errors at some points? Large errors at all points? Define "large."
that I'm not able to explain.What are your expectations of the differences in results obtained from different programs? Especially since one is closed-source.

interp as called this way retu[r]ns the natural cubic spline, right?That's what the NR book indicates.
What about spline.m?I don't have Matlab, but that's what the Matlab documentation says. I do have GNU Octave, and running your program through GNU Octave gives the same values (within a precision five or six significant decimal digits) as I got from the approximation with a C++ program using Spline_interp::interp.

The results aren't exactly the same, of course. The only "blips" that I notice are in a few values near the right end point. (Maximum relative error of GNU Octave spline is something like 7.9e-4 and with GNU Octave it's on the order of 4.1e-5.) Perhaps some details differ in the way the functions handle end point evaluation, but without knowing your requirements, I can't say whether this should be considered a "problem."

Bottom line: Tell us exactly what part of the results didn't agree with your expectations.


Regards,

Dave

Hi,

Thanks for your prompt response to my question. I get higher relative errors near end points as the first image reveals. I'm not knowledgeable much about details of cubic spline interpolation but my NAIVE expectation is to see errors in the order of double precision provided that both functions implement the same type of interpolation. Digging into the details of MATLAB spline function I found the following description.

Here is what MATLAB help for the spline functions says:

"Ordinarily, the not-a-knot end conditions are used. However, if Y contains two more values than X has entries, then the first and last value in Y are used as the endslopes for the cubic spline. "

So, appending and prepending zeros to the y vector fed to the function, I hoped to force MATLAB to implement a natural cubic spline interpolation but the relative error increased even more as seen in the second image.

I would deeply appreciate any further comments.

Here is what MATLAB help for the spline functions says:
the first and last value in Y are used as the endslopes for the cubic spline
.
.
.
So, appending and prepending zeros to the y vector fed to the function, ...

You are setting derivatives of the spline function to zero at the end points.

To do that with Spline_interp, you give it some extra arguments:

//Set values of derivative equal to zero at the end points
Spline_interp * myfunc= new Spline_interp(xx, yy, 0, 0);


With this change and by putting zero as values of the first and last elements of the augmented y-vector for the Octave spline call, the maximum absolute value of the differences between the GNU Octave spline approximation and the NR spline interpolation at your 401 points of evaluation is on the order of 1.1e-16. The plots are visually similar. That is, they "look the same" to me. (See Attachments.)



Regards,

Dave

Hi Dave,

They now "look same to me" too :). That's indeed what I was looking for. Now that we could make both functions implement natural cubic spline, I'd like to ask the following.

1-) Can you point me to the place where was I supposed to find this info. about Spline_interp instead of asking you (to implement natural cubic spline)?

2-) What type of cubic spline interp. does Spline_interp perform in the absence of these prepended and appended zeros?

Thanks,
Onur

.
Can you point me...Source code and some exposition in the text. It says that it is the natural spline.


I did a little more digging and I found that the GSL (Gnu Scientific Library) spline functions give answers very close to those that I got from the Numerical Recipes spline_interp.

The GNU Octave help file for the spline function says that it uses the "not-a-knot" end conditions and gives answers that look more like your Matlab output (differences between the Octave values and the NR values near the right-hand end point). Somehow I thought the "not-a-knot" conditions were equivalent to the natural spline. (That shows how much I know, right?)

Matlab is, of course, closed-source. I seem to recall that they have another spline function (cscvn) that they used to say was a "natural" spline, but it wasn't exactly the same as the usual definition of natural spline (I think).

Maybe you can follow through on this.

Note that the maximum absolute error of the GNU Octave interpolated values near the right end are somewhat less than those of the NR or the GSL functions (something on the order of 1e-6 versus 1e-4). Maybe they are doing something "better" than the classical natural spline?

Summary: NR functions have fairly easy to read source code and some narrative.

Source is available for GNU Octave, and, in fact, the spline function is in an m-file, but I have never had the patience or stamina to try to work through it.

Same for the C code in the source code for the GNU Scientific Library.


Regards,

Dave

Hi Dave,

Since your last post, I've discovered a few things about the end point conditions of MATLAB's and NR3's spline functions which is likely to resolve the exisiting ambiguity at all.

matlab spline: by default implements cubic runout spline interpolation which imposes the so-called "not a-knot" end conditions. That is, end points' second derivatives are extrapolated from the second derivatives of the corresponding two neighbouring points. (please see http://dmpeli.mcmaster.ca/Matlab/Math1J03/ComputerLabs/Lab8.htm) If y input to the spline function is padded with zeros, howeer, as we did in our second iteration, then it forces the first derivatives of the end points to zero which is obviously not the "natural cubic spline" I was trying to end up with which requires second derivatives to be set to zero.

nr3 spline: by default implements "natural cubic spline" (see page 123 of the book), and when forced as shown below, it sets endpoints' first derivatives to zero just like MATLAB does and that's how we could make them output the same thing within the numerical precision of double.

Spline_interp * myfunc= new Spline_interp(xx, yy, 0, 0);

In conclusion default end point behaviours of these functions are different and one can make them equal by forcing conditions on the first derivatives of the end points.

Hi Dave,

Since your last post, I've discovered a few things...

Hi!


Thanks for the update. The reason that I hang around here is that I like learning about things like this.


Regards,

Dave

Hello,
i have used your suggestions to implent a C code for the cubic spline interpolartion with the nr functions but the obtained results are different from the matlab code.

This is the c code:

int n=len;
VecDoub xx(n), yy(n);
for(int i=0;i<len;i++)
{
xx[i]=(Doub)freq_n[i]; // vectors with the values to //be be interpolated
yy[i]=(Doub)mag_log[i];
}

Spline_interp myfunc(xx,yy,0,0);

for (int i=0;i<Frame;i++){
x = ((double)i)/((double)Frame);
y = myfunc.interp(x);
H[i]=y;
}

Then in matlab I do:

yy=spline(freq_n,mag_log,x);

and yy should be equal to H... but the differences are high at the initial points.