boost/odeint
2
0
Fork 0
mirror of https://github.com/boostorg/odeint.git synced 2026-01-19 04:22:12 +00:00
Code Issues Projects Releases Wiki Activity

Enable the adaption of the maximal step size of dense output steppers. (#225)

Browse Source 
* Enable the adaption of the maximal step size of dense output steppers.

For efficient simulation of "hybrid" systems the integrator must approach
the sample points where the discrete variables change their value.

(hybrid systems = systems of ODEs which include discrete variables, beeing
internal variables of the system which only change their value at discrete
sample points)

Approaching sample points can be done by adapting the maximal integrator
step size to min(max_step_size, next_sample_point_time - current_time)
before each do_step.

To achive this in odeint for all dense output steppers the following
changes must be done (which does not change the existing API):
- make private members in bulirsch_stoer_dense_out,
  default_step_adjuster, rosenbrock4_controller protected.
- allow std::ref/boost::ref for step_adjuster in controlled_runge_kutta
  and controlled_runge_kutta and for stepper in rosenbrock4_dense_output
  by unwrapping these before use.
This allows to pass the step adjusters by reference to the dense output
steppers which than allows to change the maximal step size (in the step
adjuster) before each call to do_step.

* Added test for a reference controller in the Rosenbrock4 dense output stepper.

* Make in bulirsch_stoer_dense_out only the required m_max_dt member
protected not all.

Extend the test in rosenbrock4.cpp to test that the controller is a
reference and the maximal step size is applied.

* Fixed build with gcc-4.8
This commit is contained in:
Markus Friedrich
2018-02-13 05:13:12 +01:00
committed by Mario Mulansky
parent a393540e17
commit 5dd9519b7b
5 changed files with 65 additions and 20 deletions
Download Patch File Download Diff File Expand all files Collapse all files

7
include/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp
View File

@@ -391,6 +391,11 @@ public:
}
protected:
time_type m_max_dt;
private:
template< class StateInOut , class StateVector >
@@ -667,8 +672,6 @@ private:
default_error_checker< value_type, algebra_type , operations_type > m_error_checker;
modified_midpoint_dense_out< state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > m_midpoint;
time_type m_max_dt;
bool m_control_interpolation;
bool m_last_step_rejected;

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

@@ -163,7 +163,7 @@ public:
time_type get_max_dt() { return m_max_dt; }
private:
protected:
time_type m_max_dt;
};
@@ -411,10 +411,11 @@ public:
template< class System , class StateIn , class DerivIn , class StateOut >
controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , time_type &dt )
{
if( !m_step_adjuster.check_step_size_limit(dt) )
unwrapped_step_adjuster &step_adjuster = m_step_adjuster;
if( !step_adjuster.check_step_size_limit(dt) )
{
// given dt was above step size limit - adjust and return fail;
dt = m_step_adjuster.get_max_dt();
dt = step_adjuster.get_max_dt();
return fail;
}
@@ -428,13 +429,13 @@ public:
if( max_rel_err > 1.0 )
{
// error too big, decrease step size and reject this step
dt = m_step_adjuster.decrease_step(dt, max_rel_err, m_stepper.error_order());
dt = step_adjuster.decrease_step(dt, max_rel_err, m_stepper.error_order());
return fail;
} else
{
// otherwise, increase step size and accept
t += dt;
dt = m_step_adjuster.increase_step(dt, max_rel_err, m_stepper.stepper_order());
dt = step_adjuster.increase_step(dt, max_rel_err, m_stepper.stepper_order());
return success;
}
}
@@ -505,6 +506,7 @@ private:
stepper_type m_stepper;
error_checker_type m_error_checker;
step_adjuster_type m_step_adjuster;
typedef typename unwrap_reference< step_adjuster_type >::type unwrapped_step_adjuster;
resizer_type m_dxdt_resizer;
resizer_type m_xerr_resizer;
@@ -584,7 +586,7 @@ public:
const step_adjuster_type &step_adjuster = step_adjuster_type() ,
const stepper_type &stepper = stepper_type()
)
: m_stepper( stepper ) , m_error_checker( error_checker ) , m_step_adjuster(step_adjuster) ,
: m_stepper( stepper ) , m_error_checker( error_checker ) , m_step_adjuster(step_adjuster) ,
m_first_call( true )
{ }
@@ -751,10 +753,11 @@ public:
controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt_in , time_type &t ,
StateOut &out , DerivOut &dxdt_out , time_type &dt )
{
if( !m_step_adjuster.check_step_size_limit(dt) )
unwrapped_step_adjuster &step_adjuster = m_step_adjuster;
if( !step_adjuster.check_step_size_limit(dt) )
{
// given dt was above step size limit - adjust and return fail;
dt = m_step_adjuster.get_max_dt();
dt = step_adjuster.get_max_dt();
return fail;
}
@@ -770,12 +773,12 @@ public:
if( max_rel_err > 1.0 )
{
// error too big, decrease step size and reject this step
dt = m_step_adjuster.decrease_step(dt, max_rel_err, m_stepper.error_order());
dt = step_adjuster.decrease_step(dt, max_rel_err, m_stepper.error_order());
return fail;
}
// otherwise, increase step size and accept
t += dt;
dt = m_step_adjuster.increase_step(dt, max_rel_err, m_stepper.stepper_order());
dt = step_adjuster.increase_step(dt, max_rel_err, m_stepper.stepper_order());
return success;
}
@@ -903,6 +906,7 @@ private:
stepper_type m_stepper;
error_checker_type m_error_checker;
step_adjuster_type m_step_adjuster;
typedef typename unwrap_reference< step_adjuster_type >::type unwrapped_step_adjuster;
resizer_type m_dxdt_resizer;
resizer_type m_xerr_resizer;

