1 /* Boost numeric test of the adams-bashforth steppers test file
2
3 Copyright 2013 Karsten Ahnert
4 Copyright 2013-2015 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
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
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 boost::array< double , 2 > state_type;
37
38 // harmonic oscillator, analytic solution x[0] = sin( t )
39 struct osc
40 {
operator ()osc41 void operator()( const state_type &x , state_type &dxdt , const double t ) const
42 {
43 dxdt[0] = x[1];
44 dxdt[1] = -x[0];
45 }
46 };
47
48 BOOST_AUTO_TEST_SUITE( numeric_adams_bashforth_test )
49
50
51 /* generic test for all adams bashforth steppers */
52 template< class Stepper >
53 struct perform_adams_bashforth_test
54 {
operator ()perform_adams_bashforth_test55 void operator()( void )
56 {
57 Stepper stepper;
58 const int o = stepper.order()+1; //order of the error is order of approximation + 1
59
60 const state_type x0 = {{ 0.0 , 1.0 }};
61 state_type x1 = x0;
62 double t = 0.0;
63 double dt = 0.2;
64 // initialization, does a number of steps already to fill internal buffer, t is increased
65 stepper.initialize( osc() , x1 , t , dt );
66 double A = std::sqrt( x1[0]*x1[0] + x1[1]*x1[1] );
67 double phi = std::asin(x1[0]/A) - t;
68 // do a number of steps to fill the buffer with results from adams bashforth
69 for( size_t n=0 ; n < stepper.steps ; ++n )
70 {
71 stepper.do_step( osc() , x1 , t , dt );
72 t += dt;
73 }
74 // now we do the actual step
75 stepper.do_step( osc() , x1 , t , dt );
76 // only examine the error of the adams-bashforth step, not the initialization
77 const double f = 2.0 * std::abs( A*sin(t+dt+phi) - x1[0] ) / std::pow( dt , o ); // upper bound
78
79 std::cout << o << " , "
80 << Stepper::initializing_stepper_type::order_value+1 << " , "
81 << f << std::endl;
82
83 /* as long as we have errors above machine precision */
84 while( f*std::pow( dt , o ) > 1E-16 )
85 {
86 x1 = x0;
87 t = 0.0;
88 stepper.initialize( osc() , x1 , t , dt );
89 A = std::sqrt( x1[0]*x1[0] + x1[1]*x1[1] );
90 phi = std::asin(x1[0]/A) - t;
91 // now we do the actual step
92 stepper.do_step( osc() , x1 , t , dt );
93 // only examine the error of the adams-bashforth step, not the initialization
94 std::cout << "Testing dt=" << dt << " , " << std::abs( A*sin(t+dt+phi) - x1[0] ) << std::endl;
95 BOOST_CHECK_LT( std::abs( A*sin(t+dt+phi) - x1[0] ) , f*std::pow( dt , o ) );
96 dt *= 0.5;
97 }
98 }
99 };
100
101 typedef mpl::vector<
102 adams_bashforth< 2 , state_type > ,
103 adams_bashforth< 3 , state_type > ,
104 adams_bashforth< 4 , state_type > ,
105 adams_bashforth< 5 , state_type > ,
106 adams_bashforth< 6 , state_type > ,
107 adams_bashforth< 7 , state_type > ,
108 adams_bashforth< 8 , state_type >
109 > adams_bashforth_steppers;
110
BOOST_AUTO_TEST_CASE_TEMPLATE(adams_bashforth_test,Stepper,adams_bashforth_steppers)111 BOOST_AUTO_TEST_CASE_TEMPLATE( adams_bashforth_test , Stepper, adams_bashforth_steppers )
112 {
113 perform_adams_bashforth_test< Stepper > tester;
114 tester();
115 }
116
117 BOOST_AUTO_TEST_SUITE_END()
118