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

added overflow exception to integrate_const

Browse Source 
Following the discussion in #173, the integrate_const function now
provide a mechanisms to check for TOO_MUCH_WORK situations where too
many steps are performed without any progress (i.e. observer calls).
Naturally, this only makes sense for controlled steppers or dense
output steppers.
Also, integrate_adaptive functions do not require such functionality as
there observer calls happen at every time step.
Hence, only integrate_n_steps and integrate_times will be adapted shortly
This commit is contained in:
Mario Mulansky 2015-10-06 17:26:18 +02:00
parent 703350b223
commit 23ffb209fa
8 changed files with 369 additions and 107 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
include/boost/numeric/odeint.hpp
View File

@ -71,6 +71,7 @@
#include <boost/numeric/odeint/integrate/integrate_times.hpp>
#include <boost/numeric/odeint/integrate/observer_collection.hpp>
#include <boost/numeric/odeint/integrate/max_step_checker.hpp>
#include <boost/numeric/odeint/stepper/generation.hpp>

34
include/boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp
View File

@ -25,6 +25,7 @@
#include <boost/numeric/odeint/stepper/stepper_categories.hpp>
#include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
#include <boost/numeric/odeint/integrate/null_checker.hpp>
#include <boost/numeric/odeint/integrate/detail/integrate_const.hpp>
#include <boost/numeric/odeint/util/bind.hpp>
#include <boost/numeric/odeint/util/unwrap_reference.hpp>
@ -41,11 +42,11 @@ namespace odeint {
namespace detail {
// forward declaration
template< class Stepper , class System , class State , class Time , class Observer>
template< class Stepper , class System , class State , class Time , class Observer , class StepOverflowChecker >
size_t integrate_const(
Stepper stepper , System system , State &start_state ,
Time start_time , Time end_time , Time dt ,
Observer observer , stepper_tag );
Observer observer , StepOverflowChecker checker , stepper_tag );
/*
* integrate_adaptive for simple stepper is basically an integrate_const + some last step
@ -58,7 +59,7 @@ size_t integrate_adaptive(
)
{
size_t steps = detail::integrate_const( stepper , system , start_state , start_time ,
end_time , dt , observer , stepper_tag() );
end_time , dt , observer , null_checker() , stepper_tag() );
typename odeint::unwrap_reference< Observer >::type &obs = observer;
typename odeint::unwrap_reference< Stepper >::type &st = stepper;
@ -74,17 +75,18 @@ size_t integrate_adaptive(
/*
* classical integrate adaptive
* small helper function that calls a given checker - this is for use in integrate_const for ControlledStepper
*/
template< class Stepper , class System , class State , class Time , class Observer >
size_t integrate_adaptive(
template< class Stepper , class System , class State , class Time , class Observer , class StepOverflowChecker >
size_t integrate_adaptive_checked(
Stepper stepper , System system , State &start_state ,
Time &start_time , Time end_time , Time &dt ,
Observer observer , controlled_stepper_tag
Observer observer , StepOverflowChecker checker
)
{
typename odeint::unwrap_reference< Observer >::type &obs = observer;
typename odeint::unwrap_reference< Stepper >::type &st = stepper;
typename odeint::unwrap_reference< StepOverflowChecker >::type &chk = checker;
const size_t max_attempts = 1000;
const char *error_string = "Integrate adaptive : Maximal number of iterations reached. A step size could not be found.";
@ -102,6 +104,7 @@ size_t integrate_adaptive(
do
{
res = st.try_step( system , start_state , start_time , dt );
chk();
++trials;
}
while( ( res == fail ) && ( trials < max_attempts ) );
@ -114,6 +117,23 @@ size_t integrate_adaptive(
}
/*
* classical integrate adaptive
*/
template< class Stepper , class System , class State , class Time , class Observer >
size_t integrate_adaptive(
Stepper stepper , System system , State &start_state ,
Time &start_time , Time end_time , Time &dt ,
Observer observer , controlled_stepper_tag
)
{
// call the checked version with a null_checker
return integrate_adaptive_checked(stepper, system, start_state, start_time, end_time,
dt, observer, null_checker());
}
/*
* integrate adaptive for dense output steppers
*

36
include/boost/numeric/odeint/integrate/detail/integrate_const.hpp
View File

@ -31,24 +31,28 @@ namespace odeint {
namespace detail {
// forward declaration
template< class Stepper , class System , class State , class Time , class Observer >
size_t integrate_adaptive(
template< class Stepper , class System , class State , class Time , class Observer , class StepOverflowChecker >
size_t integrate_adaptive_checked(
Stepper stepper , System system , State &start_state ,
Time &start_time , Time end_time , Time &dt ,
Observer observer , controlled_stepper_tag
Observer observer , StepOverflowChecker checker
);
template< class Stepper , class System , class State , class Time , class Observer >
template< class Stepper , class System , class State , class Time , class Observer , class StepOverflowChecker >
size_t integrate_const(
Stepper stepper , System system , State &start_state ,
Time start_time , Time end_time , Time dt ,
Observer observer , stepper_tag
Observer observer , StepOverflowChecker /* checker */ , stepper_tag
)
{
// note that for integrate_const with a basic stepper, overflow checks are not required as exactly
// one step is performed between each observer call
typename odeint::unwrap_reference< Observer >::type &obs = observer;
typename odeint::unwrap_reference< Stepper >::type &st = stepper;
Time time = start_time;
int step = 0;
// cast time+dt explicitely in case of expression templates (e.g. multiprecision)
@ -68,15 +72,15 @@ size_t integrate_const(
template< class Stepper , class System , class State , class Time , class Observer >
template< class Stepper , class System , class State , class Time , class Observer , class StepOverflowChecker >
size_t integrate_const(
Stepper stepper , System system , State &start_state ,
Time start_time , Time end_time , Time dt ,
Observer observer , controlled_stepper_tag
Observer observer , StepOverflowChecker checker , controlled_stepper_tag
)
{
typename odeint::unwrap_reference< Observer >::type &obs = observer;
Time time = start_time;
const Time time_step = dt;
int real_steps = 0;
@ -85,8 +89,8 @@ size_t integrate_const(
while( less_eq_with_sign( static_cast<Time>(time+time_step) , end_time , dt ) )
{
obs( start_state , time );
real_steps += detail::integrate_adaptive( stepper , system , start_state , time , time+time_step , dt ,
null_observer() , controlled_stepper_tag() );
real_steps += detail::integrate_adaptive_checked( stepper , system , start_state , time , time+time_step , dt ,
null_observer() , checker );
// direct computation of the time avoids error propagation happening when using time += dt
// we need clumsy type analysis to get boost units working here
step++;
@ -98,16 +102,17 @@ size_t integrate_const(
}
template< class Stepper , class System , class State , class Time , class Observer >
template< class Stepper , class System , class State , class Time , class Observer , class StepOverflowChecker >
size_t integrate_const(
Stepper stepper , System system , State &start_state ,
Time start_time , Time end_time , Time dt ,
Observer observer , dense_output_stepper_tag
Observer observer , StepOverflowChecker checker , dense_output_stepper_tag
)
{
typename odeint::unwrap_reference< Observer >::type &obs = observer;
typename odeint::unwrap_reference< Stepper >::type &st = stepper;
typename odeint::unwrap_reference< StepOverflowChecker >::type &chk = checker;
Time time = start_time;
st.initialize( start_state , time , dt );
@ -123,6 +128,7 @@ size_t integrate_const(
{
st.calc_state( time , start_state );
obs( start_state , time );
chk.reset(); // reset checker at observation points
++obs_step;
// direct computation of the time avoids error propagation happening when using time += dt
// we need clumsy type analysis to get boost units working here
@ -136,6 +142,7 @@ size_t integrate_const(
while( less_eq_with_sign( st.current_time() , time , dt ) )
{
st.do_step( system );
chk();
++real_step;
}
}
@ -143,6 +150,7 @@ size_t integrate_const(
{ // do the last step ending exactly on the end point
st.initialize( st.current_state() , st.current_time() , end_time - st.current_time() );
st.do_step( system );
chk();
++real_step;
}

195
include/boost/numeric/odeint/integrate/integrate_const.hpp
View File

@ -8,7 +8,7 @@
[end_description]
Copyright 2011-2013 Karsten Ahnert
Copyright 2011-2012 Mario Mulansky
Copyright 2011-2015 Mario Mulansky
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or
@ -23,6 +23,7 @@
#include <boost/numeric/odeint/stepper/stepper_categories.hpp>
#include <boost/numeric/odeint/integrate/null_observer.hpp>
#include <boost/numeric/odeint/integrate/null_checker.hpp>
#include <boost/numeric/odeint/integrate/detail/integrate_const.hpp>
#include <boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp>
@ -31,91 +32,107 @@ namespace numeric {
namespace odeint {
/*
* Integrates with constant time step dt.
*/
template< class Stepper , class System , class State , class Time , class Observer >
template<class Stepper, class System, class State, class Time, class Observer, class StepOverflowChecker>
size_t integrate_const(
Stepper stepper , System system , State &start_state ,
Time start_time , Time end_time , Time dt ,
Observer observer
)
{
typedef typename odeint::unwrap_reference< Stepper >::type::stepper_category stepper_category;
Stepper stepper, System system, State &start_state,
Time start_time, Time end_time, Time dt,
Observer observer, StepOverflowChecker checker
) {
typedef typename odeint::unwrap_reference<Stepper>::type::stepper_category stepper_category;
// we want to get as fast as possible to the end
if( boost::is_same< null_observer , Observer >::value )
{
// no overflow checks needed
if (boost::is_same<null_observer, Observer>::value) {
return detail::integrate_adaptive(
stepper , system , start_state ,
start_time , end_time , dt ,
observer , stepper_category() );
stepper, system, start_state,
start_time, end_time, dt,
observer, stepper_category());
}
else {
return detail::integrate_const(stepper, system, start_state,
start_time, end_time, dt,
observer, checker, stepper_category());
}
else
{
return detail::integrate_const( stepper , system , start_state ,
start_time , end_time , dt ,
observer , stepper_category() );
}
}
/**
* \brief Second version to solve the forwarding problem,
* can be called with Boost.Range as start_state.
*/
template< class Stepper , class System , class State , class Time , class Observer >
* \brief Second version to solve the forwarding problem,
* can be called with Boost.Range as start_state.
*/
template<class Stepper, class System, class State, class Time, class Observer, class StepOverflowChecker >
size_t integrate_const(
Stepper stepper , System system , const State &start_state ,
Time start_time , Time end_time , Time dt ,
Observer observer
)
{
typedef typename odeint::unwrap_reference< Stepper >::type::stepper_category stepper_category;
Stepper stepper, System system, const State &start_state,
Time start_time, Time end_time, Time dt,
Observer observer, StepOverflowChecker checker
) {
typedef typename odeint::unwrap_reference<Stepper>::type::stepper_category stepper_category;
// we want to get as fast as possible to the end
if( boost::is_same< null_observer , Observer >::value )
{
// ToDo: forward checker to this call
if (boost::is_same<null_observer, Observer>::value) {
return detail::integrate_adaptive(
stepper , system , start_state ,
start_time , end_time , dt ,
observer , stepper_category() );
stepper, system, start_state,
start_time, end_time, dt,
observer, stepper_category());
}
else
{
return detail::integrate_const( stepper , system , start_state ,
start_time , end_time , dt ,
observer , stepper_category() );
else {
return detail::integrate_const(stepper, system, start_state,
start_time, end_time, dt,
observer, checker, stepper_category());
}
}
/**
* \brief integrate_const without observer calls
*/
template< class Stepper , class System , class State , class Time >
* \brief integrate_const without step overflow checker
*/
template<class Stepper, class System, class State, class Time, class Observer>
size_t integrate_const(
Stepper stepper , System system , State &start_state ,
Time start_time , Time end_time , Time dt
)
{
return integrate_const( stepper , system , start_state , start_time , end_time , dt , null_observer() );
Stepper stepper, System system, State &start_state,
Time start_time, Time end_time, Time dt, Observer observer
) {
return integrate_const(stepper, system, start_state, start_time, end_time, dt, observer,
null_checker());
}
/**
* \brief Second version to solve the forwarding problem,
* can be called with Boost.Range as start_state.
*/
template< class Stepper , class System , class State , class Time >
* \brief Second version to solve the forwarding problem,
* can be called with Boost.Range as start_state.
*/
template<class Stepper, class System, class State, class Time, class Observer>
size_t integrate_const(
Stepper stepper , System system , const State &start_state ,
Time start_time , Time end_time , Time dt
)
{
return integrate_const( stepper , system , start_state , start_time , end_time , dt , null_observer() );
Stepper stepper, System system, const State &start_state,
Time start_time, Time end_time, Time dt, Observer observer
) {
return integrate_const(stepper, system, start_state, start_time, end_time, dt, observer,
null_checker());
}
/**
* \brief integrate_const without observer calls
*/
template<class Stepper, class System, class State, class Time>
size_t integrate_const(
Stepper stepper, System system, State &start_state,
Time start_time, Time end_time, Time dt
) {
return integrate_const(stepper, system, start_state, start_time, end_time, dt, null_observer(),
null_checker());
}
/**
* \brief Second version to solve the forwarding problem,
* can be called with Boost.Range as start_state.
*/
template<class Stepper, class System, class State, class Time>
size_t integrate_const(
Stepper stepper, System system, const State &start_state,
Time start_time, Time end_time, Time dt
) {
return integrate_const(stepper, system, start_state, start_time, end_time, dt, null_observer(),
null_checker());
}
@ -124,30 +141,36 @@ size_t integrate_const(
/********* DOXYGEN *********/
/**
* \fn integrate_const( Stepper stepper , System system , State &start_state , Time start_time , Time end_time , Time dt , Observer observer )
* \brief Integrates the ODE with constant step size.
*
* Integrates the ODE defined by system using the given stepper.
* This method ensures that the observer is called at constant intervals dt.
* If the Stepper is a normal stepper without step size control, dt is also
* used for the numerical scheme. If a ControlledStepper is provided, the
* algorithm might reduce the step size to meet the error bounds, but it is
* ensured that the observer is always called at equidistant time points
* t0 + n*dt. If a DenseOutputStepper is used, the step size also may vary
* and the dense output is used to call the observer at equidistant time
* points.
*
* \param stepper The stepper to be used for numerical integration.
* \param system Function/Functor defining the rhs of the ODE.
* \param start_state The initial condition x0.
* \param start_time The initial time t0.
* \param end_time The final integration time tend.
* \param dt The time step between observer calls, _not_ necessarily the
* time step of the integration.
* \param observer Function/Functor called at equidistant time intervals.
* \return The number of steps performed.
*/
/**
* \fn integrate_const( Stepper stepper , System system , State &start_state , Time start_time ,
* Time end_time , Time dt , Observer observer , StepOverflowChecker checker )
* \brief Integrates the ODE with constant step size.
*
* Integrates the ODE defined by system using the given stepper.
* This method ensures that the observer is called at constant intervals dt.
* If the Stepper is a normal stepper without step size control, dt is also
* used for the numerical scheme. If a ControlledStepper is provided, the
* algorithm might reduce the step size to meet the error bounds, but it is
* ensured that the observer is always called at equidistant time points
* t0 + n*dt. If a DenseOutputStepper is used, the step size also may vary
* and the dense output is used to call the observer at equidistant time
* points.
* If a max_step_checker is provided as StepOverflowChecker, an exception is
* thrown if too many steps (default: 500) are performed without progress,
* i.e. in between observer calls. If no checker is provided, no such
* overflow check is performed.
*
* \param stepper The stepper to be used for numerical integration.
* \param system Function/Functor defining the rhs of the ODE.
* \param start_state The initial condition x0.
* \param start_time The initial time t0.
* \param end_time The final integration time tend.
* \param dt The time step between observer calls, _not_ necessarily the
* time step of the integration.
* \param observer [optional] Function/Functor called at equidistant time intervals.
* \param checker [optional] Functor to check for step count overflows.
* \return The number of steps performed.
*/
} // namespace odeint
} // namespace numeric

64
include/boost/numeric/odeint/integrate/max_step_checker.hpp Normal file
View File

@ -0,0 +1,64 @@
/*
[auto_generated]
boost/numeric/odeint/integrate/max_step_checker.hpp
[begin_description]
Throws exception if too many steps are performed.
[end_description]
Copyright 2015 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)
*/
#ifndef BOOST_NUMERIC_ODEINT_INTEGRATE_MAX_STEP_CHECKER_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_INTEGRATE_MAX_STEP_CHECKER_HPP_INCLUDED
#include <stdexcept>
#include <boost/throw_exception.hpp>
namespace boost {
namespace numeric {
namespace odeint {
/**
* \brief A class for performing overflow checks on the step count in integrate functions.
*
* Provide an instance of this class to integrate functions if you want to throw an exception if
* too many steps are performed without progress during the integrate routine.
*/
struct max_step_checker{
const int m_max_steps;
int m_steps;
/**
* \brief Construct the max_step_checker.
*/
max_step_checker(const int max_steps = 500)
: m_max_steps(max_steps)
{
reset();
}
void reset()
{
m_steps = 0;
}
void operator()(void)
{
if( m_steps++ >= m_max_steps )
BOOST_THROW_EXCEPTION( std::overflow_error( "too many steps!") );
}
};
} // namespace odeint
} // namespace numeric
} // namespace boost
#endif

37
include/boost/numeric/odeint/integrate/null_checker.hpp Normal file
View File

@ -0,0 +1,37 @@
/*
[auto_generated]
boost/numeric/odeint/integrate/null_checker.hpp
[begin_description]
null_checker
[end_description]
Copyright 2015 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)
*/
#ifndef BOOST_NUMERIC_ODEINT_INTEGRATE_NULL_CHECKER_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_INTEGRATE_NULL_CHECKER_HPP_INCLUDED
namespace boost {
namespace numeric {
namespace odeint {
struct null_checker
{
void operator()( void ) const
{ }
void reset( void ) const
{ }
};
} // namespace odeint
} // namespace numeric
} // namespace boost
#endif // BOOST_NUMERIC_ODEINT_INTEGRATE_NULL_CHECKER_HPP_INCLUDED

1
test/Jamfile.v2
View File

@ -78,6 +78,7 @@ test-suite "odeint"
[ run n_step_time_iterator.cpp ]
[ run times_iterator.cpp ]
[ run times_time_iterator.cpp ]
[ run integrate_overflow.cpp ]
[ compile unwrap_boost_reference.cpp ]
[ compile unwrap_reference.cpp : <cxxflags>-std=c++0x : unwrap_reference_C++11 ]
[ compile-fail unwrap_reference.cpp : <cxxflags>-std=c++98 : unwrap_reference_C++98 ]

108
test/integrate_overflow.cpp Normal file
View File

@ -0,0 +1,108 @@
/*
[auto_generated]
libs/numeric/odeint/test/integrate_overflow.cpp
[begin_description]
This file tests the overflow exception of the integrate functions.
[end_description]
Copyright 2015 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)
*/
#define BOOST_TEST_MODULE odeint_integrate_overflow
#include <boost/test/unit_test.hpp>
#include <utility>
#include <iostream>
#include <vector>
#include <vector>
#include <cmath>
#include <iostream>
#include <boost/numeric/odeint.hpp>
using namespace boost::numeric::odeint;
typedef double value_type;
typedef std::vector< value_type > state_type;
// the famous lorenz system as usual
void lorenz( const state_type &x , state_type &dxdt , const value_type t )
{
//const value_type sigma( 10.0 );
const value_type R( 28.0 );
const value_type b( value_type( 8.0 ) / value_type( 3.0 ) );
// first component trivial
dxdt[0] = 1.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];
}
struct push_back_time
{
std::vector< double >& m_times;
push_back_time( std::vector< double > &times )
: m_times( times ) { }
void operator()( const state_type &x , double t )
{
m_times.push_back( t );
}
};
typedef runge_kutta_dopri5<state_type> stepper_type;
BOOST_AUTO_TEST_SUITE( integrate_overflow )
BOOST_AUTO_TEST_CASE( test_integrate_const )
{
state_type x(3, 5.0);
std::vector<double> times;
// check the function signatures with normal stepper
integrate_const(stepper_type(), lorenz, x, 0.0, 10.0, 1.0);
integrate_const(stepper_type(), lorenz, x, 0.0, 10.0, 1.0, null_observer());
integrate_const(stepper_type(), lorenz, x, 0.0, 10.0, 1.0, null_observer(), null_checker());
// no exceptions expected for normal steppers
integrate_const(stepper_type(), lorenz, x, 0.0, 10.0, 1.0, null_observer(), max_step_checker());
// check exceptions for controlled stepper
integrate_const(stepper_type(), lorenz, x, 0.0, 10.0, 1.0);
integrate_const(make_controlled<stepper_type>(1E-5, 1E-5), lorenz, x, 0.0, 10.0, 1.0);
// very small error terms -> standard overflow threshold of 500 should fire an exception
BOOST_CHECK_THROW(integrate_const(make_controlled<stepper_type>(1E-15, 1E-15), lorenz, x,
0.0, 10.0, 1.0, push_back_time(times), max_step_checker()),
std::overflow_error);
// small threshold of 10 -> larger error still gives an exception
BOOST_CHECK_THROW(integrate_const(make_controlled<stepper_type>(1E-5, 1E-5), lorenz, x,
0.0, 10.0, 1.0, push_back_time(times), max_step_checker(10)),
std::overflow_error);
// check exceptions for dense output stepper
integrate_const(make_dense_output<stepper_type>(1E-5, 1E-5), lorenz, x, 0.0, 10.0, 1.0);
integrate_const(make_dense_output<stepper_type>(1E-5, 1E-5), lorenz, x, 0.0, 10.0, 1.0);
// very small error terms -> standard overflow threshold of 500 should fire an exception
BOOST_CHECK_THROW(integrate_const(make_dense_output<stepper_type>(1E-15, 1E-15), lorenz, x,
0.0, 10.0, 1.0, push_back_time(times), max_step_checker()),
std::overflow_error);
// small threshold of 10 -> larger error still gives an exception
BOOST_CHECK_THROW(integrate_const(make_dense_output<stepper_type>(1E-5, 1E-5), lorenz, x,
0.0, 10.0, 1.0, push_back_time(times), max_step_checker(10)),
std::overflow_error);
}
BOOST_AUTO_TEST_SUITE_END()