boost::math:quadrature::sinh_sinh 中的错误？

Question

我正在寻找一个用于在整个实数线上进行数字求积的库，即 (-inf,inf)，然后我找到了 boost（版本 1.70.0）。我要使用的函数是boost::math::quadrature:sinh_sinh。为了测试它，我从文档中复制了示例代码：

https://www.boost.org/doc/libs/1_70_0/libs/math/doc/html/math_toolkit/double_exponential/de_sinh_sinh.html

并想出了这个代码：

#include <iostream>
#include <boost/math/quadrature/sinh_sinh.hpp>

using namespace boost::math::quadrature;

int main()
{
    sinh_sinh<double> integrator;
    auto f = [](double x) { return exp(-x*x); };
    double error;
    double L1;
    double Q = integrator.integrate(f, &error, &L1);
    std::cout << Q << "\t" << error << "\t" << L1 << std::endl;

    int i = 0;
    std::cin >> i; // Just to make sure the console does not close automatically
}

不幸的是，这不会编译，因为在文档中，"integrate" 的第二个参数不是指向实数的指针，而是一个普通的实数。所以我不得不改变这一行：

double Q = integrator.integrate(f, &error, &L1);

进入这个：

double Q = integrator.integrate(f , boost::math::tools::root_epsilon<double>() , &error, &L1);

这个编译并给出了很好的结果。但是我很好奇我是否可以写

double Q = integrator.integrate(f);

因为除了第一个参数之外的所有参数都有默认值（因此对于我对 c++ 的理解是可选的）。不幸的是，这不会用 Visual-Studio-2013 编译。错误是：

error C2783: "T boost::math::tools::root_epsilon(void)": template-Argument für "T" konnte nicht hergeleitet werden。（英语：无法导出 "T" 的模板参数）

出现在pathTo\boost_1_70_0\boost\math\quadrature\sinh_sinh.hpp

的第33行

因为我不确定这个错误是否只与Visual-Studio有关所以想问问大家。

现在我想在我感兴趣的函数上使用工作代码，即：

auto f = [](double x) {return pow(abs(x), 3) / cosh(x); };

这个函数看起来像这样：

https://www.wolframalpha.com/input/?i=plot+abs(x)%5E3%2Fcosh(x)

求积的结果应该是大约。 23.7:

https://www.wolframalpha.com/input/?i=integrate+abs(x)%5E3%2Fcosh(x)+from+-+inf+to+inf

这个程序用这个函数编译但它崩溃了，即我从 Windows 得到了 "The program has stopped working" 消息。当我在调试模式和运行下编译时，我收到以下错误消息：

所以我的问题基本上是为什么boost::math::quadrature::sinh_sinh不能集成这个功能。它在正负无穷大时衰减为零，并且没有奇点。

是否有可能是因为我使用的是 Visual-Studio 而出现所有这些错误？

Answer 1

很遗憾，Visual studio 对你不好。在你的第二个例子中，我得到了更容易理解的错误信息：

terminate called after throwing an instance of 'boost::wrapexcept<boost::math::evaluation_error>'
  what():  Error in function boost::math::quadrature::sinh_sinh<double>::integrate: The sinh_sinh quadrature evaluated your function at a singular point, leading to the value nan.
sinh_sinh quadrature cannot handle singularities in the domain.
If you are sure your function has no singularities, please submit a bug against boost.math

我添加了一些诊断代码来帮忙：

auto f = [](double x) {
    double y = pow(abs(x), 3) / cosh(x);
    if (!std::isfinite(y)) {
        std::cout << "f(" << x << ") = " << y << "\n";
    }
    return y;
};

我得到：

f(1.79769e+308) = nan
f(-1.79769e+308) = nan
f(2.01977e+137) = nan
f(-2.01977e+137) = nan
f(7.35294e+106) = nan
f(-7.35294e+106) = nan

得知 sinh-sinh quadrature 在如此巨大的参数下评估其函数，大多数人都感到非常惊讶。它还迫使他们考虑他们通常不必考虑的事情，即：

IEEE算法不能取极限

例如，您可能知道 $x \to \infty$、$x^2/(1+x^4) \to 0$。但是在IEEE浮点运算中，对于足够大的$x$，分子和分母都溢出了，怎么办？唯一明智的解决方案是让 inf/inf 成为 nan.

在您的情况下，您知道 cosh(x) 比 pow(|x|, 3) 增长得更快，但 IEEE 不知道。因此，您需要通过 $x->\infty$ 显式地告诉函数有关限制行为的信息：

#include <iostream>
#include <cmath>
#include <boost/math/quadrature/sinh_sinh.hpp>

using namespace boost::math::quadrature;

int main()
{
    sinh_sinh<double> integrator;
    auto f = [](double x) {
        double numerator = pow(abs(x), 3);
        if (!std::isfinite(numerator)) { 
            return double(0);
        }
        return numerator / cosh(x);
    };
    double error;
    double L1;
    double tolerance = std::sqrt(std::numeric_limits<double>::epsilon());
    double Q = integrator.integrate(f, tolerance, &error, &L1);
    std::cout << Q << "\t" << error << "\t" << L1 << std::endl;
}

最后一条评论：您的被积函数是偶数，因此您可以对 [0, inf] 使用 exp_sinh 求积并将结果加倍。

boost::math:quadrature::sinh_sinh 中的错误？

Bug in boost::math:quadrature::sinh_sinh?

c++

boost

visual-studio-2013