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
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 &x, const state_type &v, state_type &a,
42 const double t ) const
43 {
44 a[0] = -x[0];
45 }
46};
47
48BOOST_AUTO_TEST_SUITE( velocity_verlet_test )
49
50BOOST_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
93BOOST_AUTO_TEST_SUITE_END()
94

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