1 /* 2 [auto_generated] 3 boost/numeric/odeint/stepper/runge_kutta4_classic.hpp 4 5 [begin_description] 6 Implementation for the classical Runge Kutta stepper. 7 [end_description] 8 9 Copyright 2010-2013 Karsten Ahnert 10 Copyright 2010-2013 Mario Mulansky 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_KUTTA4_CLASSIC_HPP_INCLUDED 20 #define BOOST_NUMERIC_ODEINT_STEPPER_RUNGE_KUTTA4_CLASSIC_HPP_INCLUDED 21 22 23 24 #include <boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp> 25 #include <boost/numeric/odeint/algebra/range_algebra.hpp> 26 #include <boost/numeric/odeint/algebra/default_operations.hpp> 27 #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp> 28 #include <boost/numeric/odeint/algebra/operations_dispatcher.hpp> 29 30 #include <boost/numeric/odeint/util/state_wrapper.hpp> 31 #include <boost/numeric/odeint/util/is_resizeable.hpp> 32 #include <boost/numeric/odeint/util/resizer.hpp> 33 34 namespace boost { 35 namespace numeric { 36 namespace odeint { 37 38 template< 39 class State , 40 class Value = double , 41 class Deriv = State , 42 class Time = Value , 43 class Algebra = typename algebra_dispatcher< State >::algebra_type , 44 class Operations = typename operations_dispatcher< State >::operations_type , 45 class Resizer = initially_resizer 46 > 47 #ifndef DOXYGEN_SKIP 48 class runge_kutta4_classic 49 : public explicit_stepper_base< 50 runge_kutta4_classic< State , Value , Deriv , Time , Algebra , Operations , Resizer > , 51 4 , State , Value , Deriv , Time , Algebra , Operations , Resizer > 52 #else 53 class runge_kutta4_classic : public explicit_stepper_base 54 #endif 55 { 56 57 public : 58 59 #ifndef DOXYGEN_SKIP 60 typedef explicit_stepper_base< 61 runge_kutta4_classic< State , Value , Deriv , Time , Algebra , Operations , Resizer > , 62 4 , State , Value , Deriv , Time , Algebra , Operations , Resizer > stepper_base_type; 63 #else 64 typedef explicit_stepper_base< runge_kutta4_classic< ... > , ... > stepper_base_type; 65 #endif 66 67 typedef typename stepper_base_type::state_type state_type; 68 typedef typename stepper_base_type::value_type value_type; 69 typedef typename stepper_base_type::deriv_type deriv_type; 70 typedef typename stepper_base_type::time_type time_type; 71 typedef typename stepper_base_type::algebra_type algebra_type; 72 typedef typename stepper_base_type::operations_type operations_type; 73 typedef typename stepper_base_type::resizer_type resizer_type; 74 75 #ifndef DOXYGEN_SKIP 76 typedef typename stepper_base_type::stepper_type stepper_type; 77 typedef typename stepper_base_type::wrapped_state_type wrapped_state_type; 78 typedef typename stepper_base_type::wrapped_deriv_type wrapped_deriv_type; 79 #endif // DOXYGEN_SKIP 80 81 82 runge_kutta4_classic(const algebra_type & algebra=algebra_type ())83 runge_kutta4_classic( const algebra_type &algebra = algebra_type() ) : stepper_base_type( algebra ) 84 { } 85 86 87 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)88 void do_step_impl( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt ) 89 { 90 // ToDo : check if size of in,dxdt,out are equal? 91 92 static const value_type val1 = static_cast< value_type >( 1 ); 93 94 m_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_impl< StateIn > , detail::ref( *this ) , detail::_1 ) ); 95 96 typename odeint::unwrap_reference< System >::type &sys = system; 97 98 const time_type dh = dt / static_cast< value_type >( 2 ); 99 const time_type th = t + dh; 100 101 // dt * dxdt = k1 102 // m_x_tmp = x + dh*dxdt 103 stepper_base_type::m_algebra.for_each3( m_x_tmp.m_v , in , dxdt , 104 typename operations_type::template scale_sum2< value_type , time_type >( val1 , dh ) ); 105 106 107 // dt * m_dxt = k2 108 sys( m_x_tmp.m_v , m_dxt.m_v , th ); 109 110 // m_x_tmp = x + dh*m_dxt 111 stepper_base_type::m_algebra.for_each3( m_x_tmp.m_v , in , m_dxt.m_v , 112 typename operations_type::template scale_sum2< value_type , time_type >( val1 , dh ) ); 113 114 115 // dt * m_dxm = k3 116 sys( m_x_tmp.m_v , m_dxm.m_v , th ); 117 //m_x_tmp = x + dt*m_dxm 118 stepper_base_type::m_algebra.for_each3( m_x_tmp.m_v , in , m_dxm.m_v , 119 typename operations_type::template scale_sum2< value_type , time_type >( val1 , dt ) ); 120 121 122 // dt * m_dxh = k4 123 sys( m_x_tmp.m_v , m_dxh.m_v , t + dt ); 124 125 //x += dt/6 * ( m_dxdt + m_dxt + val2*m_dxm ) 126 time_type dt6 = dt / static_cast< value_type >( 6 ); 127 time_type dt3 = dt / static_cast< value_type >( 3 ); 128 stepper_base_type::m_algebra.for_each6( out , in , dxdt , m_dxt.m_v , m_dxm.m_v , m_dxh.m_v , 129 typename operations_type::template scale_sum5< value_type , time_type , time_type , time_type , time_type >( 1.0 , dt6 , dt3 , dt3 , dt6 ) ); 130 131 // x += dt/6 * m_dxdt + dt/3 * m_dxt ) 132 // stepper_base_type::m_algebra.for_each4( out , in , dxdt , m_dxt.m_v , 133 // typename operations_type::template scale_sum3< value_type , time_type , time_type >( 1.0 , dt6 , dt3 ) ); 134 // // x += dt/3 * m_dxm + dt/6 * m_dxh ) 135 // stepper_base_type::m_algebra.for_each4( out , out , m_dxm.m_v , m_dxh.m_v , 136 // typename operations_type::template scale_sum3< value_type , time_type , time_type >( 1.0 , dt3 , dt6 ) ); 137 138 } 139 140 template< class StateType > adjust_size(const StateType & x)141 void adjust_size( const StateType &x ) 142 { 143 resize_impl( x ); 144 stepper_base_type::adjust_size( x ); 145 } 146 147 private: 148 149 template< class StateIn > resize_impl(const StateIn & x)150 bool resize_impl( const StateIn &x ) 151 { 152 bool resized = false; 153 resized |= adjust_size_by_resizeability( m_x_tmp , x , typename is_resizeable<state_type>::type() ); 154 resized |= adjust_size_by_resizeability( m_dxm , x , typename is_resizeable<deriv_type>::type() ); 155 resized |= adjust_size_by_resizeability( m_dxt , x , typename is_resizeable<deriv_type>::type() ); 156 resized |= adjust_size_by_resizeability( m_dxh , x , typename is_resizeable<deriv_type>::type() ); 157 return resized; 158 } 159 160 161 resizer_type m_resizer; 162 163 wrapped_deriv_type m_dxt; 164 wrapped_deriv_type m_dxm; 165 wrapped_deriv_type m_dxh; 166 wrapped_state_type m_x_tmp; 167 168 }; 169 170 171 /********* DOXYGEN *********/ 172 173 /** 174 * \class runge_kutta4_classic 175 * \brief The classical Runge-Kutta stepper of fourth order. 176 * 177 * The Runge-Kutta method of fourth order is one standard method for 178 * solving ordinary differential equations and is widely used, see also 179 * <a href="http://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods">en.wikipedia.org/wiki/Runge-Kutta_methods</a> 180 * The method is explicit and fulfills the Stepper concept. Step size control 181 * or continuous output are not provided. This class implements the method directly, hence the 182 * generic Runge-Kutta algorithm is not used. 183 * 184 * This class derives from explicit_stepper_base and inherits its interface via 185 * CRTP (current recurring template pattern). For more details see 186 * explicit_stepper_base. 187 * 188 * \tparam State The state type. 189 * \tparam Value The value type. 190 * \tparam Deriv The type representing the time derivative of the state. 191 * \tparam Time The time representing the independent variable - the time. 192 * \tparam Algebra The algebra type. 193 * \tparam Operations The operations type. 194 * \tparam Resizer The resizer policy type. 195 */ 196 197 /** 198 * \fn runge_kutta4_classic::runge_kutta4_classic( const algebra_type &algebra ) 199 * \brief Constructs the runge_kutta4_classic class. This constructor can be used as a default 200 * constructor if the algebra has a default constructor. 201 * \param algebra A copy of algebra is made and stored inside explicit_stepper_base. 202 */ 203 204 205 /** 206 * \fn runge_kutta4_classic::do_step_impl( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt ) 207 * \brief This method performs one step. The derivative `dxdt` of `in` at the time `t` is passed to the method. 208 * The result is updated out of place, hence the input is in `in` and the output in `out`. 209 * Access to this step functionality is provided by explicit_stepper_base and 210 * `do_step_impl` should not be called directly. 211 * 212 * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the 213 * Simple System concept. 214 * \param in The state of the ODE which should be solved. in is not modified in this method 215 * \param dxdt The derivative of x at t. 216 * \param t The value of the time, at which the step should be performed. 217 * \param out The result of the step is written in out. 218 * \param dt The step size. 219 */ 220 221 /** 222 * \fn runge_kutta4_classic::adjust_size( const StateType &x ) 223 * \brief Adjust the size of all temporaries in the stepper manually. 224 * \param x A state from which the size of the temporaries to be resized is deduced. 225 */ 226 227 } // odeint 228 } // numeric 229 } // boost 230 231 232 #endif // BOOST_NUMERIC_ODEINT_STEPPER_RUNGE_KUTTA4_CLASSIC_HPP_INCLUDED 233