更强大的集成器 cpp

问题描述

我试图产生一个相对论 voigt 分布，一个相对论 Breit-Wigner distribution 和一个 gaussian function 的卷积。

MWE：

double relativisticBreitwigner_pdf(double energy,double width,double mass){
  double massSquare = pow(mass,2);
  double widthSquare = pow(width,2);
  double gamma = sqrt(massSquare*(massSquare+widthSquare));
  double k = (2*sqrt(2)*mass*width*gamma)/(M_PI*sqrt(massSquare + gamma) );
  return k/(pow((pow(energy,2)-massSquare),2) + massSquare*widthSquare);
}

double gaussian_pdf(double energy,double sigma,double mass){
  return (1.0/(sigma*sqrt(2*M_PI)))*exp(-(1.0/2.0)*pow((energy-mass)/sigma,2.0));
}

double relativisticVoigt_pdf(double energy,double mass,double range=100.0){

  auto f = [&](double dummy) { return ( relativisticBreitwigner_pdf(dummy+mass,width,mass)*gaussian_pdf(energy-dummy,sigma,mass) );};
  boost::math::quadrature::tanh_sinh<double> integrator;
  return integrator.integrate(f,-range,range);
  //  return boost::math::quadrature::trapezoidal(f,range,sqrt(std::numeric_limits<double>::epsilon()),10000);    }

这个函数 relativisticVoigt_pdf(...) 正确地产生了一个相对论 voigt 分布，但是由于 integrator.integrate(f,range) 返回的值，分布中有许多不正确的“下降” ;不正确。

如果我减小积分范围，这些下降的大小/数量会更小，但相对论 voigt 分布的范围会被截断。

附上一个屏幕截图，显示黑色相对论 voigt 有这个问题（与没有这个问题的粉红色非相对论 voigt 相比，只是为了检查倾角之外的值是否有效，黑色曲线应该接近粉红色曲线，但略高于峰左侧，略低于峰右侧，正如可以看到的那样）。

我认为问题与积分方法中的舍入误差有关，因为涉及的数字很多，但这只是一个猜测。是否有更强大/更可靠的集成器在这种情况下运行良好？

tanh_sinh 积分范围=100

梯形积分范围=250

解决方法

对于将来遇到类似问题的任何人，通过将 boost::math::quadrature::trapezoidal(...) 从默认值 (sqrt(std::numeric_limits::epsilon()= 1.48e-8，在我上面的代码中我明确输入，但如果没有输入，这是默认值）到更大的值 1e-6 我设法让梯形积分器在整个范围内工作。>

我不明白为什么会这样，特别是因为下跌不会全部归零，其中一些只是比实际值略低一点。如果它们都归零，我会理解它可能只是找不到所需精度的积分，因此返回零。如果有人理解为什么更精确的容差会导致积分错误，我想知道。

更新：

虽然上述方法有所帮助，但并没有解决问题。为了解决这个问题，您可以采用梯形返回的 L1 参数，如下所示：

double error; 
double L1; 
boost::math::quadrature::trapezoidal(f,-range,range,1e-6,10000,&error,&L1);

然后检查 L1==0 或非常小，如果是，则改变范围直到不是。