1 /*
2 [auto_generated]
3 libs/numeric/odeint/test/stepper_with_ranges.cpp
4
5 [begin_description]
6 This file tests if the steppers play well with Boost.Range.
7 [end_description]
8
9 Copyright 2011-2012 Karsten Ahnert
10 Copyright 2011-2013 Mario Mulansky
11
12 Distributed under the Boost Software License, Version 1.0.
13 (See accompanying file LICENSE_1_0.txt or
14 copy at http://www.boost.org/LICENSE_1_0.txt)
15 */
16
17 // disable checked iterator warning for msvc
18 #include <boost/config.hpp>
19 #ifdef BOOST_MSVC
20 #pragma warning(disable:4996)
21 #endif
22
23 #define BOOST_TEST_MODULE odeint_stepper_with_ranges
24
25 #include <boost/test/unit_test.hpp>
26
27 #include <vector>
28 #include <utility>
29 #include <iostream>
30
31 #include <boost/array.hpp>
32 #include <boost/range.hpp>
33 #include <boost/ref.hpp>
34
35 #include <boost/numeric/odeint/stepper/euler.hpp>
36 #include <boost/numeric/odeint/stepper/runge_kutta_cash_karp54_classic.hpp>
37 #include <boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp>
38 #include <boost/numeric/odeint/stepper/controlled_runge_kutta.hpp>
39 #include <boost/numeric/odeint/stepper/symplectic_euler.hpp>
40 #include <boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp>
41
42 typedef std::vector< double > state_type;
43 typedef boost::array< double , 3 > state_type2;
44
45 /* explicitly force range algebra for this array! */
46 namespace boost { namespace numeric { namespace odeint {
47
48 template<>
49 struct algebra_dispatcher< state_type2 >
50 { typedef range_algebra algebra_type; };
51
52 } } }
53
54
55 /*
56 * The two systems are needed, since for steppers with more than
57 * one internal step it is difficult to calculate the exact result
58 *
59 * system1 is suited for euler
60 */
61 struct system1
62 {
63 template< class State , class Deriv >
operator ()system164 void operator()( const State &x_ , Deriv &dxdt_ , double t )
65 {
66 typename boost::range_iterator< const State >::type x = boost::begin( x_ );
67 typename boost::range_iterator< Deriv >::type dxdt = boost::begin( dxdt_ );
68
69 dxdt[0] = x[0];
70 dxdt[1] = 2.0;
71 dxdt[2] = 3.0;
72 }
73
74 template< class State , class Deriv >
operator ()system175 void operator()( const State &x_ , const Deriv &dxdt_ , double t )
76 {
77 typename boost::range_iterator< const State >::type x = boost::begin( x_ );
78 typename boost::range_iterator< Deriv >::type dxdt = boost::begin( dxdt_ );
79
80 dxdt[0] = x[0];
81 dxdt[1] = 2.0;
82 dxdt[2] = 3.0;
83 }
84 };
85
86 /*
87 * system2 is suited for all steppers, it allows you to calculate the result analytically.
88 */
89 struct system2
90 {
91 template< class State , class Deriv >
operator ()system292 void operator()( const State &x_ , Deriv &dxdt_ , double t )
93 {
94 typename boost::range_iterator< Deriv >::type dxdt = boost::begin( dxdt_ );
95
96 dxdt[0] = 1.0;
97 dxdt[1] = 2.0;
98 dxdt[2] = 3.0;
99 }
100
101 template< class State , class Deriv >
operator ()system2102 void operator()( const State &x_ , const Deriv &dxdt_ , double t )
103 {
104 typename boost::range_iterator< Deriv >::type dxdt = boost::begin( dxdt_ );
105
106 dxdt[0] = 1.0;
107 dxdt[1] = 2.0;
108 dxdt[2] = 3.0;
109 }
110 };
111
112
113 /*
114 * Useful for Hamiltonian systems
115 */
116 struct ham_sys
117 {
118 template< class State , class Deriv >
operator ()ham_sys119 void operator()( const State &x_ , Deriv &dxdt_ )
120 {
121 typename boost::range_iterator< Deriv >::type dxdt = boost::begin( dxdt_ );
122 dxdt[0] = 1.0;
123 dxdt[1] = 2.0;
124 dxdt[2] = 3.0;
125 }
126
127 template< class State , class Deriv >
operator ()ham_sys128 void operator()( const State &x_ , const Deriv &dxdt_ )
129 {
130 typename boost::range_iterator< Deriv >::type dxdt = boost::begin( dxdt_ );
131 dxdt[0] = 1.0;
132 dxdt[1] = 2.0;
133 dxdt[2] = 3.0;
134 }
135 };
136
137
138 struct vector_fixture
139 {
140 const static size_t dim = 6;
141 boost::array< double , dim > in;
142 boost::array< double , dim > q;
143 boost::array< double , dim > p;
144 state_type err;
145
vector_fixturevector_fixture146 vector_fixture( void )
147 : in() , err( 3 )
148 {
149 for( size_t i=0 ; i<dim ; ++i )
150 {
151 in[ i ] = q[i] = p[i] = double( i );
152 }
153 for( size_t i=0 ; i<3 ; ++i )
154 {
155 err[i] = double( i ) * 10.0;
156 }
157 }
158
~vector_fixturevector_fixture159 ~vector_fixture( void )
160 {
161 }
162 };
163
164 #define CHECK_VALUES( x , x0 , x1 , x2 , x3 , x4 , x5 ) \
165 BOOST_CHECK_CLOSE( x[0] , x0 , 1.0e-8 ); \
166 BOOST_CHECK_CLOSE( x[1] , x1 , 1.0e-8 ); \
167 BOOST_CHECK_CLOSE( x[2] , x2 , 1.0e-8 ); \
168 BOOST_CHECK_CLOSE( x[3] , x3 , 1.0e-8 ); \
169 BOOST_CHECK_CLOSE( x[4] , x4 , 1.0e-8 ); \
170 BOOST_CHECK_CLOSE( x[5] , x5 , 1.0e-8 )
171
172
173
174 BOOST_AUTO_TEST_SUITE( stepper_with_ranges )
175
BOOST_AUTO_TEST_CASE(explicit_euler_with_range_v1)176 BOOST_AUTO_TEST_CASE( explicit_euler_with_range_v1 )
177 {
178 vector_fixture f;
179 boost::numeric::odeint::euler< state_type > euler;
180 euler.do_step( system1() , std::make_pair( f.in.begin() + 1 , f.in.begin() + 4 ) , 0.1 , 0.1 );
181 CHECK_VALUES( f.in , 0.0 , 1.1 , 2.2 , 3.3 , 4.0 , 5.0 );
182 }
183
184
185
BOOST_AUTO_TEST_CASE(explicit_error_k54_with_range_v1)186 BOOST_AUTO_TEST_CASE( explicit_error_k54_with_range_v1 )
187 {
188 vector_fixture f;
189 boost::numeric::odeint::runge_kutta_cash_karp54_classic< state_type > rk54;
190 rk54.do_step( system2() , std::make_pair( f.in.begin() + 1 , f.in.begin() + 4 ) , 0.1 , 0.1 );
191 CHECK_VALUES( f.in , 0.0 , 1.1 , 2.2 , 3.3 , 4.0 , 5.0 );
192 }
193
BOOST_AUTO_TEST_CASE(explicit_error_k54_with_range_v5)194 BOOST_AUTO_TEST_CASE( explicit_error_k54_with_range_v5 )
195 {
196 vector_fixture f;
197 boost::numeric::odeint::runge_kutta_cash_karp54_classic< state_type > rk54;
198 rk54.do_step( system2() , std::make_pair( f.in.begin() + 1 , f.in.begin() + 4 ) , 0.1 , 0.1 , f.err );
199 CHECK_VALUES( f.in , 0.0 , 1.1 , 2.2 , 3.3 , 4.0 , 5.0 );
200 }
201
202
BOOST_AUTO_TEST_CASE(runge_kutta_dopri5_with_range_v1)203 BOOST_AUTO_TEST_CASE( runge_kutta_dopri5_with_range_v1 )
204 {
205 vector_fixture f;
206 boost::numeric::odeint::runge_kutta_dopri5< state_type > dopri5;
207 dopri5.do_step( system2() , std::make_pair( f.in.begin() + 1 , f.in.begin() + 4 ) , 0.1 , 0.1 );
208 CHECK_VALUES( f.in , 0.0 , 1.1 , 2.2 , 3.3 , 4.0 , 5.0 );
209 }
210
211
BOOST_AUTO_TEST_CASE(runge_kutta_dopri5_with_range_v5)212 BOOST_AUTO_TEST_CASE( runge_kutta_dopri5_with_range_v5 )
213 {
214 vector_fixture f;
215 boost::numeric::odeint::runge_kutta_dopri5< state_type > dopri5;
216 dopri5.do_step( system2() , std::make_pair( f.in.begin() + 1 , f.in.begin() + 4 ) , 0.1 , 0.1 , f.err );
217 CHECK_VALUES( f.in , 0.0 , 1.1 , 2.2 , 3.3 , 4.0 , 5.0 );
218 }
219
220
BOOST_AUTO_TEST_CASE(controlled_error_stepper_rk54)221 BOOST_AUTO_TEST_CASE( controlled_error_stepper_rk54 )
222 {
223 double t = 0.0 , dt = 0.1;
224 vector_fixture f;
225 boost::numeric::odeint::controlled_runge_kutta< boost::numeric::odeint::runge_kutta_cash_karp54_classic< state_type > > stepper;
226 stepper.try_step( system2() , std::make_pair( f.in.begin() + 1 , f.in.begin() + 4 ) , t , dt );
227 CHECK_VALUES( f.in , 0.0 , 1.1 , 2.2 , 3.3 , 4.0 , 5.0 );
228 }
229
BOOST_AUTO_TEST_CASE(controlled_error_stepper_dopri5)230 BOOST_AUTO_TEST_CASE( controlled_error_stepper_dopri5 )
231 {
232 double t = 0.0 , dt = 0.1;
233 vector_fixture f;
234 boost::numeric::odeint::controlled_runge_kutta< boost::numeric::odeint::runge_kutta_dopri5< state_type > > stepper;
235 stepper.try_step( system2() , std::make_pair( f.in.begin() + 1 , f.in.begin() + 4 ) , t , dt );
236 CHECK_VALUES( f.in , 0.0 , 1.1 , 2.2 , 3.3 , 4.0 , 5.0 );
237 }
238
239
BOOST_AUTO_TEST_CASE(symplectic_euler_coor_func)240 BOOST_AUTO_TEST_CASE( symplectic_euler_coor_func )
241 {
242 vector_fixture f;
243 boost::numeric::odeint::symplectic_euler< state_type > euler;
244 euler.do_step( ham_sys() ,
245 std::make_pair( f.q.begin() + 1 , f.q.begin() + 4 ) ,
246 std::make_pair( f.p.begin() + 3 , f.p.begin() + 6 ) ,
247 0.0 , 0.1 );
248 CHECK_VALUES( f.q , 0.0 , 1.3 , 2.4 , 3.5 , 4.0 , 5.0 );
249 CHECK_VALUES( f.p , 0.0 , 1.0 , 2.0 , 3.1 , 4.2 , 5.3 );
250 }
251
BOOST_AUTO_TEST_CASE(symplectic_euler_coor_and_mom_func)252 BOOST_AUTO_TEST_CASE( symplectic_euler_coor_and_mom_func )
253 {
254 vector_fixture f;
255 boost::numeric::odeint::symplectic_euler< state_type > euler;
256 euler.do_step( std::make_pair( ham_sys() , ham_sys() ) ,
257 std::make_pair( f.q.begin() + 1 , f.q.begin() + 4 ) ,
258 std::make_pair( f.p.begin() + 3 , f.p.begin() + 6 ) ,
259 0.0 , 0.1 );
260 CHECK_VALUES( f.q , 0.0 , 1.1 , 2.2 , 3.3 , 4.0 , 5.0 );
261 CHECK_VALUES( f.p , 0.0 , 1.0 , 2.0 , 3.1 , 4.2 , 5.3 );
262 }
263
264
BOOST_AUTO_TEST_CASE(dense_output_euler_with_ranges)265 BOOST_AUTO_TEST_CASE( dense_output_euler_with_ranges )
266 {
267 using namespace boost::numeric::odeint;
268 vector_fixture f;
269 dense_output_runge_kutta< euler< state_type > > stepper;
270 stepper.initialize( std::make_pair( f.in.begin() + 1, f.in.begin() + 4 ) , 0.0 , 0.1 );
271 stepper.do_step( system1() );
272 stepper.calc_state( 0.05 , std::make_pair( f.in.begin() + 1 ,f.in.begin() +4 ) );
273 CHECK_VALUES( f.in , 0.0 , 1.05 , 2.1 , 3.15 , 4.0 , 5.0 );
274 }
275
BOOST_AUTO_TEST_CASE(dense_output_dopri5_with_ranges)276 BOOST_AUTO_TEST_CASE( dense_output_dopri5_with_ranges )
277 {
278 using namespace boost::numeric::odeint;
279 vector_fixture f;
280 dense_output_runge_kutta<
281 controlled_runge_kutta<
282 runge_kutta_dopri5< state_type >
283 > > stepper;
284 stepper.initialize( std::make_pair( f.in.begin() + 1, f.in.begin() + 4 ) , 0.0 , 0.1 );
285 stepper.do_step( system2() );
286 stepper.calc_state( 0.05 , std::make_pair( f.in.begin() + 1 ,f.in.begin() +4 ) );
287 CHECK_VALUES( f.in , 0.0 , 1.05 , 2.1 , 3.15 , 4.0 , 5.0 );
288 }
289
290
291
292 BOOST_AUTO_TEST_SUITE_END()
293