1 | /* Boost numeric test of the adams-bashforth-moulton steppers test file |
2 | |
3 | Copyright 2013 Karsten Ahnert |
4 | Copyright 2013 Mario Mulansky |
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_adams_bashforth_moulton |
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 | |
30 | using namespace boost::unit_test; |
31 | using namespace boost::numeric::odeint; |
32 | namespace mpl = boost::mpl; |
33 | |
34 | typedef double value_type; |
35 | |
36 | typedef std::array< double , 2 > state_type; |
37 | typedef runge_kutta_fehlberg78<state_type> initializing_stepper; |
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 | BOOST_AUTO_TEST_SUITE( numeric_adams_bashforth_moulton_test ) |
50 | |
51 | |
52 | /* generic test for all adams bashforth moulton steppers */ |
53 | template< class Stepper > |
54 | struct perform_adams_bashforth_moulton_test |
55 | { |
56 | void operator()( void ) |
57 | { |
58 | Stepper stepper; |
59 | initializing_stepper init_stepper; |
60 | const int o = stepper.order()+1; //order of the error is order of approximation + 1 |
61 | |
62 | const state_type x0 = {._M_elems: { 0.0 , 1.0 }}; |
63 | state_type x1 = x0; |
64 | double t = 0.0; |
65 | double dt = 0.25; |
66 | // initialization, does a number of steps already to fill internal buffer, t is increased |
67 | // we use the rk78 as initializing stepper |
68 | stepper.initialize( boost::ref(t&: init_stepper) , osc() , x1 , t , dt ); |
69 | // do a number of steps to fill the buffer with results from adams bashforth |
70 | for( size_t n=0 ; n < stepper.steps ; ++n ) |
71 | { |
72 | stepper.do_step( osc() , x1 , t , dt ); |
73 | t += dt; |
74 | } |
75 | double A = std::sqrt( x: x1[0]*x1[0] + x1[1]*x1[1] ); |
76 | double phi = std::asin(x: x1[0]/A) - t; |
77 | // now we do the actual step |
78 | stepper.do_step( osc() , x1 , t , dt ); |
79 | // only examine the error of the adams-bashforth-moulton step, not the initialization |
80 | const double f = 2.0 * std::abs( x: A*sin(x: t+dt+phi) - x1[0] ) / std::pow( x: dt , y: o ); // upper bound |
81 | |
82 | std::cout << o << " , " << f << std::endl; |
83 | |
84 | /* as long as we have errors above machine precision */ |
85 | while( f*std::pow( x: dt , y: o ) > 1E-16 ) |
86 | { |
87 | x1 = x0; |
88 | t = 0.0; |
89 | stepper.initialize( boost::ref(t&: init_stepper) , osc() , x1 , t , dt ); |
90 | A = std::sqrt( x: x1[0]*x1[0] + x1[1]*x1[1] ); |
91 | phi = std::asin(x: x1[0]/A) - t; |
92 | // now we do the actual step |
93 | stepper.do_step( osc() , x1 , t , dt ); |
94 | // only examine the error of the adams-bashforth-moulton step, not the initialization |
95 | std::cout << "Testing dt=" << dt << " , " << std::abs( x: A*sin(x: t+dt+phi) - x1[0] ) << std::endl; |
96 | BOOST_CHECK_LT( std::abs( A*sin(t+dt+phi) - x1[0] ) , f*std::pow( dt , o ) ); |
97 | dt *= 0.5; |
98 | } |
99 | } |
100 | }; |
101 | |
102 | typedef mpl::vector< |
103 | adams_bashforth_moulton< 1 , state_type > , |
104 | adams_bashforth_moulton< 2 , state_type > , |
105 | adams_bashforth_moulton< 3 , state_type > , |
106 | adams_bashforth_moulton< 4 , state_type > , |
107 | adams_bashforth_moulton< 5 , state_type > , |
108 | adams_bashforth_moulton< 6 , state_type > , |
109 | adams_bashforth_moulton< 7 , state_type > , |
110 | adams_bashforth_moulton< 8 , state_type > |
111 | > adams_bashforth_moulton_steppers; |
112 | |
113 | BOOST_AUTO_TEST_CASE_TEMPLATE( adams_bashforth_moulton_test , Stepper, adams_bashforth_moulton_steppers ) |
114 | { |
115 | perform_adams_bashforth_moulton_test< Stepper > tester; |
116 | tester(); |
117 | } |
118 | |
119 | BOOST_AUTO_TEST_SUITE_END() |
120 | |