c++ - 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
}

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

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)”:模板参数用于“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 收到“程序已停止工作”消息。当我在 Debug模式下编译并运行它时，我收到以下错误消息:

所以我的问题基本上是为什么 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 正交在如此巨大的争论中评估他们的功能。它还迫使他们考虑通常不需要考虑的事情，即:

IEEE 算法不能取限制。

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

在你的例子中，你知道 cosh(x) 比 pow(|x|, 3) 增长得更快，但 IEEE 不知道't。所以你需要通过 $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 求积并将结果加倍。