1 /* 2 [auto_generated] 3 boost/numeric/odeint/stepper/runge_kutta_cash_karp54_classic.hpp 4 5 [begin_description] 6 Classical implementation of the Runge-Kutta Cash-Karp 5(4) method. 7 [end_description] 8 9 Copyright 2010-2013 Mario Mulansky 10 Copyright 2010-2013 Karsten Ahnert 11 Copyright 2012 Christoph Koke 12 13 Distributed under the Boost Software License, Version 1.0. 14 (See accompanying file LICENSE_1_0.txt or 15 copy at http://www.boost.org/LICENSE_1_0.txt) 16 */ 17 18 19 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_RUNGE_KUTTA_CASH_KARP54_CLASSIC_HPP_INCLUDED 20 #define BOOST_NUMERIC_ODEINT_STEPPER_RUNGE_KUTTA_CASH_KARP54_CLASSIC_HPP_INCLUDED 21 22 23 #include <boost/numeric/odeint/util/bind.hpp> 24 25 #include <boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp> 26 #include <boost/numeric/odeint/algebra/range_algebra.hpp> 27 #include <boost/numeric/odeint/algebra/default_operations.hpp> 28 #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp> 29 #include <boost/numeric/odeint/algebra/operations_dispatcher.hpp> 30 #include <boost/numeric/odeint/stepper/stepper_categories.hpp> 31 #include <boost/numeric/odeint/util/state_wrapper.hpp> 32 #include <boost/numeric/odeint/util/is_resizeable.hpp> 33 #include <boost/numeric/odeint/util/resizer.hpp> 34 35 namespace boost { 36 namespace numeric { 37 namespace odeint { 38 39 40 41 42 template< 43 class State , 44 class Value = double , 45 class Deriv = State , 46 class Time = Value , 47 class Algebra = typename algebra_dispatcher< State >::algebra_type , 48 class Operations = typename operations_dispatcher< State >::operations_type , 49 class Resizer = initially_resizer 50 > 51 #ifndef DOXYGEN_SKIP 52 class runge_kutta_cash_karp54_classic 53 : public explicit_error_stepper_base< 54 runge_kutta_cash_karp54_classic< State , Value , Deriv , Time , Algebra , Operations , Resizer > , 55 5 , 5 , 4 , State , Value , Deriv , Time , Algebra , Operations , Resizer > 56 #else 57 class runge_kutta_cash_karp54_classic : public explicit_error_stepper_base 58 #endif 59 { 60 61 62 public : 63 64 #ifndef DOXYGEN_SKIP 65 typedef explicit_error_stepper_base< 66 runge_kutta_cash_karp54_classic< State , Value , Deriv , Time , Algebra , Operations , Resizer > , 67 5 , 5 , 4 , State , Value , Deriv , Time , Algebra , Operations , Resizer > stepper_base_type; 68 #else 69 typedef explicit_error_stepper_base< runge_kutta_cash_karp54_classic< ... > , ... > stepper_base_type; 70 #endif 71 72 typedef typename stepper_base_type::state_type state_type; 73 typedef typename stepper_base_type::value_type value_type; 74 typedef typename stepper_base_type::deriv_type deriv_type; 75 typedef typename stepper_base_type::time_type time_type; 76 typedef typename stepper_base_type::algebra_type algebra_type; 77 typedef typename stepper_base_type::operations_type operations_type; 78 typedef typename stepper_base_type::resizer_type resizer_type; 79 80 #ifndef DOXYGEN_SKIP 81 typedef typename stepper_base_type::wrapped_state_type wrapped_state_type; 82 typedef typename stepper_base_type::wrapped_deriv_type wrapped_deriv_type; 83 typedef typename stepper_base_type::stepper_type stepper_type; 84 #endif 85 86 runge_kutta_cash_karp54_classic(const algebra_type & algebra=algebra_type ())87 runge_kutta_cash_karp54_classic( const algebra_type &algebra = algebra_type() ) : stepper_base_type( algebra ) 88 { } 89 90 91 92 template< class System , class StateIn , class DerivIn , class StateOut , class Err > do_step_impl(System system,const StateIn & in,const DerivIn & dxdt,time_type t,StateOut & out,time_type dt,Err & xerr)93 void do_step_impl( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt , Err &xerr ) 94 { 95 const value_type c1 = static_cast<value_type> ( 37 ) / static_cast<value_type>( 378 ); 96 const value_type c3 = static_cast<value_type> ( 250 ) / static_cast<value_type>( 621 ); 97 const value_type c4 = static_cast<value_type> ( 125 ) / static_cast<value_type>( 594 ); 98 const value_type c6 = static_cast<value_type> ( 512 ) / static_cast<value_type>( 1771 ); 99 100 const value_type dc1 = c1 - static_cast<value_type> ( 2825 ) / static_cast<value_type>( 27648 ); 101 const value_type dc3 = c3 - static_cast<value_type> ( 18575 ) / static_cast<value_type>( 48384 ); 102 const value_type dc4 = c4 - static_cast<value_type> ( 13525 ) / static_cast<value_type>( 55296 ); 103 const value_type dc5 = static_cast<value_type> ( -277 ) / static_cast<value_type>( 14336 ); 104 const value_type dc6 = c6 - static_cast<value_type> ( 1 ) / static_cast<value_type> ( 4 ); 105 106 do_step_impl( system , in , dxdt , t , out , dt ); 107 108 //error estimate 109 stepper_base_type::m_algebra.for_each6( xerr , dxdt , m_k3.m_v , m_k4.m_v , m_k5.m_v , m_k6.m_v , 110 typename operations_type::template scale_sum5< time_type , time_type , time_type , time_type , time_type >( dt*dc1 , dt*dc3 , dt*dc4 , dt*dc5 , dt*dc6 )); 111 112 } 113 114 115 116 template< class System , class StateIn , class DerivIn , class StateOut > do_step_impl(System system,const StateIn & in,const DerivIn & dxdt,time_type t,StateOut & out,time_type dt)117 void do_step_impl( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt ) 118 { 119 const value_type a2 = static_cast<value_type> ( 1 ) / static_cast<value_type> ( 5 ); 120 const value_type a3 = static_cast<value_type> ( 3 ) / static_cast<value_type> ( 10 ); 121 const value_type a4 = static_cast<value_type> ( 3 ) / static_cast<value_type> ( 5 ); 122 const value_type a5 = static_cast<value_type> ( 1 ); 123 const value_type a6 = static_cast<value_type> ( 7 ) / static_cast<value_type> ( 8 ); 124 125 const value_type b21 = static_cast<value_type> ( 1 ) / static_cast<value_type> ( 5 ); 126 const value_type b31 = static_cast<value_type> ( 3 ) / static_cast<value_type>( 40 ); 127 const value_type b32 = static_cast<value_type> ( 9 ) / static_cast<value_type>( 40 ); 128 const value_type b41 = static_cast<value_type> ( 3 ) / static_cast<value_type> ( 10 ); 129 const value_type b42 = static_cast<value_type> ( -9 ) / static_cast<value_type> ( 10 ); 130 const value_type b43 = static_cast<value_type> ( 6 ) / static_cast<value_type> ( 5 ); 131 const value_type b51 = static_cast<value_type> ( -11 ) / static_cast<value_type>( 54 ); 132 const value_type b52 = static_cast<value_type> ( 5 ) / static_cast<value_type> ( 2 ); 133 const value_type b53 = static_cast<value_type> ( -70 ) / static_cast<value_type>( 27 ); 134 const value_type b54 = static_cast<value_type> ( 35 ) / static_cast<value_type>( 27 ); 135 const value_type b61 = static_cast<value_type> ( 1631 ) / static_cast<value_type>( 55296 ); 136 const value_type b62 = static_cast<value_type> ( 175 ) / static_cast<value_type>( 512 ); 137 const value_type b63 = static_cast<value_type> ( 575 ) / static_cast<value_type>( 13824 ); 138 const value_type b64 = static_cast<value_type> ( 44275 ) / static_cast<value_type>( 110592 ); 139 const value_type b65 = static_cast<value_type> ( 253 ) / static_cast<value_type>( 4096 ); 140 141 const value_type c1 = static_cast<value_type> ( 37 ) / static_cast<value_type>( 378 ); 142 const value_type c3 = static_cast<value_type> ( 250 ) / static_cast<value_type>( 621 ); 143 const value_type c4 = static_cast<value_type> ( 125 ) / static_cast<value_type>( 594 ); 144 const value_type c6 = static_cast<value_type> ( 512 ) / static_cast<value_type>( 1771 ); 145 146 typename odeint::unwrap_reference< System >::type &sys = system; 147 148 m_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_impl<StateIn> , detail::ref( *this ) , detail::_1 ) ); 149 150 //m_x1 = x + dt*b21*dxdt 151 stepper_base_type::m_algebra.for_each3( m_x_tmp.m_v , in , dxdt , 152 typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , dt*b21 ) ); 153 154 sys( m_x_tmp.m_v , m_k2.m_v , t + dt*a2 ); 155 // m_x_tmp = x + dt*b31*dxdt + dt*b32*m_x2 156 stepper_base_type::m_algebra.for_each4( m_x_tmp.m_v , in , dxdt , m_k2.m_v , 157 typename operations_type::template scale_sum3< value_type , time_type , time_type >( 1.0 , dt*b31 , dt*b32 )); 158 159 sys( m_x_tmp.m_v , m_k3.m_v , t + dt*a3 ); 160 // m_x_tmp = x + dt * (b41*dxdt + b42*m_x2 + b43*m_x3) 161 stepper_base_type::m_algebra.for_each5( m_x_tmp.m_v , in , dxdt , m_k2.m_v , m_k3.m_v , 162 typename operations_type::template scale_sum4< value_type , time_type , time_type , time_type >( 1.0 , dt*b41 , dt*b42 , dt*b43 )); 163 164 sys( m_x_tmp.m_v, m_k4.m_v , t + dt*a4 ); 165 stepper_base_type::m_algebra.for_each6( m_x_tmp.m_v , in , dxdt , m_k2.m_v , m_k3.m_v , m_k4.m_v , 166 typename operations_type::template scale_sum5< value_type , time_type , time_type , time_type , time_type >( 1.0 , dt*b51 , dt*b52 , dt*b53 , dt*b54 )); 167 168 sys( m_x_tmp.m_v , m_k5.m_v , t + dt*a5 ); 169 stepper_base_type::m_algebra.for_each7( m_x_tmp.m_v , in , dxdt , m_k2.m_v , m_k3.m_v , m_k4.m_v , m_k5.m_v , 170 typename operations_type::template scale_sum6< value_type , time_type , time_type , time_type , time_type , time_type >( 1.0 , dt*b61 , dt*b62 , dt*b63 , dt*b64 , dt*b65 )); 171 172 sys( m_x_tmp.m_v , m_k6.m_v , t + dt*a6 ); 173 stepper_base_type::m_algebra.for_each6( out , in , dxdt , m_k3.m_v , m_k4.m_v , m_k6.m_v , 174 typename operations_type::template scale_sum5< value_type , time_type , time_type , time_type , time_type >( 1.0 , dt*c1 , dt*c3 , dt*c4 , dt*c6 )); 175 176 } 177 178 /** 179 * \brief Adjust the size of all temporaries in the stepper manually. 180 * \param x A state from which the size of the temporaries to be resized is deduced. 181 */ 182 template< class StateIn > adjust_size(const StateIn & x)183 void adjust_size( const StateIn &x ) 184 { 185 resize_impl( x ); 186 stepper_base_type::adjust_size( x ); 187 } 188 189 private: 190 191 template< class StateIn > resize_impl(const StateIn & x)192 bool resize_impl( const StateIn &x ) 193 { 194 bool resized = false; 195 resized |= adjust_size_by_resizeability( m_x_tmp , x , typename is_resizeable<state_type>::type() ); 196 resized |= adjust_size_by_resizeability( m_k2 , x , typename is_resizeable<deriv_type>::type() ); 197 resized |= adjust_size_by_resizeability( m_k3 , x , typename is_resizeable<deriv_type>::type() ); 198 resized |= adjust_size_by_resizeability( m_k4 , x , typename is_resizeable<deriv_type>::type() ); 199 resized |= adjust_size_by_resizeability( m_k5 , x , typename is_resizeable<deriv_type>::type() ); 200 resized |= adjust_size_by_resizeability( m_k6 , x , typename is_resizeable<deriv_type>::type() ); 201 return resized; 202 } 203 204 205 wrapped_state_type m_x_tmp; 206 wrapped_deriv_type m_k2, m_k3, m_k4, m_k5, m_k6; 207 resizer_type m_resizer; 208 209 }; 210 211 212 213 /************ DOXYGEN *************/ 214 215 /** 216 * \class runge_kutta_cash_karp54_classic 217 * \brief The Runge-Kutta Cash-Karp method implemented without the generic Runge-Kutta algorithm. 218 * 219 * The Runge-Kutta Cash-Karp method is one of the standard methods for 220 * solving ordinary differential equations, see 221 * <a href="http://en.wikipedia.org/wiki/Cash%E2%80%93Karp_method">en.wikipedia.org/wiki/Cash-Karp_method</a>. 222 * The method is explicit and fulfills the Error Stepper concept. Step size control 223 * is provided but continuous output is not available for this method. 224 * 225 * This class derives from explicit_error_stepper_base and inherits its interface via CRTP (current recurring 226 * template pattern). This class implements the method directly, hence the generic Runge-Kutta algorithm is not used. 227 * 228 * \tparam State The state type. 229 * \tparam Value The value type. 230 * \tparam Deriv The type representing the time derivative of the state. 231 * \tparam Time The time representing the independent variable - the time. 232 * \tparam Algebra The algebra type. 233 * \tparam Operations The operations type. 234 * \tparam Resizer The resizer policy type. 235 */ 236 237 238 /** 239 * \fn runge_kutta_cash_karp54_classic::runge_kutta_cash_karp54_classic( const algebra_type &algebra ) 240 * \brief Constructs the runge_kutta_cash_karp54_classic class. This constructor can be used as a default 241 * constructor if the algebra has a default constructor. 242 * \param algebra A copy of algebra is made and stored inside explicit_stepper_base. 243 */ 244 245 246 /** 247 * \fn runge_kutta_cash_karp54_classic::do_step_impl( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt , Err &xerr ) 248 * \brief This method performs one step. The derivative `dxdt` of `in` at the time `t` is passed to the method. 249 * 250 * The result is updated out-of-place, hence the input is in `in` and the output in `out`. Futhermore, an 251 * estimation of the error is stored in `xerr`. 252 * Access to this step functionality is provided by explicit_error_stepper_base and 253 * `do_step_impl` should not be called directly. 254 255 * 256 * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the 257 * Simple System concept. 258 * \param in The state of the ODE which should be solved. in is not modified in this method 259 * \param dxdt The derivative of x at t. 260 * \param t The value of the time, at which the step should be performed. 261 * \param out The result of the step is written in out. 262 * \param dt The step size. 263 * \param xerr The result of the error estimation is written in xerr. 264 */ 265 266 /** 267 * \fn runge_kutta_cash_karp54_classic::do_step_impl( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt ) 268 * \brief This method performs one step. The derivative `dxdt` of `in` at the time `t` is passed to the method. 269 * The result is updated out-of-place, hence the input is in `in` and the output in `out`. 270 * Access to this step functionality is provided by explicit_error_stepper_base and 271 * `do_step_impl` should not be called directly. 272 * 273 * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the 274 * Simple System concept. 275 * \param in The state of the ODE which should be solved. in is not modified in this method 276 * \param dxdt The derivative of x at t. 277 * \param t The value of the time, at which the step should be performed. 278 * \param out The result of the step is written in out. 279 * \param dt The step size. 280 */ 281 282 } // odeint 283 } // numeric 284 } // boost 285 286 287 288 289 #endif // BOOST_NUMERIC_ODEINT_STEPPER_RUNGE_KUTTA_CASH_KARP54_CLASSIC_HPP_INCLUDED 290