GSL 数值积分不正确的值

Question

我有一个积分

static double Integrand(double k , void * params)
{
 long double *p = (long double *)params;
  long double MassSquared = p[0];
  long double Temp = p[1];
  long double sign = p[2];

  long double EXP = std::exp(-std::sqrt(k*k+MassSquared/(Temp*Temp)));
  double res = 0;
  if(sign == -1)
  {
    res = k*k*std::log(1-exp(-std::sqrt(k*k+MassSquared/(Temp*Temp))));
  }
  else if(sign == +1)
  {
    res = k*k*std::log(1+exp(-std::sqrt(k*k+MassSquared/(Temp*Temp))));
  }
  return res;
};

这个在一个名为 VEffClass 的 Class 中。

我的集成例程将是

long double VEffClass::NumIntInf(long double low, double sign, long double
MassSquared, long double Temp)
{
//long double  res = 0;
double up =
std::sqrt(std::log(C_threshold)*std::log(C_threshold)- MassSquared/(Temp*Temp));
long double StepSize = (up-low)/C_IntSteps;
long double par[3] = {MassSquared, Temp, sign};


 if(Temp == 0) return 0;

gsl_integration_workspace *w = gsl_integration_workspace_alloc(5e+5);
double result, error;
gsl_function F;
F.function = &VEffClass::Integrand;
F.params = &par;

gsl_integration_qags(&F, C_threshold, up, 1e-07, 0, 5e+5, w, &result,
&error);
gsl_integration_workspace_free(w);

return result;

}

与 C_threshold = 1e-3;

如果我运行 NumIntInf(10^-3, +1 , 100, 500) 结果是 1.83038。结果 Maple 给出的是 ~ 6*10^9 。如果我使用简单的辛普森一家 1/3 方法与 Maple 结果相差约 1%。

任何人都可以指出错误吗？

感谢任何建议

Answer 1

每当您提出涉及数值结果（以及错误）的问题时，您需要提供我们可以分析并运行重现您的原始问题的简明代码。所以是的，信息丢失了。我们需要

(1) 简洁的代码（我们不需要 VEffClass！）与 GSL 集成（您部分提供了 - 但出于堆栈溢出的目的，您的代码仍然相当复杂）。

(2) 简明代码与您对辛普森规则的实施。

(3) 您的枫木计算快照。

(2) 和 (3) 非常重要，因为您的问题是由于方法 (2) 和 (3) 与 (1) 不一致，而您认为这是不正确的答案。为什么您认为 (1) 是不正确的而不是相反？你有检查数字的解析解吗？

当我在 Wolfram Mathematica 中分析具有您在评论中提供的参数值的被积函数时，我得到了与 GSL 计算一致的答案。所以我无法重现你的问题，没有 (2) + (3) 我只能猜测出了什么问题...

我附上了我的 Mathematica 笔记本的快照。

GSL 数值积分不正确的值

GSL Incorrect values with numerical Integration

numerical-integration

gsl