path: root/libs/numeric/odeint/test/stepper_with_ranges.cpp

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