2016-03-30 2 views
3

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

Я также попытался использовать итераторы (например, 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>()); 

Если изменить ЬурейеЕ для value_type удвоить компилирует и работает нормально.

+0

Просто для интереса, почему вы хотите использовать арифметику с несколькими точками? Я бы предположил, что ошибка, вызванная дискретизацией системы, намного больше, чем ошибка, вызванная ограниченной точностью удвоений. –

+0

Я просто использую систему Лоренца в качестве тестового примера. Адаптивный интегратор сделает дискретизацию достаточно хорошей, чтобы ошибка упала ниже точности машины - если используется арифметика произвольной точности. – nielsw

ответ

1

Использование адаптивных интеграторов с произвольной точностью возможно с помощью 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).

+0

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