尝试将插值函数与 C GSL 集成

Question

我是 C 语言的新手，我在类型声明和指针方面遇到了一些麻烦。

我正在尝试对数据集进行插值然后对其进行积分，这两种数值方法都是使用 C GSL 库完成的。该代码改编自 GSL 关于 monte carlo 积分和 2d 插值的手动示例。在使用数据集之前，我正在使用已知且非常简单的 2D 函数进行测试。

二维插值部分运行良好，我可以从插值函数中获取值，interpFun。

但是我不能使用这个插值函数作为积分算法的输入。

#include <stdio.h>
#include <stdlib.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_vegas.h>

int
main()
{
  const gsl_interp2d_type *T = gsl_interp2d_bilinear;
  const size_t N = 100;             /* number of points to interpolate */
  double xa[] = { 0.0, 1.0 }; /* define unit square */
  double ya[] = { 0.0, 1.0 };
  double xl[] = { 0., 0.};
  double xu[] = { 1., 1.}, xp[] = {1.0, 0.0}, zz[] = {.2,.6};
  const size_t nx = sizeof(xa) / sizeof(double); /* x grid points */
  const size_t ny = sizeof(ya) / sizeof(double); /* y grid points */
  double *za = malloc(nx * ny * sizeof(double));
  double fun( double *x, size_t dim, void *params);
  gsl_spline2d *spline = gsl_spline2d_alloc(T, nx, ny);
  gsl_interp_accel *xacc = gsl_interp_accel_alloc();
  gsl_interp_accel *yacc = gsl_interp_accel_alloc();
  size_t i, j;
  const gsl_rng_type *S; 
  gsl_rng *r; 
  double res, err;
  struct fparams { gsl_spline2d *spline; gsl_interp_accel *xacc; gsl_interp_accel *yacc; };

  /* set z grid values */
  gsl_spline2d_set(spline, za, 0, 0, fun( xl, 2, 0) );
  gsl_spline2d_set(spline, za, 0, 1, fun( xa, 2, 0) );
  gsl_spline2d_set(spline, za, 1, 1, fun( xu, 2, 0) );
  gsl_spline2d_set(spline, za, 1, 0, fun( xp, 2, 0) );

  /* initialize interpolation */
  gsl_spline2d_init(spline, xa, ya, za, nx, ny);

  struct fparams params1 = { spline, xacc, yacc };

  double interpFun( double *x, size_t dim, struct fparams p ) {
  (void)(dim);
  struct fparams *par = &p;
  return gsl_spline2d_eval(par->spline, x[0], x[1], par->xacc, par->yacc);
   };

  gsl_monte_function G = { &interpFun, 2, &params1 };
  size_t calls = 500000;
  gsl_rng_env_setup ();
  S = gsl_rng_default;
  r = gsl_rng_alloc (S);
  {
    gsl_monte_vegas_state *s = gsl_monte_vegas_alloc (2);
    gsl_monte_vegas_integrate (&G, xl, xu, 2, 10000, r, s, &res, &err);
    printf ("result = %.6f\n", res);
    printf ("sigma = %.6f\n", err);
    printf ("converging...\n");
    do
      {
       gsl_monte_vegas_integrate (&G, xl, xu, 2, calls/5, r, s, &res, &err);
       printf ("result = %.6f sigma = %.6f "
               "chisq/dof = %.1f\n", res, err, gsl_monte_vegas_chisq (s));
     }
    while (fabs (gsl_monte_vegas_chisq(s)-1.0) > 0.5);
    printf ("result = %.6f\n", res);
    printf ("sigma = %.6f\n", err);
    gsl_monte_vegas_free (s);
  }

  gsl_rng_free (r);
  gsl_spline2d_free(spline);
  gsl_interp_accel_free(xacc);
  gsl_interp_accel_free(yacc);
  free(za);

  return 0;
}

double fun( double *x, size_t dim, void *p ) {
 (void)(dim);
 (void)(p);
 return 3*x[0]+5*x[1];
}

这段代码returns错误信息：

inter_test.c: In function ‘main’:
inter_test.c:55:3: warning: initialization from incompatible pointer type [enabled by default]
   gsl_monte_function G = { &interpFun, 2, &params1 };
   ^
inter_test.c:55:3: warning: (near initialization for ‘G.f’) [enabled by default]

如果我尝试集成 fun 而不是 interpFun，则没有错误消息，但两者是同一类型。

Answer 1

您程序中定义的 interpFun() 函数的类型为 double (*)(double *, size_t, struct fparams)，但根据 the GSL docs，您需要的类型为 double (*)(double *, size_t, void *)。这是一个修复：

double interpFun( double *x, size_t dim, void *p ) {
    (void)(dim);
    struct fparams *par = (struct fparams *) p;
    return gsl_spline2d_eval(par->spline, x[0], x[1], par->xacc, par->yacc);
}

此外，C 不允许嵌套函数定义，尽管有 GCC extension for this.

尝试将插值函数与 C GSL 集成

Trying to integrate an interpolated function with C GSL

c

interpolation

numerical-integration

gsl