拆分多个变量 C++ 函数以进行正交积分

Question

我想与boost::math::quadrature::trapezoidal(g, a, b, 1e-6); 进行数值积分这里我正在积分函数g(x)。问题是我必须执行双重积分。此外，我要集成的函数中有 4 个变量。其中 2 个在积分 (m,n) 时通过，另外 2 个是积分变量 (r,z)。这是我要计算的积分：

$$ \int_0^b\int_0^af(r,z)\sin{(\frac{n\pi}{a}z)}J_0(\frac{\alpha_{0,m}}{b}r)dzdr $$

我看到了这个例子 Performing 2d numerical integration with Boost Cpp 并注意到他使用 lambda 函数将主被积函数拆分为 2。到目前为止，我已经完成了这个

double integrate(int m, int n)
{ 
  auto f1 = [](double r, double z, int m, int n) { return integrand(r,z,m,n); };
  auto f = [&](double r, m) {
         auto g = [&](double z, n) {
            return f1(r, z);
              };
         //return gauss_kronrod<double, 61>::integrate(g, 0, a, 5);
         return boost::math::quadrature::trapezoidal(g, 0, a, 1e-6);
       };
  double error;
  //double Q = gauss_kronrod<double, 15>::integrate(f, 0, b, 5, 1e-9, &error);
  double Q =  boost::math::quadrature::trapezoidal(f, 0, b, 1e-6);
  //std::cout << Q << ", error estimated at " << error <<std::endl;
  return Q;
}

函数$f(r,z)$的实现和积分的其余部分如下

double initial(double r, double z, int m, int n)
{
return std::sin(M_PI*n*z/a)*std::cyl_bessel_j(0, boost::math::cyl_bessel_j_zero(0,m)*r/b);
}
double integrand(double r,double z,int n,int m)
{
  return initial(r,z,m,n)*std::sin(M_PI*n*z/a)*std::cyl_bessel_j(0, boost::math::cyl_bessel_j_zero(0,m)*r/b);
}

通常 Initial 不需要它们和 n 个变量，但在这种情况下，我需要进行一些测试。

问题是我真的不明白如何像我的问题示例中那样拆分我的函数并执行集成，因为 boost 只接受 1 个变量函数。

请帮忙

Answer 1

与往常一样，基本思想是分两步集成。为此，您首先求解内部积分并从中生成另一个一维函数，然后再次将其传递给积分器。

每当您想将多参数函数缩减为单参数函数时，都会使用 lambda。在这种情况下，您将所有非积分变量放入 lambda 捕获中。

这是伪代码：

double integrand(double r,double z, int m,int n, double a, double b)
{
   //this is the function you want to integrate
}

double integrate(int m, int n)
{
    double a=1.0;
    double b=1.0;

    auto inner_integral = [m,n,a,b](double z)
    {
        auto f = [z,m,n,a,b](double r) { return integrand(r,z,m,n,a,b);} 
        return trapezoidal(f,0,a);
    }
    return trapezoidal(inner_integral,0,b);
};

您可能不需要写出 lambda 捕获，即带有 & 的引用捕获可能也可以工作 (auto inner_integral = [&](double z){...})。