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
31using namespace boost::unit_test;
32using namespace boost::numeric::odeint;
33namespace mpl = boost::mpl;
34
35typedef double value_type;
36typedef value_type time_type;
37typedef value_type state_type;
38
39BOOST_AUTO_TEST_SUITE( order_of_convergence_test )
40
41/* defines the simple monomial f(t) = (p+1) * t^p.*/
42struct 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 */
56template< class Stepper >
57struct 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
103typedef 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
114typedef 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
139typedef 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
171BOOST_AUTO_TEST_CASE_TEMPLATE( runge_kutta_test , Stepper, runge_kutta_steppers )
172{
173 stepper_order_test< Stepper > tester;
174 tester(10);
175}
176
177
178BOOST_AUTO_TEST_CASE_TEMPLATE( adams_bashforth_test , Stepper, ab_steppers )
179{
180 stepper_order_test< Stepper > tester;
181 tester(16);
182}
183
184
185BOOST_AUTO_TEST_CASE_TEMPLATE( adams_bashforth_moultion_test , Stepper, abm_steppers )
186{
187 stepper_order_test< Stepper > tester;
188 tester(16);
189}
190
191BOOST_AUTO_TEST_SUITE_END()
192

source code of boost/libs/numeric/odeint/test/numeric/order_quadrature_formula.cpp