diff options
Diffstat (limited to 'boost/numeric/odeint/integrate/detail')
4 files changed, 45 insertions, 36 deletions
diff --git a/boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp b/boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp index 743e57709c..7516d44087 100644 --- a/boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp +++ b/boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp @@ -7,7 +7,7 @@ [end_description] Copyright 2011-2013 Karsten Ahnert - Copyright 2011-2012 Mario Mulansky + Copyright 2011-2015 Mario Mulansky Copyright 2012 Christoph Koke Distributed under the Boost Software License, Version 1.0. @@ -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/max_step_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,7 +42,7 @@ 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 > size_t integrate_const( Stepper stepper , System system , State &start_state , Time start_time , Time end_time , Time dt , @@ -74,7 +75,7 @@ size_t integrate_adaptive( /* - * classical integrate adaptive + * integrate adaptive for controlled stepper */ template< class Stepper , class System , class State , class Time , class Observer > size_t integrate_adaptive( @@ -86,8 +87,7 @@ size_t integrate_adaptive( typename odeint::unwrap_reference< Observer >::type &obs = observer; typename odeint::unwrap_reference< Stepper >::type &st = stepper; - const size_t max_attempts = 1000; - const char *error_string = "Integrate adaptive : Maximal number of iterations reached. A step size could not be found."; + failed_step_checker fail_checker; // to throw a runtime_error if step size adjustment fails size_t count = 0; while( less_with_sign( start_time , end_time , dt ) ) { @@ -97,15 +97,14 @@ size_t integrate_adaptive( dt = end_time - start_time; } - size_t trials = 0; controlled_step_result res; do { res = st.try_step( system , start_state , start_time , dt ); - ++trials; + fail_checker(); // check number of failed steps } - while( ( res == fail ) && ( trials < max_attempts ) ); - if( trials == max_attempts ) BOOST_THROW_EXCEPTION( std::overflow_error( error_string ) ); + while( res == fail ); + fail_checker.reset(); // if we reach here, the step was successful -> reset fail checker ++count; } diff --git a/boost/numeric/odeint/integrate/detail/integrate_const.hpp b/boost/numeric/odeint/integrate/detail/integrate_const.hpp index 312564ff16..7a86b32fa6 100644 --- a/boost/numeric/odeint/integrate/detail/integrate_const.hpp +++ b/boost/numeric/odeint/integrate/detail/integrate_const.hpp @@ -6,7 +6,7 @@ integrate const implementation [end_description] - Copyright 2012 Mario Mulansky + Copyright 2012-2015 Mario Mulansky Copyright 2012 Christoph Koke Copyright 2012 Karsten Ahnert @@ -43,12 +43,13 @@ template< class Stepper , class System , class State , class Time , class Observ 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 , stepper_tag ) { + 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) @@ -72,11 +73,11 @@ template< class Stepper , class System , class State , class Time , class Observ 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 , 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 +86,10 @@ 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() ); + // integrate_adaptive_checked uses the given checker to throw if an overflow occurs + real_steps += detail::integrate_adaptive(stepper, system, start_state, time, + static_cast<Time>(time + time_step), dt, + null_observer(), controlled_stepper_tag()); // 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++; @@ -102,12 +105,12 @@ template< class Stepper , class System , class State , class Time , class Observ 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 , dense_output_stepper_tag ) { typename odeint::unwrap_reference< Observer >::type &obs = observer; typename odeint::unwrap_reference< Stepper >::type &st = stepper; - + Time time = start_time; st.initialize( start_state , time , dt ); @@ -117,7 +120,7 @@ size_t integrate_const( int obs_step( 1 ); int real_step( 0 ); - while( less_with_sign( static_cast<Time>(time+dt) , end_time , dt ) ) + while( less_eq_with_sign( static_cast<Time>(time+dt) , end_time , dt ) ) { while( less_eq_with_sign( time , st.current_time() , dt ) ) { @@ -148,6 +151,7 @@ size_t integrate_const( } // last observation, if we are still in observation interval + // might happen due to finite precision problems when computing the the time if( less_eq_with_sign( time , end_time , dt ) ) { st.calc_state( time , start_state ); diff --git a/boost/numeric/odeint/integrate/detail/integrate_n_steps.hpp b/boost/numeric/odeint/integrate/detail/integrate_n_steps.hpp index 3c1d171620..2ef490d592 100644 --- a/boost/numeric/odeint/integrate/detail/integrate_n_steps.hpp +++ b/boost/numeric/odeint/integrate/detail/integrate_n_steps.hpp @@ -6,7 +6,7 @@ integrate steps implementation [end_description] - Copyright 2012 Mario Mulansky + Copyright 2012-2015 Mario Mulansky Copyright 2012 Christoph Koke Copyright 2012 Karsten Ahnert @@ -32,10 +32,10 @@ namespace detail { // forward declaration template< class Stepper , class System , class State , class Time , class Observer > -size_t integrate_adaptive( +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, controlled_stepper_tag ); @@ -66,7 +66,7 @@ Time integrate_n_steps( /* controlled version */ -template< class Stepper , class System , class State , class Time , class Observer> +template< class Stepper , class System , class State , class Time , class Observer > Time integrate_n_steps( Stepper stepper , System system , State &start_state , Time start_time , Time dt , size_t num_of_steps , @@ -80,8 +80,9 @@ Time integrate_n_steps( for( size_t step = 0; step < num_of_steps ; ++step ) { obs( start_state , time ); - detail::integrate_adaptive( stepper , system , start_state , time , static_cast<Time>(time+time_step) , dt , - null_observer() , controlled_stepper_tag() ); + // integrate_adaptive_checked uses the given checker to throw if an overflow occurs + detail::integrate_adaptive(stepper, system, start_state, time, static_cast<Time>(time + time_step), dt, + null_observer(), controlled_stepper_tag()); // direct computation of the time avoids error propagation happening when using time += dt // we need clumsy type analysis to get boost units working here time = start_time + static_cast< typename unit_value_type<Time>::type >(step+1) * time_step; @@ -93,7 +94,7 @@ Time integrate_n_steps( /* dense output version */ -template< class Stepper , class System , class State , class Time , class Observer> +template< class Stepper , class System , class State , class Time , class Observer > Time integrate_n_steps( Stepper stepper , System system , State &start_state , Time start_time , Time dt , size_t num_of_steps , @@ -135,6 +136,7 @@ Time integrate_n_steps( } } + // make sure we really end exactly where we should end while( st.current_time() < end_time ) { if( less_with_sign( end_time , @@ -144,7 +146,7 @@ Time integrate_n_steps( st.do_step( system ); } - // observation at end point, only if we ended exactly on the end-point (or above due to finite precision) + // observation at final point obs( st.current_state() , end_time ); return time; diff --git a/boost/numeric/odeint/integrate/detail/integrate_times.hpp b/boost/numeric/odeint/integrate/detail/integrate_times.hpp index d5446ba590..2e27990412 100644 --- a/boost/numeric/odeint/integrate/detail/integrate_times.hpp +++ b/boost/numeric/odeint/integrate/detail/integrate_times.hpp @@ -6,7 +6,7 @@ Default integrate times implementation. [end_description] - Copyright 2011-2012 Mario Mulansky + Copyright 2011-2015 Mario Mulansky Copyright 2012 Karsten Ahnert Copyright 2012 Christoph Koke @@ -26,6 +26,7 @@ #include <boost/numeric/odeint/util/unwrap_reference.hpp> #include <boost/numeric/odeint/stepper/controlled_step_result.hpp> #include <boost/numeric/odeint/util/detail/less_with_sign.hpp> +#include <boost/numeric/odeint/integrate/max_step_checker.hpp> namespace boost { @@ -38,15 +39,18 @@ namespace detail { /* * integrate_times for simple stepper */ -template< class Stepper , class System , class State , class TimeIterator , class Time , class Observer > +template<class Stepper, class System, class State, class TimeIterator, class Time, class Observer> size_t integrate_times( Stepper stepper , System system , State &start_state , TimeIterator start_time , TimeIterator end_time , Time dt , Observer observer , stepper_tag ) { - typename odeint::unwrap_reference< Observer >::type &obs = observer; - typename odeint::unwrap_reference< Stepper >::type &st = stepper; + typedef typename odeint::unwrap_reference< Stepper >::type stepper_type; + typedef typename odeint::unwrap_reference< Observer >::type observer_type; + + stepper_type &st = stepper; + observer_type &obs = observer; typedef typename unit_value_type<Time>::type time_type; size_t steps = 0; @@ -82,12 +86,10 @@ size_t integrate_times( typename odeint::unwrap_reference< Stepper >::type &st = stepper; typedef typename unit_value_type<Time>::type time_type; - const size_t max_attempts = 1000; - const char *error_string = "Integrate adaptive : Maximal number of iterations reached. A step size could not be found."; + failed_step_checker fail_checker; // to throw a runtime_error if step size adjustment fails size_t steps = 0; while( true ) { - size_t fail_steps = 0; Time current_time = *start_time++; obs( start_state , current_time ); if( start_time == end_time ) @@ -99,15 +101,16 @@ size_t integrate_times( if( st.try_step( system , start_state , current_time , current_dt ) == success ) { ++steps; + // successful step -> reset the fail counter, see #173 + fail_checker.reset(); // continue with the original step size if dt was reduced due to observation dt = max_abs( dt , current_dt ); } else { - ++fail_steps; + fail_checker(); // check for possible overflow of failed steps in step size adjustment dt = current_dt; } - if( fail_steps == max_attempts ) BOOST_THROW_EXCEPTION( std::overflow_error( error_string )); } } return steps; @@ -125,6 +128,7 @@ size_t integrate_times( { typename odeint::unwrap_reference< Observer >::type &obs = observer; typename odeint::unwrap_reference< Stepper >::type &st = stepper; + typedef typename unit_value_type<Time>::type time_type; if( start_time == end_time ) |