1 | /* Boost numeric test of the runge kutta 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_runge_kutta |
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 &x, const state_type &v, state_type &a, |
42 | const double t ) const |
43 | { |
44 | a[0] = -x[0]; |
45 | } |
46 | }; |
47 | |
48 | BOOST_AUTO_TEST_SUITE( velocity_verlet_test ) |
49 | |
50 | BOOST_AUTO_TEST_CASE( numeric_velocity_verlet_test ) |
51 | { |
52 | |
53 | velocity_verlet<state_type> stepper; |
54 | const int steps = 10; |
55 | // order of the error is order of approximation + 1 |
56 | const int o = stepper.order() + 1; |
57 | |
58 | const state_type x0 = {._M_elems: { 0.0 }}; |
59 | const state_type v0 = {._M_elems: { 1.0 }}; |
60 | state_type x = x0; |
61 | state_type v = v0; |
62 | const double t = 0.0; |
63 | /* do a first step with dt=0.1 to get an estimate on the prefactor of the error dx = f * dt^(order+1) */ |
64 | double dt = 0.5; |
65 | for ( int step = 0; step < steps; ++step ) |
66 | { |
67 | stepper.do_step( |
68 | system: osc(), x: std::make_pair( x: boost::ref( t&: x ), y: boost::ref( t&: v ) ), t, dt ); |
69 | } |
70 | const double f = steps * std::abs( x: sin( x: steps * dt ) - x[0] ) / |
71 | std::pow( x: dt, y: o ); // upper bound |
72 | |
73 | std::cout << o << " , " << f << std::endl; |
74 | |
75 | /* as long as we have errors above machine precision */ |
76 | while( f*std::pow( x: dt , y: o ) > 1E-16 ) |
77 | { |
78 | x = x0; |
79 | v = v0; |
80 | stepper.reset(); |
81 | for ( int step = 0; step < steps; ++step ) |
82 | { |
83 | stepper.do_step( system: osc() , x: std::make_pair(x: boost::ref(t&: x), y: boost::ref(t&: v)) , t , dt ); |
84 | } |
85 | std::cout << "Testing dt=" << dt << std::endl; |
86 | BOOST_CHECK_LT( std::abs( sin( steps * dt ) - x[0] ), |
87 | f * std::pow( dt, o ) ); |
88 | dt *= 0.5; |
89 | } |
90 | }; |
91 | |
92 | |
93 | BOOST_AUTO_TEST_SUITE_END() |
94 | |