1 /*
2 [auto_generated]
3 libs/numeric/odeint/test/integrate.cpp
4
5 [begin_description]
6 This file tests the integrate function and its variants.
7 [end_description]
8
9 Copyright 2011-2012 Mario Mulansky
10 Copyright 2011-2012 Karsten Ahnert
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
18 #define BOOST_TEST_MODULE odeint_integrate_functions
19
20 #include <boost/type_traits.hpp>
21
22 #include <vector>
23 #include <cmath>
24 #include <iostream>
25
26 #include <boost/numeric/odeint/config.hpp>
27
28 #include <boost/array.hpp>
29 #include <boost/ref.hpp>
30 #include <boost/iterator/counting_iterator.hpp>
31
32 #include <boost/test/unit_test.hpp>
33
34 #include <boost/mpl/vector.hpp>
35
36 // nearly everything from odeint is used in these tests
37 #ifndef ODEINT_INTEGRATE_ITERATOR
38 #include <boost/numeric/odeint/integrate/integrate_const.hpp>
39 #include <boost/numeric/odeint/integrate/integrate_adaptive.hpp>
40 #include <boost/numeric/odeint/integrate/integrate_times.hpp>
41 #include <boost/numeric/odeint/integrate/integrate_n_steps.hpp>
42 #else
43 #include <boost/numeric/odeint/iterator/integrate/integrate_const.hpp>
44 #include <boost/numeric/odeint/iterator/integrate/integrate_adaptive.hpp>
45 #include <boost/numeric/odeint/iterator/integrate/integrate_times.hpp>
46 #include <boost/numeric/odeint/iterator/integrate/integrate_n_steps.hpp>
47 #endif
48 #include <boost/numeric/odeint/stepper/euler.hpp>
49 #include <boost/numeric/odeint/stepper/modified_midpoint.hpp>
50 #include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
51 #include <boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp>
52 #include <boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp>
53 #include <boost/numeric/odeint/stepper/runge_kutta_fehlberg78.hpp>
54 #include <boost/numeric/odeint/stepper/controlled_runge_kutta.hpp>
55 #include <boost/numeric/odeint/stepper/bulirsch_stoer.hpp>
56 #include <boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp>
57 #include <boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp>
58
59 #include <boost/numeric/odeint/stepper/adaptive_adams_bashforth_moulton.hpp>
60 #include <boost/numeric/odeint/stepper/controlled_adams_bashforth_moulton.hpp>
61
62 #include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
63
64 using namespace boost::unit_test;
65 using namespace boost::numeric::odeint;
66 namespace mpl = boost::mpl;
67
68
69 typedef double value_type;
70 typedef std::vector< value_type > state_type;
71
lorenz(const state_type & x,state_type & dxdt,const value_type t)72 void lorenz( const state_type &x , state_type &dxdt , const value_type t )
73 {
74 //const value_type sigma( 10.0 );
75 const value_type R( 28.0 );
76 const value_type b( value_type( 8.0 ) / value_type( 3.0 ) );
77
78 // first component trivial
79 dxdt[0] = 1.0; //sigma * ( x[1] - x[0] );
80 dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
81 dxdt[2] = -b * x[2] + x[0] * x[1];
82 }
83
84 struct push_back_time
85 {
86 std::vector< double >& m_times;
87
88 state_type& m_x;
89
push_back_timepush_back_time90 push_back_time( std::vector< double > × , state_type &x )
91 : m_times( times ) , m_x( x ) { }
92
operator ()push_back_time93 void operator()( const state_type &x , double t )
94 {
95 m_times.push_back( t );
96 boost::numeric::odeint::copy( x , m_x );
97 }
98 };
99
100 template< class Stepper >
101 struct perform_integrate_const_test
102 {
operator ()perform_integrate_const_test103 void operator()( const value_type t_end , const value_type dt )
104 {
105 // std::cout << "Testing integrate_const with " << typeid( Stepper ).name() << std::endl;
106
107 state_type x( 3 , 10.0 ) , x_end( 3 );
108
109 std::vector< value_type > times;
110
111 size_t steps = integrate_const( Stepper() , lorenz , x , 0.0 , t_end ,
112 dt , push_back_time( times , x_end ) );
113
114 std::cout.precision(16);
115 std::cout << t_end << " (" << dt << "), " << steps << " , " << times.size() << " , " << 10.0+dt*steps << "=" << x_end[0] << std::endl;
116
117 std::cout << static_cast<int>(floor(t_end/dt)) << " , " << t_end/dt << std::endl;
118
119 BOOST_CHECK_EQUAL( static_cast<int>(times.size()) , static_cast<int>(floor(t_end/dt))+1 );
120
121 for( size_t i=0 ; i<times.size() ; ++i )
122 {
123 //std::cout << i << " , " << times[i] << " , " << static_cast< value_type >(i)*dt << std::endl;
124 // check if observer was called at times 0,1,2,...
125 BOOST_CHECK_SMALL( times[i] - static_cast< value_type >(i)*dt , (i+1) * 2E-16 );
126 }
127
128 // check first, trivial, component
129 BOOST_CHECK_SMALL( (10.0 + times[times.size()-1]) - x_end[0] , 1E-6 ); // precision of steppers: 1E-6
130 //BOOST_CHECK_EQUAL( x[1] , x_end[1] );
131 //BOOST_CHECK_EQUAL( x[2] , x_end[2] );
132 }
133 };
134
135 template< class Stepper >
136 struct perform_integrate_adaptive_test
137 {
operator ()perform_integrate_adaptive_test138 void operator()( const value_type t_end = 10.0 , const value_type dt = 0.03 )
139 {
140 // std::cout << "Testing integrate_adaptive with " << typeid( Stepper ).name() << std::endl;
141
142 state_type x( 3 , 10.0 ) , x_end( 3 );
143
144 std::vector< value_type > times;
145
146 size_t steps = integrate_adaptive( Stepper() , lorenz , x , 0.0 , t_end ,
147 dt , push_back_time( times , x_end ) );
148
149 // std::cout << t_end << " , " << steps << " , " << times.size() << " , " << dt << " , " << 10.0+t_end << "=" << x_end[0] << std::endl;
150
151 BOOST_CHECK_EQUAL( times.size() , steps+1 );
152
153 BOOST_CHECK_SMALL( times[0] - 0.0 , 2E-16 );
154 BOOST_CHECK_SMALL( times[times.size()-1] - t_end , times.size() * 2E-16 );
155
156 // check first, trivial, component
157 BOOST_CHECK_SMALL( (10.0 + t_end) - x_end[0] , 1E-6 ); // precision of steppers: 1E-6
158 // BOOST_CHECK_EQUAL( x[1] , x_end[1] );
159 // BOOST_CHECK_EQUAL( x[2] , x_end[2] );
160 }
161 };
162
163
164 template< class Stepper >
165 struct perform_integrate_times_test
166 {
operator ()perform_integrate_times_test167 void operator()( const int n = 10 , const int dn=1 , const value_type dt = 0.03 )
168 {
169 // std::cout << "Testing integrate_times with " << typeid( Stepper ).name() << std::endl;
170
171 state_type x( 3 ) , x_end( 3 );
172 x[0] = x[1] = x[2] = 10.0;
173
174 std::vector< double > times;
175
176 std::vector< double > obs_times( abs(n) );
177 for( int i=0 ; boost::numeric::odeint::detail::less_with_sign( static_cast<double>(i) ,
178 static_cast<double>(obs_times.size()) ,
179 dt ) ; i+=dn )
180 {
181 obs_times[i] = i;
182 }
183 // simple stepper
184 integrate_times( Stepper() , lorenz , x , obs_times.begin() , obs_times.end() ,
185 dt , push_back_time( times , x_end ) );
186
187 BOOST_CHECK_EQUAL( static_cast<int>(times.size()) , abs(n) );
188
189 for( size_t i=0 ; i<times.size() ; ++i )
190 // check if observer was called at times 0,1,2,...
191 BOOST_CHECK_EQUAL( times[i] , static_cast<double>(i) );
192
193 // check first, trivial, component
194 BOOST_CHECK_SMALL( (10.0 + 1.0*times[times.size()-1]) - x_end[0] , 1E-6 ); // precision of steppers: 1E-6
195 // BOOST_CHECK_EQUAL( x[1] , x_end[1] );
196 // BOOST_CHECK_EQUAL( x[2] , x_end[2] );
197 }
198 };
199
200 template< class Stepper >
201 struct perform_integrate_n_steps_test
202 {
operator ()perform_integrate_n_steps_test203 void operator()( const int n = 200 , const value_type dt = 0.01 )
204 {
205 // std::cout << "Testing integrate_n_steps with " << typeid( Stepper ).name() << ". ";
206 // std::cout << "dt=" << dt << std::endl;
207
208 state_type x( 3 ) , x_end( 3 );
209 x[0] = x[1] = x[2] = 10.0;
210
211 std::vector< double > times;
212
213 // simple stepper
214 value_type end_time = integrate_n_steps( Stepper() , lorenz , x , 0.0 , dt , n , push_back_time( times , x_end ) );
215
216 BOOST_CHECK_SMALL( end_time - n*dt , 2E-16 );
217 BOOST_CHECK_EQUAL( static_cast<int>(times.size()) , n+1 );
218
219 for( size_t i=0 ; i<times.size() ; ++i )
220 {
221 // check if observer was called at times 0,1,2,...
222 BOOST_CHECK_SMALL(times[i] - static_cast< value_type >(i) * dt, 2E-16);
223 }
224 // check first, trivial, component
225 BOOST_CHECK_SMALL( (10.0 + end_time) - x_end[0] , 1E-6 ); // precision of steppers: 1E-6
226 // BOOST_CHECK_EQUAL( x[1] , x_end[1] );
227 // BOOST_CHECK_EQUAL( x[2] , x_end[2] );
228
229 }
230 };
231
232
233
234 class stepper_methods : public mpl::vector<
235 euler< state_type > ,
236 modified_midpoint< state_type > ,
237 runge_kutta4< state_type > ,
238 runge_kutta_cash_karp54< state_type > ,
239 runge_kutta_dopri5< state_type > ,
240 runge_kutta_fehlberg78< state_type > ,
241 controlled_runge_kutta< runge_kutta_cash_karp54< state_type > > ,
242 controlled_runge_kutta< runge_kutta_dopri5< state_type > > ,
243 controlled_runge_kutta< runge_kutta_fehlberg78< state_type > > ,
244 bulirsch_stoer< state_type > ,
245 dense_output_runge_kutta< controlled_runge_kutta< runge_kutta_dopri5< state_type > > >,
246 adaptive_adams_bashforth_moulton<3, state_type>,
247 adaptive_adams_bashforth_moulton<5, state_type>,
248 adaptive_adams_bashforth_moulton<7, state_type>,
249 controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<3, state_type> >,
250 controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<5, state_type> >
251 //bulirsch_stoer_dense_out< state_type >
252 > { };
253
254
255
256 BOOST_AUTO_TEST_SUITE( integrate_test )
257
BOOST_AUTO_TEST_CASE_TEMPLATE(integrate_const_test_case,Stepper,stepper_methods)258 BOOST_AUTO_TEST_CASE_TEMPLATE( integrate_const_test_case , Stepper, stepper_methods )
259 {
260 perform_integrate_const_test< Stepper > tester;
261 tester( 1.005 , 0.01 );
262 tester( 1.0 , 0.01 );
263 tester( 1.1 , 0.01 );
264 tester( -1.005 , -0.01 );
265 }
266
267
BOOST_AUTO_TEST_CASE_TEMPLATE(integrate_adaptive_test_case,Stepper,stepper_methods)268 BOOST_AUTO_TEST_CASE_TEMPLATE( integrate_adaptive_test_case , Stepper, stepper_methods )
269 {
270 perform_integrate_adaptive_test< Stepper > tester;
271 tester( 1.005 , 0.01 );
272 tester( 1.0 , 0.01 );
273 tester( 1.1 , 0.01 );
274 tester( -1.005 , -0.01 );
275 }
276
277
BOOST_AUTO_TEST_CASE_TEMPLATE(integrate_times_test_case,Stepper,stepper_methods)278 BOOST_AUTO_TEST_CASE_TEMPLATE( integrate_times_test_case , Stepper, stepper_methods )
279 {
280 perform_integrate_times_test< Stepper > tester;
281 tester();
282 //tester( -10 , -0.01 );
283 }
284
BOOST_AUTO_TEST_CASE_TEMPLATE(integrate_n_steps_test_case,Stepper,stepper_methods)285 BOOST_AUTO_TEST_CASE_TEMPLATE( integrate_n_steps_test_case , Stepper, stepper_methods )
286 {
287 if(!boost::is_same<Stepper, controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<3, state_type> > >::value &&
288 !boost::is_same<Stepper, controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<5, state_type> > >::value)
289 {
290 perform_integrate_n_steps_test< Stepper > tester;
291 tester();
292 tester( 200 , 0.01 );
293 tester( 200 , 0.01 );
294 tester( 200 , 0.01 );
295 tester( 200 , -0.01 );
296 }
297 }
298
299 BOOST_AUTO_TEST_SUITE_END()
300