| /* |
| [auto_generated] |
| libs/numeric/odeint/test/step_size_limitation.cpp |
| |
| [begin_description] |
| Tests the step size limitation functionality |
| [end_description] |
| |
| Copyright 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) |
| */ |
| |
| #define BOOST_TEST_MODULE odeint_integrate_times |
| |
| #include <boost/test/unit_test.hpp> |
| |
| #include <utility> |
| #include <iostream> |
| #include <vector> |
| |
| #include <boost/numeric/odeint.hpp> |
| |
| using namespace boost::unit_test; |
| using namespace boost::numeric::odeint; |
| |
| typedef double value_type; |
| typedef std::vector< value_type > state_type; |
| |
| /*********************************************** |
| * first part of the tests: explicit methods |
| *********************************************** |
| */ |
| |
| void damped_osc( const state_type &x , state_type &dxdt , const value_type t ) |
| { |
| const value_type gam( 0.1); |
| |
| dxdt[0] = x[1]; |
| dxdt[1] = -x[0] - gam*x[1]; |
| } |
| |
| |
| struct push_back_time |
| { |
| std::vector< double >& m_times; |
| |
| push_back_time( std::vector< double > × ) |
| : m_times( times ) { } |
| |
| template<typename State> |
| void operator()( const State &x , double t ) |
| { |
| m_times.push_back( t ); |
| } |
| }; |
| |
| BOOST_AUTO_TEST_SUITE( step_size_limitation_test ) |
| |
| BOOST_AUTO_TEST_CASE( test_step_adjuster ) |
| { |
| // first use adjuster without step size limitation |
| default_step_adjuster<double, double> step_adjuster; |
| const double dt = 0.1; |
| double dt_new = step_adjuster.decrease_step(dt, 1.5, 2); |
| BOOST_CHECK(dt_new < dt*2.0/3.0); |
| |
| dt_new = step_adjuster.increase_step(dt, 0.8, 1); |
| // for errors > 0.5 no increase is performed |
| BOOST_CHECK(dt_new == dt); |
| |
| dt_new = step_adjuster.increase_step(dt, 0.4, 1); |
| // smaller errors should lead to step size increase |
| std::cout << dt_new << std::endl; |
| BOOST_CHECK(dt_new > dt); |
| |
| |
| // now test with step size limitation max_dt = 0.1 |
| default_step_adjuster<double, double> |
| limited_adjuster(dt); |
| |
| dt_new = limited_adjuster.decrease_step(dt, 1.5, 2); |
| // decreasing works as before |
| BOOST_CHECK(dt_new < dt*2.0/3.0); |
| |
| dt_new = limited_adjuster.decrease_step(2*dt, 1.1, 2); |
| // decreasing a large step size should give max_dt |
| BOOST_CHECK(dt_new == dt); |
| |
| dt_new = limited_adjuster.increase_step(dt, 0.8, 1); |
| // for errors > 0.5 no increase is performed, still valid |
| BOOST_CHECK(dt_new == dt); |
| |
| dt_new = limited_adjuster.increase_step(dt, 0.4, 1); |
| // but even for smaller errors, we should at most get 0.1 |
| BOOST_CHECK_EQUAL(dt_new, dt); |
| |
| dt_new = limited_adjuster.increase_step(0.9*dt, 0.1, 1); |
| std::cout << dt_new << std::endl; |
| // check that we don't increase beyond max_dt |
| BOOST_CHECK(dt_new == dt); |
| } |
| |
| |
| template<class Stepper> |
| void test_explicit_stepper(Stepper stepper, const double max_dt) |
| { |
| state_type x( 2 ); |
| x[0] = x[1] = 10.0; |
| const size_t steps = 100; |
| |
| std::vector<double> times; |
| |
| integrate_adaptive(stepper, damped_osc, x, 0.0, steps*max_dt, max_dt, push_back_time(times)); |
| |
| BOOST_CHECK_EQUAL(times.size(), steps+1); |
| // check that dt remains at exactly max_dt |
| for( size_t i=0 ; i<times.size() ; ++i ) |
| // check if observer was called at times 0,1,2,... |
| BOOST_CHECK_SMALL( times[i] - static_cast<double>(i)*max_dt , 1E-15); |
| times.clear(); |
| |
| // this should also work when we provide some bigger initial dt |
| integrate_adaptive(stepper, damped_osc, x, 0.0, steps*max_dt, 10*max_dt, push_back_time(times)); |
| |
| BOOST_CHECK_EQUAL(times.size(), steps+1); |
| // check that dt remains at exactly max_dt |
| for( size_t i=0 ; i<times.size() ; ++i ) |
| // check if observer was called at times 0,1,2,... |
| BOOST_CHECK_SMALL( times[i] - static_cast<double>(i)*max_dt , 1E-15); |
| times.clear(); |
| } |
| |
| |
| BOOST_AUTO_TEST_CASE( test_controlled ) |
| { |
| const double max_dt = 0.01; |
| |
| test_explicit_stepper(make_controlled(1E-2, 1E-2, max_dt, |
| runge_kutta_dopri5<state_type>()), |
| max_dt); |
| test_explicit_stepper(make_controlled(1E-2, 1E-2, -max_dt, |
| runge_kutta_dopri5<state_type>()), |
| -max_dt); |
| |
| test_explicit_stepper(make_controlled(1E-2, 1E-2, max_dt, |
| runge_kutta_cash_karp54<state_type>()), |
| max_dt); |
| test_explicit_stepper(make_controlled(1E-2, 1E-2, -max_dt, |
| runge_kutta_cash_karp54<state_type>()), |
| -max_dt); |
| |
| test_explicit_stepper(bulirsch_stoer<state_type>(1E-2, 1E-2, 1.0, 1.0, max_dt), |
| max_dt); |
| test_explicit_stepper(bulirsch_stoer<state_type>(1E-2, 1E-2, 1.0, 1.0, -max_dt), |
| -max_dt); |
| } |
| |
| |
| BOOST_AUTO_TEST_CASE( test_dense_out ) |
| { |
| const double max_dt = 0.01; |
| test_explicit_stepper(make_dense_output(1E-2, 1E-2, max_dt, |
| runge_kutta_dopri5<state_type>()), |
| max_dt); |
| test_explicit_stepper(make_dense_output(1E-2, 1E-2, -max_dt, |
| runge_kutta_dopri5<state_type>()), |
| -max_dt); |
| |
| test_explicit_stepper(bulirsch_stoer_dense_out<state_type>(1E-2, 1E-2, 1, 1, max_dt), |
| max_dt); |
| |
| test_explicit_stepper(bulirsch_stoer_dense_out<state_type>(1E-2, 1E-2, 1, 1, -max_dt), |
| -max_dt); |
| } |
| |
| |
| /*********************************************** |
| * second part of the tests: implicit Rosenbrock |
| *********************************************** |
| */ |
| |
| typedef boost::numeric::ublas::vector< value_type > vector_type; |
| typedef boost::numeric::ublas::matrix< value_type > matrix_type; |
| |
| |
| // harmonic oscillator, analytic solution x[0] = sin( t ) |
| struct osc_rhs |
| { |
| void operator()( const vector_type &x , vector_type &dxdt , const value_type &t ) const |
| { |
| dxdt( 0 ) = x( 1 ); |
| dxdt( 1 ) = -x( 0 ); |
| } |
| }; |
| |
| struct osc_jacobi |
| { |
| void operator()( const vector_type &x , matrix_type &jacobi , const value_type &t , vector_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; |
| } |
| }; |
| |
| |
| template<class Stepper> |
| void test_rosenbrock_stepper(Stepper stepper, const double max_dt) |
| { |
| vector_type x( 2 ); |
| x(0) = x(1) = 10.0; |
| const size_t steps = 100; |
| |
| std::vector<double> times; |
| |
| integrate_adaptive(stepper, |
| std::make_pair(osc_rhs(), osc_jacobi()), |
| x, 0.0, steps*max_dt, max_dt, push_back_time(times)); |
| |
| BOOST_CHECK_EQUAL(times.size(), steps+1); |
| // check that dt remains at exactly max_dt |
| for( size_t i=0 ; i<times.size() ; ++i ) |
| // check if observer was called at times 0,1,2,... |
| BOOST_CHECK_SMALL( times[i] - static_cast<double>(i)*max_dt , 1E-15); |
| times.clear(); |
| |
| // this should also work when we provide some bigger initial dt |
| integrate_adaptive(stepper, |
| std::make_pair(osc_rhs(), osc_jacobi()), |
| x, 0.0, steps*max_dt, 10*max_dt, push_back_time(times)); |
| |
| BOOST_CHECK_EQUAL(times.size(), steps+1); |
| // check that dt remains at exactly max_dt |
| for( size_t i=0 ; i<times.size() ; ++i ) |
| // check if observer was called at times 0,1,2,... |
| BOOST_CHECK_SMALL( times[i] - static_cast<double>(i)*max_dt , 1E-15); |
| times.clear(); |
| } |
| |
| |
| BOOST_AUTO_TEST_CASE( test_controlled_rosenbrock ) |
| { |
| const double max_dt = 0.01; |
| |
| test_rosenbrock_stepper(make_controlled(1E-2, 1E-2, max_dt, rosenbrock4<value_type>()), |
| max_dt); |
| test_rosenbrock_stepper(make_controlled(1E-2, 1E-2, -max_dt, rosenbrock4<value_type>()), |
| -max_dt); |
| } |
| |
| |
| BOOST_AUTO_TEST_CASE( test_dense_out_rosenbrock ) |
| { |
| const double max_dt = 0.01; |
| |
| test_rosenbrock_stepper(make_dense_output(1E-2, 1E-2, max_dt, rosenbrock4<value_type>()), |
| max_dt); |
| test_rosenbrock_stepper(make_dense_output(1E-2, 1E-2, -max_dt, rosenbrock4<value_type>()), |
| -max_dt); |
| } |
| |
| BOOST_AUTO_TEST_SUITE_END() |