c++ - gsl_interp_linear 的插值错误

Question

我需要对修改后的贝塞尔函数进行插值，为此我使用了函数 gsl_spline_eval()。我得到一个 'gsl: interp.c:145: ERROR: interpolation error' 的代码如下:

 #define MAXPOINTS 10000 
 int main()
{
    double nu;
    gsl_function F;
    F.function = &sync;
    F.params =0;
    size_t size;    
    double table[2][MAXPOINTS];
    size = read_table(table[0],table[1]);   
    gsl_interp_accel *acc      =      gsl_interp_accel_alloc ();
    gsl_spline *spline         =      gsl_spline_alloc   (gsl_interp_linear, 5401); 
    gsl_spline_init(spline,table[0],table[1],size);
    cout <<" nu:   "<< nu<< "  GSL_FN_EVAL(&F,nu):   "<<nu*GSL_FN_EVAL(&F,nu)<<endl;

 return 0;
}

其中 sync() 使用 gsl 函数 gsl_integration_qag() 返回函数 sync_kern() 的积​​分。函数read_table()定义如下:

size_t read_table(double *xa, double *ya)
{
    size_t datapoints;
    datapoints = 5401;
    double x[] = {
        1.388794e-11,
        1.395756e-11,
        ...       };
    double y[] = {
        5.166810e-04,
        5.175428e-04,
        ...
                 };
    int i = 0;
    while(i<datapoints){ 
        xa[i]=log(x[i]);
        ya[i]=log(y[i]);
        i++;
        }

return datapoints;  
}

sync_kern() 函数如下:

double sync_kern(double gamma, void *params)
{
    struct func_params *part= (struct func_params *)params;
    double result;
    double P,x, nuc, nu_0;
    double nu = *(double *)params;                
    gamma = exp(gamma);
    nu_0    =   (3*E_COULOMB*B*sqrt(2/3.0))/(4.0*M_PI*M*C);
    nuc     =   nu_0*sq(gamma);     
    x   =   nu/nuc;

    gsl_interp_accel *acc      =      gsl_interp_accel_alloc ();
    gsl_spline *spline     =      gsl_spline_alloc (gsl_interp_linear, 5401);


     /* double table[2][MAXPOINTS];
    size = read_table(table[0],table[1]);   
    gsl_interp_accel *acc      =      gsl_interp_accel_alloc ();
        gsl_spline *spline    =      gsl_spline_alloc (gsl_interp_interp, 5401);        
    gsl_spline_init(spline,table[0],table[1],size);  */

    P   =   gsl_spline_eval(spline,log(x),acc); 
    result  =   exp(P)*f(gamma)*gamma;
    return(result);
} 

当我包含这些行时

 /* double table[2][MAXPOINTS];
    size = read_table(table[0],table[1]);   
    gsl_interp_accel *acc      =      gsl_interp_accel_alloc ();
        gsl_spline *spline    =      gsl_spline_alloc (gsl_interp_interp, 5401);        
    gsl_spline_init(spline,table[0],table[1],size);  */

在函数 sync_kern() 中而不是在 main() 然后它起作用:我得到了 GSL_FN_EVAL 的正确值(&F,nu) 但这当然需要太多时间...

我希望我的消息不会太长...我不是 C++ 的专家，我已经查看了问题的所有方面，但我仍然不明白它来自哪里...有人吗有什么想法吗？请不要犹豫，问我是否需要任何进一步的信息。非常感谢您的帮助!

也许我应该准确地说我的代码是从主管的代码中获取的。略有不同:

int main(int argc,char *argv[])
{
    FILE *fp;
    struct func_params params;
    size_t size;
    int i;
    /* initializing parameters */
    initialize(argc, argv, &params);

    /* Interpolation of  Synchrotronfunction F(x) */

        double table[2][MAXPOINTS];
        size = read_table(table[0],table[1]);       /*Tabelle fuer Interpolation einlesen aus create_table.c*/
        params.acc = gsl_interp_accel_alloc();
        params.spline = gsl_spline_alloc(gsl_interp_linear, size);
        gsl_spline_init(params.spline,table[0],table[1],size);
...
return (0);
}   


double
sync_kern(double gamma, void *params)
{
    struct func_params *part= (struct func_params *)params;
    double result;
    double P,nu_c,x;
    gamma = exp(gamma);
    nu_c    =   part->nu_0*sq(gamma);
    x   =   part->nu_s/nu_c;
    P   =   gsl_spline_eval(part->spline,log(x),part->acc); /*Aus Interpolation: Pointer auf spline */

    result  =   exp(P)*f(gamma)*gamma;
    return(result);
}

其中 spline 和 acc 在 file.h 中定义:

struct func_params
    {   
    gsl_spline *spline;
    gsl_interp_accel *acc;
    ...
    }
void initialize(int argc, char *argv[], void *params)
{
struct func_params *part= (struct func_params *)params; 
/* Default values for parameters: */
    part->gmin=1.0e7/M_E_EV;
    part->gamma=1.0e10/M_E_EV;
    part->gmax=2.0e18/M_E_EV;   /*Integrationgrenzen fuer gamma festlegen*/
    part->nu=1e-9/H_EV; 
...
}

我在我的函数 sync_kern() 中定义了 spline 和 acc 并且他在一个结构中。有人知道错误可能来自哪里吗？

Answer 1

The function read_table() gives the arrays 'x[]' and 'y[]' used for the interpolation

double table[2][MAXPOINTS];
size = read_table(table[0],table[1]);

您定义了一个 double 组，但不对其进行初始化。它包含随机值。它很可能会产生荒谬的结果。

您将两个值传递给 read_table()。如果你想让函数改变它们，这个函数的定义应该是

size_t read_table(double& val1, double& val2);

检查您的代码。有些变量没有声明(nu_0、nuc、acc)

编辑

在您编辑问题后，我看到了 read_table() 是如何定义的。它使用两个指针。您只有一个数组，二维，固定大小。由于这个固定大小，可以使用 table[0] 和 table[1] 这两个指针来传递数组。

double table[2][N] 等同于double** table。所以传递 table[1] 就是传递一个指向 array[N] 的指针。

gsl 错误是由于 gsl:interp.c 中的这段代码造成的 https://github.com/ampl/gsl/blob/master/interpolation/interp.c

  if (x < interp->xmin || x > interp->xmax)
    {
      GSL_ERROR_VAL("interpolation error", GSL_EDOM, GSL_NAN);
    }

看来您正在尝试超出限制的插值(称为“外推”)