gsl Interpolation Question

[Elias Posoukidis]

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