2
include/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp
View File

@@ -200,7 +200,7 @@ public:
private:
protected:
template< class StateIn >
bool resize_m_xerr( const StateIn &x )

19
include/boost/numeric/odeint/stepper/rosenbrock4_dense_output.hpp
View File

@@ -41,7 +41,8 @@ class rosenbrock4_dense_output
public:
typedef ControlledStepper controlled_stepper_type;
typedef typename controlled_stepper_type::stepper_type stepper_type;
typedef typename unwrap_reference< controlled_stepper_type >::type unwrapped_controlled_stepper_type;
typedef typename unwrapped_controlled_stepper_type::stepper_type stepper_type;
typedef typename stepper_type::value_type value_type;
typedef typename stepper_type::state_type state_type;
typedef typename stepper_type::wrapped_state_type wrapped_state_type;
@@ -54,7 +55,7 @@ public:
typedef rosenbrock4_dense_output< ControlledStepper > dense_output_stepper_type;
rosenbrock4_dense_output( const controlled_stepper_type &stepper = controlled_stepper_type() )
: m_stepper( stepper ) ,
: m_stepper( stepper ) ,
m_x1() , m_x2() ,
m_current_state_x1( true ) ,
m_t() , m_t_old() , m_dt()
@@ -75,16 +76,17 @@ public:
template< class System >
std::pair< time_type , time_type > do_step( System system )
{
unwrapped_controlled_stepper_type &stepper = m_stepper;
failed_step_checker fail_checker; // to throw a runtime_error if step size adjustment fails
controlled_step_result res = fail;
m_t_old = m_t;
do
{
res = m_stepper.try_step( system , get_current_state() , m_t , get_old_state() , m_dt );
res = stepper.try_step( system , get_current_state() , m_t , get_old_state() , m_dt );
fail_checker(); // check for overflow of failed steps
}
while( res == fail );
m_stepper.stepper().prepare_dense_output();
stepper.stepper().prepare_dense_output();
this->toggle_current_state();
return std::make_pair( m_t_old , m_t );
}
@@ -96,20 +98,23 @@ public:
template< class StateOut >
void calc_state( time_type t , StateOut &x )
{
m_stepper.stepper().calc_state( t , x , get_old_state() , m_t_old , get_current_state() , m_t );
unwrapped_controlled_stepper_type &stepper = m_stepper;
stepper.stepper().calc_state( t , x , get_old_state() , m_t_old , get_current_state() , m_t );
}
template< class StateOut >
void calc_state( time_type t , const StateOut &x )
{
m_stepper.stepper().calc_state( t , x , get_old_state() , m_t_old , get_current_state() , m_t );
unwrapped_controlled_stepper_type &stepper = m_stepper;
stepper.stepper().calc_state( t , x , get_old_state() , m_t_old , get_current_state() , m_t );
}
template< class StateType >
void adjust_size( const StateType &x )
{
m_stepper.adjust_size( x );
unwrapped_controlled_stepper_type &stepper = m_stepper;
stepper.adjust_size( x );
resize_impl( x );
}

33
test/rosenbrock4.cpp
View File

@@ -127,6 +127,39 @@ BOOST_AUTO_TEST_CASE( test_rosenbrock4_dense_output )
stepper.calc_state( 0.5 * ( tr.first + tr.second ) , x );
}
class rosenbrock4_controller_max_dt_adaptable : public rosenbrock4_controller< rosenbrock4< value_type > >
{
public:
void set_max_dt(value_type max_dt)
{
m_max_dt = max_dt;
}
};
BOOST_AUTO_TEST_CASE( test_rosenbrock4_dense_output_ref )
{
typedef rosenbrock4_dense_output< boost::reference_wrapper< rosenbrock4_controller_max_dt_adaptable > > stepper_type;
rosenbrock4_controller_max_dt_adaptable c_stepper;
stepper_type stepper( boost::ref( c_stepper ) );
typedef stepper_type::state_type state_type;
typedef stepper_type::value_type stepper_value_type;
typedef stepper_type::deriv_type deriv_type;
typedef stepper_type::time_type time_type;
state_type x( 2 );
x( 0 ) = 0.0 ; x(1) = 1.0;
stepper.initialize( x , 0.0 , 0.1 );
std::pair< value_type , value_type > tr = stepper.do_step( std::make_pair( sys() , jacobi() ) );
stepper.calc_state( 0.5 * ( tr.first + tr.second ) , x );
// adapt the maximal step size to a small value (smaller than the step size) and make a step
const double max_dt = 1e-8;
c_stepper.set_max_dt(max_dt);
stepper.do_step( std::make_pair( sys() , jacobi() ) );
// check if the step was done with the requested maximal step size
BOOST_CHECK_CLOSE(max_dt, stepper.current_time_step(), 1e-14);
}
BOOST_AUTO_TEST_CASE( test_rosenbrock4_copy_dense_output )
{
typedef rosenbrock4_controller< rosenbrock4< value_type > > controlled_stepper_type;