diff --git a/include/boost/numeric/odeint/stepper/adaptive_adams_bashforth_moulton.hpp b/include/boost/numeric/odeint/stepper/adaptive_adams_bashforth_moulton.hpp index 7856ddad..5ba319c1 100644 --- a/include/boost/numeric/odeint/stepper/adaptive_adams_bashforth_moulton.hpp +++ b/include/boost/numeric/odeint/stepper/adaptive_adams_bashforth_moulton.hpp @@ -49,7 +49,6 @@ class adaptive_adams_bashforth_moulton: public algebra_stepper_base< Algebra , O typedef state_wrapper< state_type > wrapped_state_type; typedef state_wrapper< deriv_type > wrapped_deriv_type; - typedef boost::array< wrapped_state_type , 3 > error_storage_type; typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type; typedef typename algebra_stepper_base_type::algebra_type algebra_type; @@ -74,14 +73,14 @@ class adaptive_adams_bashforth_moulton: public algebra_stepper_base< Algebra , O { m_xnew_resizer.adjust_size( inOut , detail::bind( &stepper_type::template resize_xnew_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); - do_step(system, inOut, t, m_xnew.m_v, dt, m_xerr[1].m_v); + do_step(system, inOut, t, m_xnew.m_v, dt, m_xerr.m_v); boost::numeric::odeint::copy( m_xnew.m_v , inOut); }; template< class System > void do_step(System system, const state_type &in, time_type t, state_type &out, time_type dt ) { - do_step(system, in, t, out, dt, m_xerr); + do_step(system, in, t, out, dt, m_xerr.m_v); }; template< class System > @@ -89,16 +88,14 @@ class adaptive_adams_bashforth_moulton: public algebra_stepper_base< Algebra , O { m_xnew_resizer.adjust_size( inOut , detail::bind( &stepper_type::template resize_xnew_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); - do_step(system, inOut, t, m_xnew.m_v, dt, m_xerr[1].m_v); + do_step(system, inOut, t, m_xnew.m_v, dt, xerr); boost::numeric::odeint::copy( m_xnew.m_v , inOut); - boost::numeric::odeint::copy( m_xerr[1].m_v , xerr); }; template< class System > void do_step(System system, const state_type &in, time_type t, state_type &out, time_type dt , state_type &xerr) { - do_step_impl(system, in, t, out, dt, m_xerr); - boost::numeric::odeint::copy( m_xerr[1].m_v , xerr); + do_step_impl(system, in, t, out, dt, xerr); system(out, m_dxdt.m_v, t+dt); m_coeff.do_step(m_dxdt.m_v); @@ -142,7 +139,7 @@ class adaptive_adams_bashforth_moulton: public algebra_stepper_base< Algebra , O }; template< class System > - void do_step_impl(System system, const state_type & in, time_type t, state_type & out, time_type &dt, error_storage_type &xerr) + void do_step_impl(System system, const state_type & in, time_type t, state_type & out, time_type &dt, state_type &xerr) { size_t eO = m_coeff.m_eo; @@ -150,7 +147,7 @@ class adaptive_adams_bashforth_moulton: public algebra_stepper_base< Algebra , O m_dxdt_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); m_coeff.predict(t, dt); - if (eO == 1) + if (m_coeff.m_steps_init == 1) { system(in, m_dxdt.m_v, t); m_coeff.do_step(m_dxdt.m_v, 1); @@ -170,22 +167,8 @@ class adaptive_adams_bashforth_moulton: public algebra_stepper_base< Algebra , O typename Operations::template scale_sum2(1.0, dt*m_coeff.g[eO])); // error for current order - this->m_algebra.for_each2(xerr[1].m_v, m_coeff.phi[0][eO].m_v, - typename Operations::template scale_sum1(dt*(m_coeff.g[eO]-m_coeff.g[eO-1]))); - - // error for order below - if (eO > 1) - { - this->m_algebra.for_each2(xerr[0].m_v, m_coeff.phi[0][eO-1].m_v, - typename Operations::template scale_sum1(dt*(m_coeff.g[eO-1]-m_coeff.g[eO-2]))); - } - - // error for order above - if(m_coeff.m_steps_init > eO) - { - this->m_algebra.for_each2(xerr[2].m_v, m_coeff.phi[0][eO+1].m_v, - typename Operations::template scale_sum1(dt*(m_coeff.g[eO+1]-m_coeff.g[eO]))); - } + this->m_algebra.for_each2(xerr, m_coeff.phi[0][eO].m_v, + typename Operations::template scale_sum1(dt*(m_coeff.g[eO]))); }; const coeff_type& coeff() const { return m_coeff; }; @@ -208,13 +191,7 @@ class adaptive_adams_bashforth_moulton: public algebra_stepper_base< Algebra , O template< class StateType > bool resize_xerr_impl( const StateType &x ) { - bool resized( false ); - - for(size_t i=0; i<3; ++i) - { - resized |= adjust_size_by_resizeability( m_xerr[i], x, typename is_resizeable::type() ); - } - return resized; + return adjust_size_by_resizeability( m_xerr, x, typename is_resizeable::type() ); }; coeff_type m_coeff; @@ -225,7 +202,7 @@ class adaptive_adams_bashforth_moulton: public algebra_stepper_base< Algebra , O wrapped_deriv_type m_dxdt; wrapped_state_type m_xnew; - error_storage_type m_xerr; + wrapped_state_type m_xerr; }; } // odeint diff --git a/include/boost/numeric/odeint/stepper/controlled_adams_bashforth_moulton.hpp b/include/boost/numeric/odeint/stepper/controlled_adams_bashforth_moulton.hpp index adc537f5..9372fd77 100644 --- a/include/boost/numeric/odeint/stepper/controlled_adams_bashforth_moulton.hpp +++ b/include/boost/numeric/odeint/stepper/controlled_adams_bashforth_moulton.hpp @@ -39,41 +39,58 @@ class default_order_adjuster : m_algebra( algebra ) {}; - void adjust_order(size_t &order, const boost::array & xerr) + size_t adjust_order(size_t order, size_t init, boost::array &xerr) { using std::abs; - value_type errm = abs(m_algebra.norm_inf(xerr[0].m_v)); - value_type errc = abs(m_algebra.norm_inf(xerr[1].m_v)); - value_type errp = abs(m_algebra.norm_inf(xerr[2].m_v)); + value_type errc = abs(m_algebra.norm_inf(xerr[2].m_v)); - if(order == 1) - { - if(errp < errc) - { - order ++; - } + value_type errm1 = 3*errc; + value_type errm2 = 3*errc; + + if(order > 2) + { + errm2 = abs(m_algebra.norm_inf(xerr[0].m_v)); } - else if(order == MaxOrder) + if(order >= 2) { - if(errm < errc) - { - order --; - } + errm1 = abs(m_algebra.norm_inf(xerr[1].m_v)); } - else + + size_t o_new = order; + + if(order == 2 && errm1 <= 0.5*errc) + { + o_new = order - 1; + } + else if(order > 2 && errm2 < errc && errm1 < errc) + { + o_new = order - 1; + } + + if(init < order) + { + return order+1; + } + else if(o_new == order - 1) + { + return order-1; + } + else if(order <= MaxOrder) { - if(errm < errc && - errm < errp) + value_type errp = abs(m_algebra.norm_inf(xerr[3].m_v)); + + if(order > 1 && errm1 < errc && errp) { - order--; + return order-1; } - else if(errp < errm && - errp < errc) + else if(order < MaxOrder && errp < (0.5-0.25*order/MaxOrder) * errc) { - order++; + return order+1; } } + + return order; }; private: algebra_type m_algebra; @@ -118,7 +135,7 @@ class controlled_adams_bashforth_moulton typedef typename stepper_type::wrapped_state_type wrapped_state_type; typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type; - typedef boost::array< wrapped_state_type , 3 > error_storage_type; + typedef boost::array< wrapped_state_type , 4 > error_storage_type; typedef typename stepper_type::coeff_type coeff_type; typedef controlled_adams_bashforth_moulton< ErrorStepper , StepAdjuster , OrderAdjuster , Resizer > controlled_stepper_type; @@ -162,39 +179,54 @@ class controlled_adams_bashforth_moulton m_xerr_resizer.adjust_size( in , detail::bind( &controlled_stepper_type::template resize_xerr_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); m_dxdt_resizer.adjust_size( in , detail::bind( &controlled_stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); - m_stepper.do_step_impl(system, in, t, out, dt, m_xerr); + m_stepper.do_step_impl(system, in, t, out, dt, m_xerr[2].m_v); coeff_type &coeff = m_stepper.coeff(); time_type dtPrev = dt; - size_t prevOrder = coeff.m_eo; - - if(coeff.m_steps_init > coeff.m_eo) - { - m_order_adjuster.adjust_order(coeff.m_eo, m_xerr); - dt = m_step_adjuster.adjust_stepsize(coeff.m_eo, dt, m_xerr[1 + coeff.m_eo - prevOrder].m_v, in, m_stepper.dxdt()); - } - else - { - if(coeff.m_eo < order_value) - { - coeff.m_eo ++; - } - dt = m_step_adjuster.adjust_stepsize(coeff.m_eo, dt, m_xerr[1].m_v, in, m_stepper.dxdt()); - } + dt = m_step_adjuster.adjust_stepsize(coeff.m_eo, dt, m_xerr[2].m_v, out, m_stepper.dxdt() ); if(dt / dtPrev >= step_adjuster_type::threshold()) { - system(out, m_dxdt.m_v, t+dt); + system(out, m_dxdt.m_v, t+dtPrev); + coeff.do_step(m_dxdt.m_v); coeff.confirm(); t += dtPrev; + + size_t eo = coeff.m_eo; + + // estimate errors for next step + double factor = 1; + algebra_type m_algebra; + + m_algebra.for_each2(m_xerr[2].m_v, coeff.phi[1][eo].m_v, + typename operations_type::template scale_sum1(factor*dt*(coeff.gs[eo]))); + + if(eo > 1) + { + m_algebra.for_each2(m_xerr[1].m_v, coeff.phi[1][eo-1].m_v, + typename operations_type::template scale_sum1(factor*dt*(coeff.gs[eo-1]))); + } + if(eo > 2) + { + m_algebra.for_each2(m_xerr[0].m_v, coeff.phi[1][eo-2].m_v, + typename operations_type::template scale_sum1(factor*dt*(coeff.gs[eo-2]))); + } + if(eo < order_value && coeff.m_eo < coeff.m_steps_init-1) + { + m_algebra.for_each2(m_xerr[3].m_v, coeff.phi[1][eo+1].m_v, + typename operations_type::template scale_sum1(factor*dt*(coeff.gs[eo+1]))); + } + + // adjust order + coeff.m_eo = m_order_adjuster.adjust_order(coeff.m_eo, coeff.m_steps_init-1, m_xerr); + return success; } else { - coeff.m_eo = prevOrder; return fail; } }; @@ -212,7 +244,7 @@ class controlled_adams_bashforth_moulton { bool resized( false ); - for(size_t i=0; i<3; ++i) + for(size_t i=0; i::type() ); } diff --git a/include/boost/numeric/odeint/stepper/detail/adaptive_adams_coefficients.hpp b/include/boost/numeric/odeint/stepper/detail/adaptive_adams_coefficients.hpp index 1a24d857..a68c151f 100644 --- a/include/boost/numeric/odeint/stepper/detail/adaptive_adams_coefficients.hpp +++ b/include/boost/numeric/odeint/stepper/detail/adaptive_adams_coefficients.hpp @@ -42,7 +42,6 @@ class adaptive_adams_coefficients typedef Time time_type; typedef state_wrapper< deriv_type > wrapped_deriv_type; - typedef rotating_buffer< wrapped_deriv_type , steps+1 > step_storage_type; // +1 for moulton typedef rotating_buffer< time_type , steps+1 > time_storage_type; typedef Algebra algebra_type; @@ -58,13 +57,30 @@ class adaptive_adams_coefficients { for (size_t i=0; im_algebra.for_each3(phi[o][i].m_v, phi[o][i-1].m_v, phi[o+1][i-1].m_v, - typename Operations::template scale_sum2(1.0, -beta[o][i-1])); - } + if (o == 0) + { + this->m_algebra.for_each3(phi[o][i].m_v, phi[o][i-1].m_v, phi[o+1][i-1].m_v, + typename Operations::template scale_sum2(1.0, -beta[o][i-1])); + } + else + { + this->m_algebra.for_each2(phi[o][i].m_v, phi[o][i-1].m_v, + typename Operations::template scale_sum1(1.0)); + } + } }; void confirm() @@ -124,7 +149,7 @@ class adaptive_adams_coefficients } }; - void reset() { m_eo = 1; }; + void reset() { m_eo = 1; m_steps_init = 1; }; size_t m_eo; size_t m_steps_init; @@ -132,12 +157,12 @@ class adaptive_adams_coefficients rotating_buffer, 2> beta; // beta[0] = beta(n) rotating_buffer, 3> phi; // phi[0] = phi(n+1) boost::array g; + boost::array gs; private: template< class StateType > bool resize_phi_impl( const StateType &x ) { - bool resized( false ); for(size_t i=0; i<(order_value + 2); ++i) @@ -152,7 +177,8 @@ class adaptive_adams_coefficients size_t m_ns; time_storage_type m_time_storage; - boost::array, order_value + 2> c; + static const size_t c_size = order_value + 2; + boost::array c; algebra_type m_algebra;