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