hi all, wanted to turn the cubic spline function from 3.3 into an inline version that uses 0-pad (mixing 1 and 0-pad code is in my opinion terrible style) in C. so i determined that i had to subtract 1 from every index to the vector u. i'm not getting the result i'd hoped tho, maybe i missed something.



#declare SIZE 32

int main() {

int i; double sig, p; double x[SIZE], y[SIZE], u[SIZE - 1], y2[SIZE - 1];

// fill up with interesting data
for (i = 0; i < SIZE; i++) {
x[i] = i;
y[i] = pow(i, 2) + i;
}

// left boundary conditions: for these, 0 is fine for the moment
u[0] = 0.0; y2[0] = 0.0;

// decomposition loop
for (i = 2; i < SIZE - 1; i++) {
sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]);
p = sig * y2[i - 1] + 2.0; y2[i] = (sig - 1.0) / p;
u[i-1] = (y[i+1]-y[i]) / (x[i+1]-x[i]) - (y[i]-y[i-1]) / (x[i]-x[i-1]);
u[i-1] = (6.0 * u[i-1] / (x[i+1] - x[i-1]) - sig * u[i-2]) / p;
}

// right boundary conditions
u[SIZE - 2] = 0.0; y2[SIZE - 2] = 0.0;

// backsubstitution loop
for (i = SIZE - 1; i >= 1; i--) { // stays saem
y2[i] = y2[i] * y2[i + 1] + u[i - 1];
}
}


anybody have any ideas? y2 does not look anything like y. are there any other tricks to turn 1-pad code in the C edition to 0-pad C?

...y2 does not look anything like y....

1. The NR spline() function does not actually interpolate anything. It gives coefficients, corresponding to the second derivatives, that can be used by splint() to perform the actual interpolation. This is explained in section 3.3.

2. Instead of my trying to figure what your program is doing and trying to see whether it is wrong, I offer the following as an example that uses the standard Numerical Recipes Version 2 functions. Note that for zero-based 1-D arrays in the main() function, it's just a matter of sending adjusted pointer values to the one-based NR functions.

Here's a test program using your function. I compiled with GNU gcc.


/*
* Example of cubic spline interpolation
* using spline() and splint()
*
* davekw7x
*/

#include <stdio.h>
#define NRANSI
#include "nr.h"

#define NPOINTS 32

/* Function to be interpolated */
float f(float x)
{
return x*x+x;
}

/* Second derivative of f(x) */
float fpp(float x)
{
return 2.0;
}

int main(void)
{
int i;
float x, y, yp1, ypn;

float xa[NPOINTS];
float ya[NPOINTS];
float y2[NPOINTS];

for (i = 0; i < NPOINTS; i++) {
xa[i] = i;
ya[i] = f(xa[i]);
}

yp1 = fpp(xa[0]); /* Second derivative at first point */
ypn = fpp(xa[NPOINTS-1]); /* Second derivative at last point */

/*
* Call spline to get second derivatives.
*
* Since the Numerical Recipes function assumes
* that arrays are 1-based, adjust the pointer
* values accordingly
*/
spline(xa - 1, ya - 1, NPOINTS, yp1, ypn, y2 - 1);

/*
* Call splint for interpolations. The same
* considerations apply (one-based arrays).
*/
printf("%9s %13s %17s\n", "x", "f(x)", "from splint");
for (i = 0; i < NPOINTS; i++) {
x = i;
splint(xa - 1, ya - 1, y2 - 1, NPOINTS, x, &y);
printf("%13.5e %13.5e %13.5e\n", x, f(x), y);
}
printf("\n");


printf("%9s %13s %17s\n", "x", "f(x)", "from splint");
for (i = 0; i < 10; i++) {
x = 20.0 + i/10.0;
splint(xa - 1, ya - 1, y2 - 1, NPOINTS, x, &y);
printf("%13.5e %13.5e %13.5e\n", x, f(x), y);
}
return 0;
}


Output:

x f(x) from splint
0.00000e+00 0.00000e+00 0.00000e+00
1.00000e+00 2.00000e+00 2.00000e+00
2.00000e+00 6.00000e+00 6.00000e+00
3.00000e+00 1.20000e+01 1.20000e+01
4.00000e+00 2.00000e+01 2.00000e+01
5.00000e+00 3.00000e+01 3.00000e+01
6.00000e+00 4.20000e+01 4.20000e+01
7.00000e+00 5.60000e+01 5.60000e+01
8.00000e+00 7.20000e+01 7.20000e+01
9.00000e+00 9.00000e+01 9.00000e+01
1.00000e+01 1.10000e+02 1.10000e+02
1.10000e+01 1.32000e+02 1.32000e+02
1.20000e+01 1.56000e+02 1.56000e+02
1.30000e+01 1.82000e+02 1.82000e+02
1.40000e+01 2.10000e+02 2.10000e+02
1.50000e+01 2.40000e+02 2.40000e+02
1.60000e+01 2.72000e+02 2.72000e+02
1.70000e+01 3.06000e+02 3.06000e+02
1.80000e+01 3.42000e+02 3.42000e+02
1.90000e+01 3.80000e+02 3.80000e+02
2.00000e+01 4.20000e+02 4.20000e+02
2.10000e+01 4.62000e+02 4.62000e+02
2.20000e+01 5.06000e+02 5.06000e+02
2.30000e+01 5.52000e+02 5.52000e+02
2.40000e+01 6.00000e+02 6.00000e+02
2.50000e+01 6.50000e+02 6.50000e+02
2.60000e+01 7.02000e+02 7.02000e+02
2.70000e+01 7.56000e+02 7.56000e+02
2.80000e+01 8.12000e+02 8.12000e+02
2.90000e+01 8.70000e+02 8.70000e+02
3.00000e+01 9.30000e+02 9.30000e+02
3.10000e+01 9.92000e+02 9.92000e+02

x f(x) from splint
2.00000e+01 4.20000e+02 4.20000e+02
2.01000e+01 4.24110e+02 4.24110e+02
2.02000e+01 4.28240e+02 4.28240e+02
2.03000e+01 4.32390e+02 4.32390e+02
2.04000e+01 4.36560e+02 4.36560e+02
2.05000e+01 4.40750e+02 4.40750e+02
2.06000e+01 4.44960e+02 4.44960e+02
2.07000e+01 4.49190e+02 4.49190e+02
2.08000e+01 4.53440e+02 4.53440e+02
2.09000e+01 4.57710e+02 4.57710e+02


The first table shows values for the given data points. They must agree.

The second table shows some interior points in the interval between 20 and 21. Interpolated values agree, as expected for this function.

My point: Once you have a reference design with a given function, then try your code with the same function. If the results are different, then look at the changes that you made from the function(s) as given in the book. For starters, you can try printing out the values from NR spline() and compare with your code (so far). Then implement the actual interpolation Ã* la splint(). Stuff like that.

Regards,

Dave