odeint/stepper/controlled_adams_bashforth_moulton.hpp

0001 /*
0002  boost/numeric/odeint/stepper/detail/controlled_adams_bashforth_moulton.hpp
0003
0004  [begin_description]
0005  Implemetation of an controlled adams bashforth moulton stepper.
0006  [end_description]
0007
0008  Copyright 2017 Valentin Noah Hartmann
0009
0010  Distributed under the Boost Software License, Version 1.0.
0011  (See accompanying file LICENSE_1_0.txt or
0012  copy at http://www.boost.org/LICENSE_1_0.txt)
0013 */
0014
0015 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLED_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED
0016 #define BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLED_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED
0017
0018 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
0019 #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
0020
0021 #include <boost/numeric/odeint/stepper/adaptive_adams_bashforth_moulton.hpp>
0022 #include <boost/numeric/odeint/stepper/detail/pid_step_adjuster.hpp>
0023
0024 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
0025 #include <boost/numeric/odeint/util/is_resizeable.hpp>
0026 #include <boost/numeric/odeint/util/resizer.hpp>
0027
0028 #include <boost/numeric/odeint/util/copy.hpp>
0029 #include <boost/numeric/odeint/util/bind.hpp>
0030
0031 #include <iostream>
0032
0033 namespace boost {
0034 namespace numeric {
0035 namespace odeint {
0036
0037 template<
0038 size_t MaxOrder,
0039 class State,
0040 class Value = double,
0041 class Algebra = typename algebra_dispatcher< State >::algebra_type
0042 >
0043 class default_order_adjuster
0044 {
0045 public:
0046     typedef State state_type;
0047     typedef Value value_type;
0048     typedef state_wrapper< state_type > wrapped_state_type;
0049
0050     typedef Algebra algebra_type;
0051
0052     default_order_adjuster( const algebra_type &algebra = algebra_type() )
0053     : m_algebra( algebra )
0054     {};
0055
0056     size_t adjust_order(size_t order, size_t init, boost::array<wrapped_state_type, 4> &xerr)
0057     {
0058         using std::abs;
0059
0060         value_type errc = abs(m_algebra.norm_inf(xerr[2].m_v));
0061
0062         value_type errm1 = 3*errc;
0063         value_type errm2 = 3*errc;
0064
0065         if(order > 2)
0066         {
0067             errm2 = abs(m_algebra.norm_inf(xerr[0].m_v));
0068         }
0069         if(order >= 2)
0070         {
0071             errm1 = abs(m_algebra.norm_inf(xerr[1].m_v));
0072         }
0073
0074         size_t o_new = order;
0075
0076         if(order == 2 && errm1 <= 0.5*errc)
0077         {
0078             o_new = order - 1;
0079         }
0080         else if(order > 2 && errm2 < errc && errm1 < errc)
0081         {
0082             o_new = order - 1;
0083         }
0084
0085         if(init < order)
0086         {
0087             return order+1;
0088         }
0089         else if(o_new == order - 1)
0090         {
0091             return order-1;
0092         }
0093         else if(order <= MaxOrder)
0094         {
0095             value_type errp = abs(m_algebra.norm_inf(xerr[3].m_v));
0096
0097             if(order > 1 && errm1 < errc && errp)
0098             {
0099                 return order-1;
0100             }
0101             else if(order < MaxOrder && errp < (0.5-0.25*order/MaxOrder) * errc)
0102             {
0103                 return order+1;
0104             }
0105         }
0106
0107         return order;
0108     };
0109 private:
0110     algebra_type m_algebra;
0111 };
0112
0113 template<
0114 class ErrorStepper,
0115 class StepAdjuster = detail::pid_step_adjuster< typename ErrorStepper::state_type,
0116     typename ErrorStepper::value_type,
0117     typename ErrorStepper::deriv_type,
0118     typename ErrorStepper::time_type,
0119     typename ErrorStepper::algebra_type,
0120     typename ErrorStepper::operations_type,
0121     detail::H211PI
0122     >,
0123 class OrderAdjuster = default_order_adjuster< ErrorStepper::order_value,
0124     typename ErrorStepper::state_type,
0125     typename ErrorStepper::value_type,
0126     typename ErrorStepper::algebra_type
0127     >,
0128 class Resizer = initially_resizer
0129 >
0130 class controlled_adams_bashforth_moulton
0131 {
0132 public:
0133     typedef ErrorStepper stepper_type;
0134
0135     static const typename stepper_type::order_type order_value = stepper_type::order_value;
0136
0137     typedef typename stepper_type::state_type state_type;
0138     typedef typename stepper_type::value_type value_type;
0139     typedef typename stepper_type::deriv_type deriv_type;
0140     typedef typename stepper_type::time_type time_type;
0141
0142     typedef typename stepper_type::algebra_type algebra_type;
0143     typedef typename stepper_type::operations_type operations_type;
0144     typedef Resizer resizer_type;
0145
0146     typedef StepAdjuster step_adjuster_type;
0147     typedef OrderAdjuster order_adjuster_type;
0148     typedef controlled_stepper_tag stepper_category;
0149
0150     typedef typename stepper_type::wrapped_state_type wrapped_state_type;
0151     typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type;
0152     typedef boost::array< wrapped_state_type , 4 > error_storage_type;
0153
0154     typedef typename stepper_type::coeff_type coeff_type;
0155     typedef controlled_adams_bashforth_moulton< ErrorStepper , StepAdjuster , OrderAdjuster , Resizer > controlled_stepper_type;
0156
0157     controlled_adams_bashforth_moulton(step_adjuster_type step_adjuster = step_adjuster_type())
0158     :m_stepper(),
0159     m_dxdt_resizer(), m_xerr_resizer(), m_xnew_resizer(),
0160     m_step_adjuster( step_adjuster ), m_order_adjuster()
0161     {};
0162
0163     template< class ExplicitStepper, class System >
0164     void initialize(ExplicitStepper stepper, System system, state_type &inOut, time_type &t, time_type dt)
0165     {
0166         m_stepper.initialize(stepper, system, inOut, t, dt);
0167     };
0168
0169     template< class System >
0170     void initialize(System system, state_type &inOut, time_type &t, time_type dt)
0171     {
0172         m_stepper.initialize(system, inOut, t, dt);
0173     };
0174
0175     template< class ExplicitStepper, class System >
0176     void initialize_controlled(ExplicitStepper stepper, System system, state_type &inOut, time_type &t, time_type &dt)
0177     {
0178         reset();
0179         coeff_type &coeff = m_stepper.coeff();
0180
0181         m_dxdt_resizer.adjust_size( inOut , detail::bind( &controlled_stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
0182
0183         controlled_step_result res = fail;
0184
0185         for( size_t i=0 ; i<order_value; ++i )
0186         {
0187             do
0188             {
0189                 res = stepper.try_step( system, inOut, t, dt );
0190             }
0191             while(res != success);
0192
0193             system( inOut , m_dxdt.m_v , t );
0194
0195             coeff.predict(t-dt, dt);
0196             coeff.do_step(m_dxdt.m_v);
0197             coeff.confirm();
0198
0199             if(coeff.m_eo < order_value)
0200             {
0201                 ++coeff.m_eo;
0202             }
0203         }
0204     }
0205
0206     template< class System >
0207     controlled_step_result try_step(System system, state_type & inOut, time_type &t, time_type &dt)
0208     {
0209         m_xnew_resizer.adjust_size( inOut , detail::bind( &controlled_stepper_type::template resize_xnew_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
0210
0211         controlled_step_result res = try_step(system, inOut, t, m_xnew.m_v, dt);
0212
0213         if(res == success)
0214         {
0215             boost::numeric::odeint::copy( m_xnew.m_v , inOut);
0216         }
0217
0218         return res;
0219     };
0220
0221     template< class System >
0222     controlled_step_result try_step(System system, const state_type & in, time_type &t, state_type & out, time_type &dt)
0223     {
0224         m_xerr_resizer.adjust_size( in , detail::bind( &controlled_stepper_type::template resize_xerr_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
0225         m_dxdt_resizer.adjust_size( in , detail::bind( &controlled_stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
0226
0227         m_stepper.do_step_impl(system, in, t, out, dt, m_xerr[2].m_v);
0228
0229         coeff_type &coeff = m_stepper.coeff();
0230
0231         time_type dtPrev = dt;
0232         dt = m_step_adjuster.adjust_stepsize(coeff.m_eo, dt, m_xerr[2].m_v, out, m_stepper.dxdt() );
0233
0234         if(dt / dtPrev >= step_adjuster_type::threshold())
0235         {
0236             system(out, m_dxdt.m_v, t+dtPrev);
0237
0238             coeff.do_step(m_dxdt.m_v);
0239             coeff.confirm();
0240
0241             t += dtPrev;
0242
0243             size_t eo = coeff.m_eo;
0244
0245             // estimate errors for next step
0246             double factor = 1;
0247             algebra_type m_algebra;
0248
0249             m_algebra.for_each2(m_xerr[2].m_v, coeff.phi[1][eo].m_v,
0250                 typename operations_type::template scale_sum1<double>(factor*dt*(coeff.gs[eo])));
0251
0252             if(eo > 1)
0253             {
0254                 m_algebra.for_each2(m_xerr[1].m_v, coeff.phi[1][eo-1].m_v,
0255                     typename operations_type::template scale_sum1<double>(factor*dt*(coeff.gs[eo-1])));
0256             }
0257             if(eo > 2)
0258             {
0259                 m_algebra.for_each2(m_xerr[0].m_v, coeff.phi[1][eo-2].m_v,
0260                     typename operations_type::template scale_sum1<double>(factor*dt*(coeff.gs[eo-2])));
0261             }
0262             if(eo < order_value && coeff.m_eo < coeff.m_steps_init-1)
0263             {
0264                 m_algebra.for_each2(m_xerr[3].m_v, coeff.phi[1][eo+1].m_v,
0265                     typename operations_type::template scale_sum1<double>(factor*dt*(coeff.gs[eo+1])));
0266             }
0267
0268             // adjust order
0269             coeff.m_eo = m_order_adjuster.adjust_order(coeff.m_eo, coeff.m_steps_init-1, m_xerr);
0270
0271             return success;
0272         }
0273         else
0274         {
0275             return fail;
0276         }
0277     };
0278
0279     void reset() { m_stepper.reset(); };
0280
0281 private:
0282     template< class StateType >
0283     bool resize_dxdt_impl( const StateType &x )
0284     {
0285         return adjust_size_by_resizeability( m_dxdt, x, typename is_resizeable<deriv_type>::type() );
0286     };
0287     template< class StateType >
0288     bool resize_xerr_impl( const StateType &x )
0289     {
0290         bool resized( false );
0291
0292         for(size_t i=0; i<m_xerr.size(); ++i)
0293         {
0294             resized |= adjust_size_by_resizeability( m_xerr[i], x, typename is_resizeable<state_type>::type() );
0295         }
0296         return resized;
0297     };
0298     template< class StateType >
0299     bool resize_xnew_impl( const StateType &x )
0300     {
0301         return adjust_size_by_resizeability( m_xnew, x, typename is_resizeable<state_type>::type() );
0302     };
0303
0304     stepper_type m_stepper;
0305
0306     wrapped_deriv_type m_dxdt;
0307     error_storage_type m_xerr;
0308     wrapped_state_type m_xnew;
0309
0310     resizer_type m_dxdt_resizer;
0311     resizer_type m_xerr_resizer;
0312     resizer_type m_xnew_resizer;
0313
0314     step_adjuster_type m_step_adjuster;
0315     order_adjuster_type m_order_adjuster;
0316 };
0317
0318 } // odeint
0319 } // numeric
0320 } // boost
0321
0322 #endif