| /* 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/test/unit_test.hpp> |
| |
| #include <boost/mpl/vector.hpp> |
| |
| #include <boost/numeric/odeint.hpp> |
| |
| #include <boost/numeric/ublas/vector.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() |