ODEint: адаптивная интеграция с произвольной точностью

Может ли ODEint использовать адаптивные процедуры интеграции с произвольной арифметикой точности? Например, я хотел бы использовать библиотеки мультиточности Boost с функцией Integrated_adaptive() с управляемым степпером. В документации ODEint приводятся примеры использования арифметики произвольной точности для Integrated_const(), но я не могу модифицировать их для использования адаптивного интегратора.

Я также пытался использовать итераторы (например, make_adaptive_time_iterator...), но столкнулся с похожими проблемами. Для конкретности это простой код, который я ищу для работы:

#include <iostream>

//[ mp_lorenz_defs
#include <boost/numeric/odeint.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>

using namespace std;
using namespace boost::numeric::odeint;

typedef boost::multiprecision::cpp_dec_float_50 value_type;
//typedef double value_type;

typedef boost::array< value_type , 3 > state_type;
//]

//[ mp_lorenz_rhs
struct lorenz
{
    void operator()( const state_type &x , state_type &dxdt , value_type t ) const
    {
        const value_type sigma( 10 );
        const value_type R( 28 );
        const value_type b( value_type( 8 ) / value_type( 3 ) );

        dxdt[0] = sigma * ( x[1] - x[0] );
        dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
        dxdt[2] = -b * x[2] + x[0] * x[1];
    }
};
//]

int main( int argc , char **argv )
{
    //[ mp_lorenz_int
    state_type x = {{ value_type( 10.0 ) , value_type( 10.0 ) , value_type( 10.0 ) }};

    auto stepper = make_controlled( 1.0e-16 , 1.0e-16 , runge_kutta_cash_karp54< state_type >() );


    cout.precision(50);
    integrate_adaptive( stepper ,
           lorenz() , x , value_type( 0.0 ) , value_type( 0.1 ) , value_type( value_type( 1.0 ) / value_type( 2000.0 ) ) );
    //]
    cout << x[0] << endl;

    return 0;
}

Компиляция возвращает ошибку:

Lorenz_mp2.cpp:52:19: error: no matching function for call to 'make_controlled'
        auto stepper = make_controlled( value_type(1.0e-16) , value_type(1.0e-16) , runge_kutta_cash_karp54< state_type >() );

Если я изменю typedef для value_type на double, он скомпилируется и запустится нормально.


c++ odeint
person nielsw    schedule 30.03.2016    source источник
comment
Просто для интереса, почему вы хотите использовать здесь многоточную арифметику? Я бы предположил, что ошибка, вызванная дискретизацией системы, намного больше, чем ошибка, вызванная ограниченной точностью двойников.   -  person Frank Puffer    schedule 30.03.2016
comment
Я просто использую систему Лоренца в качестве теста. Адаптивный интегратор сделает дискретизацию достаточно тонкой, чтобы позволить ошибке упасть ниже машинной точности, если используется арифметика с произвольной точностью.   -  person nielsw    schedule 30.03.2016

Ответы (1)


arrow_upward
1
arrow_downward

Использование адаптивных интеграторов с произвольной точностью возможно с помощью odeint. Ваш код почти правильный, вы только забыли также настроить value_type (используется для внутренних констант шагового двигателя и для переменной «время» t) как произвольный тип точности. Если вы вернетесь к документации (http://headmyshoulder.github.io/odeint-v2/doc/boost_numeric_odeint/tutorial/using_arbitrary_precision_floating_point_types.html), вы увидите, что это делается вторым аргументом шаблона для шагового двигателя. Таким образом, правильное определение шагового двигателя в make_controlled должно быть:

runge_kutta_cash_karp54< state_type , value_type >()

При этом он должен компилироваться и работать с произвольной точностью.

Независимо от этой проблемы, я хотел бы добавить, что использование преобразования из типа double в тип multi_prec потенциально проблематично. Например, 0.1 нельзя точно представить как двойное число (http://www.exploringbinary.com/why-0-point-1-does-not-exist-in-floating-point/). Следовательно, вы получите значение multi_prec, которое не точно равно 0,1. Я бы посоветовал всегда преобразовывать целые числа и, например, выражать 0,1 как value_type(1)/value_type(10).

person mariomulansky    schedule 30.03.2016
comment
Спасибо! Ваш комментарий о преобразовании из типа double в тип multi_prec также был очень полезен, так как это сразу вызвало проблемы после компиляции кода. Исправление заставило все работать красиво. - person nielsw; 30.03.2016