293 lines
8.8 KiB
C++
293 lines
8.8 KiB
C++
/*
|
|
[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()
|