1 | /* Boost numeric test for orders of quadrature formulas test file |
2 | |
3 | Copyright 2015 Gregor de Cillia |
4 | Copyright 2015 Mario Mulansky <mario.mulansky@gmx.net> |
5 | |
6 | Distributed under the Boost Software License, Version 1.0. |
7 | (See accompanying file LICENSE_1_0.txt or |
8 | copy at http://www.boost.org/LICENSE_1_0.txt) |
9 | */ |
10 | |
11 | // disable checked iterator warning for msvc |
12 | #include <boost/config.hpp> |
13 | |
14 | #ifdef BOOST_MSVC |
15 | #pragma warning(disable:4996) |
16 | #endif |
17 | |
18 | #define BOOST_TEST_MODULE order_quadrature_formula |
19 | |
20 | #include <iostream> |
21 | #include <cmath> |
22 | |
23 | #include <boost/test/unit_test.hpp> |
24 | |
25 | #include <boost/mpl/vector.hpp> |
26 | |
27 | #include <boost/numeric/odeint.hpp> |
28 | |
29 | #include <boost/numeric/ublas/vector.hpp> |
30 | |
31 | using namespace boost::unit_test; |
32 | using namespace boost::numeric::odeint; |
33 | namespace mpl = boost::mpl; |
34 | |
35 | typedef double value_type; |
36 | typedef value_type time_type; |
37 | typedef value_type state_type; |
38 | |
39 | BOOST_AUTO_TEST_SUITE( order_of_convergence_test ) |
40 | |
41 | /* defines the simple monomial f(t) = (p+1) * t^p.*/ |
42 | struct monomial |
43 | { |
44 | int power; |
45 | |
46 | monomial(int p = 0) : power( p ){}; |
47 | |
48 | void operator()( const state_type &x , state_type &dxdt , const time_type t ) |
49 | { |
50 | dxdt = ( 1.0 + power ) * pow( x: t, y: power ); |
51 | } |
52 | }; |
53 | |
54 | |
55 | /* generic test for all steppers that support integrate_const */ |
56 | template< class Stepper > |
57 | struct stepper_order_test |
58 | { |
59 | void operator()( int steps = 1 ) |
60 | { |
61 | const int estimated_order = estimate_order( steps ); |
62 | const int defined_order = Stepper::order_value; |
63 | |
64 | std::cout << boost::format( "%-20i%-20i\n" ) |
65 | % estimated_order % defined_order; |
66 | |
67 | BOOST_REQUIRE_EQUAL( estimated_order, defined_order ); |
68 | } |
69 | |
70 | /* |
71 | the order of the stepper is estimated by trying to solve the ODE |
72 | x'(t) = (p+1) * t^p |
73 | until the errors are too big to be justified by finite precision. |
74 | the first value p for which the problem is *not* solved within the |
75 | finite precision tolerance is the estimate for the order of the scheme. |
76 | */ |
77 | int estimate_order( int steps ) |
78 | { |
79 | const double dt = 1.0/steps; |
80 | const double tolerance = steps*1E-15; |
81 | int p; |
82 | for( p = 0; true; p++ ) |
83 | { |
84 | // begin with x'(t) = t^0 = 1 |
85 | // => x (t) = t |
86 | // then use x'(t) = 2*t^1 |
87 | // => x (t) = t^2 |
88 | // ... |
89 | state_type x = 0.0; |
90 | |
91 | double t = integrate_n_steps( Stepper(), monomial( p ), x, 0.0, dt, |
92 | steps ); |
93 | if( fabs( x: x - pow( x: t, y: ( 1.0 + p ) ) ) > tolerance ) |
94 | break; |
95 | } |
96 | // the smallest power p for which the test failed is the estimated order, |
97 | // as the solution for this power is x(t) = t^{p+1} |
98 | return p; |
99 | } |
100 | }; |
101 | |
102 | |
103 | typedef mpl::vector< |
104 | euler< state_type > , |
105 | modified_midpoint< state_type > , |
106 | runge_kutta4< state_type > , |
107 | runge_kutta4_classic< state_type > , |
108 | runge_kutta_cash_karp54_classic< state_type > , |
109 | runge_kutta_cash_karp54< state_type > , |
110 | runge_kutta_dopri5< state_type > , |
111 | runge_kutta_fehlberg78< state_type > |
112 | > runge_kutta_steppers; |
113 | |
114 | typedef mpl::vector< |
115 | adams_bashforth< 2, state_type, double, state_type, double, |
116 | vector_space_algebra, default_operations, |
117 | initially_resizer, runge_kutta_fehlberg78< state_type > >, |
118 | adams_bashforth< 3, state_type, double, state_type, double, |
119 | vector_space_algebra, default_operations, |
120 | initially_resizer, runge_kutta_fehlberg78< state_type > >, |
121 | adams_bashforth< 4, state_type, double, state_type, double, |
122 | vector_space_algebra, default_operations, |
123 | initially_resizer, runge_kutta_fehlberg78< state_type > >, |
124 | adams_bashforth< 5, state_type, double, state_type, double, |
125 | vector_space_algebra, default_operations, |
126 | initially_resizer, runge_kutta_fehlberg78< state_type > >, |
127 | adams_bashforth< 6, state_type, double, state_type, double, |
128 | vector_space_algebra, default_operations, |
129 | initially_resizer, runge_kutta_fehlberg78< state_type > >, |
130 | adams_bashforth< 7, state_type, double, state_type, double, |
131 | vector_space_algebra, default_operations, |
132 | initially_resizer, runge_kutta_fehlberg78< state_type > >, |
133 | adams_bashforth< 8, state_type, double, state_type, double, |
134 | vector_space_algebra, default_operations, |
135 | initially_resizer, runge_kutta_fehlberg78< state_type > > |
136 | > ab_steppers; |
137 | |
138 | |
139 | typedef mpl::vector< |
140 | adams_bashforth_moulton< 2, state_type, double, state_type, double, |
141 | vector_space_algebra, default_operations, |
142 | initially_resizer, |
143 | runge_kutta_fehlberg78< state_type > >, |
144 | adams_bashforth_moulton< 3, state_type, double, state_type, double, |
145 | vector_space_algebra, default_operations, |
146 | initially_resizer, |
147 | runge_kutta_fehlberg78< state_type > >, |
148 | adams_bashforth_moulton< 4, state_type, double, state_type, double, |
149 | vector_space_algebra, default_operations, |
150 | initially_resizer, |
151 | runge_kutta_fehlberg78< state_type > >, |
152 | adams_bashforth_moulton< 5, state_type, double, state_type, double, |
153 | vector_space_algebra, default_operations, |
154 | initially_resizer, |
155 | runge_kutta_fehlberg78< state_type > >, |
156 | adams_bashforth_moulton< 6, state_type, double, state_type, double, |
157 | vector_space_algebra, default_operations, |
158 | initially_resizer, |
159 | runge_kutta_fehlberg78< state_type > >, |
160 | adams_bashforth_moulton< 7, state_type, double, state_type, double, |
161 | vector_space_algebra, default_operations, |
162 | initially_resizer, |
163 | runge_kutta_fehlberg78< state_type > >, |
164 | adams_bashforth_moulton< 8, state_type, double, state_type, double, |
165 | vector_space_algebra, default_operations, |
166 | initially_resizer, |
167 | runge_kutta_fehlberg78< state_type > > |
168 | > abm_steppers; |
169 | |
170 | |
171 | BOOST_AUTO_TEST_CASE_TEMPLATE( runge_kutta_test , Stepper, runge_kutta_steppers ) |
172 | { |
173 | stepper_order_test< Stepper > tester; |
174 | tester(10); |
175 | } |
176 | |
177 | |
178 | BOOST_AUTO_TEST_CASE_TEMPLATE( adams_bashforth_test , Stepper, ab_steppers ) |
179 | { |
180 | stepper_order_test< Stepper > tester; |
181 | tester(16); |
182 | } |
183 | |
184 | |
185 | BOOST_AUTO_TEST_CASE_TEMPLATE( adams_bashforth_moultion_test , Stepper, abm_steppers ) |
186 | { |
187 | stepper_order_test< Stepper > tester; |
188 | tester(16); |
189 | } |
190 | |
191 | BOOST_AUTO_TEST_SUITE_END() |
192 | |