boost::math:quadrature::sinh_sinh中的Bug？

Question

我在整个实行(-inf，inf)上搜索一个数字求积库，我找到了boost (版本1.70.0)。我想使用的函数是boost::math::quadrature:sinh_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编译。错误是：

错误C2783："T boost::math::tools::root_epsilon(void)"：模板-参数für "T“konnte nicht hergeleitet werden。(英文:它无法导出"T“的模板-参数)

发生在pathTo\boost_1_70_0\boost\math\quadrature\sinh_sinh.hpp的第33行

由于我不确定此错误是否仅与有关，所以我想问所有的人。

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

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收到了“程序停止工作”的消息。在调试模式下编译并运行它时，会收到以下错误消息：

​

​

所以我的问题基本上是为什么boost::math::quadrature::sinh_sinh不能集成这个函数。对于正无穷大和负无穷远，它衰变为零，而且它没有奇点。

是否可能所有这些错误都是因为我使用Visual而发生的？

Answer 1

不幸的是，Visual对你不友好。在第二个例子中，我得到了更容易理解的错误消息：

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

大多数人都很惊讶地得知，辛辛-辛正交在这么大的争论中评估了他们的功能。它也迫使他们去思考他们通常不需要考虑的事情，即：

IEEE算法不能接受极限。

例如，您可能知道作为$x \to \infty$，$x^2/(1+x^4) \to 0$。但是在IEEE浮点算法中，对于足够大的$x$，分子溢出和分母溢出，可以做些什么呢？唯一明智的解决方案就是使inf/inf成为nan。

在您的例子中，您知道cosh(x)的增长速度比pow(|x|, 3)快，但是IEEE没有。因此，您需要将限制行为的函数显式地告诉函数$x->\infty$ via：

#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;
}

最后一个注释:您的integrand是偶数，所以您可以在[0, inf]上使用[0, inf]求积并将结果加倍。