odeint/test/rosenbrock4.cpp
Markus Friedrich 5dd9519b7b Enable the adaption of the maximal step size of dense output steppers. (#225)
* 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
2018-02-12 20:13:12 -08:00

174 lines
5.4 KiB
C++

/*
[auto_generated]
libs/numeric/odeint/test/rosenbrock4.cpp
[begin_description]
This file tests the Rosenbrock 4 stepper and its controller and dense output stepper.
[end_description]
Copyright 2011-2012 Karsten Ahnert
Copyright 2011 Mario Mulansky
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or
copy at http://www.boost.org/LICENSE_1_0.txt)
*/
// disable checked iterator warning for msvc
#include <boost/config.hpp>
#ifdef BOOST_MSVC
#pragma warning(disable:4996)
#endif
#define BOOST_TEST_MODULE odeint_rosenbrock4
#include <utility>
#include <iostream>
#include <boost/test/unit_test.hpp>
#include <boost/numeric/odeint/stepper/rosenbrock4.hpp>
#include <boost/numeric/odeint/stepper/rosenbrock4_controller.hpp>
#include <boost/numeric/odeint/stepper/rosenbrock4_dense_output.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
using namespace boost::unit_test;
using namespace boost::numeric::odeint;
typedef double value_type;
typedef boost::numeric::ublas::vector< value_type > state_type;
typedef boost::numeric::ublas::matrix< value_type > matrix_type;
struct sys
{
void operator()( const state_type &x , state_type &dxdt , const value_type &t ) const
{
dxdt( 0 ) = x( 0 ) + 2 * x( 1 );
dxdt( 1 ) = x( 1 );
}
};
struct jacobi
{
void operator()( const state_type &x , matrix_type &jacobi , const value_type &t , state_type &dfdt ) const
{
jacobi( 0 , 0 ) = 1;
jacobi( 0 , 1 ) = 2;
jacobi( 1 , 0 ) = 0;
jacobi( 1 , 1 ) = 1;
dfdt( 0 ) = 0.0;
dfdt( 1 ) = 0.0;
}
};
BOOST_AUTO_TEST_SUITE( rosenbrock4_test )
BOOST_AUTO_TEST_CASE( test_rosenbrock4_stepper )
{
typedef rosenbrock4< value_type > stepper_type;
stepper_type 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 ) , xerr( 2 );
x(0) = 0.0; x(1) = 1.0;
stepper.do_step( std::make_pair( sys() , jacobi() ) , x , 0.0 , 0.1 , xerr );
stepper.do_step( std::make_pair( sys() , jacobi() ) , x , 0.0 , 0.1 );
// using std::abs;
// value_type eps = 1E-12;
//
// // compare with analytic solution of above system
// BOOST_CHECK_SMALL( abs( x(0) - 20.0/81.0 ) , eps );
// BOOST_CHECK_SMALL( abs( x(1) - 10.0/9.0 ) , eps );
}
BOOST_AUTO_TEST_CASE( test_rosenbrock4_controller )
{
typedef rosenbrock4_controller< rosenbrock4< value_type > > stepper_type;
stepper_type 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;
value_type t = 0.0 , dt = 0.01;
stepper.try_step( std::make_pair( sys() , jacobi() ) , x , t , dt );
}
BOOST_AUTO_TEST_CASE( test_rosenbrock4_dense_output )
{
typedef rosenbrock4_dense_output< rosenbrock4_controller< rosenbrock4< value_type > > > stepper_type;
typedef rosenbrock4_controller< rosenbrock4< value_type > > controlled_stepper_type;
controlled_stepper_type c_stepper;
stepper_type stepper( 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 );
}
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;
typedef rosenbrock4_dense_output< controlled_stepper_type > stepper_type;
controlled_stepper_type c_stepper;
stepper_type stepper( c_stepper );
stepper_type stepper2( stepper );
}
BOOST_AUTO_TEST_SUITE_END()