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