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 <array>
23
24#include <boost/test/unit_test.hpp>
25
26#include <boost/mpl/vector.hpp>
27
28#include <boost/numeric/odeint.hpp>
29
30using namespace boost::unit_test;
31using namespace boost::numeric::odeint;
32namespace mpl = boost::mpl;
33
34typedef double value_type;
35
36typedef std::array< double , 2 > state_type;
37
38// harmonic oscillator, analytic solution x[0] = sin( t )
39struct osc
40{
41 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
48BOOST_AUTO_TEST_SUITE( numeric_adams_bashforth_test )
49
50
51/* generic test for all adams bashforth steppers */
52template< class Stepper >
53struct perform_adams_bashforth_test
54{
55 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 = {._M_elems: { 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( x: x1[0]*x1[0] + x1[1]*x1[1] );
67 double phi = std::asin(x: 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( x: A*sin(x: t+dt+phi) - x1[0] ) / std::pow( x: dt , y: 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( x: dt , y: o ) > 1E-16 )
85 {
86 x1 = x0;
87 t = 0.0;
88 stepper.initialize( osc() , x1 , t , dt );
89 A = std::sqrt( x: x1[0]*x1[0] + x1[1]*x1[1] );
90 phi = std::asin(x: 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( x: A*sin(x: 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
101typedef 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
111BOOST_AUTO_TEST_CASE_TEMPLATE( adams_bashforth_test , Stepper, adams_bashforth_steppers )
112{
113 perform_adams_bashforth_test< Stepper > tester;
114 tester();
115}
116
117BOOST_AUTO_TEST_SUITE_END()
118

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