Squashed 'third_party/boostorg/odeint/' content from commit 6ff2719
Change-Id: If4892e29c1a5e6cf3a7aa51486a2725c251b0c7d
git-subtree-dir: third_party/boostorg/odeint
git-subtree-split: 6ff2719b6907b86596c3d43e88c1bcfdf29df560
diff --git a/test/numeric/Jamfile.v2 b/test/numeric/Jamfile.v2
new file mode 100644
index 0000000..1915a32
--- /dev/null
+++ b/test/numeric/Jamfile.v2
@@ -0,0 +1,35 @@
+# Copyright 2012 Karsten Ahnert
+# Copyright 2012 Mario Mulansky
+# Distributed under the Boost Software License, Version 1.0.
+# (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+# bring in rules for testing
+
+
+import testing ;
+
+use-project boost : $(BOOST_ROOT) ;
+
+project
+ : requirements
+ <library>/boost/test//boost_unit_test_framework
+ <define>BOOST_ALL_NO_LIB=1
+ <include>../../include
+ <link>static
+ <toolset>clang:<cxxflags>-Wno-unused-variable
+
+# <cxxflags>-D_SCL_SECURE_NO_WARNINGS
+ ;
+
+test-suite "odeint"
+ :
+ [ run runge_kutta.cpp ]
+ [ run symplectic.cpp ]
+ [ run rosenbrock.cpp ]
+ [ run adams_bashforth.cpp ]
+ [ run adams_bashforth_moulton.cpp ]
+ [ run abm_time_dependent.cpp ]
+ [ run order_quadrature_formula.cpp ]
+ [ run velocity_verlet.cpp ]
+ : <testing.launcher>valgrind
+ ;
diff --git a/test/numeric/abm_time_dependent.cpp b/test/numeric/abm_time_dependent.cpp
new file mode 100644
index 0000000..019a17d
--- /dev/null
+++ b/test/numeric/abm_time_dependent.cpp
@@ -0,0 +1,85 @@
+/* Boost numeric test of the adams-bashforth-moulton steppers test file
+
+ Copyright 2013 Karsten Ahnert
+ Copyright 2013-2015 Mario Mulansky
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+*/
+
+// disable checked iterator warning for msvc
+#include <boost/config.hpp>
+#ifdef BOOST_MSVC
+ #pragma warning(disable:4996)
+#endif
+
+#define BOOST_TEST_MODULE numeric_adams_bashforth_moulton
+
+#include <iostream>
+#include <cmath>
+
+#include <boost/test/unit_test.hpp>
+
+#include <boost/mpl/vector.hpp>
+
+#include <boost/numeric/odeint.hpp>
+
+using namespace boost::unit_test;
+using namespace boost::numeric::odeint;
+namespace mpl = boost::mpl;
+
+typedef double value_type;
+
+typedef value_type state_type;
+
+
+// simple time-dependent rhs, analytic solution x = 0.5*t^2
+struct simple_rhs
+{
+ void operator()( const state_type& x , state_type &dxdt , const double t ) const
+ {
+ dxdt = t;
+ }
+};
+
+BOOST_AUTO_TEST_SUITE( numeric_abm_time_dependent_test )
+
+
+/* generic test for all adams bashforth moulton steppers */
+template< class Stepper >
+struct perform_abm_time_dependent_test
+{
+ void operator()( void )
+ {
+ Stepper stepper;
+ const int o = stepper.order()+1; //order of the error is order of approximation + 1
+
+ const state_type x0 = 0.0;
+ state_type x1 = x0;
+ double t = 0.0;
+ double dt = 0.1;
+ const int steps = 10;
+
+ integrate_n_steps( boost::ref(stepper) , simple_rhs(), x1 , t , dt , steps );
+ BOOST_CHECK_LT( std::abs( 0.5 - x1 ) , std::pow( dt , o ) );
+ }
+};
+
+typedef mpl::vector<
+ adams_bashforth_moulton< 2 , state_type > ,
+ adams_bashforth_moulton< 3 , state_type > ,
+ adams_bashforth_moulton< 4 , state_type > ,
+ adams_bashforth_moulton< 5 , state_type > ,
+ adams_bashforth_moulton< 6 , state_type > ,
+ adams_bashforth_moulton< 7 , state_type > ,
+ adams_bashforth_moulton< 8 , state_type >
+ > adams_bashforth_moulton_steppers;
+
+BOOST_AUTO_TEST_CASE_TEMPLATE( abm_time_dependent_test , Stepper, adams_bashforth_moulton_steppers )
+{
+ perform_abm_time_dependent_test< Stepper > tester;
+ tester();
+}
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/test/numeric/adams_bashforth.cpp b/test/numeric/adams_bashforth.cpp
new file mode 100644
index 0000000..ca8a6ce
--- /dev/null
+++ b/test/numeric/adams_bashforth.cpp
@@ -0,0 +1,117 @@
+/* Boost numeric test of the adams-bashforth steppers test file
+
+ Copyright 2013 Karsten Ahnert
+ Copyright 2013-2015 Mario Mulansky
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+*/
+
+// disable checked iterator warning for msvc
+#include <boost/config.hpp>
+#ifdef BOOST_MSVC
+ #pragma warning(disable:4996)
+#endif
+
+#define BOOST_TEST_MODULE numeric_adams_bashforth
+
+#include <iostream>
+#include <cmath>
+
+#include <boost/array.hpp>
+
+#include <boost/test/unit_test.hpp>
+
+#include <boost/mpl/vector.hpp>
+
+#include <boost/numeric/odeint.hpp>
+
+using namespace boost::unit_test;
+using namespace boost::numeric::odeint;
+namespace mpl = boost::mpl;
+
+typedef double value_type;
+
+typedef boost::array< double , 2 > state_type;
+
+// harmonic oscillator, analytic solution x[0] = sin( t )
+struct osc
+{
+ void operator()( const state_type &x , state_type &dxdt , const double t ) const
+ {
+ dxdt[0] = x[1];
+ dxdt[1] = -x[0];
+ }
+};
+
+BOOST_AUTO_TEST_SUITE( numeric_adams_bashforth_test )
+
+
+/* generic test for all adams bashforth steppers */
+template< class Stepper >
+struct perform_adams_bashforth_test
+{
+ void operator()( void )
+ {
+ Stepper stepper;
+ const int o = stepper.order()+1; //order of the error is order of approximation + 1
+
+ const state_type x0 = {{ 0.0 , 1.0 }};
+ state_type x1 = x0;
+ double t = 0.0;
+ double dt = 0.2;
+ // initialization, does a number of steps already to fill internal buffer, t is increased
+ stepper.initialize( osc() , x1 , t , dt );
+ double A = std::sqrt( x1[0]*x1[0] + x1[1]*x1[1] );
+ double phi = std::asin(x1[0]/A) - t;
+ // do a number of steps to fill the buffer with results from adams bashforth
+ for( size_t n=0 ; n < stepper.steps ; ++n )
+ {
+ stepper.do_step( osc() , x1 , t , dt );
+ t += dt;
+ }
+ // now we do the actual step
+ stepper.do_step( osc() , x1 , t , dt );
+ // only examine the error of the adams-bashforth step, not the initialization
+ const double f = 2.0 * std::abs( A*sin(t+dt+phi) - x1[0] ) / std::pow( dt , o ); // upper bound
+
+ std::cout << o << " , "
+ << Stepper::initializing_stepper_type::order_value+1 << " , "
+ << f << std::endl;
+
+ /* as long as we have errors above machine precision */
+ while( f*std::pow( dt , o ) > 1E-16 )
+ {
+ x1 = x0;
+ t = 0.0;
+ stepper.initialize( osc() , x1 , t , dt );
+ A = std::sqrt( x1[0]*x1[0] + x1[1]*x1[1] );
+ phi = std::asin(x1[0]/A) - t;
+ // now we do the actual step
+ stepper.do_step( osc() , x1 , t , dt );
+ // only examine the error of the adams-bashforth step, not the initialization
+ std::cout << "Testing dt=" << dt << " , " << std::abs( A*sin(t+dt+phi) - x1[0] ) << std::endl;
+ BOOST_CHECK_LT( std::abs( A*sin(t+dt+phi) - x1[0] ) , f*std::pow( dt , o ) );
+ dt *= 0.5;
+ }
+ }
+};
+
+typedef mpl::vector<
+ adams_bashforth< 2 , state_type > ,
+ adams_bashforth< 3 , state_type > ,
+ adams_bashforth< 4 , state_type > ,
+ adams_bashforth< 5 , state_type > ,
+ adams_bashforth< 6 , state_type > ,
+ adams_bashforth< 7 , state_type > ,
+ adams_bashforth< 8 , state_type >
+ > adams_bashforth_steppers;
+
+BOOST_AUTO_TEST_CASE_TEMPLATE( adams_bashforth_test , Stepper, adams_bashforth_steppers )
+{
+ perform_adams_bashforth_test< Stepper > tester;
+ tester();
+}
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/test/numeric/adams_bashforth_moulton.cpp b/test/numeric/adams_bashforth_moulton.cpp
new file mode 100644
index 0000000..dea5718
--- /dev/null
+++ b/test/numeric/adams_bashforth_moulton.cpp
@@ -0,0 +1,119 @@
+/* Boost numeric test of the adams-bashforth-moulton steppers test file
+
+ Copyright 2013 Karsten Ahnert
+ Copyright 2013 Mario Mulansky
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+*/
+
+// disable checked iterator warning for msvc
+#include <boost/config.hpp>
+#ifdef BOOST_MSVC
+ #pragma warning(disable:4996)
+#endif
+
+#define BOOST_TEST_MODULE numeric_adams_bashforth_moulton
+
+#include <iostream>
+#include <cmath>
+
+#include <boost/array.hpp>
+
+#include <boost/test/unit_test.hpp>
+
+#include <boost/mpl/vector.hpp>
+
+#include <boost/numeric/odeint.hpp>
+
+using namespace boost::unit_test;
+using namespace boost::numeric::odeint;
+namespace mpl = boost::mpl;
+
+typedef double value_type;
+
+typedef boost::array< double , 2 > state_type;
+typedef runge_kutta_fehlberg78<state_type> initializing_stepper;
+
+// harmonic oscillator, analytic solution x[0] = sin( t )
+struct osc
+{
+ void operator()( const state_type &x , state_type &dxdt , const double t ) const
+ {
+ dxdt[0] = x[1];
+ dxdt[1] = -x[0];
+ }
+};
+
+BOOST_AUTO_TEST_SUITE( numeric_adams_bashforth_moulton_test )
+
+
+/* generic test for all adams bashforth moulton steppers */
+template< class Stepper >
+struct perform_adams_bashforth_moulton_test
+{
+ void operator()( void )
+ {
+ Stepper stepper;
+ initializing_stepper init_stepper;
+ const int o = stepper.order()+1; //order of the error is order of approximation + 1
+
+ const state_type x0 = {{ 0.0 , 1.0 }};
+ state_type x1 = x0;
+ double t = 0.0;
+ double dt = 0.25;
+ // initialization, does a number of steps already to fill internal buffer, t is increased
+ // we use the rk78 as initializing stepper
+ stepper.initialize( boost::ref(init_stepper) , osc() , x1 , t , dt );
+ // do a number of steps to fill the buffer with results from adams bashforth
+ for( size_t n=0 ; n < stepper.steps ; ++n )
+ {
+ stepper.do_step( osc() , x1 , t , dt );
+ t += dt;
+ }
+ double A = std::sqrt( x1[0]*x1[0] + x1[1]*x1[1] );
+ double phi = std::asin(x1[0]/A) - t;
+ // now we do the actual step
+ stepper.do_step( osc() , x1 , t , dt );
+ // only examine the error of the adams-bashforth-moulton step, not the initialization
+ const double f = 2.0 * std::abs( A*sin(t+dt+phi) - x1[0] ) / std::pow( dt , o ); // upper bound
+
+ std::cout << o << " , " << f << std::endl;
+
+ /* as long as we have errors above machine precision */
+ while( f*std::pow( dt , o ) > 1E-16 )
+ {
+ x1 = x0;
+ t = 0.0;
+ stepper.initialize( boost::ref(init_stepper) , osc() , x1 , t , dt );
+ A = std::sqrt( x1[0]*x1[0] + x1[1]*x1[1] );
+ phi = std::asin(x1[0]/A) - t;
+ // now we do the actual step
+ stepper.do_step( osc() , x1 , t , dt );
+ // only examine the error of the adams-bashforth-moulton step, not the initialization
+ std::cout << "Testing dt=" << dt << " , " << std::abs( A*sin(t+dt+phi) - x1[0] ) << std::endl;
+ BOOST_CHECK_LT( std::abs( A*sin(t+dt+phi) - x1[0] ) , f*std::pow( dt , o ) );
+ dt *= 0.5;
+ }
+ }
+};
+
+typedef mpl::vector<
+ adams_bashforth_moulton< 1 , state_type > ,
+ adams_bashforth_moulton< 2 , state_type > ,
+ adams_bashforth_moulton< 3 , state_type > ,
+ adams_bashforth_moulton< 4 , state_type > ,
+ adams_bashforth_moulton< 5 , state_type > ,
+ adams_bashforth_moulton< 6 , state_type > ,
+ adams_bashforth_moulton< 7 , state_type > ,
+ adams_bashforth_moulton< 8 , state_type >
+ > adams_bashforth_moulton_steppers;
+
+BOOST_AUTO_TEST_CASE_TEMPLATE( adams_bashforth_moulton_test , Stepper, adams_bashforth_moulton_steppers )
+{
+ perform_adams_bashforth_moulton_test< Stepper > tester;
+ tester();
+}
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/test/numeric/order_quadrature_formula.cpp b/test/numeric/order_quadrature_formula.cpp
new file mode 100644
index 0000000..cefab4d
--- /dev/null
+++ b/test/numeric/order_quadrature_formula.cpp
@@ -0,0 +1,195 @@
+/* Boost numeric test for orders of quadrature formulas test file
+
+ Copyright 2015 Gregor de Cillia
+ Copyright 2015 Mario Mulansky <mario.mulansky@gmx.net>
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+*/
+
+// disable checked iterator warning for msvc
+#include <boost/config.hpp>
+
+#ifdef BOOST_MSVC
+ #pragma warning(disable:4996)
+#endif
+
+#define BOOST_TEST_MODULE order_quadrature_formula
+
+#include <iostream>
+#include <cmath>
+
+#include "boost/format.hpp"
+
+#include <boost/test/unit_test.hpp>
+
+#include <boost/mpl/vector.hpp>
+
+#include <boost/numeric/odeint.hpp>
+
+#include <boost/numeric/ublas/vector.hpp>
+
+#include <boost/format.hpp>
+
+using namespace boost::unit_test;
+using namespace boost::numeric::odeint;
+namespace mpl = boost::mpl;
+
+typedef double value_type;
+typedef value_type time_type;
+typedef value_type state_type;
+
+BOOST_AUTO_TEST_SUITE( order_of_convergence_test )
+
+/* defines the simple monomial f(t) = (p+1) * t^p.*/
+struct monomial
+{
+ int power;
+
+ monomial(int p = 0) : power( p ){};
+
+ void operator()( const state_type &x , state_type &dxdt , const time_type t )
+ {
+ dxdt = ( 1.0 + power ) * pow( t, power );
+ }
+};
+
+
+/* generic test for all steppers that support integrate_const */
+template< class Stepper >
+struct stepper_order_test
+{
+ void operator()( int steps = 1 )
+ {
+ const int estimated_order = estimate_order( steps );
+ const int defined_order = Stepper::order_value;
+
+ std::cout << boost::format( "%-20i%-20i\n" )
+ % estimated_order % defined_order;
+
+ BOOST_REQUIRE_EQUAL( estimated_order, defined_order );
+ }
+
+ /*
+ the order of the stepper is estimated by trying to solve the ODE
+ x'(t) = (p+1) * t^p
+ until the errors are too big to be justified by finite precision.
+ the first value p for which the problem is *not* solved within the
+ finite precision tolerance is the estimate for the order of the scheme.
+ */
+ int estimate_order( int steps )
+ {
+ const double dt = 1.0/steps;
+ const double tolerance = steps*1E-15;
+ int p;
+ for( p = 0; true; p++ )
+ {
+ // begin with x'(t) = t^0 = 1
+ // => x (t) = t
+ // then use x'(t) = 2*t^1
+ // => x (t) = t^2
+ // ...
+ state_type x = 0.0;
+
+ double t = integrate_n_steps( Stepper(), monomial( p ), x, 0.0, dt,
+ steps );
+ if( fabs( x - pow( t, ( 1.0 + p ) ) ) > tolerance )
+ break;
+ }
+ // the smallest power p for which the test failed is the estimated order,
+ // as the solution for this power is x(t) = t^{p+1}
+ return p;
+ }
+};
+
+
+typedef mpl::vector<
+ euler< state_type > ,
+ modified_midpoint< state_type > ,
+ runge_kutta4< state_type > ,
+ runge_kutta4_classic< state_type > ,
+ runge_kutta_cash_karp54_classic< state_type > ,
+ runge_kutta_cash_karp54< state_type > ,
+ runge_kutta_dopri5< state_type > ,
+ runge_kutta_fehlberg78< state_type >
+ > runge_kutta_steppers;
+
+typedef mpl::vector<
+ adams_bashforth< 2, state_type, double, state_type, double,
+ vector_space_algebra, default_operations,
+ initially_resizer, runge_kutta_fehlberg78< state_type > >,
+ adams_bashforth< 3, state_type, double, state_type, double,
+ vector_space_algebra, default_operations,
+ initially_resizer, runge_kutta_fehlberg78< state_type > >,
+ adams_bashforth< 4, state_type, double, state_type, double,
+ vector_space_algebra, default_operations,
+ initially_resizer, runge_kutta_fehlberg78< state_type > >,
+ adams_bashforth< 5, state_type, double, state_type, double,
+ vector_space_algebra, default_operations,
+ initially_resizer, runge_kutta_fehlberg78< state_type > >,
+ adams_bashforth< 6, state_type, double, state_type, double,
+ vector_space_algebra, default_operations,
+ initially_resizer, runge_kutta_fehlberg78< state_type > >,
+ adams_bashforth< 7, state_type, double, state_type, double,
+ vector_space_algebra, default_operations,
+ initially_resizer, runge_kutta_fehlberg78< state_type > >,
+ adams_bashforth< 8, state_type, double, state_type, double,
+ vector_space_algebra, default_operations,
+ initially_resizer, runge_kutta_fehlberg78< state_type > >
+ > ab_steppers;
+
+
+typedef mpl::vector<
+ adams_bashforth_moulton< 2, state_type, double, state_type, double,
+ vector_space_algebra, default_operations,
+ initially_resizer,
+ runge_kutta_fehlberg78< state_type > >,
+ adams_bashforth_moulton< 3, state_type, double, state_type, double,
+ vector_space_algebra, default_operations,
+ initially_resizer,
+ runge_kutta_fehlberg78< state_type > >,
+ adams_bashforth_moulton< 4, state_type, double, state_type, double,
+ vector_space_algebra, default_operations,
+ initially_resizer,
+ runge_kutta_fehlberg78< state_type > >,
+ adams_bashforth_moulton< 5, state_type, double, state_type, double,
+ vector_space_algebra, default_operations,
+ initially_resizer,
+ runge_kutta_fehlberg78< state_type > >,
+ adams_bashforth_moulton< 6, state_type, double, state_type, double,
+ vector_space_algebra, default_operations,
+ initially_resizer,
+ runge_kutta_fehlberg78< state_type > >,
+ adams_bashforth_moulton< 7, state_type, double, state_type, double,
+ vector_space_algebra, default_operations,
+ initially_resizer,
+ runge_kutta_fehlberg78< state_type > >,
+ adams_bashforth_moulton< 8, state_type, double, state_type, double,
+ vector_space_algebra, default_operations,
+ initially_resizer,
+ runge_kutta_fehlberg78< state_type > >
+ > abm_steppers;
+
+
+BOOST_AUTO_TEST_CASE_TEMPLATE( runge_kutta_test , Stepper, runge_kutta_steppers )
+{
+ stepper_order_test< Stepper > tester;
+ tester(10);
+}
+
+
+BOOST_AUTO_TEST_CASE_TEMPLATE( adams_bashforth_test , Stepper, ab_steppers )
+{
+ stepper_order_test< Stepper > tester;
+ tester(16);
+}
+
+
+BOOST_AUTO_TEST_CASE_TEMPLATE( adams_bashforth_moultion_test , Stepper, abm_steppers )
+{
+ stepper_order_test< Stepper > tester;
+ tester(16);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/test/numeric/rosenbrock.cpp b/test/numeric/rosenbrock.cpp
new file mode 100644
index 0000000..2b6a267
--- /dev/null
+++ b/test/numeric/rosenbrock.cpp
@@ -0,0 +1,89 @@
+/* Boost numeric test of the rosenbrock4 stepper test file
+
+ Copyright 2012 Karsten Ahnert
+ Copyright 2012 Mario Mulansky
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+*/
+
+// disable checked iterator warning for msvc
+#include <boost/config.hpp>
+#ifdef BOOST_MSVC
+ #pragma warning(disable:4996)
+#endif
+
+#define BOOST_TEST_MODULE numeric_rosenbrock
+
+#include <iostream>
+#include <cmath>
+
+#include <boost/array.hpp>
+
+#include <boost/test/unit_test.hpp>
+
+#include <boost/numeric/odeint/stepper/rosenbrock4.hpp>
+#include <boost/numeric/ublas/vector.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+
+using namespace boost::unit_test;
+using namespace boost::numeric::odeint;
+
+typedef double value_type;
+typedef boost::numeric::ublas::vector< value_type > state_type;
+typedef boost::numeric::ublas::matrix< value_type > matrix_type;
+
+// harmonic oscillator, analytic solution x[0] = sin( t )
+struct sys
+{
+ void operator()( const state_type &x , state_type &dxdt , const value_type &t ) const
+ {
+ dxdt( 0 ) = x( 1 );
+ dxdt( 1 ) = -x( 0 );
+ }
+};
+
+struct jacobi
+{
+ void operator()( const state_type &x , matrix_type &jacobi , const value_type &t , state_type &dfdt ) const
+ {
+ jacobi( 0 , 0 ) = 0;
+ jacobi( 0 , 1 ) = 1;
+ jacobi( 1 , 0 ) = -1;
+ jacobi( 1 , 1 ) = 0;
+ dfdt( 0 ) = 0.0;
+ dfdt( 1 ) = 0.0;
+ }
+};
+
+
+BOOST_AUTO_TEST_SUITE( numeric_rosenbrock4 )
+
+BOOST_AUTO_TEST_CASE( rosenbrock4_numeric_test )
+{
+ typedef rosenbrock4< value_type > stepper_type;
+ stepper_type stepper;
+
+ const int o = stepper.order()+1;
+
+ state_type x0( 2 ) , x1( 2 );
+ x0(0) = 0.0; x0(1) = 1.0;
+
+ double dt = 0.5;
+
+ stepper.do_step( std::make_pair( sys() , jacobi() ) , x0 , 0.0 , x1 , dt );
+ const double f = 2.0 * std::abs( std::sin(dt) - x1(0) ) / std::pow( dt , o );
+
+ std::cout << o << " , " << f << std::endl;
+
+ while( f*std::pow( dt , o ) > 1E-16 )
+ {
+ stepper.do_step( std::make_pair( sys() , jacobi() ) , x0 , 0.0 , x1 , dt );
+ std::cout << "Testing dt=" << dt << std::endl;
+ BOOST_CHECK_SMALL( std::abs( std::sin(dt) - x1(0) ) , f*std::pow( dt , o ) );
+ dt *= 0.5;
+ }
+}
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/test/numeric/runge_kutta.cpp b/test/numeric/runge_kutta.cpp
new file mode 100644
index 0000000..a71701b
--- /dev/null
+++ b/test/numeric/runge_kutta.cpp
@@ -0,0 +1,177 @@
+/* Boost numeric test of the runge kutta steppers test file
+
+ Copyright 2012 Mario Mulansky
+ Copyright 2012 Karsten Ahnert
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+*/
+
+// disable checked iterator warning for msvc
+#include <boost/config.hpp>
+#ifdef BOOST_MSVC
+ #pragma warning(disable:4996)
+#endif
+
+#define BOOST_TEST_MODULE numeric_runge_kutta
+
+#include <iostream>
+#include <cmath>
+
+#include <boost/array.hpp>
+
+#include <boost/test/unit_test.hpp>
+
+#include <boost/mpl/vector.hpp>
+
+#include <boost/numeric/odeint.hpp>
+#include <boost/numeric/odeint/stepper/extrapolation_stepper.hpp>
+
+using namespace boost::unit_test;
+using namespace boost::numeric::odeint;
+namespace mpl = boost::mpl;
+
+typedef double value_type;
+
+typedef boost::array< double , 2 > state_type;
+
+// harmonic oscillator, analytic solution x[0] = sin( t )
+struct osc
+{
+ void operator()( const state_type &x , state_type &dxdt , const double t ) const
+ {
+ dxdt[0] = x[1];
+ dxdt[1] = -x[0];
+ }
+};
+
+/* reset dispatcher */
+template< class StepperCategory >
+struct resetter
+{
+ template< class Stepper >
+ static void reset( Stepper &stepper ) { }
+};
+
+template< >
+struct resetter< explicit_error_stepper_fsal_tag >
+{
+ template< class Stepper >
+ static void reset( Stepper &stepper )
+ { stepper.reset(); }
+};
+
+
+BOOST_AUTO_TEST_SUITE( numeric_runge_kutta_test )
+
+
+/* generic test for all runge kutta steppers */
+template< class Stepper >
+struct perform_runge_kutta_test
+{
+ void operator()( void )
+ {
+
+ Stepper stepper;
+ const int o = stepper.order()+1; //order of the error is order of approximation + 1
+
+ const state_type x0 = {{ 0.0 , 1.0 }};
+ state_type x1;
+ const double t = 0.0;
+ /* do a first step with dt=0.1 to get an estimate on the prefactor of the error dx = f * dt^(order+1) */
+ double dt = 0.5;
+ stepper.do_step( osc() , x0 , t , x1 , dt );
+ const double f = 2.0 * std::abs( sin(dt) - x1[0] ) / std::pow( dt , o ); // upper bound
+
+ std::cout << o << " , " << f << std::endl;
+
+ /* as long as we have errors above machine precision */
+ while( f*std::pow( dt , o ) > 1E-16 )
+ {
+ // reset stepper which require resetting (fsal steppers)
+ resetter< typename Stepper::stepper_category >::reset( stepper );
+
+ stepper.do_step( osc() , x0 , t , x1 , dt );
+ std::cout << "Testing dt=" << dt << std::endl;
+ BOOST_CHECK_LT( std::abs( sin(dt) - x1[0] ) , f*std::pow( dt , o ) );
+ dt *= 0.5;
+ }
+ }
+};
+
+
+/* generic error test for all runge kutta steppers */
+template< class Stepper >
+struct perform_runge_kutta_error_test
+{
+ void operator()( void )
+ {
+ Stepper stepper;
+ const int o = stepper.error_order()+1; //order of the error is order of approximation + 1
+
+ const state_type x0 = {{ 0.0 , 1.0 }};
+ state_type x1 , x_err;
+ const double t = 0.0;
+ /* do a first step with dt=0.1 to get an estimate on the prefactor of the error dx = f * dt^(order+1) */
+ double dt = 0.5;
+ stepper.do_step( osc() , x0 , t , x1 , dt , x_err );
+ const double f = 2.0 * std::abs( x_err[0] ) / std::pow( dt , o );
+
+ std::cout << o << " , " << f << " , " << x0[0] << std::endl;
+
+ /* as long as we have errors above machine precision */
+ while( f*std::pow( dt , o ) > 1E-16 )
+ {
+ // reset stepper which require resetting (fsal steppers)
+ resetter< typename Stepper::stepper_category >::reset( stepper );
+
+ stepper.do_step( osc() , x0 , t , x1 , dt , x_err );
+ std::cout << "Testing dt=" << dt << ": " << x_err[1] << std::endl;
+ BOOST_CHECK_SMALL( std::abs( x_err[0] ) , f*std::pow( dt , o ) );
+ dt *= 0.5;
+ }
+ }
+};
+
+
+typedef mpl::vector<
+ euler< state_type > ,
+ modified_midpoint< state_type > ,
+ runge_kutta4< state_type > ,
+ runge_kutta4_classic< state_type > ,
+ runge_kutta_cash_karp54_classic< state_type > ,
+ runge_kutta_cash_karp54< state_type > ,
+ runge_kutta_dopri5< state_type > ,
+ runge_kutta_fehlberg78< state_type > ,
+ extrapolation_stepper< 4, state_type > ,
+ extrapolation_stepper< 6, state_type > ,
+ extrapolation_stepper< 8, state_type > ,
+ extrapolation_stepper< 10, state_type >
+ > runge_kutta_steppers;
+
+BOOST_AUTO_TEST_CASE_TEMPLATE( runge_kutta_test , Stepper, runge_kutta_steppers )
+{
+ perform_runge_kutta_test< Stepper > tester;
+ tester();
+}
+
+
+typedef mpl::vector<
+ runge_kutta_cash_karp54_classic< state_type > ,
+ runge_kutta_cash_karp54< state_type > ,
+ runge_kutta_dopri5< state_type > ,
+ runge_kutta_fehlberg78< state_type > ,
+ extrapolation_stepper< 4, state_type > ,
+ extrapolation_stepper< 6, state_type > ,
+ extrapolation_stepper< 8, state_type > ,
+ extrapolation_stepper< 10, state_type >
+ > runge_kutta_error_steppers;
+
+BOOST_AUTO_TEST_CASE_TEMPLATE( runge_kutta_error_test , Stepper, runge_kutta_error_steppers )
+{
+ perform_runge_kutta_error_test< Stepper > tester;
+ tester();
+}
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/test/numeric/symplectic.cpp b/test/numeric/symplectic.cpp
new file mode 100644
index 0000000..6d3f75b
--- /dev/null
+++ b/test/numeric/symplectic.cpp
@@ -0,0 +1,95 @@
+/* Boost numeric test of the symplectic steppers test file
+
+ Copyright 2012 Mario Mulansky
+ Copyright 2012 Karsten Ahnert
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+*/
+
+// disable checked iterator warning for msvc
+#include <boost/config.hpp>
+#ifdef BOOST_MSVC
+ #pragma warning(disable:4996)
+#endif
+
+#define BOOST_TEST_MODULE numeric_symplectic
+
+#include <iostream>
+#include <cmath>
+
+#include <boost/array.hpp>
+
+#include <boost/test/unit_test.hpp>
+
+#include <boost/mpl/vector.hpp>
+
+#include <boost/numeric/odeint.hpp>
+
+using namespace boost::unit_test;
+using namespace boost::numeric::odeint;
+namespace mpl = boost::mpl;
+
+typedef double value_type;
+
+typedef boost::array< double ,1 > state_type;
+
+// harmonic oscillator, analytic solution x[0] = sin( t )
+struct osc
+{
+ void operator()( const state_type &q , state_type &dpdt ) const
+ {
+ dpdt[0] = -q[0];
+ }
+};
+
+BOOST_AUTO_TEST_SUITE( numeric_symplectic_test )
+
+
+/* generic test for all symplectic steppers */
+template< class Stepper >
+struct perform_symplectic_test
+{
+ void operator()( void )
+ {
+
+ Stepper stepper;
+ const int o = stepper.order()+1; //order of the error is order of approximation + 1
+
+ const state_type q0 = {{ 0.0 }};
+ const state_type p0 = {{ 1.0 }};
+ state_type q1,p1;
+ std::pair< state_type , state_type >x1( q1 , p1 );
+ const double t = 0.0;
+ /* do a first step with dt=0.1 to get an estimate on the prefactor of the error dx = f * dt^(order+1) */
+ double dt = 0.5;
+ stepper.do_step( osc() , std::make_pair( q0 , p0 ) , t , x1 , dt );
+ const double f = 2.0 * std::abs( sin(dt) - x1.first[0] ) / std::pow( dt , o );
+
+ std::cout << o << " , " << f << std::endl;
+
+ /* as long as we have errors above machine precision */
+ while( f*std::pow( dt , o ) > 1E-16 )
+ {
+ stepper.do_step( osc() , std::make_pair( q0 , p0 ) , t , x1 , dt );
+ std::cout << "Testing dt=" << dt << std::endl;
+ BOOST_CHECK_SMALL( std::abs( sin(dt) - x1.first[0] ) , f*std::pow( dt , o ) );
+ dt *= 0.5;
+ }
+ }
+};
+
+
+typedef mpl::vector<
+ symplectic_euler< state_type > ,
+ symplectic_rkn_sb3a_mclachlan< state_type >
+ > symplectic_steppers;
+
+BOOST_AUTO_TEST_CASE_TEMPLATE( symplectic_test , Stepper, symplectic_steppers )
+{
+ perform_symplectic_test< Stepper > tester;
+ tester();
+}
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/test/numeric/velocity_verlet.cpp b/test/numeric/velocity_verlet.cpp
new file mode 100644
index 0000000..5ac7780
--- /dev/null
+++ b/test/numeric/velocity_verlet.cpp
@@ -0,0 +1,93 @@
+/* Boost numeric test of the runge kutta steppers test file
+
+ Copyright 2012 Mario Mulansky
+ Copyright 2012 Karsten Ahnert
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+*/
+
+// disable checked iterator warning for msvc
+#include <boost/config.hpp>
+#ifdef BOOST_MSVC
+ #pragma warning(disable:4996)
+#endif
+
+#define BOOST_TEST_MODULE numeric_runge_kutta
+
+#include <iostream>
+#include <cmath>
+
+#include <boost/array.hpp>
+
+#include <boost/test/unit_test.hpp>
+
+#include <boost/mpl/vector.hpp>
+
+#include <boost/numeric/odeint.hpp>
+
+using namespace boost::unit_test;
+using namespace boost::numeric::odeint;
+namespace mpl = boost::mpl;
+
+typedef double value_type;
+
+typedef boost::array< double , 1 > state_type;
+
+// harmonic oscillator, analytic solution x[0] = sin( t )
+struct osc
+{
+ void operator()( const state_type &x, const state_type &v, state_type &a,
+ const double t ) const
+ {
+ a[0] = -x[0];
+ }
+};
+
+BOOST_AUTO_TEST_SUITE( velocity_verlet_test )
+
+BOOST_AUTO_TEST_CASE( numeric_velocity_verlet_test )
+{
+
+ velocity_verlet<state_type> stepper;
+ const int steps = 10;
+ // order of the error is order of approximation + 1
+ const int o = stepper.order() + 1;
+
+ const state_type x0 = {{ 0.0 }};
+ const state_type v0 = {{ 1.0 }};
+ state_type x = x0;
+ state_type v = v0;
+ const double t = 0.0;
+ /* do a first step with dt=0.1 to get an estimate on the prefactor of the error dx = f * dt^(order+1) */
+ double dt = 0.5;
+ for ( int step = 0; step < steps; ++step )
+ {
+ stepper.do_step(
+ osc(), std::make_pair( boost::ref( x ), boost::ref( v ) ), t, dt );
+ }
+ const double f = steps * std::abs( sin( steps * dt ) - x[0] ) /
+ std::pow( dt, o ); // upper bound
+
+ std::cout << o << " , " << f << std::endl;
+
+ /* as long as we have errors above machine precision */
+ while( f*std::pow( dt , o ) > 1E-16 )
+ {
+ x = x0;
+ v = v0;
+ stepper.reset();
+ for ( int step = 0; step < steps; ++step )
+ {
+ stepper.do_step( osc() , std::make_pair(boost::ref(x), boost::ref(v)) , t , dt );
+ }
+ std::cout << "Testing dt=" << dt << std::endl;
+ BOOST_CHECK_LT( std::abs( sin( steps * dt ) - x[0] ),
+ f * std::pow( dt, o ) );
+ dt *= 0.5;
+ }
+};
+
+
+BOOST_AUTO_TEST_SUITE_END()