1 | /* Boost numeric test of the runge kutta steppers test file |
2 | |
3 | Copyright 2012 Mario Mulansky |
4 | Copyright 2012 Karsten Ahnert |
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 | #ifdef BOOST_MSVC |
14 | #pragma warning(disable:4996) |
15 | #endif |
16 | |
17 | #define BOOST_TEST_MODULE numeric_runge_kutta |
18 | |
19 | #include <iostream> |
20 | #include <cmath> |
21 | |
22 | #include <array> |
23 | |
24 | #include <boost/test/unit_test.hpp> |
25 | |
26 | #include <boost/mpl/vector.hpp> |
27 | |
28 | #include <boost/numeric/odeint.hpp> |
29 | #include <boost/numeric/odeint/stepper/extrapolation_stepper.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 | |
37 | typedef std::array< double , 2 > state_type; |
38 | |
39 | // harmonic oscillator, analytic solution x[0] = sin( t ) |
40 | struct osc |
41 | { |
42 | void operator()( const state_type &x , state_type &dxdt , const double t ) const |
43 | { |
44 | dxdt[0] = x[1]; |
45 | dxdt[1] = -x[0]; |
46 | } |
47 | }; |
48 | |
49 | /* reset dispatcher */ |
50 | template< class StepperCategory > |
51 | struct resetter |
52 | { |
53 | template< class Stepper > |
54 | static void reset( Stepper &stepper ) { } |
55 | }; |
56 | |
57 | template< > |
58 | struct resetter< explicit_error_stepper_fsal_tag > |
59 | { |
60 | template< class Stepper > |
61 | static void reset( Stepper &stepper ) |
62 | { stepper.reset(); } |
63 | }; |
64 | |
65 | |
66 | BOOST_AUTO_TEST_SUITE( numeric_runge_kutta_test ) |
67 | |
68 | |
69 | /* generic test for all runge kutta steppers */ |
70 | template< class Stepper > |
71 | struct perform_runge_kutta_test |
72 | { |
73 | void operator()( void ) |
74 | { |
75 | |
76 | Stepper stepper; |
77 | const int o = stepper.order()+1; //order of the error is order of approximation + 1 |
78 | |
79 | const state_type x0 = {._M_elems: { 0.0 , 1.0 }}; |
80 | state_type x1; |
81 | const double t = 0.0; |
82 | /* do a first step with dt=0.1 to get an estimate on the prefactor of the error dx = f * dt^(order+1) */ |
83 | double dt = 0.5; |
84 | stepper.do_step( osc() , x0 , t , x1 , dt ); |
85 | const double f = 2.0 * std::abs( x: sin(x: dt) - x1[0] ) / std::pow( x: dt , y: o ); // upper bound |
86 | |
87 | std::cout << o << " , " << f << std::endl; |
88 | |
89 | /* as long as we have errors above machine precision */ |
90 | while( f*std::pow( x: dt , y: o ) > 1E-16 ) |
91 | { |
92 | // reset stepper which require resetting (fsal steppers) |
93 | resetter< typename Stepper::stepper_category >::reset( stepper ); |
94 | |
95 | stepper.do_step( osc() , x0 , t , x1 , dt ); |
96 | std::cout << "Testing dt=" << dt << std::endl; |
97 | BOOST_CHECK_LT( std::abs( sin(dt) - x1[0] ) , f*std::pow( dt , o ) ); |
98 | dt *= 0.5; |
99 | } |
100 | } |
101 | }; |
102 | |
103 | |
104 | /* generic error test for all runge kutta steppers */ |
105 | template< class Stepper > |
106 | struct perform_runge_kutta_error_test |
107 | { |
108 | void operator()( void ) |
109 | { |
110 | Stepper stepper; |
111 | const int o = stepper.error_order()+1; //order of the error is order of approximation + 1 |
112 | |
113 | const state_type x0 = {._M_elems: { 0.0 , 1.0 }}; |
114 | state_type x1 , x_err; |
115 | const double t = 0.0; |
116 | /* do a first step with dt=0.1 to get an estimate on the prefactor of the error dx = f * dt^(order+1) */ |
117 | double dt = 0.5; |
118 | stepper.do_step( osc() , x0 , t , x1 , dt , x_err ); |
119 | const double f = 2.0 * std::abs( x: x_err[0] ) / std::pow( x: dt , y: o ); |
120 | |
121 | std::cout << o << " , " << f << " , " << x0[0] << std::endl; |
122 | |
123 | /* as long as we have errors above machine precision */ |
124 | while( f*std::pow( x: dt , y: o ) > 1E-16 ) |
125 | { |
126 | // reset stepper which require resetting (fsal steppers) |
127 | resetter< typename Stepper::stepper_category >::reset( stepper ); |
128 | |
129 | stepper.do_step( osc() , x0 , t , x1 , dt , x_err ); |
130 | std::cout << "Testing dt=" << dt << ": " << x_err[1] << std::endl; |
131 | BOOST_CHECK_SMALL( std::abs( x_err[0] ) , f*std::pow( dt , o ) ); |
132 | dt *= 0.5; |
133 | } |
134 | } |
135 | }; |
136 | |
137 | |
138 | typedef mpl::vector< |
139 | euler< state_type > , |
140 | modified_midpoint< state_type > , |
141 | runge_kutta4< state_type > , |
142 | runge_kutta4_classic< state_type > , |
143 | runge_kutta_cash_karp54_classic< state_type > , |
144 | runge_kutta_cash_karp54< state_type > , |
145 | runge_kutta_dopri5< state_type > , |
146 | runge_kutta_fehlberg78< state_type > , |
147 | extrapolation_stepper< 4, state_type > , |
148 | extrapolation_stepper< 6, state_type > , |
149 | extrapolation_stepper< 8, state_type > , |
150 | extrapolation_stepper< 10, state_type > |
151 | > runge_kutta_steppers; |
152 | |
153 | BOOST_AUTO_TEST_CASE_TEMPLATE( runge_kutta_test , Stepper, runge_kutta_steppers ) |
154 | { |
155 | perform_runge_kutta_test< Stepper > tester; |
156 | tester(); |
157 | } |
158 | |
159 | |
160 | typedef mpl::vector< |
161 | runge_kutta_cash_karp54_classic< state_type > , |
162 | runge_kutta_cash_karp54< state_type > , |
163 | runge_kutta_dopri5< state_type > , |
164 | runge_kutta_fehlberg78< state_type > , |
165 | extrapolation_stepper< 4, state_type > , |
166 | extrapolation_stepper< 6, state_type > , |
167 | extrapolation_stepper< 8, state_type > , |
168 | extrapolation_stepper< 10, state_type > |
169 | > runge_kutta_error_steppers; |
170 | |
171 | BOOST_AUTO_TEST_CASE_TEMPLATE( runge_kutta_error_test , Stepper, runge_kutta_error_steppers ) |
172 | { |
173 | perform_runge_kutta_error_test< Stepper > tester; |
174 | tester(); |
175 | } |
176 | |
177 | BOOST_AUTO_TEST_SUITE_END() |
178 | |