1 /* 2 boost/numeric/odeint/stepper/detail/adaptive_adams_bashforth_moulton.hpp 3 4 [begin_description] 5 Implemetation of an adaptive adams bashforth moulton stepper. 6 Used as the stepper for the controlled adams bashforth moulton stepper. 7 [end_description] 8 9 Copyright 2017 Valentin Noah Hartmann 10 11 Distributed under the Boost Software License, Version 1.0. 12 (See accompanying file LICENSE_1_0.txt or 13 copy at http://www.boost.org/LICENSE_1_0.txt) 14 */ 15 16 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_ADAPTIVE_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED 17 #define BOOST_NUMERIC_ODEINT_STEPPER_ADAPTIVE_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED 18 19 #include <boost/numeric/odeint/stepper/detail/adaptive_adams_coefficients.hpp> 20 21 #include <boost/numeric/odeint/util/unwrap_reference.hpp> 22 #include <boost/numeric/odeint/util/bind.hpp> 23 #include <boost/numeric/odeint/util/copy.hpp> 24 25 #include <boost/numeric/odeint/algebra/default_operations.hpp> 26 #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp> 27 #include <boost/numeric/odeint/algebra/operations_dispatcher.hpp> 28 29 #include <boost/numeric/odeint/util/state_wrapper.hpp> 30 #include <boost/numeric/odeint/util/is_resizeable.hpp> 31 #include <boost/numeric/odeint/util/resizer.hpp> 32 33 #include <boost/numeric/odeint/stepper/stepper_categories.hpp> 34 35 #include <boost/numeric/odeint/stepper/base/algebra_stepper_base.hpp> 36 #include <boost/numeric/odeint/stepper/detail/rotating_buffer.hpp> 37 38 namespace boost { 39 namespace numeric { 40 namespace odeint { 41 42 template< 43 size_t Steps, 44 class State, 45 class Value = double, 46 class Deriv = State, 47 class Time = Value, 48 class Algebra = typename algebra_dispatcher< State >::algebra_type , 49 class Operations = typename operations_dispatcher< State >::operations_type, 50 class Resizer = initially_resizer 51 > 52 class adaptive_adams_bashforth_moulton: public algebra_stepper_base< Algebra , Operations > 53 { 54 public: 55 static const size_t steps = Steps; 56 57 typedef unsigned short order_type; 58 static const order_type order_value = steps; 59 60 typedef State state_type; 61 typedef Value value_type; 62 typedef Deriv deriv_type; 63 typedef Time time_type; 64 65 typedef state_wrapper< state_type > wrapped_state_type; 66 typedef state_wrapper< deriv_type > wrapped_deriv_type; 67 68 typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type; 69 typedef typename algebra_stepper_base_type::algebra_type algebra_type; 70 typedef typename algebra_stepper_base_type::operations_type operations_type; 71 typedef Resizer resizer_type; 72 typedef error_stepper_tag stepper_category; 73 74 typedef detail::adaptive_adams_coefficients< Steps , Deriv , Value , Time , Algebra , Operations , Resizer > coeff_type; 75 typedef adaptive_adams_bashforth_moulton< Steps , State , Value , Deriv , Time , Algebra , Operations , Resizer > stepper_type; 76 order() const77 order_type order() const { return order_value; }; stepper_order() const78 order_type stepper_order() const { return order_value + 1; }; error_order() const79 order_type error_order() const { return order_value; }; 80 adaptive_adams_bashforth_moulton(const algebra_type & algebra=algebra_type ())81 adaptive_adams_bashforth_moulton( const algebra_type &algebra = algebra_type() ) 82 :algebra_stepper_base_type( algebra ), m_coeff(), 83 m_dxdt_resizer(), m_xnew_resizer(), m_xerr_resizer() 84 {}; 85 86 template< class System > do_step(System system,state_type & inOut,time_type t,time_type dt)87 void do_step(System system, state_type &inOut, time_type t, time_type dt ) 88 { 89 m_xnew_resizer.adjust_size( inOut , detail::bind( &stepper_type::template resize_xnew_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); 90 91 do_step(system, inOut, t, m_xnew.m_v, dt, m_xerr.m_v); 92 boost::numeric::odeint::copy( m_xnew.m_v , inOut); 93 }; 94 95 template< class System > do_step(System system,const state_type & in,time_type t,state_type & out,time_type dt)96 void do_step(System system, const state_type &in, time_type t, state_type &out, time_type dt ) 97 { 98 do_step(system, in, t, out, dt, m_xerr.m_v); 99 }; 100 101 template< class System > do_step(System system,state_type & inOut,time_type t,time_type dt,state_type & xerr)102 void do_step(System system, state_type &inOut, time_type t, time_type dt, state_type &xerr) 103 { 104 m_xnew_resizer.adjust_size( inOut , detail::bind( &stepper_type::template resize_xnew_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); 105 106 do_step(system, inOut, t, m_xnew.m_v, dt, xerr); 107 boost::numeric::odeint::copy( m_xnew.m_v , inOut); 108 }; 109 110 template< class System > do_step(System system,const state_type & in,time_type t,state_type & out,time_type dt,state_type & xerr)111 void do_step(System system, const state_type &in, time_type t, state_type &out, time_type dt , state_type &xerr) 112 { 113 do_step_impl(system, in, t, out, dt, xerr); 114 115 system(out, m_dxdt.m_v, t+dt); 116 m_coeff.do_step(m_dxdt.m_v); 117 m_coeff.confirm(); 118 119 if(m_coeff.m_eo < order_value) 120 { 121 m_coeff.m_eo ++; 122 } 123 }; 124 125 template< class ExplicitStepper, class System > initialize(ExplicitStepper stepper,System system,state_type & inOut,time_type & t,time_type dt)126 void initialize(ExplicitStepper stepper, System system, state_type &inOut, time_type &t, time_type dt) 127 { 128 reset(); 129 dt = dt/static_cast< time_type >(order_value); 130 131 m_dxdt_resizer.adjust_size( inOut , detail::bind( &stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); 132 133 system( inOut , m_dxdt.m_v , t ); 134 for( size_t i=0 ; i<order_value; ++i ) 135 { 136 stepper.do_step_dxdt_impl( system, inOut, m_dxdt.m_v, t, dt ); 137 138 system( inOut , m_dxdt.m_v , t + dt); 139 140 m_coeff.predict(t, dt); 141 m_coeff.do_step(m_dxdt.m_v); 142 m_coeff.confirm(); 143 144 t += dt; 145 146 if(m_coeff.m_eo < order_value) 147 { 148 ++m_coeff.m_eo; 149 } 150 } 151 }; 152 153 template< class System > initialize(System system,state_type & inOut,time_type & t,time_type dt)154 void initialize(System system, state_type &inOut, time_type &t, time_type dt) 155 { 156 reset(); 157 dt = dt/static_cast< time_type >(order_value); 158 159 for(size_t i=0; i<order_value; ++i) 160 { 161 this->do_step(system, inOut, t, dt); 162 t += dt; 163 }; 164 }; 165 166 template< class System > do_step_impl(System system,const state_type & in,time_type t,state_type & out,time_type & dt,state_type & xerr)167 void do_step_impl(System system, const state_type & in, time_type t, state_type & out, time_type &dt, state_type &xerr) 168 { 169 size_t eO = m_coeff.m_eo; 170 171 m_xerr_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_xerr_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); 172 m_dxdt_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); 173 174 m_coeff.predict(t, dt); 175 if (m_coeff.m_steps_init == 1) 176 { 177 system(in, m_dxdt.m_v, t); 178 m_coeff.do_step(m_dxdt.m_v, 1); 179 } 180 181 boost::numeric::odeint::copy( in , out ); 182 for(size_t i=0; i<eO; ++i) 183 { 184 this->m_algebra.for_each3(out, out, m_coeff.phi[1][i].m_v, 185 typename Operations::template scale_sum2<double, double>(1.0, dt*m_coeff.g[i]*m_coeff.beta[0][i])); 186 } 187 188 system(out, m_dxdt.m_v, t+dt); 189 m_coeff.do_step(m_dxdt.m_v); 190 191 this->m_algebra.for_each3(out, out, m_coeff.phi[0][eO].m_v, 192 typename Operations::template scale_sum2<double, double>(1.0, dt*m_coeff.g[eO])); 193 194 // error for current order 195 this->m_algebra.for_each2(xerr, m_coeff.phi[0][eO].m_v, 196 typename Operations::template scale_sum1<double>(dt*(m_coeff.g[eO]))); 197 }; 198 coeff() const199 const coeff_type& coeff() const { return m_coeff; }; coeff()200 coeff_type & coeff() { return m_coeff; }; 201 reset()202 void reset() { m_coeff.reset(); }; dxdt() const203 const deriv_type & dxdt() const { return m_dxdt.m_v; }; 204 205 private: 206 template< class StateType > resize_dxdt_impl(const StateType & x)207 bool resize_dxdt_impl( const StateType &x ) 208 { 209 return adjust_size_by_resizeability( m_dxdt, x, typename is_resizeable<deriv_type>::type() ); 210 }; 211 template< class StateType > resize_xnew_impl(const StateType & x)212 bool resize_xnew_impl( const StateType &x ) 213 { 214 return adjust_size_by_resizeability( m_xnew, x, typename is_resizeable<state_type>::type() ); 215 }; 216 template< class StateType > resize_xerr_impl(const StateType & x)217 bool resize_xerr_impl( const StateType &x ) 218 { 219 return adjust_size_by_resizeability( m_xerr, x, typename is_resizeable<state_type>::type() ); 220 }; 221 222 coeff_type m_coeff; 223 224 resizer_type m_dxdt_resizer; 225 resizer_type m_xnew_resizer; 226 resizer_type m_xerr_resizer; 227 228 wrapped_deriv_type m_dxdt; 229 wrapped_state_type m_xnew; 230 wrapped_state_type m_xerr; 231 }; 232 233 } // odeint 234 } // numeric 235 } // boost 236 237 #endif