gsl akima样条的插值错误答案

【问题标题】：Interpolation error with gsl akima splinegsl akima样条的插值错误
【发布时间】：2014-12-09 17:20:54
【问题描述】：

我收到一个带有以下代码的“gsl: interp.c:150: ERROR: interpolation error”。一些谷歌搜索说，当您尝试使用 interp 函数进行推断时会发生此错误，但我看不出这里是如何发生的。帮助将不胜感激。谢谢。

函数randomground() 只是返回一个随机数（双精度）。

#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;
}

【问题讨论】：

我应该补充一点，当我在我的 Windows 机器上运行它时没有错误，但是当我在实验室的 Linux 服务器上运行它时，我得到了 interp 错误
浮点精度错误可能是问题所在。向 tx 的第一个和最后一个元素添加一个小偏移以增加范围。这不会影响您的最终结果，因为您的计算所需的精度可能远低于双精度。
@ViniciusMiranda，抱歉，我不明白您所说的“在 tx 的第一个和最后一个元素中添加一个小偏移”是什么意思

标签： c interpolation gsl

【解决方案1】：

通过向 tx 和 theight 添加一个额外的元素来解决问题。我猜这是你要求我做的，@ViniciusMiranda。代码现在读取

                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;
            }

我仍然不明白为什么需要此修复程序。

【讨论】：

问题的根源是您要求插值为您提供边界处的高度但浮点数错误（ 0.0 == 0.0 通常是错误的，因为双打的第 16 位截断！）可以给你外推错误。这就是为什么我建议你增加范围。数组中的第一个元素是 0，您开始要求插值为零。测试此错误是否存在的一种快速技巧是在边界 tx[0] = 0 - 1e-14 和 tx[size-1] = last_element*(1+1e-14) 处添加一个小偏移。

【解决方案2】：

Akima 样条是一种局部的子样条插值。对于函数f(x)，如果您尝试获取x_i <= x <= x_i+1 的x 的值，Akima 样条需要(x_j, f_j) 对j = i-2, i-1, i, i+1, i+2, i+3.

【讨论】：

【解决方案3】：

我在我的代码zigzag.sourceforge.net 中发现了同样的错误，并通过注释gsl library source 中的行，编译并重新安装来解决。

在 1.14 版本之前，interp.c 中有第 150 行

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

我评论是为了解决这个问题！

【讨论】：