diff --git a/include/boost/numeric/odeint/stepper/generation/generation_rosenbrock4.hpp b/include/boost/numeric/odeint/stepper/generation/generation_rosenbrock4.hpp index d5322c8d..366cb5e9 100644 --- a/include/boost/numeric/odeint/stepper/generation/generation_rosenbrock4.hpp +++ b/include/boost/numeric/odeint/stepper/generation/generation_rosenbrock4.hpp @@ -54,12 +54,19 @@ struct dense_output_factory< Stepper , rosenbrock4_dense_output< rosenbrock4_con typedef Stepper stepper_type; typedef rosenbrock4_controller< stepper_type > controller_type; typedef typename stepper_type::value_type value_type; + typedef typename stepper_type::time_type time_type; typedef rosenbrock4_dense_output< controller_type > dense_output_type; dense_output_type operator()( value_type abs_error , value_type rel_error , const stepper_type &stepper ) { return dense_output_type( controller_type( abs_error , rel_error , stepper ) ); } + + dense_output_type operator()( value_type abs_error , value_type rel_error , + time_type max_dt, const stepper_type &stepper ) + { + return dense_output_type( controller_type( abs_error , rel_error , max_dt , stepper ) ); + } }; diff --git a/include/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp b/include/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp index df4e6c48..0e5edd32 100644 --- a/include/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp +++ b/include/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp @@ -27,6 +27,7 @@ #include #include +#include #include @@ -55,12 +56,20 @@ public: typedef rosenbrock4_controller< Stepper > controller_type; - rosenbrock4_controller( value_type atol = 1.0e-6 , value_type rtol = 1.0e-6 , const stepper_type &stepper = stepper_type() ) - : m_stepper( stepper ) , m_atol( atol ) , m_rtol( rtol ) , - m_first_step( true ) , m_err_old( 0.0 ) , m_dt_old( 0.0 ) , - m_last_rejected( false ) + rosenbrock4_controller( value_type atol = 1.0e-6 , value_type rtol = 1.0e-6 , + const stepper_type &stepper = stepper_type() ) + : m_stepper( stepper ) , m_atol( atol ) , m_rtol( rtol ) , + m_max_dt( static_cast(0) ) , + m_first_step( true ) , m_err_old( 0.0 ) , m_dt_old( 0.0 ) , + m_last_rejected( false ) { } + rosenbrock4_controller( value_type atol, value_type rtol, time_type max_dt, + const stepper_type &stepper = stepper_type() ) + : m_stepper( stepper ) , m_atol( atol ) , m_rtol( rtol ) , m_max_dt( max_dt ) , + m_first_step( true ) , m_err_old( 0.0 ) , m_dt_old( 0.0 ) , + m_last_rejected( false ) + { } value_type error( const state_type &x , const state_type &xold , const state_type &xerr ) { @@ -104,6 +113,14 @@ public: boost::numeric::odeint::controlled_step_result try_step( System sys , const state_type &x , time_type &t , state_type &xout , time_type &dt ) { + if( m_max_dt != static_cast(0) && detail::less_with_sign(m_max_dt, dt, dt) ) + { + // given step size is bigger then max_dt + // set limit and return fail + dt = m_max_dt; + return fail; + } + BOOST_USING_STD_MIN(); BOOST_USING_STD_MAX(); using std::pow; @@ -142,7 +159,11 @@ public: min BOOST_PREVENT_MACRO_SUBSTITUTION ( dt_new , dt ) : max BOOST_PREVENT_MACRO_SUBSTITUTION ( dt_new , dt ) ); t += dt; - dt = dt_new; + // limit step size to max_dt + if( m_max_dt != static_cast(0) ) + { + dt = detail::min_abs(m_max_dt, dt_new); + } m_last_rejected = false; return success; } @@ -198,6 +219,7 @@ private: wrapped_state_type m_xerr; wrapped_state_type m_xnew; value_type m_atol , m_rtol; + time_type m_max_dt; bool m_first_step; value_type m_err_old , m_dt_old; bool m_last_rejected; diff --git a/test/step_size_limitation.cpp b/test/step_size_limitation.cpp index dfacee91..2d8e4725 100644 --- a/test/step_size_limitation.cpp +++ b/test/step_size_limitation.cpp @@ -29,6 +29,10 @@ using namespace boost::numeric::odeint; typedef double value_type; typedef std::vector< value_type > state_type; +/*********************************************** + * first part of the tests: explicit methods + *********************************************** + */ void damped_osc( const state_type &x , state_type &dxdt , const value_type t ) { @@ -46,7 +50,8 @@ struct push_back_time push_back_time( std::vector< double > × ) : m_times( times ) { } - void operator()( const state_type &x , double t ) + template + void operator()( const State &x , double t ) { m_times.push_back( t ); } @@ -90,7 +95,7 @@ BOOST_AUTO_TEST_CASE( test_step_adjuster ) dt_new = limited_adjuster.increase_step(dt, 0.4, 1); // but even for smaller errors, we should at most get 0.1 - BOOST_CHECK(dt_new == dt); + BOOST_CHECK_EQUAL(dt_new, dt); dt_new = limited_adjuster.increase_step(0.9*dt, 0.1, 1); std::cout << dt_new << std::endl; @@ -99,27 +104,28 @@ BOOST_AUTO_TEST_CASE( test_step_adjuster ) } -template -void test_controlled_stepper(ControlledStepper stepper, const double max_dt) +template +void test_explicit_stepper(Stepper stepper, const double max_dt) { state_type x( 2 ); x[0] = x[1] = 10.0; + const size_t steps = 100; std::vector times; - integrate_adaptive(stepper, damped_osc, x, 0.0, 100*max_dt, max_dt, push_back_time(times)); + integrate_adaptive(stepper, damped_osc, x, 0.0, steps*max_dt, max_dt, push_back_time(times)); + + BOOST_CHECK_EQUAL(times.size(), steps+1); // 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(i)*max_dt , 1E-15); times.clear(); - return; // this should also work when we provide some bigger initial dt - integrate_adaptive(stepper, damped_osc, x, 0.0, 100*max_dt, 10*max_dt, push_back_time(times)); + integrate_adaptive(stepper, damped_osc, x, 0.0, steps*max_dt, 10*max_dt, push_back_time(times)); + + BOOST_CHECK_EQUAL(times.size(), steps+1); // check that dt remains at exactly max_dt for( size_t i=0 ; i()), max_dt); - test_controlled_stepper(make_controlled(1E-2, 1E-2, -max_dt, + test_explicit_stepper(make_controlled(1E-2, 1E-2, -max_dt, runge_kutta_dopri5()), -max_dt); - std::cout << "DOPRI5 FINISHED" << std::endl; - test_controlled_stepper(make_controlled(1E-2, 1E-2, max_dt, + test_explicit_stepper(make_controlled(1E-2, 1E-2, max_dt, runge_kutta_cash_karp54()), max_dt); - test_controlled_stepper(make_controlled(1E-2, 1E-2, -max_dt, + test_explicit_stepper(make_controlled(1E-2, 1E-2, -max_dt, runge_kutta_cash_karp54()), -max_dt); - std::cout << "RK-CK54 FINISHED" << std::endl; - test_controlled_stepper(bulirsch_stoer(1E-2, 1E-2, 1.0, 1.0, max_dt), + test_explicit_stepper(bulirsch_stoer(1E-2, 1E-2, 1.0, 1.0, max_dt), max_dt); - test_controlled_stepper(bulirsch_stoer(1E-2, 1E-2, 1.0, 1.0, -max_dt), + test_explicit_stepper(bulirsch_stoer(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 ) { - state_type x( 2 ); - x[0] = x[1] = 10.0; const double max_dt = 0.01; + test_explicit_stepper(make_dense_output(1E-2, 1E-2, max_dt, + runge_kutta_dopri5()), + max_dt); + test_explicit_stepper(make_dense_output(1E-2, 1E-2, -max_dt, + runge_kutta_dopri5()), + -max_dt); + + test_explicit_stepper(bulirsch_stoer_dense_out(1E-2, 1E-2, 1, 1, max_dt), + max_dt); + + test_explicit_stepper(bulirsch_stoer_dense_out(1E-2, 1E-2, 1, 1, -max_dt), + -max_dt); +} + + +/*********************************************** + * second part of the tests: implicit Rosenbrock + *********************************************** + */ + +typedef boost::numeric::ublas::vector< value_type > vector_type; +typedef boost::numeric::ublas::matrix< value_type > matrix_type; + + +// harmonic oscillator, analytic solution x[0] = sin( t ) +struct osc_rhs +{ + void operator()( const vector_type &x , vector_type &dxdt , const value_type &t ) const + { + dxdt( 0 ) = x( 1 ); + dxdt( 1 ) = -x( 0 ); + } +}; + +struct osc_jacobi +{ + void operator()( const vector_type &x , matrix_type &jacobi , const value_type &t , vector_type &dfdt ) const + { + jacobi( 0 , 0 ) = 0; + jacobi( 0 , 1 ) = 1; + jacobi( 1 , 0 ) = -1; + jacobi( 1 , 1 ) = 0; + dfdt( 0 ) = 0.0; + dfdt( 1 ) = 0.0; + } +}; + + +template +void test_rosenbrock_stepper(Stepper stepper, const double max_dt) +{ + vector_type x( 2 ); + x(0) = x(1) = 10.0; + const size_t steps = 100; std::vector times; - integrate_adaptive(make_dense_output(1E-2, 1E-2, max_dt, runge_kutta_dopri5()), - damped_osc, x, 0.0, 1.0, max_dt, push_back_time(times)); + integrate_adaptive(stepper, + std::make_pair(osc_rhs(), osc_jacobi()), + x, 0.0, steps*max_dt, max_dt, push_back_time(times)); + + BOOST_CHECK_EQUAL(times.size(), steps+1); // check that dt remains at exactly max_dt for( size_t i=0 ; i()), - damped_osc, x, 0.0, 1.0, 0.1, push_back_time(times)); + integrate_adaptive(stepper, + std::make_pair(osc_rhs(), osc_jacobi()), + x, 0.0, steps*max_dt, 10*max_dt, push_back_time(times)); + + BOOST_CHECK_EQUAL(times.size(), steps+1); // check that dt remains at exactly max_dt for( size_t i=0 ; i()), + max_dt); + test_rosenbrock_stepper(make_controlled(1E-2, 1E-2, -max_dt, rosenbrock4()), + -max_dt); +} + + +BOOST_AUTO_TEST_CASE( test_dense_out_rosenbrock ) +{ + const double max_dt = 0.01; + + test_rosenbrock_stepper(make_dense_output(1E-2, 1E-2, max_dt, rosenbrock4()), + max_dt); + test_rosenbrock_stepper(make_dense_output(1E-2, 1E-2, -max_dt, rosenbrock4()), + -max_dt); +} + BOOST_AUTO_TEST_SUITE_END()