| /* |
| boost/numeric/odeint/stepper/detail/adaptive_adams_bashforth_moulton.hpp |
| |
| [begin_description] |
| Implemetation of an adaptive adams bashforth moulton stepper. |
| Used as the stepper for the 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_ADAPTIVE_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED |
| #define BOOST_NUMERIC_ODEINT_STEPPER_ADAPTIVE_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED |
| |
| #include <boost/numeric/odeint/stepper/detail/adaptive_adams_coefficients.hpp> |
| |
| #include <boost/numeric/odeint/util/unwrap_reference.hpp> |
| #include <boost/numeric/odeint/util/bind.hpp> |
| #include <boost/numeric/odeint/util/copy.hpp> |
| |
| #include <boost/numeric/odeint/algebra/default_operations.hpp> |
| #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp> |
| #include <boost/numeric/odeint/algebra/operations_dispatcher.hpp> |
| |
| #include <boost/numeric/odeint/util/state_wrapper.hpp> |
| #include <boost/numeric/odeint/util/is_resizeable.hpp> |
| #include <boost/numeric/odeint/util/resizer.hpp> |
| |
| #include <boost/numeric/odeint/stepper/stepper_categories.hpp> |
| |
| #include <boost/numeric/odeint/stepper/base/algebra_stepper_base.hpp> |
| #include <boost/numeric/odeint/stepper/detail/rotating_buffer.hpp> |
| |
| namespace boost { |
| namespace numeric { |
| namespace odeint { |
| |
| template< |
| size_t Steps, |
| class State, |
| class Value = double, |
| class Deriv = State, |
| class Time = Value, |
| class Algebra = typename algebra_dispatcher< State >::algebra_type , |
| class Operations = typename operations_dispatcher< State >::operations_type, |
| class Resizer = initially_resizer |
| > |
| class adaptive_adams_bashforth_moulton: public algebra_stepper_base< Algebra , Operations > |
| { |
| public: |
| static const size_t steps = Steps; |
| |
| typedef unsigned short order_type; |
| static const order_type order_value = steps; |
| |
| typedef State state_type; |
| typedef Value value_type; |
| typedef Deriv deriv_type; |
| typedef Time time_type; |
| |
| typedef state_wrapper< state_type > wrapped_state_type; |
| typedef state_wrapper< deriv_type > wrapped_deriv_type; |
| |
| typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type; |
| typedef typename algebra_stepper_base_type::algebra_type algebra_type; |
| typedef typename algebra_stepper_base_type::operations_type operations_type; |
| typedef Resizer resizer_type; |
| typedef error_stepper_tag stepper_category; |
| |
| typedef detail::adaptive_adams_coefficients< Steps , Deriv , Value , Time , Algebra , Operations , Resizer > coeff_type; |
| typedef adaptive_adams_bashforth_moulton< Steps , State , Value , Deriv , Time , Algebra , Operations , Resizer > stepper_type; |
| |
| order_type order() const { return order_value; }; |
| order_type stepper_order() const { return order_value + 1; }; |
| order_type error_order() const { return order_value; }; |
| |
| adaptive_adams_bashforth_moulton( const algebra_type &algebra = algebra_type() ) |
| :algebra_stepper_base_type( algebra ), m_coeff(), |
| m_dxdt_resizer(), m_xnew_resizer(), m_xerr_resizer() |
| {}; |
| |
| template< class System > |
| void do_step(System system, state_type &inOut, time_type t, time_type dt ) |
| { |
| 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.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.m_v); |
| }; |
| |
| template< class System > |
| void do_step(System system, state_type &inOut, time_type t, time_type dt, state_type &xerr) |
| { |
| 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, xerr); |
| 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 , state_type &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); |
| m_coeff.confirm(); |
| |
| if(m_coeff.m_eo < order_value) |
| { |
| m_coeff.m_eo ++; |
| } |
| }; |
| |
| template< class ExplicitStepper, class System > |
| void initialize(ExplicitStepper stepper, System system, state_type &inOut, time_type &t, time_type dt) |
| { |
| reset(); |
| dt = dt/static_cast< time_type >(order_value); |
| |
| m_dxdt_resizer.adjust_size( inOut , detail::bind( &stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); |
| |
| system( inOut , m_dxdt.m_v , t ); |
| for( size_t i=0 ; i<order_value; ++i ) |
| { |
| stepper.do_step_dxdt_impl( system, inOut, m_dxdt.m_v, t, dt ); |
| |
| system( inOut , m_dxdt.m_v , t + dt); |
| |
| m_coeff.predict(t, dt); |
| m_coeff.do_step(m_dxdt.m_v); |
| m_coeff.confirm(); |
| |
| t += dt; |
| |
| if(m_coeff.m_eo < order_value) |
| { |
| ++m_coeff.m_eo; |
| } |
| } |
| }; |
| |
| template< class System > |
| void initialize(System system, state_type &inOut, time_type &t, time_type dt) |
| { |
| reset(); |
| dt = dt/static_cast< time_type >(order_value); |
| |
| for(size_t i=0; i<order_value; ++i) |
| { |
| this->do_step(system, inOut, t, dt); |
| t += dt; |
| }; |
| }; |
| |
| template< class System > |
| 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; |
| |
| m_xerr_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_xerr_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); |
| 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 (m_coeff.m_steps_init == 1) |
| { |
| system(in, m_dxdt.m_v, t); |
| m_coeff.do_step(m_dxdt.m_v, 1); |
| } |
| |
| boost::numeric::odeint::copy( in , out ); |
| for(size_t i=0; i<eO; ++i) |
| { |
| this->m_algebra.for_each3(out, out, m_coeff.phi[1][i].m_v, |
| typename Operations::template scale_sum2<double, double>(1.0, dt*m_coeff.g[i]*m_coeff.beta[0][i])); |
| } |
| |
| system(out, m_dxdt.m_v, t+dt); |
| m_coeff.do_step(m_dxdt.m_v); |
| |
| this->m_algebra.for_each3(out, out, m_coeff.phi[0][eO].m_v, |
| typename Operations::template scale_sum2<double, double>(1.0, dt*m_coeff.g[eO])); |
| |
| // error for current order |
| this->m_algebra.for_each2(xerr, m_coeff.phi[0][eO].m_v, |
| typename Operations::template scale_sum1<double>(dt*(m_coeff.g[eO]))); |
| }; |
| |
| const coeff_type& coeff() const { return m_coeff; }; |
| coeff_type & coeff() { return m_coeff; }; |
| |
| void reset() { m_coeff.reset(); }; |
| const deriv_type & dxdt() const { return m_dxdt.m_v; }; |
| |
| 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_xnew_impl( const StateType &x ) |
| { |
| return adjust_size_by_resizeability( m_xnew, x, typename is_resizeable<state_type>::type() ); |
| }; |
| template< class StateType > |
| bool resize_xerr_impl( const StateType &x ) |
| { |
| return adjust_size_by_resizeability( m_xerr, x, typename is_resizeable<state_type>::type() ); |
| }; |
| |
| coeff_type m_coeff; |
| |
| resizer_type m_dxdt_resizer; |
| resizer_type m_xnew_resizer; |
| resizer_type m_xerr_resizer; |
| |
| wrapped_deriv_type m_dxdt; |
| wrapped_state_type m_xnew; |
| wrapped_state_type m_xerr; |
| }; |
| |
| } // odeint |
| } // numeric |
| } // boost |
| |
| #endif |