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

add max_dt to rosenbrock steppers

Browse Source 
now all controlled and dense out steppers have step size limitation
functionality.
This commit is contained in:
Mario Mulansky 2015-11-01 14:21:02 -06:00
parent 936af851bb
commit 60cd09818d
3 changed files with 142 additions and 31 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
include/boost/numeric/odeint/stepper/generation/generation_rosenbrock4.hpp
View File

@ -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 ) );
}
};

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

@ -27,6 +27,7 @@
#include <boost/numeric/odeint/util/copy.hpp>
#include <boost/numeric/odeint/util/is_resizeable.hpp>
#include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
#include <boost/numeric/odeint/stepper/rosenbrock4.hpp>
@ -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<time_type>(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<time_type>(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<time_type>(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;

134
test/step_size_limitation.cpp
View File

@ -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 > &times )
: m_times( times ) { }
void operator()( const state_type &x , double t )
template<typename State>
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<class ControlledStepper>
void test_controlled_stepper(ControlledStepper stepper, const double max_dt)
template<class Stepper>
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<double> 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<times.size() ; ++i )
// check if observer was called at times 0,1,2,...
BOOST_CHECK_SMALL( times[i] - static_cast<double>(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<times.size() ; ++i )
// check if observer was called at times 0,1,2,...
@ -132,41 +138,92 @@ BOOST_AUTO_TEST_CASE( test_controlled )
{
const double max_dt = 0.01;
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<state_type>()),
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<state_type>()),
-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<state_type>()),
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<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_explicit_stepper(bulirsch_stoer<state_type>(1E-2, 1E-2, 1.0, 1.0, max_dt),
max_dt);
test_controlled_stepper(bulirsch_stoer<state_type>(1E-2, 1E-2, 1.0, 1.0, -max_dt),
test_explicit_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 )
{
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<state_type>()),
max_dt);
test_explicit_stepper(make_dense_output(1E-2, 1E-2, -max_dt,
runge_kutta_dopri5<state_type>()),
-max_dt);
test_explicit_stepper(bulirsch_stoer_dense_out<state_type>(1E-2, 1E-2, 1, 1, max_dt),
max_dt);
test_explicit_stepper(bulirsch_stoer_dense_out<state_type>(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<class Stepper>
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<double> times;
integrate_adaptive(make_dense_output(1E-2, 1E-2, max_dt, runge_kutta_dopri5<state_type>()),
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<times.size() ; ++i )
// check if observer was called at times 0,1,2,...
@ -174,8 +231,11 @@ BOOST_AUTO_TEST_CASE( test_dense_out )
times.clear();
// 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>()),
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<times.size() ; ++i )
// check if observer was called at times 0,1,2,...
@ -183,4 +243,26 @@ BOOST_AUTO_TEST_CASE( test_dense_out )
times.clear();
}
BOOST_AUTO_TEST_CASE( test_controlled_rosenbrock )
{
const double max_dt = 0.01;
test_rosenbrock_stepper(make_controlled(1E-2, 1E-2, max_dt, rosenbrock4<value_type>()),
max_dt);
test_rosenbrock_stepper(make_controlled(1E-2, 1E-2, -max_dt, rosenbrock4<value_type>()),
-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<value_type>()),
max_dt);
test_rosenbrock_stepper(make_dense_output(1E-2, 1E-2, -max_dt, rosenbrock4<value_type>()),
-max_dt);
}
BOOST_AUTO_TEST_SUITE_END()