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
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;
37typedef runge_kutta_fehlberg78<state_type> initializing_stepper;
38
39// harmonic oscillator, analytic solution x[0] = sin( t )
40struct 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
49BOOST_AUTO_TEST_SUITE( numeric_adams_bashforth_moulton_test )
50
51
52/* generic test for all adams bashforth moulton steppers */
53template< class Stepper >
54struct 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
102typedef 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
113BOOST_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
119BOOST_AUTO_TEST_SUITE_END()
120

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