123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268 |
- /*
- [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()
|