1 /*
2  [auto_generated]
3  libs/numeric/odeint/test/integrate_times.cpp
4 
5  [begin_description]
6  This file tests the integrate_times function and its variants.
7  [end_description]
8 
9  Copyright 2011-2012 Karsten Ahnert
10  Copyright 2011-2012 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 #define BOOST_TEST_MODULE odeint_integrate_times
18 
19 #include <boost/test/unit_test.hpp>
20 
21 #include <utility>
22 #include <iostream>
23 #include <vector>
24 
25 #include <boost/ref.hpp>
26 #include <boost/iterator/counting_iterator.hpp>
27 
28 #ifndef ODEINT_INTEGRATE_ITERATOR
29 #include <boost/numeric/odeint/integrate/integrate_times.hpp>
30 #include <boost/numeric/odeint/integrate/integrate_adaptive.hpp>
31 #else
32 #include <boost/numeric/odeint/iterator/integrate/integrate_times.hpp>
33 #include <boost/numeric/odeint/iterator/integrate/integrate_adaptive.hpp>
34 #endif
35 #include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
36 #include <boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp>
37 #include <boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp>
38 #include <boost/numeric/odeint/stepper/controlled_runge_kutta.hpp>
39 #include <boost/numeric/odeint/stepper/bulirsch_stoer.hpp>
40 #include <boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp>
41 #include <boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp>
42 
43 using namespace boost::unit_test;
44 using namespace boost::numeric::odeint;
45 
46 typedef double value_type;
47 typedef std::vector< value_type > state_type;
48 
49 
lorenz(const state_type & x,state_type & dxdt,const value_type t)50 void lorenz( const state_type &x , state_type &dxdt , const value_type t )
51 {
52     BOOST_CHECK( t >= 0.0 );
53 
54     const value_type sigma( 10.0 );
55     const value_type R( 28.0 );
56     const value_type b( value_type( 8.0 ) / value_type( 3.0 ) );
57 
58     dxdt[0] = sigma * ( x[1] - x[0] );
59     dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
60     dxdt[2] = -b * x[2] + x[0] * x[1];
61 }
62 
63 struct push_back_time
64 {
65     std::vector< double >& m_times;
66 
push_back_timepush_back_time67     push_back_time( std::vector< double > &times )
68     :  m_times( times ) { }
69 
operator ()push_back_time70     void operator()( const state_type &x , double t )
71     {
72         m_times.push_back( t );
73     }
74 };
75 
76 BOOST_AUTO_TEST_SUITE( integrate_times_test )
77 
BOOST_AUTO_TEST_CASE(test_integrate_times)78 BOOST_AUTO_TEST_CASE( test_integrate_times )
79 {
80 
81     state_type x( 3 );
82     x[0] = x[1] = x[2] = 10.0;
83 
84     const value_type dt = 0.03;
85 
86     std::vector< double > times;
87 
88     std::cout << "test rk4 stepper" << std::endl;
89 
90     // simple stepper
91     integrate_times( runge_kutta4< state_type >() , lorenz , x , boost::counting_iterator<int>(0) , boost::counting_iterator<int>(10) ,
92                 dt , push_back_time( times ) );
93 
94     for( int i=0 ; i<10 ; ++i )
95         // check if observer was called at times 0,1,2,...
96         BOOST_CHECK_EQUAL( times[i] , static_cast<double>(i) );
97     times.clear();
98 
99     std::cout << "test dopri5 stepper" << std::endl;
100 
101     // controlled stepper
102     integrate_times( controlled_runge_kutta< runge_kutta_dopri5< state_type > >() , lorenz , x ,
103                 boost::counting_iterator<int>(0) , boost::counting_iterator<int>(10) ,
104                 dt , push_back_time( times ) );
105 
106     for( int i=0 ; i<10 ; ++i )
107         // check if observer was called at times 0,1,2,...
108         BOOST_CHECK_EQUAL( times[i] , static_cast<double>(i) );
109     times.clear();
110 
111     std::cout << "test BS stepper" << std::endl;
112 
113     //another controlled stepper
114     integrate_times( bulirsch_stoer< state_type >() , lorenz , x ,
115                 boost::counting_iterator<int>(0) , boost::counting_iterator<int>(10) ,
116                 dt , push_back_time( times ) );
117 
118     for( int i=0 ; i<10 ; ++i )
119         // check if observer was called at times 0,1,2,...
120         BOOST_CHECK_EQUAL( times[i] , static_cast<double>(i) );
121     times.clear();
122 
123     std::cout << "test dense_out stepper" << std::endl;
124 
125     // dense output stepper
126     integrate_times( dense_output_runge_kutta< controlled_runge_kutta< runge_kutta_dopri5< state_type > > >() ,
127                      lorenz , x ,
128                      boost::counting_iterator<int>(0) , boost::counting_iterator<int>(10) ,
129                      dt , push_back_time( times ) );
130 
131     for( int i=0 ; i<10 ; ++i )
132         // check if observer was called at times 0,1,2,...
133         BOOST_CHECK_EQUAL( times[i] , static_cast<double>(i) );
134 
135     std::cout << "test BS_do stepper" << std::endl;
136 
137     integrate_times( bulirsch_stoer_dense_out< state_type >() , lorenz , x ,
138                 boost::counting_iterator<int>(0) , boost::counting_iterator<int>(10) ,
139                 dt , push_back_time( times ) );
140 
141     for( int i=0 ; i<10 ; ++i )
142         // check if observer was called at times 0,1,2,...
143         BOOST_CHECK_EQUAL( times[i] , static_cast<double>(i) );
144 
145 }
146 
147 
BOOST_AUTO_TEST_CASE(test_integrate_times_ranges)148 BOOST_AUTO_TEST_CASE( test_integrate_times_ranges )
149 {
150 
151     state_type x( 3 );
152     x[0] = x[1] = x[2] = 10.0;
153 
154     const value_type dt = 0.03;
155 
156     std::vector< double > times;
157 
158     // simple stepper
159     integrate_times( runge_kutta4< state_type >() , lorenz , x ,
160                 std::make_pair( boost::counting_iterator<int>(0) , boost::counting_iterator<int>(10) ) ,
161                 dt , push_back_time( times ) );
162 
163     for( int i=0 ; i<10 ; ++i )
164         // check if observer was called at times 0,1,2,...
165         BOOST_CHECK_EQUAL( times[i] , static_cast<double>(i) );
166     times.clear();
167 
168     // controlled stepper
169     integrate_times( controlled_runge_kutta< runge_kutta_dopri5< state_type > >() , lorenz , x ,
170                 std::make_pair( boost::counting_iterator<int>(0) , boost::counting_iterator<int>(10) ) ,
171                 dt , push_back_time( times ) );
172 
173     for( int i=0 ; i<10 ; ++i )
174         // check if observer was called at times 0,1,2,...
175         BOOST_CHECK_EQUAL( times[i] , static_cast<double>(i) );
176     times.clear();
177 
178     //another controlled stepper
179     integrate_times( bulirsch_stoer< state_type >() , lorenz , x ,
180                 std::make_pair( boost::counting_iterator<int>(0) , boost::counting_iterator<int>(10) ) ,
181                 dt , push_back_time( times ) );
182 
183     for( int i=0 ; i<10 ; ++i )
184         // check if observer was called at times 0,1,2,...
185         BOOST_CHECK_EQUAL( times[i] , static_cast<double>(i) );
186     times.clear();
187 
188 
189     // dense output stepper
190     integrate_times( bulirsch_stoer_dense_out< state_type >() , lorenz , x ,
191                 std::make_pair( boost::counting_iterator<int>(0) , boost::counting_iterator<int>(10) ) ,
192                 dt , push_back_time( times ) );
193 
194     for( int i=0 ; i<10 ; ++i )
195         // check if observer was called at times 0,1,2,...
196         BOOST_CHECK_EQUAL( times[i] , static_cast<double>(i) );
197 
198 }
199 
BOOST_AUTO_TEST_CASE(test_integrate_times_overshoot)200 BOOST_AUTO_TEST_CASE( test_integrate_times_overshoot )
201 {
202     state_type x( 3 );
203     x[0] = x[1] = x[2] = 10.0;
204     double dt = -0.1;
205 
206     std::vector<double> times( 10 );
207     for( int i=0 ; i<10 ; ++i )
208             times[i] = 1.0-i*1.0/9.0;
209 
210     std::cout << "test rk4 stepper" << std::endl;
211     // simple stepper
212     std::vector<double> obs_times;
213     int steps = integrate_times( runge_kutta4< state_type >() , lorenz , x ,
214                                  times.begin() , times.end() ,
215                                  dt , push_back_time( obs_times ) );
216 // different behavior for the iterator based integrate implementaton
217 #ifndef ODEINT_INTEGRATE_ITERATOR
218     BOOST_CHECK_EQUAL( steps , 18 ); // we really need 18 steps because dt and
219                                      // the difference of the observation times
220                                      // are so out of sync
221 #else
222     // iterator based implementation can only return the number of iteration steps
223     BOOST_CHECK_EQUAL( steps , 9 );
224 #endif
225     for( int i=0 ; i<10 ; ++i )
226         BOOST_CHECK_EQUAL( times[i] , obs_times[i] );
227 
228     std::cout << "test rk_ck stepper" << std::endl;
229     // controlled stepper
230     obs_times.clear();
231     integrate_times( controlled_runge_kutta< runge_kutta_cash_karp54< state_type > >() , lorenz , x ,
232                      times.begin() , times.end() ,
233                      dt , push_back_time( obs_times ) );
234     for( int i=0 ; i<10 ; ++i )
235         BOOST_CHECK_EQUAL( times[i] , obs_times[i] );
236 
237     std::cout << "test dopri5 stepper" << std::endl;
238     // controlled stepper
239     obs_times.clear();
240     integrate_times( dense_output_runge_kutta< controlled_runge_kutta< runge_kutta_dopri5< state_type > > >() , lorenz , x ,
241                      times.begin() , times.end() ,
242                      dt , push_back_time( obs_times ) );
243     for( int i=0 ; i<10 ; ++i )
244         BOOST_CHECK_EQUAL( times[i] , obs_times[i] );
245 }
246 
247 BOOST_AUTO_TEST_SUITE_END()
248