1 | /* Boost numeric test of the symplectic 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_symplectic |
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 ,1 > state_type; |
37 | |
38 | // harmonic oscillator, analytic solution x[0] = sin( t ) |
39 | struct osc |
40 | { |
41 | void operator()( const state_type &q , state_type &dpdt ) const |
42 | { |
43 | dpdt[0] = -q[0]; |
44 | } |
45 | }; |
46 | |
47 | BOOST_AUTO_TEST_SUITE( numeric_symplectic_test ) |
48 | |
49 | |
50 | /* generic test for all symplectic steppers */ |
51 | template< class Stepper > |
52 | struct perform_symplectic_test |
53 | { |
54 | void operator()( void ) |
55 | { |
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 q0 = {._M_elems: { 0.0 }}; |
61 | const state_type p0 = {._M_elems: { 1.0 }}; |
62 | state_type q1,p1; |
63 | std::pair< state_type , state_type >x1( q1 , p1 ); |
64 | const double t = 0.0; |
65 | /* do a first step with dt=0.1 to get an estimate on the prefactor of the error dx = f * dt^(order+1) */ |
66 | double dt = 0.5; |
67 | stepper.do_step( osc() , std::make_pair( x: q0 , y: p0 ) , t , x1 , dt ); |
68 | const double f = 2.0 * std::abs( x: sin(x: dt) - x1.first[0] ) / std::pow( x: dt , y: o ); |
69 | |
70 | std::cout << o << " , " << f << std::endl; |
71 | |
72 | /* as long as we have errors above machine precision */ |
73 | while( f*std::pow( x: dt , y: o ) > 1E-16 ) |
74 | { |
75 | stepper.do_step( osc() , std::make_pair( x: q0 , y: p0 ) , t , x1 , dt ); |
76 | std::cout << "Testing dt=" << dt << std::endl; |
77 | BOOST_CHECK_SMALL( std::abs( sin(dt) - x1.first[0] ) , f*std::pow( dt , o ) ); |
78 | dt *= 0.5; |
79 | } |
80 | } |
81 | }; |
82 | |
83 | |
84 | typedef mpl::vector< |
85 | symplectic_euler< state_type > , |
86 | symplectic_rkn_sb3a_mclachlan< state_type > |
87 | > symplectic_steppers; |
88 | |
89 | BOOST_AUTO_TEST_CASE_TEMPLATE( symplectic_test , Stepper, symplectic_steppers ) |
90 | { |
91 | perform_symplectic_test< Stepper > tester; |
92 | tester(); |
93 | } |
94 | |
95 | BOOST_AUTO_TEST_SUITE_END() |
96 | |