123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292 |
- /*
- [auto_generated]
- libs/numeric/odeint/test/stepper_with_ranges.cpp
- [begin_description]
- This file tests if the steppers play well with Boost.Range.
- [end_description]
- Copyright 2011-2012 Karsten Ahnert
- Copyright 2011-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 odeint_stepper_with_ranges
- #include <boost/test/unit_test.hpp>
- #include <vector>
- #include <utility>
- #include <iostream>
- #include <boost/array.hpp>
- #include <boost/range.hpp>
- #include <boost/ref.hpp>
- #include <boost/numeric/odeint/stepper/euler.hpp>
- #include <boost/numeric/odeint/stepper/runge_kutta_cash_karp54_classic.hpp>
- #include <boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp>
- #include <boost/numeric/odeint/stepper/controlled_runge_kutta.hpp>
- #include <boost/numeric/odeint/stepper/symplectic_euler.hpp>
- #include <boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp>
- typedef std::vector< double > state_type;
- typedef boost::array< double , 3 > state_type2;
- /* explicitly force range algebra for this array! */
- namespace boost { namespace numeric { namespace odeint {
- template<>
- struct algebra_dispatcher< state_type2 >
- { typedef range_algebra algebra_type; };
- } } }
- /*
- * The two systems are needed, since for steppers with more than
- * one internal step it is difficult to calculate the exact result
- *
- * system1 is suited for euler
- */
- struct system1
- {
- template< class State , class Deriv >
- void operator()( const State &x_ , Deriv &dxdt_ , double t )
- {
- typename boost::range_iterator< const State >::type x = boost::begin( x_ );
- typename boost::range_iterator< Deriv >::type dxdt = boost::begin( dxdt_ );
- dxdt[0] = x[0];
- dxdt[1] = 2.0;
- dxdt[2] = 3.0;
- }
- template< class State , class Deriv >
- void operator()( const State &x_ , const Deriv &dxdt_ , double t )
- {
- typename boost::range_iterator< const State >::type x = boost::begin( x_ );
- typename boost::range_iterator< Deriv >::type dxdt = boost::begin( dxdt_ );
- dxdt[0] = x[0];
- dxdt[1] = 2.0;
- dxdt[2] = 3.0;
- }
- };
- /*
- * system2 is suited for all steppers, it allows you to calculate the result analytically.
- */
- struct system2
- {
- template< class State , class Deriv >
- void operator()( const State &x_ , Deriv &dxdt_ , double t )
- {
- typename boost::range_iterator< Deriv >::type dxdt = boost::begin( dxdt_ );
- dxdt[0] = 1.0;
- dxdt[1] = 2.0;
- dxdt[2] = 3.0;
- }
- template< class State , class Deriv >
- void operator()( const State &x_ , const Deriv &dxdt_ , double t )
- {
- typename boost::range_iterator< Deriv >::type dxdt = boost::begin( dxdt_ );
- dxdt[0] = 1.0;
- dxdt[1] = 2.0;
- dxdt[2] = 3.0;
- }
- };
- /*
- * Useful for Hamiltonian systems
- */
- struct ham_sys
- {
- template< class State , class Deriv >
- void operator()( const State &x_ , Deriv &dxdt_ )
- {
- typename boost::range_iterator< Deriv >::type dxdt = boost::begin( dxdt_ );
- dxdt[0] = 1.0;
- dxdt[1] = 2.0;
- dxdt[2] = 3.0;
- }
- template< class State , class Deriv >
- void operator()( const State &x_ , const Deriv &dxdt_ )
- {
- typename boost::range_iterator< Deriv >::type dxdt = boost::begin( dxdt_ );
- dxdt[0] = 1.0;
- dxdt[1] = 2.0;
- dxdt[2] = 3.0;
- }
- };
- struct vector_fixture
- {
- const static size_t dim = 6;
- boost::array< double , dim > in;
- boost::array< double , dim > q;
- boost::array< double , dim > p;
- state_type err;
- vector_fixture( void )
- : in() , err( 3 )
- {
- for( size_t i=0 ; i<dim ; ++i )
- {
- in[ i ] = q[i] = p[i] = double( i );
- }
- for( size_t i=0 ; i<3 ; ++i )
- {
- err[i] = double( i ) * 10.0;
- }
- }
- ~vector_fixture( void )
- {
- }
- };
- #define CHECK_VALUES( x , x0 , x1 , x2 , x3 , x4 , x5 ) \
- BOOST_CHECK_CLOSE( x[0] , x0 , 1.0e-8 ); \
- BOOST_CHECK_CLOSE( x[1] , x1 , 1.0e-8 ); \
- BOOST_CHECK_CLOSE( x[2] , x2 , 1.0e-8 ); \
- BOOST_CHECK_CLOSE( x[3] , x3 , 1.0e-8 ); \
- BOOST_CHECK_CLOSE( x[4] , x4 , 1.0e-8 ); \
- BOOST_CHECK_CLOSE( x[5] , x5 , 1.0e-8 )
- BOOST_AUTO_TEST_SUITE( stepper_with_ranges )
- BOOST_AUTO_TEST_CASE( explicit_euler_with_range_v1 )
- {
- vector_fixture f;
- boost::numeric::odeint::euler< state_type > euler;
- euler.do_step( system1() , std::make_pair( f.in.begin() + 1 , f.in.begin() + 4 ) , 0.1 , 0.1 );
- CHECK_VALUES( f.in , 0.0 , 1.1 , 2.2 , 3.3 , 4.0 , 5.0 );
- }
- BOOST_AUTO_TEST_CASE( explicit_error_k54_with_range_v1 )
- {
- vector_fixture f;
- boost::numeric::odeint::runge_kutta_cash_karp54_classic< state_type > rk54;
- rk54.do_step( system2() , std::make_pair( f.in.begin() + 1 , f.in.begin() + 4 ) , 0.1 , 0.1 );
- CHECK_VALUES( f.in , 0.0 , 1.1 , 2.2 , 3.3 , 4.0 , 5.0 );
- }
- BOOST_AUTO_TEST_CASE( explicit_error_k54_with_range_v5 )
- {
- vector_fixture f;
- boost::numeric::odeint::runge_kutta_cash_karp54_classic< state_type > rk54;
- rk54.do_step( system2() , std::make_pair( f.in.begin() + 1 , f.in.begin() + 4 ) , 0.1 , 0.1 , f.err );
- CHECK_VALUES( f.in , 0.0 , 1.1 , 2.2 , 3.3 , 4.0 , 5.0 );
- }
- BOOST_AUTO_TEST_CASE( runge_kutta_dopri5_with_range_v1 )
- {
- vector_fixture f;
- boost::numeric::odeint::runge_kutta_dopri5< state_type > dopri5;
- dopri5.do_step( system2() , std::make_pair( f.in.begin() + 1 , f.in.begin() + 4 ) , 0.1 , 0.1 );
- CHECK_VALUES( f.in , 0.0 , 1.1 , 2.2 , 3.3 , 4.0 , 5.0 );
- }
- BOOST_AUTO_TEST_CASE( runge_kutta_dopri5_with_range_v5 )
- {
- vector_fixture f;
- boost::numeric::odeint::runge_kutta_dopri5< state_type > dopri5;
- dopri5.do_step( system2() , std::make_pair( f.in.begin() + 1 , f.in.begin() + 4 ) , 0.1 , 0.1 , f.err );
- CHECK_VALUES( f.in , 0.0 , 1.1 , 2.2 , 3.3 , 4.0 , 5.0 );
- }
- BOOST_AUTO_TEST_CASE( controlled_error_stepper_rk54 )
- {
- double t = 0.0 , dt = 0.1;
- vector_fixture f;
- boost::numeric::odeint::controlled_runge_kutta< boost::numeric::odeint::runge_kutta_cash_karp54_classic< state_type > > stepper;
- stepper.try_step( system2() , std::make_pair( f.in.begin() + 1 , f.in.begin() + 4 ) , t , dt );
- CHECK_VALUES( f.in , 0.0 , 1.1 , 2.2 , 3.3 , 4.0 , 5.0 );
- }
- BOOST_AUTO_TEST_CASE( controlled_error_stepper_dopri5 )
- {
- double t = 0.0 , dt = 0.1;
- vector_fixture f;
- boost::numeric::odeint::controlled_runge_kutta< boost::numeric::odeint::runge_kutta_dopri5< state_type > > stepper;
- stepper.try_step( system2() , std::make_pair( f.in.begin() + 1 , f.in.begin() + 4 ) , t , dt );
- CHECK_VALUES( f.in , 0.0 , 1.1 , 2.2 , 3.3 , 4.0 , 5.0 );
- }
- BOOST_AUTO_TEST_CASE( symplectic_euler_coor_func )
- {
- vector_fixture f;
- boost::numeric::odeint::symplectic_euler< state_type > euler;
- euler.do_step( ham_sys() ,
- std::make_pair( f.q.begin() + 1 , f.q.begin() + 4 ) ,
- std::make_pair( f.p.begin() + 3 , f.p.begin() + 6 ) ,
- 0.0 , 0.1 );
- CHECK_VALUES( f.q , 0.0 , 1.3 , 2.4 , 3.5 , 4.0 , 5.0 );
- CHECK_VALUES( f.p , 0.0 , 1.0 , 2.0 , 3.1 , 4.2 , 5.3 );
- }
- BOOST_AUTO_TEST_CASE( symplectic_euler_coor_and_mom_func )
- {
- vector_fixture f;
- boost::numeric::odeint::symplectic_euler< state_type > euler;
- euler.do_step( std::make_pair( ham_sys() , ham_sys() ) ,
- std::make_pair( f.q.begin() + 1 , f.q.begin() + 4 ) ,
- std::make_pair( f.p.begin() + 3 , f.p.begin() + 6 ) ,
- 0.0 , 0.1 );
- CHECK_VALUES( f.q , 0.0 , 1.1 , 2.2 , 3.3 , 4.0 , 5.0 );
- CHECK_VALUES( f.p , 0.0 , 1.0 , 2.0 , 3.1 , 4.2 , 5.3 );
- }
- BOOST_AUTO_TEST_CASE( dense_output_euler_with_ranges )
- {
- using namespace boost::numeric::odeint;
- vector_fixture f;
- dense_output_runge_kutta< euler< state_type > > stepper;
- stepper.initialize( std::make_pair( f.in.begin() + 1, f.in.begin() + 4 ) , 0.0 , 0.1 );
- stepper.do_step( system1() );
- stepper.calc_state( 0.05 , std::make_pair( f.in.begin() + 1 ,f.in.begin() +4 ) );
- CHECK_VALUES( f.in , 0.0 , 1.05 , 2.1 , 3.15 , 4.0 , 5.0 );
- }
- BOOST_AUTO_TEST_CASE( dense_output_dopri5_with_ranges )
- {
- using namespace boost::numeric::odeint;
- vector_fixture f;
- dense_output_runge_kutta<
- controlled_runge_kutta<
- runge_kutta_dopri5< state_type >
- > > stepper;
- stepper.initialize( std::make_pair( f.in.begin() + 1, f.in.begin() + 4 ) , 0.0 , 0.1 );
- stepper.do_step( system2() );
- stepper.calc_state( 0.05 , std::make_pair( f.in.begin() + 1 ,f.in.begin() +4 ) );
- CHECK_VALUES( f.in , 0.0 , 1.05 , 2.1 , 3.15 , 4.0 , 5.0 );
- }
- BOOST_AUTO_TEST_SUITE_END()
|