]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 | ||
b32b8144 FG |
23 | #include "boost/format.hpp" |
24 | ||
7c673cae FG |
25 | #include <boost/test/unit_test.hpp> |
26 | ||
27 | #include <boost/mpl/vector.hpp> | |
28 | ||
29 | #include <boost/numeric/odeint.hpp> | |
30 | ||
31 | #include <boost/numeric/ublas/vector.hpp> | |
32 | ||
b32b8144 FG |
33 | #include <boost/format.hpp> |
34 | ||
7c673cae FG |
35 | using namespace boost::unit_test; |
36 | using namespace boost::numeric::odeint; | |
37 | namespace mpl = boost::mpl; | |
38 | ||
39 | typedef double value_type; | |
40 | typedef value_type time_type; | |
41 | typedef value_type state_type; | |
42 | ||
43 | BOOST_AUTO_TEST_SUITE( order_of_convergence_test ) | |
44 | ||
45 | /* defines the simple monomial f(t) = (p+1) * t^p.*/ | |
46 | struct monomial | |
47 | { | |
48 | int power; | |
49 | ||
50 | monomial(int p = 0) : power( p ){}; | |
51 | ||
52 | void operator()( const state_type &x , state_type &dxdt , const time_type t ) | |
53 | { | |
54 | dxdt = ( 1.0 + power ) * pow( t, power ); | |
55 | } | |
56 | }; | |
57 | ||
58 | ||
59 | /* generic test for all steppers that support integrate_const */ | |
60 | template< class Stepper > | |
61 | struct stepper_order_test | |
62 | { | |
63 | void operator()( int steps = 1 ) | |
64 | { | |
65 | const int estimated_order = estimate_order( steps ); | |
66 | const int defined_order = Stepper::order_value; | |
67 | ||
68 | std::cout << boost::format( "%-20i%-20i\n" ) | |
69 | % estimated_order % defined_order; | |
70 | ||
71 | BOOST_REQUIRE_EQUAL( estimated_order, defined_order ); | |
72 | } | |
73 | ||
74 | /* | |
75 | the order of the stepper is estimated by trying to solve the ODE | |
76 | x'(t) = (p+1) * t^p | |
77 | until the errors are too big to be justified by finite precision. | |
78 | the first value p for which the problem is *not* solved within the | |
79 | finite precision tolerance is the estimate for the order of the scheme. | |
80 | */ | |
81 | int estimate_order( int steps ) | |
82 | { | |
83 | const double dt = 1.0/steps; | |
84 | const double tolerance = steps*1E-15; | |
85 | int p; | |
86 | for( p = 0; true; p++ ) | |
87 | { | |
88 | // begin with x'(t) = t^0 = 1 | |
89 | // => x (t) = t | |
90 | // then use x'(t) = 2*t^1 | |
91 | // => x (t) = t^2 | |
92 | // ... | |
93 | state_type x = 0.0; | |
94 | ||
95 | double t = integrate_n_steps( Stepper(), monomial( p ), x, 0.0, dt, | |
96 | steps ); | |
97 | if( fabs( x - pow( t, ( 1.0 + p ) ) ) > tolerance ) | |
98 | break; | |
99 | } | |
100 | // the smallest power p for which the test failed is the estimated order, | |
101 | // as the solution for this power is x(t) = t^{p+1} | |
102 | return p; | |
103 | } | |
104 | }; | |
105 | ||
106 | ||
107 | typedef mpl::vector< | |
108 | euler< state_type > , | |
109 | modified_midpoint< state_type > , | |
110 | runge_kutta4< state_type > , | |
111 | runge_kutta4_classic< state_type > , | |
112 | runge_kutta_cash_karp54_classic< state_type > , | |
113 | runge_kutta_cash_karp54< state_type > , | |
114 | runge_kutta_dopri5< state_type > , | |
115 | runge_kutta_fehlberg78< state_type > | |
116 | > runge_kutta_steppers; | |
117 | ||
118 | typedef mpl::vector< | |
119 | adams_bashforth< 2, state_type, double, state_type, double, | |
120 | vector_space_algebra, default_operations, | |
121 | initially_resizer, runge_kutta_fehlberg78< state_type > >, | |
122 | adams_bashforth< 3, state_type, double, state_type, double, | |
123 | vector_space_algebra, default_operations, | |
124 | initially_resizer, runge_kutta_fehlberg78< state_type > >, | |
125 | adams_bashforth< 4, state_type, double, state_type, double, | |
126 | vector_space_algebra, default_operations, | |
127 | initially_resizer, runge_kutta_fehlberg78< state_type > >, | |
128 | adams_bashforth< 5, state_type, double, state_type, double, | |
129 | vector_space_algebra, default_operations, | |
130 | initially_resizer, runge_kutta_fehlberg78< state_type > >, | |
131 | adams_bashforth< 6, state_type, double, state_type, double, | |
132 | vector_space_algebra, default_operations, | |
133 | initially_resizer, runge_kutta_fehlberg78< state_type > >, | |
134 | adams_bashforth< 7, state_type, double, state_type, double, | |
135 | vector_space_algebra, default_operations, | |
136 | initially_resizer, runge_kutta_fehlberg78< state_type > >, | |
137 | adams_bashforth< 8, state_type, double, state_type, double, | |
138 | vector_space_algebra, default_operations, | |
139 | initially_resizer, runge_kutta_fehlberg78< state_type > > | |
140 | > ab_steppers; | |
141 | ||
142 | ||
143 | typedef mpl::vector< | |
144 | adams_bashforth_moulton< 2, 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< 3, 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< 4, 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< 5, 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< 6, 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< 7, state_type, double, state_type, double, | |
165 | vector_space_algebra, default_operations, | |
166 | initially_resizer, | |
167 | runge_kutta_fehlberg78< state_type > >, | |
168 | adams_bashforth_moulton< 8, state_type, double, state_type, double, | |
169 | vector_space_algebra, default_operations, | |
170 | initially_resizer, | |
171 | runge_kutta_fehlberg78< state_type > > | |
172 | > abm_steppers; | |
173 | ||
174 | ||
175 | BOOST_AUTO_TEST_CASE_TEMPLATE( runge_kutta_test , Stepper, runge_kutta_steppers ) | |
176 | { | |
177 | stepper_order_test< Stepper > tester; | |
178 | tester(10); | |
179 | } | |
180 | ||
181 | ||
182 | BOOST_AUTO_TEST_CASE_TEMPLATE( adams_bashforth_test , Stepper, ab_steppers ) | |
183 | { | |
184 | stepper_order_test< Stepper > tester; | |
185 | tester(16); | |
186 | } | |
187 | ||
188 | ||
189 | BOOST_AUTO_TEST_CASE_TEMPLATE( adams_bashforth_moultion_test , Stepper, abm_steppers ) | |
190 | { | |
191 | stepper_order_test< Stepper > tester; | |
192 | tester(16); | |
193 | } | |
194 | ||
195 | BOOST_AUTO_TEST_SUITE_END() |