diff options
Diffstat (limited to 'boost/numeric/odeint')
7 files changed, 39 insertions, 18 deletions
diff --git a/boost/numeric/odeint/integrate/integrate_n_steps.hpp b/boost/numeric/odeint/integrate/integrate_n_steps.hpp index 7f3a49bddc..5cc8aa0e7b 100644 --- a/boost/numeric/odeint/integrate/integrate_n_steps.hpp +++ b/boost/numeric/odeint/integrate/integrate_n_steps.hpp @@ -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/detail/integrate_n_steps.hpp> +#include <boost/numeric/odeint/integrate/check_adapter.hpp> namespace boost { namespace numeric { diff --git a/boost/numeric/odeint/integrate/max_step_checker.hpp b/boost/numeric/odeint/integrate/max_step_checker.hpp index 6808a57bdb..654c95bc7d 100644 --- a/boost/numeric/odeint/integrate/max_step_checker.hpp +++ b/boost/numeric/odeint/integrate/max_step_checker.hpp @@ -69,7 +69,7 @@ public: if( m_steps++ >= m_max_steps ) { char error_msg[200]; - sprintf(error_msg, "Max number of iterations exceeded (%d).", m_max_steps); + std::sprintf(error_msg, "Max number of iterations exceeded (%d).", m_max_steps); BOOST_THROW_EXCEPTION( no_progress_error(error_msg) ); } } @@ -101,7 +101,7 @@ public: if( m_steps++ >= m_max_steps ) { char error_msg[200]; - sprintf(error_msg, "Max number of iterations exceeded (%d). A new step size was not found.", m_max_steps); + std::sprintf(error_msg, "Max number of iterations exceeded (%d). A new step size was not found.", m_max_steps); BOOST_THROW_EXCEPTION( step_adjustment_error(error_msg) ); } } @@ -111,4 +111,4 @@ public: } // namespace numeric } // namespace boost -#endif
\ No newline at end of file +#endif diff --git a/boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp b/boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp index 2f7cc4c6fb..f3edce1989 100644 --- a/boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp +++ b/boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp @@ -150,6 +150,12 @@ public : } + void reset(void) + { + m_adams_bashforth.reset(); + } + + private: @@ -175,7 +181,7 @@ private: { m_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) ); m_adams_bashforth.do_step( system , in , t , m_x.m_v , dt ); - m_adams_moulton.do_step( system , in , m_x.m_v , t , out , dt , m_adams_bashforth.step_storage() ); + m_adams_moulton.do_step( system , in , m_x.m_v , t+dt , out , dt , m_adams_bashforth.step_storage() ); } else { @@ -293,6 +299,11 @@ private: * \param dt The step size. */ + /** + * \fn adams_bashforth_moulton::reset( void ) + * \brief Resets the internal buffers of the stepper. + */ + } // odeint } // numeric diff --git a/boost/numeric/odeint/stepper/bulirsch_stoer.hpp b/boost/numeric/odeint/stepper/bulirsch_stoer.hpp index 4b908333ba..02c37492b9 100644 --- a/boost/numeric/odeint/stepper/bulirsch_stoer.hpp +++ b/boost/numeric/odeint/stepper/bulirsch_stoer.hpp @@ -98,6 +98,7 @@ public: m_interval_sequence( m_k_max+1 ) , m_coeff( m_k_max+1 ) , m_cost( m_k_max+1 ) , + m_facmin_table( m_k_max+1 ) , m_table( m_k_max ) , STEPFAC1( 0.65 ) , STEPFAC2( 0.94 ) , STEPFAC3( 0.02 ) , STEPFAC4( 4.0 ) , KFAC1( 0.8 ) , KFAC2( 0.9 ) { @@ -112,21 +113,14 @@ public: else m_cost[i] = m_cost[i-1] + m_interval_sequence[i]; m_coeff[i].resize(i); + m_facmin_table[i] = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , static_cast< value_type >(1) / static_cast< value_type >( 2*i+1 ) ); for( size_t k = 0 ; k < i ; ++k ) { const value_type r = static_cast< value_type >( m_interval_sequence[i] ) / static_cast< value_type >( m_interval_sequence[k] ); m_coeff[i][k] = 1.0 / ( r*r - static_cast< value_type >( 1.0 ) ); // coefficients for extrapolation } - - // crude estimate of optimal order - - m_current_k_opt = 4; - /* no calculation because log10 might not exist for value_type! - const value_type logfact( -log10( max BOOST_PREVENT_MACRO_SUBSTITUTION( eps_rel , static_cast< value_type >(1.0E-12) ) ) * 0.6 + 0.5 ); - m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>( 1 ) , min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>( m_k_max-1 ) , logfact )); - */ } - + reset(); } @@ -282,7 +276,7 @@ public: { m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<int>(m_k_max-1) , static_cast<int>(m_current_k_opt)+1 ); new_h = h_opt[k]; - new_h *= m_cost[m_current_k_opt]/m_cost[k]; + new_h *= static_cast<value_type>(m_cost[m_current_k_opt])/static_cast<value_type>(m_cost[k]); } else new_h = h_opt[m_current_k_opt]; break; @@ -344,6 +338,12 @@ public: { m_first = true; m_last_step_rejected = false; + // crude estimate of optimal order + m_current_k_opt = 4; + /* no calculation because log10 might not exist for value_type! + const value_type logfact( -log10( max BOOST_PREVENT_MACRO_SUBSTITUTION( eps_rel , static_cast< value_type >(1.0E-12) ) ) * 0.6 + 0.5 ); + m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>( 1 ) , min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>( m_k_max-1 ) , logfact )); + */ } @@ -417,7 +417,7 @@ private: BOOST_USING_STD_MAX(); using std::pow; value_type expo( 1.0/(2*k+1) ); - value_type facmin = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , expo ); + value_type facmin = m_facmin_table[k]; value_type fac; if (error == 0.0) fac=1.0/facmin; @@ -506,6 +506,7 @@ private: int_vector m_interval_sequence; // stores the successive interval counts value_matrix m_coeff; int_vector m_cost; // costs for interval count + value_vector m_facmin_table; // for precomputed facmin to save pow calls state_table_type m_table; // sequence of states for extrapolation diff --git a/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp b/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp index d876ca3d36..6a1eed15fb 100644 --- a/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp +++ b/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp @@ -110,6 +110,7 @@ public: m_interval_sequence( m_k_max+1 ) , m_coeff( m_k_max+1 ) , m_cost( m_k_max+1 ) , + m_facmin_table( m_k_max+1 ) , m_table( m_k_max ) , m_mp_states( m_k_max+1 ) , m_derivs( m_k_max+1 ) , @@ -125,9 +126,13 @@ public: m_interval_sequence[i] = 2 + 4*i; // 2 6 10 14 ... m_derivs[i].resize( m_interval_sequence[i] ); if( i == 0 ) + { m_cost[i] = m_interval_sequence[i]; - else + } else + { m_cost[i] = m_cost[i-1] + m_interval_sequence[i]; + } + m_facmin_table[i] = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , static_cast< value_type >(1) / static_cast< value_type >( 2*i+1 ) ); m_coeff[i].resize(i); for( size_t k = 0 ; k < i ; ++k ) { @@ -429,7 +434,7 @@ private: using std::pow; value_type expo = static_cast<value_type>(1)/(m_interval_sequence[k-1]); - value_type facmin = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , expo ); + value_type facmin = m_facmin_table[k]; value_type fac; if (error == 0.0) fac = static_cast<value_type>(1)/facmin; @@ -692,6 +697,7 @@ private: int_vector m_interval_sequence; // stores the successive interval counts value_matrix m_coeff; int_vector m_cost; // costs for interval count + value_vector m_facmin_table; // for precomputed facmin to save pow calls state_vector_type m_table; // sequence of states for extrapolation diff --git a/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp b/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp index 8ae627fe1c..aac2b02dca 100644 --- a/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp +++ b/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp @@ -143,7 +143,7 @@ public: error = max BOOST_PREVENT_MACRO_SUBSTITUTION ( static_cast<value_type>( pow( static_cast<value_type>(5.0) , -static_cast<value_type>(stepper_order) ) ) , error); - time_type dt_old = dt; + // time_type dt_old = dt; unused variable warning //error too small - increase dt and keep the evolution and limit scaling factor to 5.0 dt *= static_cast<value_type>(9)/static_cast<value_type>(10) * pow(error, static_cast<value_type>(-1) / stepper_order); diff --git a/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp b/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp index 0e5edd32c0..61d6e51191 100644 --- a/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp +++ b/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp @@ -163,6 +163,8 @@ public: if( m_max_dt != static_cast<time_type>(0) ) { dt = detail::min_abs(m_max_dt, dt_new); + } else { + dt = dt_new; } m_last_rejected = false; return success; |