1 /* 2 boost/numeric/odeint/stepper/detail/controlled_adams_bashforth_moulton.hpp 3 4 [begin_description] 5 Implemetation of an controlled adams bashforth moulton stepper. 6 [end_description] 7 8 Copyright 2017 Valentin Noah Hartmann 9 10 Distributed under the Boost Software License, Version 1.0. 11 (See accompanying file LICENSE_1_0.txt or 12 copy at http://www.boost.org/LICENSE_1_0.txt) 13 */ 14 15 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLED_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED 16 #define BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLED_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED 17 18 #include <boost/numeric/odeint/stepper/stepper_categories.hpp> 19 #include <boost/numeric/odeint/stepper/controlled_step_result.hpp> 20 21 #include <boost/numeric/odeint/stepper/adaptive_adams_bashforth_moulton.hpp> 22 #include <boost/numeric/odeint/stepper/detail/pid_step_adjuster.hpp> 23 24 #include <boost/numeric/odeint/util/unwrap_reference.hpp> 25 #include <boost/numeric/odeint/util/is_resizeable.hpp> 26 #include <boost/numeric/odeint/util/resizer.hpp> 27 28 #include <boost/numeric/odeint/util/copy.hpp> 29 #include <boost/numeric/odeint/util/bind.hpp> 30 31 #include <iostream> 32 33 namespace boost { 34 namespace numeric { 35 namespace odeint { 36 37 template< 38 size_t MaxOrder, 39 class State, 40 class Value = double, 41 class Algebra = typename algebra_dispatcher< State >::algebra_type 42 > 43 class default_order_adjuster 44 { 45 public: 46 typedef State state_type; 47 typedef Value value_type; 48 typedef state_wrapper< state_type > wrapped_state_type; 49 50 typedef Algebra algebra_type; 51 default_order_adjuster(const algebra_type & algebra=algebra_type ())52 default_order_adjuster( const algebra_type &algebra = algebra_type() ) 53 : m_algebra( algebra ) 54 {}; 55 adjust_order(size_t order,size_t init,boost::array<wrapped_state_type,4> & xerr)56 size_t adjust_order(size_t order, size_t init, boost::array<wrapped_state_type, 4> &xerr) 57 { 58 using std::abs; 59 60 value_type errc = abs(m_algebra.norm_inf(xerr[2].m_v)); 61 62 value_type errm1 = 3*errc; 63 value_type errm2 = 3*errc; 64 65 if(order > 2) 66 { 67 errm2 = abs(m_algebra.norm_inf(xerr[0].m_v)); 68 } 69 if(order >= 2) 70 { 71 errm1 = abs(m_algebra.norm_inf(xerr[1].m_v)); 72 } 73 74 size_t o_new = order; 75 76 if(order == 2 && errm1 <= 0.5*errc) 77 { 78 o_new = order - 1; 79 } 80 else if(order > 2 && errm2 < errc && errm1 < errc) 81 { 82 o_new = order - 1; 83 } 84 85 if(init < order) 86 { 87 return order+1; 88 } 89 else if(o_new == order - 1) 90 { 91 return order-1; 92 } 93 else if(order <= MaxOrder) 94 { 95 value_type errp = abs(m_algebra.norm_inf(xerr[3].m_v)); 96 97 if(order > 1 && errm1 < errc && errp) 98 { 99 return order-1; 100 } 101 else if(order < MaxOrder && errp < (0.5-0.25*order/MaxOrder) * errc) 102 { 103 return order+1; 104 } 105 } 106 107 return order; 108 }; 109 private: 110 algebra_type m_algebra; 111 }; 112 113 template< 114 class ErrorStepper, 115 class StepAdjuster = detail::pid_step_adjuster< typename ErrorStepper::state_type, 116 typename ErrorStepper::value_type, 117 typename ErrorStepper::deriv_type, 118 typename ErrorStepper::time_type, 119 typename ErrorStepper::algebra_type, 120 typename ErrorStepper::operations_type, 121 detail::H211PI 122 >, 123 class OrderAdjuster = default_order_adjuster< ErrorStepper::order_value, 124 typename ErrorStepper::state_type, 125 typename ErrorStepper::value_type, 126 typename ErrorStepper::algebra_type 127 >, 128 class Resizer = initially_resizer 129 > 130 class controlled_adams_bashforth_moulton 131 { 132 public: 133 typedef ErrorStepper stepper_type; 134 135 static const typename stepper_type::order_type order_value = stepper_type::order_value; 136 137 typedef typename stepper_type::state_type state_type; 138 typedef typename stepper_type::value_type value_type; 139 typedef typename stepper_type::deriv_type deriv_type; 140 typedef typename stepper_type::time_type time_type; 141 142 typedef typename stepper_type::algebra_type algebra_type; 143 typedef typename stepper_type::operations_type operations_type; 144 typedef Resizer resizer_type; 145 146 typedef StepAdjuster step_adjuster_type; 147 typedef OrderAdjuster order_adjuster_type; 148 typedef controlled_stepper_tag stepper_category; 149 150 typedef typename stepper_type::wrapped_state_type wrapped_state_type; 151 typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type; 152 typedef boost::array< wrapped_state_type , 4 > error_storage_type; 153 154 typedef typename stepper_type::coeff_type coeff_type; 155 typedef controlled_adams_bashforth_moulton< ErrorStepper , StepAdjuster , OrderAdjuster , Resizer > controlled_stepper_type; 156 controlled_adams_bashforth_moulton(step_adjuster_type step_adjuster=step_adjuster_type ())157 controlled_adams_bashforth_moulton(step_adjuster_type step_adjuster = step_adjuster_type()) 158 :m_stepper(), 159 m_dxdt_resizer(), m_xerr_resizer(), m_xnew_resizer(), 160 m_step_adjuster( step_adjuster ), m_order_adjuster() 161 {}; 162 163 template< class ExplicitStepper, class System > initialize(ExplicitStepper stepper,System system,state_type & inOut,time_type & t,time_type dt)164 void initialize(ExplicitStepper stepper, System system, state_type &inOut, time_type &t, time_type dt) 165 { 166 m_stepper.initialize(stepper, system, inOut, t, dt); 167 }; 168 169 template< class System > initialize(System system,state_type & inOut,time_type & t,time_type dt)170 void initialize(System system, state_type &inOut, time_type &t, time_type dt) 171 { 172 m_stepper.initialize(system, inOut, t, dt); 173 }; 174 175 template< class ExplicitStepper, class System > initialize_controlled(ExplicitStepper stepper,System system,state_type & inOut,time_type & t,time_type & dt)176 void initialize_controlled(ExplicitStepper stepper, System system, state_type &inOut, time_type &t, time_type &dt) 177 { 178 reset(); 179 coeff_type &coeff = m_stepper.coeff(); 180 181 m_dxdt_resizer.adjust_size( inOut , detail::bind( &controlled_stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); 182 183 controlled_step_result res = fail; 184 185 for( size_t i=0 ; i<order_value; ++i ) 186 { 187 do 188 { 189 res = stepper.try_step( system, inOut, t, dt ); 190 } 191 while(res != success); 192 193 system( inOut , m_dxdt.m_v , t ); 194 195 coeff.predict(t-dt, dt); 196 coeff.do_step(m_dxdt.m_v); 197 coeff.confirm(); 198 199 if(coeff.m_eo < order_value) 200 { 201 ++coeff.m_eo; 202 } 203 } 204 } 205 206 template< class System > try_step(System system,state_type & inOut,time_type & t,time_type & dt)207 controlled_step_result try_step(System system, state_type & inOut, time_type &t, time_type &dt) 208 { 209 m_xnew_resizer.adjust_size( inOut , detail::bind( &controlled_stepper_type::template resize_xnew_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); 210 211 controlled_step_result res = try_step(system, inOut, t, m_xnew.m_v, dt); 212 213 if(res == success) 214 { 215 boost::numeric::odeint::copy( m_xnew.m_v , inOut); 216 } 217 218 return res; 219 }; 220 221 template< class System > try_step(System system,const state_type & in,time_type & t,state_type & out,time_type & dt)222 controlled_step_result try_step(System system, const state_type & in, time_type &t, state_type & out, time_type &dt) 223 { 224 m_xerr_resizer.adjust_size( in , detail::bind( &controlled_stepper_type::template resize_xerr_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); 225 m_dxdt_resizer.adjust_size( in , detail::bind( &controlled_stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); 226 227 m_stepper.do_step_impl(system, in, t, out, dt, m_xerr[2].m_v); 228 229 coeff_type &coeff = m_stepper.coeff(); 230 231 time_type dtPrev = dt; 232 dt = m_step_adjuster.adjust_stepsize(coeff.m_eo, dt, m_xerr[2].m_v, out, m_stepper.dxdt() ); 233 234 if(dt / dtPrev >= step_adjuster_type::threshold()) 235 { 236 system(out, m_dxdt.m_v, t+dtPrev); 237 238 coeff.do_step(m_dxdt.m_v); 239 coeff.confirm(); 240 241 t += dtPrev; 242 243 size_t eo = coeff.m_eo; 244 245 // estimate errors for next step 246 double factor = 1; 247 algebra_type m_algebra; 248 249 m_algebra.for_each2(m_xerr[2].m_v, coeff.phi[1][eo].m_v, 250 typename operations_type::template scale_sum1<double>(factor*dt*(coeff.gs[eo]))); 251 252 if(eo > 1) 253 { 254 m_algebra.for_each2(m_xerr[1].m_v, coeff.phi[1][eo-1].m_v, 255 typename operations_type::template scale_sum1<double>(factor*dt*(coeff.gs[eo-1]))); 256 } 257 if(eo > 2) 258 { 259 m_algebra.for_each2(m_xerr[0].m_v, coeff.phi[1][eo-2].m_v, 260 typename operations_type::template scale_sum1<double>(factor*dt*(coeff.gs[eo-2]))); 261 } 262 if(eo < order_value && coeff.m_eo < coeff.m_steps_init-1) 263 { 264 m_algebra.for_each2(m_xerr[3].m_v, coeff.phi[1][eo+1].m_v, 265 typename operations_type::template scale_sum1<double>(factor*dt*(coeff.gs[eo+1]))); 266 } 267 268 // adjust order 269 coeff.m_eo = m_order_adjuster.adjust_order(coeff.m_eo, coeff.m_steps_init-1, m_xerr); 270 271 return success; 272 } 273 else 274 { 275 return fail; 276 } 277 }; 278 reset()279 void reset() { m_stepper.reset(); }; 280 281 private: 282 template< class StateType > resize_dxdt_impl(const StateType & x)283 bool resize_dxdt_impl( const StateType &x ) 284 { 285 return adjust_size_by_resizeability( m_dxdt, x, typename is_resizeable<deriv_type>::type() ); 286 }; 287 template< class StateType > resize_xerr_impl(const StateType & x)288 bool resize_xerr_impl( const StateType &x ) 289 { 290 bool resized( false ); 291 292 for(size_t i=0; i<m_xerr.size(); ++i) 293 { 294 resized |= adjust_size_by_resizeability( m_xerr[i], x, typename is_resizeable<state_type>::type() ); 295 } 296 return resized; 297 }; 298 template< class StateType > resize_xnew_impl(const StateType & x)299 bool resize_xnew_impl( const StateType &x ) 300 { 301 return adjust_size_by_resizeability( m_xnew, x, typename is_resizeable<state_type>::type() ); 302 }; 303 304 stepper_type m_stepper; 305 306 wrapped_deriv_type m_dxdt; 307 error_storage_type m_xerr; 308 wrapped_state_type m_xnew; 309 310 resizer_type m_dxdt_resizer; 311 resizer_type m_xerr_resizer; 312 resizer_type m_xnew_resizer; 313 314 step_adjuster_type m_step_adjuster; 315 order_adjuster_type m_order_adjuster; 316 }; 317 318 } // odeint 319 } // numeric 320 } // boost 321 322 #endif 323