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 <boost/array.hpp>
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 boost::array< double , 2 > state_type;
38
39 // harmonic oscillator, analytic solution x[0] = sin( t )
40 struct osc
41 {
operator ()osc42 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 >
resetresetter54 static void reset( Stepper &stepper ) { }
55 };
56
57 template< >
58 struct resetter< explicit_error_stepper_fsal_tag >
59 {
60 template< class Stepper >
resetresetter61 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 {
operator ()perform_runge_kutta_test73 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 = {{ 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( sin(dt) - x1[0] ) / std::pow( dt , 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( dt , 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 {
operator ()perform_runge_kutta_error_test108 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 = {{ 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_err[0] ) / std::pow( dt , 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( dt , 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
BOOST_AUTO_TEST_CASE_TEMPLATE(runge_kutta_test,Stepper,runge_kutta_steppers)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
BOOST_AUTO_TEST_CASE_TEMPLATE(runge_kutta_error_test,Stepper,runge_kutta_error_steppers)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