odeint/stepper/adaptive_adams_bashforth_moulton.hpp

0001 /*
0002  boost/numeric/odeint/stepper/detail/adaptive_adams_bashforth_moulton.hpp
0003
0004  [begin_description]
0005  Implemetation of an adaptive adams bashforth moulton stepper.
0006  Used as the stepper for the controlled adams bashforth moulton stepper.
0007  [end_description]
0008
0009  Copyright 2017 Valentin Noah Hartmann
0010
0011  Distributed under the Boost Software License, Version 1.0.
0012  (See accompanying file LICENSE_1_0.txt or
0013  copy at http://www.boost.org/LICENSE_1_0.txt)
0014  */
0015
0016 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_ADAPTIVE_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED
0017 #define BOOST_NUMERIC_ODEINT_STEPPER_ADAPTIVE_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED
0018
0019 #include <boost/numeric/odeint/stepper/detail/adaptive_adams_coefficients.hpp>
0020
0021 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
0022 #include <boost/numeric/odeint/util/bind.hpp>
0023 #include <boost/numeric/odeint/util/copy.hpp>
0024
0025 #include <boost/numeric/odeint/algebra/default_operations.hpp>
0026 #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp>
0027 #include <boost/numeric/odeint/algebra/operations_dispatcher.hpp>
0028
0029 #include <boost/numeric/odeint/util/state_wrapper.hpp>
0030 #include <boost/numeric/odeint/util/is_resizeable.hpp>
0031 #include <boost/numeric/odeint/util/resizer.hpp>
0032
0033 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
0034
0035 #include <boost/numeric/odeint/stepper/base/algebra_stepper_base.hpp>
0036 #include <boost/numeric/odeint/stepper/detail/rotating_buffer.hpp>
0037
0038 namespace boost {
0039 namespace numeric {
0040 namespace odeint {
0041
0042 template<
0043 size_t Steps,
0044 class State,
0045 class Value = double,
0046 class Deriv = State,
0047 class Time = Value,
0048 class Algebra = typename algebra_dispatcher< State >::algebra_type ,
0049 class Operations = typename operations_dispatcher< State >::operations_type,
0050 class Resizer = initially_resizer
0051 >
0052 class adaptive_adams_bashforth_moulton: public algebra_stepper_base< Algebra , Operations >
0053 {
0054 public:
0055     static const size_t steps = Steps;
0056
0057     typedef unsigned short order_type;
0058     static const order_type order_value = steps;
0059
0060     typedef State state_type;
0061     typedef Value value_type;
0062     typedef Deriv deriv_type;
0063     typedef Time time_type;
0064
0065     typedef state_wrapper< state_type > wrapped_state_type;
0066     typedef state_wrapper< deriv_type > wrapped_deriv_type;
0067
0068     typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type;
0069     typedef typename algebra_stepper_base_type::algebra_type algebra_type;
0070     typedef typename algebra_stepper_base_type::operations_type operations_type;
0071     typedef Resizer resizer_type;
0072     typedef error_stepper_tag stepper_category;
0073
0074     typedef detail::adaptive_adams_coefficients< Steps , Deriv , Value , Time , Algebra , Operations , Resizer > coeff_type;
0075     typedef adaptive_adams_bashforth_moulton< Steps , State , Value , Deriv , Time , Algebra , Operations , Resizer > stepper_type;
0076
0077     order_type order() const { return order_value; };
0078     order_type stepper_order() const { return order_value + 1; };
0079     order_type error_order() const { return order_value; };
0080
0081     adaptive_adams_bashforth_moulton( const algebra_type &algebra = algebra_type() )
0082     :algebra_stepper_base_type( algebra ), m_coeff(),
0083     m_dxdt_resizer(), m_xnew_resizer(), m_xerr_resizer()
0084     {};
0085
0086     template< class System >
0087     void do_step(System system, state_type &inOut, time_type t, time_type dt )
0088     {
0089         m_xnew_resizer.adjust_size( inOut , detail::bind( &stepper_type::template resize_xnew_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
0090
0091         do_step(system, inOut, t, m_xnew.m_v, dt, m_xerr.m_v);
0092         boost::numeric::odeint::copy( m_xnew.m_v , inOut);
0093     };
0094
0095     template< class System >
0096     void do_step(System system, const state_type &in, time_type t, state_type &out, time_type dt )
0097     {
0098         do_step(system, in, t, out, dt, m_xerr.m_v);
0099     };
0100
0101     template< class System >
0102     void do_step(System system, state_type &inOut, time_type t, time_type dt, state_type &xerr)
0103     {
0104         m_xnew_resizer.adjust_size( inOut , detail::bind( &stepper_type::template resize_xnew_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
0105
0106         do_step(system, inOut, t, m_xnew.m_v, dt, xerr);
0107         boost::numeric::odeint::copy( m_xnew.m_v , inOut);
0108     };
0109
0110     template< class System >
0111     void do_step(System system, const state_type &in, time_type t, state_type &out, time_type dt , state_type &xerr)
0112     {
0113         do_step_impl(system, in, t, out, dt, xerr);
0114
0115         system(out, m_dxdt.m_v, t+dt);
0116         m_coeff.do_step(m_dxdt.m_v);
0117         m_coeff.confirm();
0118
0119         if(m_coeff.m_eo < order_value)
0120         {
0121             m_coeff.m_eo ++;
0122         }
0123     };
0124
0125     template< class ExplicitStepper, class System >
0126     void initialize(ExplicitStepper stepper, System system, state_type &inOut, time_type &t, time_type dt)
0127     {
0128         reset();
0129         dt = dt/static_cast< time_type >(order_value);
0130
0131         m_dxdt_resizer.adjust_size( inOut , detail::bind( &stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
0132
0133         system( inOut , m_dxdt.m_v , t );
0134         for( size_t i=0 ; i<order_value; ++i )
0135         {
0136             stepper.do_step_dxdt_impl( system, inOut, m_dxdt.m_v, t, dt );
0137
0138             system( inOut , m_dxdt.m_v , t + dt);
0139
0140             m_coeff.predict(t, dt);
0141             m_coeff.do_step(m_dxdt.m_v);
0142             m_coeff.confirm();
0143
0144             t += dt;
0145
0146             if(m_coeff.m_eo < order_value)
0147             {
0148                 ++m_coeff.m_eo;
0149             }
0150         }
0151     };
0152
0153     template< class System >
0154     void initialize(System system, state_type &inOut, time_type &t, time_type dt)
0155     {
0156         reset();
0157         dt = dt/static_cast< time_type >(order_value);
0158
0159         for(size_t i=0; i<order_value; ++i)
0160         {
0161             this->do_step(system, inOut, t, dt);
0162             t += dt;
0163         };
0164     };
0165
0166     template< class System >
0167     void do_step_impl(System system, const state_type & in, time_type t, state_type & out, time_type &dt, state_type &xerr)
0168     {
0169         size_t eO = m_coeff.m_eo;
0170
0171         m_xerr_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_xerr_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
0172         m_dxdt_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
0173
0174         m_coeff.predict(t, dt);
0175         if (m_coeff.m_steps_init == 1)
0176         {
0177             system(in, m_dxdt.m_v, t);
0178             m_coeff.do_step(m_dxdt.m_v, 1);
0179         }
0180
0181         boost::numeric::odeint::copy( in , out );
0182         for(size_t i=0; i<eO; ++i)
0183         {
0184             this->m_algebra.for_each3(out, out, m_coeff.phi[1][i].m_v,
0185                 typename Operations::template scale_sum2<double, double>(1.0, dt*m_coeff.g[i]*m_coeff.beta[0][i]));
0186         }
0187
0188         system(out, m_dxdt.m_v, t+dt);
0189         m_coeff.do_step(m_dxdt.m_v);
0190
0191         this->m_algebra.for_each3(out, out, m_coeff.phi[0][eO].m_v,
0192             typename Operations::template scale_sum2<double, double>(1.0, dt*m_coeff.g[eO]));
0193
0194         // error for current order
0195         this->m_algebra.for_each2(xerr, m_coeff.phi[0][eO].m_v,
0196             typename Operations::template scale_sum1<double>(dt*(m_coeff.g[eO])));
0197     };
0198
0199     const coeff_type& coeff() const { return m_coeff; };
0200     coeff_type & coeff() { return m_coeff; };
0201
0202     void reset() { m_coeff.reset(); };
0203     const deriv_type & dxdt() const { return m_dxdt.m_v; };
0204
0205 private:
0206     template< class StateType >
0207     bool resize_dxdt_impl( const StateType &x )
0208     {
0209         return adjust_size_by_resizeability( m_dxdt, x, typename is_resizeable<deriv_type>::type() );
0210     };
0211     template< class StateType >
0212     bool resize_xnew_impl( const StateType &x )
0213     {
0214         return adjust_size_by_resizeability( m_xnew, x, typename is_resizeable<state_type>::type() );
0215     };
0216     template< class StateType >
0217     bool resize_xerr_impl( const StateType &x )
0218     {
0219         return adjust_size_by_resizeability( m_xerr, x, typename is_resizeable<state_type>::type() );
0220     };
0221
0222     coeff_type m_coeff;
0223
0224     resizer_type m_dxdt_resizer;
0225     resizer_type m_xnew_resizer;
0226     resizer_type m_xerr_resizer;
0227
0228     wrapped_deriv_type m_dxdt;
0229     wrapped_state_type m_xnew;
0230     wrapped_state_type m_xerr;
0231 };
0232
0233 } // odeint
0234 } // numeric
0235 } // boost
0236
0237 #endif