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
30using namespace boost::unit_test;
31using namespace boost::numeric::odeint;
32namespace mpl = boost::mpl;
33
34typedef double value_type;
35
36typedef std::array< double ,1 > state_type;
37
38// harmonic oscillator, analytic solution x[0] = sin( t )
39struct osc
40{
41 void operator()( const state_type &q , state_type &dpdt ) const
42 {
43 dpdt[0] = -q[0];
44 }
45};
46
47BOOST_AUTO_TEST_SUITE( numeric_symplectic_test )
48
49
50/* generic test for all symplectic steppers */
51template< class Stepper >
52struct 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
84typedef mpl::vector<
85 symplectic_euler< state_type > ,
86 symplectic_rkn_sb3a_mclachlan< state_type >
87 > symplectic_steppers;
88
89BOOST_AUTO_TEST_CASE_TEMPLATE( symplectic_test , Stepper, symplectic_steppers )
90{
91 perform_symplectic_test< Stepper > tester;
92 tester();
93}
94
95BOOST_AUTO_TEST_SUITE_END()
96

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