c - 使用 GSL 在 C 中进行高效的双重集成

Question

考虑二重积分

I = int int [(a^k)*b] da db

我们要为 [0,1] 之间的 a 和 [0,1] 之间的 b 集成code> 和 k 是一些常量。我正在使用 GSL 数值积分库，但遇到内存分配问题。

我的代码如下

#include <stdlib.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_integration.h>

double innerIntegrand(double a, void *params) {
    double *cast_params = (double *) params;
    double b = params[0];
    double k = params[1];

    return pow(a,k)*b;
}

然后我可以计算给定 b 的内部积分(以获得外部被积函数)，如下所示

double outerIntegrand(double b, void *params) {
    // params = {holder for double b, k}
    double *cast_params = (double *) params;
    cast_params[0] = b;

    // Allocate integration workspace
    gsl_integration_workspace *giw = gsl_integration_workspace_alloc(100);

    // Create GSL function
    gsl_function F;
    F.function = &innerIntegrand;
    F.params = params;

    // Initialise values to put the result in
    double result;
    double abserror;

    // Perform integration
    gsl_integration_qag(&F, 0, 1, 0.001, 0.001, 100, 1, giw, &result, &abserror);

    // Free the integration workspace
    gsl_integration_workspace_free(giw);

    // Return result
    return result
}

但是请注意，我必须在函数内分配和释放集成工作区。这意味着在评估最终的集成函数时会进行多次

double Integral(double k) {
    // Create params
    double *params = malloc(2*sizeof(double));
    params[1] = k;

    // Allocate integration workspace
    gsl_integration_workspace *giw = gsl_integration_workspace_alloc(100);

    // Create GSL function
    gsl_function F;
    F.function = &outerIntegrand;
    F.params = params;

    // Initialise values to put the result in
    double result;
    double abserror;

    // Perform integration
    gsl_integration_qag(&F, 0, 1, 0.001, 0.001, 100, 1, giw, &result, &abserror);

    // Free the integration workspace
    gsl_integration_workspace_free(giw);


    // Free memory
    free(params);

    // Return result
    return result
}

理想情况下，我想要的是两个全局 gsl_integration_workspace 变量，一个用于 outerIntegrand 中的积分，另一个用于 Integral 中的积分。但是，当我尝试将它们声明为全局值时，我收到了一个initializer element is not constant错误。

谁能找到一种无需重复内存分配和释放即可完成此二重积分的方法？我在想我们也可以通过 params 参数传递工作区，尽管它随后开始变得非常困惑。

Answer 1

我设法用 C++ 构建了一个外观不错的程序，用于基于 GSL 的双重集成，以干净的方式避免了重复分配。我用这个众所周知的函数来玩:

f(x,y)=exp(-x*x-y*y)

在整个平面上对其进行积分(结果 pi 可以通过切换到极坐标轻松获得)。修改它并通过 lambda 捕获添加参数是微不足道的。

#include <iostream>
#include <gsl/gsl_integration.h>

// Simple RAII wrapper 
class IntegrationWorkspace {
  gsl_integration_workspace * wsp;

  public:
  IntegrationWorkspace(const size_t n=1000):
    wsp(gsl_integration_workspace_alloc(n)) {}
  ~IntegrationWorkspace() { gsl_integration_workspace_free(wsp); }

  operator gsl_integration_workspace*() { return wsp; }
};

// Build gsl_function from lambda
template <typename F>
class gsl_function_pp: public gsl_function {
  const F func;
  static double invoke(double x, void *params) {
    return static_cast<gsl_function_pp*>(params)->func(x);
  }
  public:
  gsl_function_pp(const F& f) : func(f) {
    function = &gsl_function_pp::invoke; //inherited from gsl_function
    params   = this;                     //inherited from gsl_function
  }
  operator gsl_function*(){return this;}
};

// Helper function for template construction
template <typename F>
gsl_function_pp<F> make_gsl_function(const F& func) {
  return gsl_function_pp<F>(func);
}

int main() {
  double epsabs = 1e-8;
  double epsrel = 1e-8;
  size_t limit = 100;
  double result, abserr, inner_result, inner_abserr;

  IntegrationWorkspace wsp1(limit);
  IntegrationWorkspace wsp2(limit);

  auto outer = make_gsl_function( [&](double x) {
    auto inner = make_gsl_function( [&](double y) {return exp(-x*x-y*y);} );
    gsl_integration_qagi(inner, epsabs, epsrel, limit, wsp1,
                         &inner_result, &inner_abserr);
    return inner_result;
  } );
  gsl_integration_qagi(outer, epsabs, epsrel, limit, wsp2, &result, &abserr);

  std::cout << result << std::endl;
}