Interpolation error with gsl akima spline

Question

I'm getting a 'gsl: interp.c:150: ERROR: interpolation error' with the following code. Some googling says that this error occurs when you try to extrapolate using the interp function but I don't see how that is happening here. Help would be greatly appreciated. Thanks.

The function randomground() just returns a random number (double).

#define NSTEPS 100   

int main()
{
           int j, q, space = 1, refine = 100;
           double xi = 0.0, tx[2*NSTEPS] = {0}, theight[2*NSTEPS] = {0};

           double terrain[(int) (2*NSTEPS*100)] = {0};
           double terrainsl[(int) (2*NSTEPS*100)] = {0};

            for (j = 0; j < 2*NSTEPS; j++)
            {
                tx[j] = (double) j*space;
                theight[j] = randomground();
            }

            gsl_interp_accel *acc = gsl_interp_accel_alloc();
            gsl_spline *spline = gsl_spline_alloc(gsl_interp_akima, 2*NSTEPS);
            gsl_spline_init(spline, tx, theight, 2*NSTEPS);

            for (q = 0; q< 2*NSTEPS*100; q++)
            {
                terrain[q] = gsl_spline_eval(spline,xi,acc);
                terrainsl[q] = gsl_spline_eval_deriv(spline,xi,acc);
                xi = xi+(double) space/refine;
            }
return 0;
}

I should add that when I run this on my Windows machine there is no error but when I run it on the Linux server in the lab, I get the interp error — NKD, Dec 09 '14 at 17:39
Float precision error may be the problem. Add a small shift to tx's first and last elements to increase the range. This won't affect your final result because the required precision in your calculation is probably much lower than double precision. — Vivian Miranda, Dec 09 '14 at 18:19
@ViniciusMiranda, sorry I don't get what you mean by 'Add a small shift to tx's first and last elements' — NKD, Dec 10 '14 at 07:43

score 1 · Answer 1 · answered Dec 10 '14 at 09:22

Fixed the problem by adding an extra element to tx and theight. I'm guessing this is what you had asked me to do, @ViniciusMiranda. The code now reads

                double tx[2*NSTEPS+1] = {0}, theight[2*NSTEPS+1] = {0};
                double terrain[(int) (2*NSTEPS*100)] = {0};
                double terrainsl[(int) (2*NSTEPS*100)] = {0};

            for (j = 0; j < 2*NSTEPS+1; j++)
            {
                tx[j] = (double) j*space;
                theight[j] = randomground();
            }

            gsl_interp_accel *acc = gsl_interp_accel_alloc();
            gsl_spline *spline = gsl_spline_alloc(gsl_interp_akima, 2*NSTEPS+1);
            gsl_spline_init(spline, tx, theight, 2*NSTEPS+1);

            for (q = 0; q< 2*NSTEPS*100; q++)
            {
                terrain[q] = gsl_spline_eval(spline,xi,acc);
                terrainsl[q] = gsl_spline_eval_deriv(spline,xi,acc);
                xi = xi+(double) space/refine;
            }

I still don't understand why this fix was required though.

The source of the problem is that you asked the interpolation to provide you the height at the boundary but float number errors ( 0.0 == 0.0 is usually false because of truncation in the 16th digit for doubles!) can give you extrapolation errors. That is why I proposed you to increase the range. The first element in array is 0 and you start asking the interpolation at zero. One fast hack to test the presence of this error is to add a small shift at the boundaries tx[0] = 0 - 1e-14 and tx[size-1] = last_element*(1+1e-14). — Vivian Miranda, Dec 10 '14 at 16:38

score 0 · Answer 2 · answered Mar 29 '15 at 17:54

0

Akima spline is a local, sub-spline interpolation. For a function f(x), if you are trying to obtain the value of x for x_i <= x <= x_i+1, Akima spline needs the pairs (x_j, f_j) for j = i-2, i-1, i, i+1, i+2, i+3.

answered Mar 29 '15 at 17:54

namu

205
2
8

score 0 · Answer 3 · edited Sep 07 '19 at 10:29

0

I found the same error in my code zigzag.sourceforge.net and resolved by commenting the line in the gsl library source, compiling it and reinstalling.

Before version 1.14 there is in interp.c the line 150

// GSL_ERROR_VAL("interpolation error", GSL_EDOM, GSL_NAN);

that I comment to solve the break !

edited Sep 07 '19 at 10:29

Fahim Al Mahmud Ashik

1,083
10
26

answered Sep 07 '19 at 00:07

Janos Erzinger

1

Interpolation error with gsl akima spline

3 Answers3