Ошибка в boost :: math: quadrature :: sinh_sinh? - PullRequest
Новая
       22

Ошибка в boost :: math: quadrature :: sinh_sinh?

1 голос
/ 02 июня 2019

Я искал библиотеку для числовой квадратуры по всей реальной строке, т. Е. (-Inf, inf), и нашел повышение (версия 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);

, потому что все arugments, кроме первого, имеют значения по умолчанию (и, следовательно, являются необязательными для моего понимания c ++).К сожалению, это не скомпилируется с Visual-Studio-2013.Ошибка:

ошибка C2783: "T boost :: math :: tools :: root_epsilon (void)": template-Argument für "T" konnte nicht hergeleitet werden.(по-английски: он не смог получить аргумент шаблона для «T»)

Получается в строке 33 pathTo \ boost_1_70_0 \ boost \ math \ quadrature \ sinh_sinh.hpp

AsЯ не уверен, что эта ошибка связана только с 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.Когда я компилирую в режиме отладки и запускаю его, я получаю следующее сообщение об ошибке:

enter image description here

Так что мой вопрос в основном, почему boost :: math :: quadrature:: sinh_sinh не может интегрировать эту функцию.Он уменьшается до нуля для плюс и минус бесконечность и не имеет особенностей.

Возможно ли, что все эти ошибки происходят из-за того, что я использую Visual-Studio?

1 Ответ

1 голос
/ 02 июня 2019

К сожалению, 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

Большинство людей очень удивлены, узнав, что квадратурный синх-синх оценивает их функцию при таком огромном споре. Это также заставляет их думать о вещах, которые им обычно не нужны, а именно:

Арифметика IEEE не может принимать ограничения.

Например, вы можете знать, что как $ x \ to \ infty $, $ x ^ 2 / (1 + x ^ 4) \ to 0 $. Но в арифметике IEEE с плавающей запятой при достаточно больших $ x $ переполнение числителя и знаменателя и что можно сделать? Единственное разумное решение - просто сделать инф / инф нан.

В вашем случае, вы знаете, что 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;
}

Один последний комментарий: ваша интеграция четная, поэтому вы можете использовать exp_sinh квадратуры над [0, inf] и удвоить результат.

...