test/step_size_limitation.cpp

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