ODEint:具有任意精度的自适应积分
ODEint: adaptive integration with arbitrary precision
ODEint 是否可以使用具有任意精度算术的自适应 积分例程?例如,我想将带有 integrate_adaptive() 函数的 Boost 多精度库与受控步进器一起使用。 ODEint 文档提供了对 integrate_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;
}
编译这个returns错误:
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 的 typedef 更改为 double 它编译并且 运行 没问题。
odeint 可以使用任意精度的自适应积分器。您的代码几乎是正确的,您只是忘了将 value_type
(用于步进器的内部常量,以及 "time" 变量 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/)。因此,您最终会得到一个不完全为 0.1 的 multi_prec 值。我建议始终从整数转换,例如将 0.1 表示为 value_type(1)/value_type(10)
.
ODEint 是否可以使用具有任意精度算术的自适应 积分例程?例如,我想将带有 integrate_adaptive() 函数的 Boost 多精度库与受控步进器一起使用。 ODEint 文档提供了对 integrate_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;
}
编译这个returns错误:
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 的 typedef 更改为 double 它编译并且 运行 没问题。
odeint 可以使用任意精度的自适应积分器。您的代码几乎是正确的,您只是忘了将 value_type
(用于步进器的内部常量,以及 "time" 变量 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/)。因此,您最终会得到一个不完全为 0.1 的 multi_prec 值。我建议始终从整数转换,例如将 0.1 表示为 value_type(1)/value_type(10)
.