Hi,
while i was interpolation some test-data i got the following results.
3.736 0.0298362 -2.66867
3.737 0.0271675 -2.66867
3.738 0.0244988 -2.66868
3.739 0.0218302 -2.66865
3.74 0.0191615 -2.66875
3.741 0.0164928 -2.6684
3.742 0.0138242 -2.66969
3.743 0.0111555 -2.66489
3.744 0.0084868 -2.68281
3.745 0.00581813 -2.61591
3.746 0.00314946 -2.8656
3.747 0.000480783 -1.93374
3.748 -0.000604339 -0.660811
--------------------------------------
3.736 0.0298362 -2.66867
3.737 0.0271675 -2.66867
3.738 0.0244988 -2.66867
3.739 0.0218302 -2.66867
3.74 0.0191615 -2.66867
3.741 0.0164928 -2.66867
3.742 0.0138242 -2.66867
3.743 0.0111555 -2.66867
3.744 0.0084868 -2.66867
3.745 0.00581813 -2.66867
3.746 0.00314946 -2.66867
3.747 0.000480783 -2.66867
3.748 -0.000604339 -1.08512
The data where generated with the following C++ program. I used the
usual gsl-functions with spline interpolation and got the first part of
the data. For the second part of the data i used internal gsl-functions
from the file "cspline.c". The second part of the data contains the
expected result (the third column should give a constant value) while
the first part contains almost identical results. For 3.745 we get 3.745
0.00581813 -2.61591 which makes no sense. I included these internal
functions in the program, so that everyone can reproduce the results. As
far as i know these internal function which produced the expected
results are used by those produces the wrong ones. What is going wrong?
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <stdio.h>
namespace gslcode {
typedef struct {
double *c;
double *g;
double *diag;
double *offdiag;
} cspline_state_t;
/* function for common coefficient determination
*/
inline void coeff_calc(const double c_array[], double dy, double dx,
size_t index, double *b, double *c, double *d) {
const double c_i = c_array[index];
const double c_ip1 = c_array[index + 1];
*b = (dy / dx) - dx * (c_ip1 + 2.0 * c_i) / 3.0;
*c = c_i;
*d = (c_ip1 - c_i) / (3.0 * dx);
}
int cspline_eval(const void *vstate, const double x_array[],
const double y_array[], size_t size, double x,
gsl_interp_accel *a, double *y) {
const cspline_state_t *state = (const cspline_state_t *)vstate;
double x_lo, x_hi;
double dx;
size_t index;
if (a != 0) {
index = gsl_interp_accel_find(a, x_array, size, x);
} else {
index = gsl_interp_bsearch(x_array, x, 0, size - 1);
}
/* evaluate */
x_hi = x_array[index + 1];
x_lo = x_array[index];
dx = x_hi - x_lo;
if (dx > 0.0) {
const double y_lo = y_array[index];
const double y_hi = y_array[index + 1];
const double dy = y_hi - y_lo;
double delx = x - x_lo;
double b_i, c_i, d_i;
coeff_calc(state->c, dy, dx, index, &b_i, &c_i, &d_i);
*y = y_lo + delx * (b_i + delx * (c_i + delx * d_i));
return GSL_SUCCESS;
} else {
*y = 0.0;
return GSL_EINVAL;
}
}
static int cspline_eval_deriv(const void *vstate, const double x_array[],
const double y_array[], size_t size,
double x,
gsl_interp_accel *a, double *dydx) {
const cspline_state_t *state = (const cspline_state_t *)vstate;
double x_lo, x_hi;
double dx;
size_t index;
if (a != 0) {
index = gsl_interp_accel_find(a, x_array, size, x);
} else {
index = gsl_interp_bsearch(x_array, x, 0, size - 1);
}
/* evaluate */
x_hi = x_array[index + 1];
x_lo = x_array[index];
dx = x_hi - x_lo;
if (dx > 0.0) {
const double y_lo = y_array[index];
const double y_hi = y_array[index + 1];
const double dy = y_hi - y_lo;
double delx = x - x_lo;
double b_i, c_i, d_i;
coeff_calc(state->c, dy, dx, index, &b_i, &c_i, &d_i);
*dydx = b_i + delx * (2.0 * c_i + 3.0 * d_i * delx);
return GSL_SUCCESS;
} else {
*dydx = 0.0;
return GSL_EINVAL;
}
}
} // namespace gslcode
int main() {
const size_t N = 13;
const double t[] = {3.736, 3.737, 3.738, 3.739, 3.740, 3.741, 3.742,
3.743, 3.744, 3.745, 3.746, 3.747, 3.748};
const double v[] = {0.02983619, 0.027167517, 0.024498843, 0.02183017,
0.019161497, 0.016492823, 0.01382415, 0.011155476,
0.008486803, 0.00581813, 0.003149456, 0.000480783,
-0.000604339};
gsl_interp_accel *acc = gsl_interp_accel_alloc();
gsl_interp *spline = gsl_interp_alloc(gsl_interp_cspline, N);
gsl_interp_init(spline, t, v, N);
for (double xi = 3.736; xi < 3.748; xi += 0.001) {
double val1 = gsl_interp_eval(spline, t, v, xi, acc);
double val2 = gsl_interp_eval_deriv(spline, t, v, xi, acc);
printf("%g %g %g\n", xi, val1, val2);
}
printf("--------------------------------------\n");
for (double xi = 3.736; xi < 3.748; xi += 0.001) {
double val1, val2;
gslcode::cspline_eval(spline, t, v, N, xi, acc, &val1);
gslcode::cspline_eval_deriv(spline, t, v, N, xi, acc, &val2);
printf("%g %g %g\n", xi, val1, val2);
}
gsl_interp_free(spline);
gsl_interp_accel_free(acc);
}
Elias
