1 /*
2 [auto_generated]
3 libs/numeric/odeint/test/step_size_limitation.cpp
4
5 [begin_description]
6 Tests the step size limitation functionality
7 [end_description]
8
9 Copyright 2015 Mario Mulansky
10
11 Distributed under the Boost Software License, Version 1.0.
12 (See accompanying file LICENSE_1_0.txt or
13 copy at http://www.boost.org/LICENSE_1_0.txt)
14 */
15
16 #define BOOST_TEST_MODULE odeint_integrate_times
17
18 #include <boost/test/unit_test.hpp>
19
20 #include <utility>
21 #include <iostream>
22 #include <vector>
23
24 #include <boost/numeric/odeint.hpp>
25
26 using namespace boost::unit_test;
27 using namespace boost::numeric::odeint;
28
29 typedef double value_type;
30 typedef std::vector< value_type > state_type;
31
32 /***********************************************
33 * first part of the tests: explicit methods
34 ***********************************************
35 */
36
damped_osc(const state_type & x,state_type & dxdt,const value_type t)37 void damped_osc( const state_type &x , state_type &dxdt , const value_type t )
38 {
39 const value_type gam( 0.1);
40
41 dxdt[0] = x[1];
42 dxdt[1] = -x[0] - gam*x[1];
43 }
44
45
46 struct push_back_time
47 {
48 std::vector< double >& m_times;
49
push_back_timepush_back_time50 push_back_time( std::vector< double > × )
51 : m_times( times ) { }
52
53 template<typename State>
operator ()push_back_time54 void operator()( const State &x , double t )
55 {
56 m_times.push_back( t );
57 }
58 };
59
60 BOOST_AUTO_TEST_SUITE( step_size_limitation_test )
61
BOOST_AUTO_TEST_CASE(test_step_adjuster)62 BOOST_AUTO_TEST_CASE( test_step_adjuster )
63 {
64 // first use adjuster without step size limitation
65 default_step_adjuster<double, double> step_adjuster;
66 const double dt = 0.1;
67 double dt_new = step_adjuster.decrease_step(dt, 1.5, 2);
68 BOOST_CHECK(dt_new < dt*2.0/3.0);
69
70 dt_new = step_adjuster.increase_step(dt, 0.8, 1);
71 // for errors > 0.5 no increase is performed
72 BOOST_CHECK(dt_new == dt);
73
74 dt_new = step_adjuster.increase_step(dt, 0.4, 1);
75 // smaller errors should lead to step size increase
76 std::cout << dt_new << std::endl;
77 BOOST_CHECK(dt_new > dt);
78
79
80 // now test with step size limitation max_dt = 0.1
81 default_step_adjuster<double, double>
82 limited_adjuster(dt);
83
84 dt_new = limited_adjuster.decrease_step(dt, 1.5, 2);
85 // decreasing works as before
86 BOOST_CHECK(dt_new < dt*2.0/3.0);
87
88 dt_new = limited_adjuster.decrease_step(2*dt, 1.1, 2);
89 // decreasing a large step size should give max_dt
90 BOOST_CHECK(dt_new == dt);
91
92 dt_new = limited_adjuster.increase_step(dt, 0.8, 1);
93 // for errors > 0.5 no increase is performed, still valid
94 BOOST_CHECK(dt_new == dt);
95
96 dt_new = limited_adjuster.increase_step(dt, 0.4, 1);
97 // but even for smaller errors, we should at most get 0.1
98 BOOST_CHECK_EQUAL(dt_new, dt);
99
100 dt_new = limited_adjuster.increase_step(0.9*dt, 0.1, 1);
101 std::cout << dt_new << std::endl;
102 // check that we don't increase beyond max_dt
103 BOOST_CHECK(dt_new == dt);
104 }
105
106
107 template<class Stepper>
test_explicit_stepper(Stepper stepper,const double max_dt)108 void test_explicit_stepper(Stepper stepper, const double max_dt)
109 {
110 state_type x( 2 );
111 x[0] = x[1] = 10.0;
112 const size_t steps = 100;
113
114 std::vector<double> times;
115
116 integrate_adaptive(stepper, damped_osc, x, 0.0, steps*max_dt, max_dt, push_back_time(times));
117
118 BOOST_CHECK_EQUAL(times.size(), steps+1);
119 // check that dt remains at exactly max_dt
120 for( size_t i=0 ; i<times.size() ; ++i )
121 // check if observer was called at times 0,1,2,...
122 BOOST_CHECK_SMALL( times[i] - static_cast<double>(i)*max_dt , 1E-15);
123 times.clear();
124
125 // this should also work when we provide some bigger initial dt
126 integrate_adaptive(stepper, damped_osc, x, 0.0, steps*max_dt, 10*max_dt, push_back_time(times));
127
128 BOOST_CHECK_EQUAL(times.size(), steps+1);
129 // check that dt remains at exactly max_dt
130 for( size_t i=0 ; i<times.size() ; ++i )
131 // check if observer was called at times 0,1,2,...
132 BOOST_CHECK_SMALL( times[i] - static_cast<double>(i)*max_dt , 1E-15);
133 times.clear();
134 }
135
136
BOOST_AUTO_TEST_CASE(test_controlled)137 BOOST_AUTO_TEST_CASE( test_controlled )
138 {
139 const double max_dt = 0.01;
140
141 test_explicit_stepper(make_controlled(1E-2, 1E-2, max_dt,
142 runge_kutta_dopri5<state_type>()),
143 max_dt);
144 test_explicit_stepper(make_controlled(1E-2, 1E-2, -max_dt,
145 runge_kutta_dopri5<state_type>()),
146 -max_dt);
147
148 test_explicit_stepper(make_controlled(1E-2, 1E-2, max_dt,
149 runge_kutta_cash_karp54<state_type>()),
150 max_dt);
151 test_explicit_stepper(make_controlled(1E-2, 1E-2, -max_dt,
152 runge_kutta_cash_karp54<state_type>()),
153 -max_dt);
154
155 test_explicit_stepper(bulirsch_stoer<state_type>(1E-2, 1E-2, 1.0, 1.0, max_dt),
156 max_dt);
157 test_explicit_stepper(bulirsch_stoer<state_type>(1E-2, 1E-2, 1.0, 1.0, -max_dt),
158 -max_dt);
159 }
160
161
BOOST_AUTO_TEST_CASE(test_dense_out)162 BOOST_AUTO_TEST_CASE( test_dense_out )
163 {
164 const double max_dt = 0.01;
165 test_explicit_stepper(make_dense_output(1E-2, 1E-2, max_dt,
166 runge_kutta_dopri5<state_type>()),
167 max_dt);
168 test_explicit_stepper(make_dense_output(1E-2, 1E-2, -max_dt,
169 runge_kutta_dopri5<state_type>()),
170 -max_dt);
171
172 test_explicit_stepper(bulirsch_stoer_dense_out<state_type>(1E-2, 1E-2, 1, 1, max_dt),
173 max_dt);
174
175 test_explicit_stepper(bulirsch_stoer_dense_out<state_type>(1E-2, 1E-2, 1, 1, -max_dt),
176 -max_dt);
177 }
178
179
180 /***********************************************
181 * second part of the tests: implicit Rosenbrock
182 ***********************************************
183 */
184
185 typedef boost::numeric::ublas::vector< value_type > vector_type;
186 typedef boost::numeric::ublas::matrix< value_type > matrix_type;
187
188
189 // harmonic oscillator, analytic solution x[0] = sin( t )
190 struct osc_rhs
191 {
operator ()osc_rhs192 void operator()( const vector_type &x , vector_type &dxdt , const value_type &t ) const
193 {
194 dxdt( 0 ) = x( 1 );
195 dxdt( 1 ) = -x( 0 );
196 }
197 };
198
199 struct osc_jacobi
200 {
operator ()osc_jacobi201 void operator()( const vector_type &x , matrix_type &jacobi , const value_type &t , vector_type &dfdt ) const
202 {
203 jacobi( 0 , 0 ) = 0;
204 jacobi( 0 , 1 ) = 1;
205 jacobi( 1 , 0 ) = -1;
206 jacobi( 1 , 1 ) = 0;
207 dfdt( 0 ) = 0.0;
208 dfdt( 1 ) = 0.0;
209 }
210 };
211
212
213 template<class Stepper>
test_rosenbrock_stepper(Stepper stepper,const double max_dt)214 void test_rosenbrock_stepper(Stepper stepper, const double max_dt)
215 {
216 vector_type x( 2 );
217 x(0) = x(1) = 10.0;
218 const size_t steps = 100;
219
220 std::vector<double> times;
221
222 integrate_adaptive(stepper,
223 std::make_pair(osc_rhs(), osc_jacobi()),
224 x, 0.0, steps*max_dt, max_dt, push_back_time(times));
225
226 BOOST_CHECK_EQUAL(times.size(), steps+1);
227 // check that dt remains at exactly max_dt
228 for( size_t i=0 ; i<times.size() ; ++i )
229 // check if observer was called at times 0,1,2,...
230 BOOST_CHECK_SMALL( times[i] - static_cast<double>(i)*max_dt , 1E-15);
231 times.clear();
232
233 // this should also work when we provide some bigger initial dt
234 integrate_adaptive(stepper,
235 std::make_pair(osc_rhs(), osc_jacobi()),
236 x, 0.0, steps*max_dt, 10*max_dt, push_back_time(times));
237
238 BOOST_CHECK_EQUAL(times.size(), steps+1);
239 // check that dt remains at exactly max_dt
240 for( size_t i=0 ; i<times.size() ; ++i )
241 // check if observer was called at times 0,1,2,...
242 BOOST_CHECK_SMALL( times[i] - static_cast<double>(i)*max_dt , 1E-15);
243 times.clear();
244 }
245
246
BOOST_AUTO_TEST_CASE(test_controlled_rosenbrock)247 BOOST_AUTO_TEST_CASE( test_controlled_rosenbrock )
248 {
249 const double max_dt = 0.01;
250
251 test_rosenbrock_stepper(make_controlled(1E-2, 1E-2, max_dt, rosenbrock4<value_type>()),
252 max_dt);
253 test_rosenbrock_stepper(make_controlled(1E-2, 1E-2, -max_dt, rosenbrock4<value_type>()),
254 -max_dt);
255 }
256
257
BOOST_AUTO_TEST_CASE(test_dense_out_rosenbrock)258 BOOST_AUTO_TEST_CASE( test_dense_out_rosenbrock )
259 {
260 const double max_dt = 0.01;
261
262 test_rosenbrock_stepper(make_dense_output(1E-2, 1E-2, max_dt, rosenbrock4<value_type>()),
263 max_dt);
264 test_rosenbrock_stepper(make_dense_output(1E-2, 1E-2, -max_dt, rosenbrock4<value_type>()),
265 -max_dt);
266 }
267
268 BOOST_AUTO_TEST_SUITE_END()
269