| /* |
| boost/numeric/odeint/stepper/detail/controlled_adams_bashforth_moulton.hpp |
| |
| [begin_description] |
| Implemetation of an controlled adams bashforth moulton stepper. |
| [end_description] |
| |
| Copyright 2017 Valentin Noah Hartmann |
| |
| 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_STEPPER_CONTROLLED_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED |
| #define BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLED_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED |
| |
| #include <boost/numeric/odeint/stepper/stepper_categories.hpp> |
| #include <boost/numeric/odeint/stepper/controlled_step_result.hpp> |
| |
| #include <boost/numeric/odeint/stepper/adaptive_adams_bashforth_moulton.hpp> |
| #include <boost/numeric/odeint/stepper/detail/pid_step_adjuster.hpp> |
| |
| #include <boost/numeric/odeint/util/unwrap_reference.hpp> |
| #include <boost/numeric/odeint/util/is_resizeable.hpp> |
| #include <boost/numeric/odeint/util/resizer.hpp> |
| |
| #include <boost/numeric/odeint/util/copy.hpp> |
| #include <boost/numeric/odeint/util/bind.hpp> |
| |
| #include <iostream> |
| |
| namespace boost { |
| namespace numeric { |
| namespace odeint { |
| |
| template< |
| size_t MaxOrder, |
| class State, |
| class Value = double, |
| class Algebra = typename algebra_dispatcher< State >::algebra_type |
| > |
| class default_order_adjuster |
| { |
| public: |
| typedef State state_type; |
| typedef Value value_type; |
| typedef state_wrapper< state_type > wrapped_state_type; |
| |
| typedef Algebra algebra_type; |
| |
| default_order_adjuster( const algebra_type &algebra = algebra_type() ) |
| : m_algebra( algebra ) |
| {}; |
| |
| size_t adjust_order(size_t order, size_t init, boost::array<wrapped_state_type, 4> &xerr) |
| { |
| using std::abs; |
| |
| value_type errc = abs(m_algebra.norm_inf(xerr[2].m_v)); |
| |
| value_type errm1 = 3*errc; |
| value_type errm2 = 3*errc; |
| |
| if(order > 2) |
| { |
| errm2 = abs(m_algebra.norm_inf(xerr[0].m_v)); |
| } |
| if(order >= 2) |
| { |
| errm1 = abs(m_algebra.norm_inf(xerr[1].m_v)); |
| } |
| |
| 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) |
| { |
| value_type errp = abs(m_algebra.norm_inf(xerr[3].m_v)); |
| |
| if(order > 1 && errm1 < errc && errp) |
| { |
| return order-1; |
| } |
| else if(order < MaxOrder && errp < (0.5-0.25*order/MaxOrder) * errc) |
| { |
| return order+1; |
| } |
| } |
| |
| return order; |
| }; |
| private: |
| algebra_type m_algebra; |
| }; |
| |
| template< |
| class ErrorStepper, |
| class StepAdjuster = detail::pid_step_adjuster< typename ErrorStepper::state_type, |
| typename ErrorStepper::value_type, |
| typename ErrorStepper::deriv_type, |
| typename ErrorStepper::time_type, |
| typename ErrorStepper::algebra_type, |
| typename ErrorStepper::operations_type, |
| detail::H211PI |
| >, |
| class OrderAdjuster = default_order_adjuster< ErrorStepper::order_value, |
| typename ErrorStepper::state_type, |
| typename ErrorStepper::value_type, |
| typename ErrorStepper::algebra_type |
| >, |
| class Resizer = initially_resizer |
| > |
| class controlled_adams_bashforth_moulton |
| { |
| public: |
| typedef ErrorStepper stepper_type; |
| |
| static const typename stepper_type::order_type order_value = stepper_type::order_value; |
| |
| typedef typename stepper_type::state_type state_type; |
| typedef typename stepper_type::value_type value_type; |
| typedef typename stepper_type::deriv_type deriv_type; |
| typedef typename stepper_type::time_type time_type; |
| |
| typedef typename stepper_type::algebra_type algebra_type; |
| typedef typename stepper_type::operations_type operations_type; |
| typedef Resizer resizer_type; |
| |
| typedef StepAdjuster step_adjuster_type; |
| typedef OrderAdjuster order_adjuster_type; |
| typedef controlled_stepper_tag stepper_category; |
| |
| 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 , 4 > error_storage_type; |
| |
| typedef typename stepper_type::coeff_type coeff_type; |
| typedef controlled_adams_bashforth_moulton< ErrorStepper , StepAdjuster , OrderAdjuster , Resizer > controlled_stepper_type; |
| |
| controlled_adams_bashforth_moulton(step_adjuster_type step_adjuster = step_adjuster_type()) |
| :m_stepper(), |
| m_dxdt_resizer(), m_xerr_resizer(), m_xnew_resizer(), |
| m_step_adjuster( step_adjuster ), m_order_adjuster() |
| {}; |
| |
| template< class ExplicitStepper, class System > |
| void initialize(ExplicitStepper stepper, System system, state_type &inOut, time_type &t, time_type dt) |
| { |
| m_stepper.initialize(stepper, system, inOut, t, dt); |
| }; |
| |
| template< class System > |
| void initialize(System system, state_type &inOut, time_type &t, time_type dt) |
| { |
| m_stepper.initialize(system, inOut, t, dt); |
| }; |
| |
| template< class ExplicitStepper, class System > |
| void initialize_controlled(ExplicitStepper stepper, System system, state_type &inOut, time_type &t, time_type &dt) |
| { |
| reset(); |
| coeff_type &coeff = m_stepper.coeff(); |
| |
| m_dxdt_resizer.adjust_size( inOut , detail::bind( &controlled_stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); |
| |
| controlled_step_result res = fail; |
| |
| for( size_t i=0 ; i<order_value; ++i ) |
| { |
| do |
| { |
| res = stepper.try_step( system, inOut, t, dt ); |
| } |
| while(res != success); |
| |
| system( inOut , m_dxdt.m_v , t ); |
| |
| coeff.predict(t-dt, dt); |
| coeff.do_step(m_dxdt.m_v); |
| coeff.confirm(); |
| |
| if(coeff.m_eo < order_value) |
| { |
| ++coeff.m_eo; |
| } |
| } |
| } |
| |
| template< class System > |
| controlled_step_result try_step(System system, state_type & inOut, time_type &t, time_type &dt) |
| { |
| m_xnew_resizer.adjust_size( inOut , detail::bind( &controlled_stepper_type::template resize_xnew_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); |
| |
| controlled_step_result res = try_step(system, inOut, t, m_xnew.m_v, dt); |
| |
| if(res == success) |
| { |
| boost::numeric::odeint::copy( m_xnew.m_v , inOut); |
| } |
| |
| return res; |
| }; |
| |
| template< class System > |
| controlled_step_result try_step(System system, const state_type & in, time_type &t, state_type & out, time_type &dt) |
| { |
| 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[2].m_v); |
| |
| coeff_type &coeff = m_stepper.coeff(); |
| |
| time_type dtPrev = dt; |
| 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+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<double>(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<double>(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<double>(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<double>(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 |
| { |
| return fail; |
| } |
| }; |
| |
| void reset() { m_stepper.reset(); }; |
| |
| private: |
| template< class StateType > |
| bool resize_dxdt_impl( const StateType &x ) |
| { |
| return adjust_size_by_resizeability( m_dxdt, x, typename is_resizeable<deriv_type>::type() ); |
| }; |
| template< class StateType > |
| bool resize_xerr_impl( const StateType &x ) |
| { |
| bool resized( false ); |
| |
| for(size_t i=0; i<m_xerr.size(); ++i) |
| { |
| resized |= adjust_size_by_resizeability( m_xerr[i], x, typename is_resizeable<state_type>::type() ); |
| } |
| return resized; |
| }; |
| template< class StateType > |
| bool resize_xnew_impl( const StateType &x ) |
| { |
| return adjust_size_by_resizeability( m_xnew, x, typename is_resizeable<state_type>::type() ); |
| }; |
| |
| stepper_type m_stepper; |
| |
| wrapped_deriv_type m_dxdt; |
| error_storage_type m_xerr; |
| wrapped_state_type m_xnew; |
| |
| resizer_type m_dxdt_resizer; |
| resizer_type m_xerr_resizer; |
| resizer_type m_xnew_resizer; |
| |
| step_adjuster_type m_step_adjuster; |
| order_adjuster_type m_order_adjuster; |
| }; |
| |
| } // odeint |
| } // numeric |
| } // boost |
| |
| #endif |
