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