boost_mirror/odeint
1
0
Fork 0
mirror of https://github.com/boostorg/odeint.git synced 2025-05-11 05:24:01 +00:00
Code Issues Packages Projects Releases Wiki Activity

test step limiter with negative dt + bugfix

Browse Source 
Bugfix in error computation with negative dt.
This commit is contained in:
Mario Mulansky 2015-10-31 14:22:29 -06:00
parent 3d87ce360e
commit 80da40f6b5
2 changed files with 31 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
include/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp
View File

@ -79,9 +79,10 @@ public:
template< class State , class Deriv , class Err, class Time > template< class State , class Deriv , class Err, class Time >
value_type error( algebra_type &algebra , const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const value_type error( algebra_type &algebra , const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const
{ {
using std::abs;
// this overwrites x_err ! // this overwrites x_err !
algebra.for_each3( x_err , x_old , dxdt_old , algebra.for_each3( x_err , x_old , dxdt_old ,
typename operations_type::template rel_error< value_type >( m_eps_abs , m_eps_rel , m_a_x , m_a_dxdt * get_unit_value( dt ) ) ); typename operations_type::template rel_error< value_type >( m_eps_abs , m_eps_rel , m_a_x , m_a_dxdt * abs(get_unit_value( dt )) ) );
// value_type res = algebra.reduce( x_err , // value_type res = algebra.reduce( x_err ,
// typename operations_type::template maximum< value_type >() , static_cast< value_type >( 0 ) ); // typename operations_type::template maximum< value_type >() , static_cast< value_type >( 0 ) );

46
test/step_size_limitation.cpp
View File

@ -30,17 +30,12 @@ typedef double value_type;
typedef std::vector< value_type > state_type; typedef std::vector< value_type > state_type;
void lorenz( const state_type &x , state_type &dxdt , const value_type t ) void damped_osc( const state_type &x , state_type &dxdt , const value_type t )
{ {
BOOST_CHECK( t >= 0.0 ); const value_type gam( 0.1);
const value_type sigma( 10.0 ); dxdt[0] = x[1];
const value_type R( 28.0 ); dxdt[1] = -x[0] - gam*x[1];
const value_type b( value_type( 8.0 ) / value_type( 3.0 ) );
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];
} }
@ -107,20 +102,24 @@ BOOST_AUTO_TEST_CASE( test_step_adjuster )
template<class ControlledStepper> template<class ControlledStepper>
void test_controlled_stepper(ControlledStepper stepper, const double max_dt) void test_controlled_stepper(ControlledStepper stepper, const double max_dt)
{ {
state_type x( 3 ); state_type x( 2 );
x[0] = x[1] = x[2] = 10.0; x[0] = x[1] = 10.0;
std::vector<double> times; std::vector<double> times;
integrate_adaptive(stepper, lorenz, x, 0.0, 1.0, max_dt, push_back_time(times)); integrate_adaptive(stepper, damped_osc, x, 0.0, 100*max_dt, max_dt, push_back_time(times));
// check that dt remains at exactly max_dt // check that dt remains at exactly max_dt
std::cout << "STEP: " << times.size() << " , END: " << times[times.size()-1] << std::endl;
for( size_t i=0 ; i<times.size() ; ++i ) for( size_t i=0 ; i<times.size() ; ++i )
// check if observer was called at times 0,1,2,... // check if observer was called at times 0,1,2,...
BOOST_CHECK_SMALL( times[i] - static_cast<double>(i)*max_dt , 1E-15); BOOST_CHECK_SMALL( times[i] - static_cast<double>(i)*max_dt , 1E-15);
times.clear(); times.clear();
return;
// this should also work when we provide some bigger initial dt // this should also work when we provide some bigger initial dt
integrate_adaptive(stepper, lorenz, x, 0.0, 1.0, 0.1, push_back_time(times)); integrate_adaptive(stepper, damped_osc, x, 0.0, 100*max_dt, 10*max_dt, push_back_time(times));
// check that dt remains at exactly max_dt // check that dt remains at exactly max_dt
for( size_t i=0 ; i<times.size() ; ++i ) for( size_t i=0 ; i<times.size() ; ++i )
// check if observer was called at times 0,1,2,... // check if observer was called at times 0,1,2,...
@ -136,25 +135,38 @@ BOOST_AUTO_TEST_CASE( test_controlled )
test_controlled_stepper(make_controlled(1E-2, 1E-2, max_dt, test_controlled_stepper(make_controlled(1E-2, 1E-2, max_dt,
runge_kutta_dopri5<state_type>()), runge_kutta_dopri5<state_type>()),
max_dt); max_dt);
test_controlled_stepper(make_controlled(1E-2, 1E-2, -max_dt,
runge_kutta_dopri5<state_type>()),
-max_dt);
std::cout << "DOPRI5 FINISHED" << std::endl;
test_controlled_stepper(make_controlled(1E-2, 1E-2, max_dt, test_controlled_stepper(make_controlled(1E-2, 1E-2, max_dt,
runge_kutta_cash_karp54<state_type>()), runge_kutta_cash_karp54<state_type>()),
max_dt); max_dt);
test_controlled_stepper(make_controlled(1E-2, 1E-2, -max_dt,
runge_kutta_cash_karp54<state_type>()),
-max_dt);
std::cout << "RK-CK54 FINISHED" << std::endl;
test_controlled_stepper(bulirsch_stoer<state_type>(1E-2, 1E-2, 1.0, 1.0, max_dt), test_controlled_stepper(bulirsch_stoer<state_type>(1E-2, 1E-2, 1.0, 1.0, max_dt),
max_dt); max_dt);
test_controlled_stepper(bulirsch_stoer<state_type>(1E-2, 1E-2, 1.0, 1.0, -max_dt),
-max_dt);
std::cout << "BS FINISHED" << std::endl;
} }
BOOST_AUTO_TEST_CASE( test_dense_out ) BOOST_AUTO_TEST_CASE( test_dense_out )
{ {
state_type x( 3 ); state_type x( 2 );
x[0] = x[1] = x[2] = 10.0; x[0] = x[1] = 10.0;
const double max_dt = 0.01; const double max_dt = 0.01;
std::vector<double> times; std::vector<double> times;
integrate_adaptive(make_dense_output(1E-2, 1E-2, max_dt, runge_kutta_dopri5<state_type>()), integrate_adaptive(make_dense_output(1E-2, 1E-2, max_dt, runge_kutta_dopri5<state_type>()),
lorenz, x, 0.0, 1.0, max_dt, push_back_time(times)); damped_osc, x, 0.0, 1.0, max_dt, push_back_time(times));
// check that dt remains at exactly max_dt // check that dt remains at exactly max_dt
for( size_t i=0 ; i<times.size() ; ++i ) for( size_t i=0 ; i<times.size() ; ++i )
// check if observer was called at times 0,1,2,... // check if observer was called at times 0,1,2,...
@ -163,7 +175,7 @@ BOOST_AUTO_TEST_CASE( test_dense_out )
// this should also work when we provide some bigger initial dt // this should also work when we provide some bigger initial dt
integrate_adaptive(make_dense_output(1E-2, 1E-2, max_dt, runge_kutta_dopri5<state_type>()), integrate_adaptive(make_dense_output(1E-2, 1E-2, max_dt, runge_kutta_dopri5<state_type>()),
lorenz, x, 0.0, 1.0, 0.1, push_back_time(times)); damped_osc, x, 0.0, 1.0, 0.1, push_back_time(times));
// check that dt remains at exactly max_dt // check that dt remains at exactly max_dt
for( size_t i=0 ; i<times.size() ; ++i ) for( size_t i=0 ; i<times.size() ; ++i )
// check if observer was called at times 0,1,2,... // check if observer was called at times 0,1,2,...