boost_mirror/odeint
1
0
Fork 0
mirror of https://github.com/boostorg/odeint.git synced 2025-05-09 23:24:01 +00:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'master' of github.com:headmyshoulder/odeint-v2 into develop

Browse Source
This commit is contained in:
Mario Mulansky 2019-04-20 11:49:37 -07:00
commit fc43a2ef44
36 changed files with 1741 additions and 142 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
.travis.yml
View File

@ -19,7 +19,7 @@ addons:
env:
- CXXSTD=""
- CXXSTD="cxxflags='-std=c++0x'"
- CXXSTD="-std=c++0x"
# For now disable gcc on osx as g++4.8 is not yet available
matrix:
@ -28,21 +28,19 @@ matrix:
compiler: gcc
before_install:
# 1.55: http://sourceforge.net/projects/boost/files/boost/1.55.0/boost_1_55_0.tar.bz2/download
# 1.57: http://sourceforge.net/projects/boost/files/boost/1.57.0/boost_1_57_0.tar.bz2/download
# 1.58: http://downloads.sourceforge.net/project/boost/boost/1.58.0/boost_1_58_0.tar.bz2\?r\=\&ts\=1435589970\&use_mirror\=garr
- wget http://downloads.sourceforge.net/project/boost/boost/1.58.0/boost_1_58_0.tar.bz2\?r\=\&ts\=1435589970\&use_mirror\=garr -O /tmp/boost.tar.bz2
# 1.63: https://sourceforge.net/projects/boost/files/boost/1.63.0/boost_1_63_0.tar.bz2/download?use_mirror=superb-dca2
- wget https://sourceforge.net/projects/boost/files/boost/1.63.0/boost_1_63_0.tar.bz2/download?use_mirror=superb-dca2 -O /tmp/boost.tar.bz2
- tar jxf /tmp/boost.tar.bz2
- mv boost_1_58_0 $PWD/boost-trunk
# patch the boost build system - not neccessary with 1.58 anymore
# - patch $PWD/boost-trunk/tools/build/v2/build/toolset.jam toolset.jam.patch
- mv boost_1_63_0 $PWD/boost-trunk
- export BOOST_ROOT="$PWD/boost-trunk"
- cd $BOOST_ROOT
- ./bootstrap.sh
- cd $TRAVIS_BUILD_DIR
- if [ "$TRAVIS_OS_NAME" = "osx" ] && [ "$CC" = "gcc" ]; then export CC=gcc-4.8; fi
- echo $CC
- $CC --version
script:
- $BOOST_ROOT/b2 toolset=$CC $CXXSTD
- $BOOST_ROOT/b2 toolset=$CC cxxflags='$CXXSTD' -j2

2
README.md
View File

@ -1,3 +1,3 @@
[![Bitdeli Badge](https://d2weczhvl823v0.cloudfront.net/headmyshoulder/odeint-v2/trend.png)](https://bitdeli.com/free "Bitdeli Badge")
[![Build Status](https://travis-ci.org/headmyshoulder/odeint-v2.svg?branch=master)](https://travis-ci.org/headmyshoulder/odeint-v2)
odeint is a highly flexible library for solving ordinary differential equations.

3
include/boost/numeric/odeint.hpp
View File

@ -51,6 +51,9 @@
#include <boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp>
#include <boost/numeric/odeint/stepper/adaptive_adams_bashforth_moulton.hpp>
#include <boost/numeric/odeint/stepper/controlled_adams_bashforth_moulton.hpp>
#include <boost/numeric/odeint/stepper/implicit_euler.hpp>
#include <boost/numeric/odeint/stepper/rosenbrock4.hpp>
#include <boost/numeric/odeint/stepper/rosenbrock4_controller.hpp>

18
include/boost/numeric/odeint/algebra/multi_array_algebra.hpp
View File

@ -19,6 +19,8 @@
#define BOOST_NUMERIC_ODEINT_ALGEBRA_MULTI_ARRAY_ALGEBRA_HPP_DEFINED
#include <boost/multi_array.hpp>
#include <boost/numeric/odeint/algebra/detail/for_each.hpp>
#include <boost/numeric/odeint/algebra/detail/norm_inf.hpp>
#include <boost/numeric/odeint/algebra/norm_result_type.hpp>
@ -121,7 +123,7 @@ struct multi_array_algebra
{
detail::for_each15( s1.data() , s1.data() + s1.num_elements() , s2.data() , s3.data() , s4.data() , s5.data() , s6.data() , s7.data() , s8.data() , s9.data() , s10.data() , s11.data() , s12.data() , s13.data() , s14.data() , s15.data() , op );
}
template< typename S >
static typename norm_result_type<S>::type norm_inf( const S &s )
{
@ -129,15 +131,11 @@ struct multi_array_algebra
}
};
// template< class T , size_t N >
// struct algebra_dispatcher< boost::array< T , N > >
// {
// typedef array_algebra algebra_type;
// };
template< class T , size_t N >
struct algebra_dispatcher< boost::multi_array< T , N > >
{
typedef multi_array_algebra algebra_type;
};
} // namespace odeint

45
include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp vendored
View File

@ -25,42 +25,36 @@
// from odeint
// (that is, it lets odeint treat the eigen matrices correctly, knowing
// how to add, multiply, compute the norm, etc)
namespace Eigen {
template<typename D>
inline const
typename Eigen::CwiseUnaryOp<
typename Eigen::internal::scalar_add_op<
typename Eigen::internal::traits<D>::Scalar>,
const D >
typename Eigen::CwiseBinaryOp<
internal::scalar_sum_op<typename internal::traits<D>::Scalar>,
typename DenseBase<D>::ConstantReturnType,
const D>
operator+(const typename Eigen::MatrixBase<D> &m,
const typename Eigen::internal::traits<D>::Scalar &s) {
return Eigen::CwiseUnaryOp<
typename Eigen::internal::scalar_add_op<
typename Eigen::internal::traits<D>::Scalar>,
const D >(m.derived(),Eigen::internal::scalar_add_op<
typename Eigen::internal::traits<D>::Scalar>(s));
return CwiseBinaryOp<
internal::scalar_sum_op<typename internal::traits<D>::Scalar>,
typename DenseBase<D>::ConstantReturnType,
const D>(DenseBase<D>::Constant(m.rows(), m.cols(), s), m.derived());
}
template<typename D>
inline const
typename Eigen::CwiseUnaryOp<
typename Eigen::internal::scalar_add_op<
typename Eigen::internal::traits<D>::Scalar>,
const D >
typename Eigen::CwiseBinaryOp<
internal::scalar_sum_op<typename internal::traits<D>::Scalar>,
typename DenseBase<D>::ConstantReturnType,
const D>
operator+(const typename Eigen::internal::traits<D>::Scalar &s,
const typename Eigen::MatrixBase<D> &m) {
return Eigen::CwiseUnaryOp<
typename Eigen::internal::scalar_add_op<
typename Eigen::internal::traits<D>::Scalar>,
const D >(m.derived(),Eigen::internal::scalar_add_op<
typename Eigen::internal::traits<D>::Scalar>(s));
const typename Eigen::MatrixBase<D> &m) {
return CwiseBinaryOp<
internal::scalar_sum_op<typename internal::traits<D>::Scalar>,
typename DenseBase<D>::ConstantReturnType,
const D>(DenseBase<D>::Constant(m.rows(), m.cols(), s), m.derived());
}
template<typename D1,typename D2>
inline const
typename Eigen::CwiseBinaryOp<
@ -82,14 +76,9 @@ abs( const Eigen::MatrixBase< D > &m ) {
return m.cwiseAbs();
}
} // end Eigen namespace
namespace boost {
namespace numeric {
namespace odeint {

1
include/boost/numeric/odeint/integrate/integrate_n_steps.hpp
View File

@ -23,6 +23,7 @@
#include <boost/numeric/odeint/stepper/stepper_categories.hpp>
#include <boost/numeric/odeint/integrate/null_observer.hpp>
#include <boost/numeric/odeint/integrate/detail/integrate_n_steps.hpp>
#include <boost/numeric/odeint/integrate/check_adapter.hpp>
namespace boost {
namespace numeric {

13
include/boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp
View File

@ -150,6 +150,12 @@ public :
}
void reset(void)
{
m_adams_bashforth.reset();
}
private:
@ -175,7 +181,7 @@ private:
{
m_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) );
m_adams_bashforth.do_step( system , in , t , m_x.m_v , dt );
m_adams_moulton.do_step( system , in , m_x.m_v , t , out , dt , m_adams_bashforth.step_storage() );
m_adams_moulton.do_step( system , in , m_x.m_v , t+dt , out , dt , m_adams_bashforth.step_storage() );
}
else
{
@ -293,6 +299,11 @@ private:
* \param dt The step size.
*/
/**
* \fn adams_bashforth_moulton::reset( void )
* \brief Resets the internal buffers of the stepper.
*/
} // odeint
} // numeric

237
include/boost/numeric/odeint/stepper/adaptive_adams_bashforth_moulton.hpp Normal file
View File

@ -0,0 +1,237 @@
/*
boost/numeric/odeint/stepper/detail/adaptive_adams_bashforth_moulton.hpp
[begin_description]
Implemetation of an adaptive adams bashforth moulton stepper.
Used as the stepper for the controlled adams bashforth moulton stepper.
[end_description]
Copyright 2017 Valentin Noah Hartmann
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)
*/
#ifndef BOOST_NUMERIC_ODEINT_STEPPER_ADAPTIVE_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_STEPPER_ADAPTIVE_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED
#include <boost/numeric/odeint/stepper/detail/adaptive_adams_coefficients.hpp>
#include <boost/numeric/odeint/util/unwrap_reference.hpp>
#include <boost/numeric/odeint/util/bind.hpp>
#include <boost/numeric/odeint/util/copy.hpp>
#include <boost/numeric/odeint/algebra/default_operations.hpp>
#include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp>
#include <boost/numeric/odeint/algebra/operations_dispatcher.hpp>
#include <boost/numeric/odeint/util/state_wrapper.hpp>
#include <boost/numeric/odeint/util/is_resizeable.hpp>
#include <boost/numeric/odeint/util/resizer.hpp>
#include <boost/numeric/odeint/stepper/stepper_categories.hpp>
#include <boost/numeric/odeint/stepper/base/algebra_stepper_base.hpp>
#include <boost/numeric/odeint/stepper/detail/rotating_buffer.hpp>
namespace boost {
namespace numeric {
namespace odeint {
template<
size_t Steps,
class State,
class Value = double,
class Deriv = State,
class Time = Value,
class Algebra = typename algebra_dispatcher< State >::algebra_type ,
class Operations = typename operations_dispatcher< State >::operations_type,
class Resizer = initially_resizer
>
class adaptive_adams_bashforth_moulton: public algebra_stepper_base< Algebra , Operations >
{
public:
static const size_t steps = Steps;
typedef unsigned short order_type;
static const order_type order_value = steps;
typedef State state_type;
typedef Value value_type;
typedef Deriv deriv_type;
typedef Time time_type;
typedef state_wrapper< state_type > wrapped_state_type;
typedef state_wrapper< deriv_type > wrapped_deriv_type;
typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type;
typedef typename algebra_stepper_base_type::algebra_type algebra_type;
typedef typename algebra_stepper_base_type::operations_type operations_type;
typedef Resizer resizer_type;
typedef error_stepper_tag stepper_category;
typedef detail::adaptive_adams_coefficients< Steps , Deriv , Value , Time , Algebra , Operations , Resizer > coeff_type;
typedef adaptive_adams_bashforth_moulton< Steps , State , Value , Deriv , Time , Algebra , Operations , Resizer > stepper_type;
order_type order() const { return order_value; };
order_type stepper_order() const { return order_value + 1; };
order_type error_order() const { return order_value; };
adaptive_adams_bashforth_moulton( const algebra_type &algebra = algebra_type() )
:algebra_stepper_base_type( algebra ), m_coeff(),
m_dxdt_resizer(), m_xnew_resizer(), m_xerr_resizer()
{};
template< class System >
void do_step(System system, state_type &inOut, time_type t, time_type dt )
{
m_xnew_resizer.adjust_size( inOut , detail::bind( &stepper_type::template resize_xnew_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
do_step(system, inOut, t, m_xnew.m_v, dt, m_xerr.m_v);
boost::numeric::odeint::copy( m_xnew.m_v , inOut);
};
template< class System >
void do_step(System system, const state_type &in, time_type t, state_type &out, time_type dt )
{
do_step(system, in, t, out, dt, m_xerr.m_v);
};
template< class System >
void do_step(System system, state_type &inOut, time_type t, time_type dt, state_type &xerr)
{
m_xnew_resizer.adjust_size( inOut , detail::bind( &stepper_type::template resize_xnew_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
do_step(system, inOut, t, m_xnew.m_v, dt, xerr);
boost::numeric::odeint::copy( m_xnew.m_v , inOut);
};
template< class System >
void do_step(System system, const state_type &in, time_type t, state_type &out, time_type dt , state_type &xerr)
{
do_step_impl(system, in, t, out, dt, xerr);
system(out, m_dxdt.m_v, t+dt);
m_coeff.do_step(m_dxdt.m_v);
m_coeff.confirm();
if(m_coeff.m_eo < order_value)
{
m_coeff.m_eo ++;
}
};
template< class ExplicitStepper, class System >
void initialize(ExplicitStepper stepper, System system, state_type &inOut, time_type &t, time_type dt)
{
reset();
dt = dt/static_cast< time_type >(order_value);
m_dxdt_resizer.adjust_size( inOut , detail::bind( &stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
system( inOut , m_dxdt.m_v , t );
for( size_t i=0 ; i<order_value; ++i )
{
stepper.do_step_dxdt_impl( system, inOut, m_dxdt.m_v, t, dt );
system( inOut , m_dxdt.m_v , t + dt);
m_coeff.predict(t, dt);
m_coeff.do_step(m_dxdt.m_v);
m_coeff.confirm();
t += dt;
if(m_coeff.m_eo < order_value)
{
++m_coeff.m_eo;
}
}
};
template< class System >
void initialize(System system, state_type &inOut, time_type &t, time_type dt)
{
reset();
dt = dt/static_cast< time_type >(order_value);
for(size_t i=0; i<order_value; ++i)
{
this->do_step(system, inOut, t, dt);
t += dt;
};
};
template< class System >
void do_step_impl(System system, const state_type & in, time_type t, state_type & out, time_type &dt, state_type &xerr)
{
size_t eO = m_coeff.m_eo;
m_xerr_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_xerr_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
m_dxdt_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
m_coeff.predict(t, dt);
if (m_coeff.m_steps_init == 1)
{
system(in, m_dxdt.m_v, t);
m_coeff.do_step(m_dxdt.m_v, 1);
}
boost::numeric::odeint::copy( in , out );
for(size_t i=0; i<eO; ++i)
{
this->m_algebra.for_each3(out, out, m_coeff.phi[1][i].m_v,
typename Operations::template scale_sum2<double, double>(1.0, dt*m_coeff.g[i]*m_coeff.beta[0][i]));
}
system(out, m_dxdt.m_v, t+dt);
m_coeff.do_step(m_dxdt.m_v);
this->m_algebra.for_each3(out, out, m_coeff.phi[0][eO].m_v,
typename Operations::template scale_sum2<double, double>(1.0, dt*m_coeff.g[eO]));
// error for current order
this->m_algebra.for_each2(xerr, m_coeff.phi[0][eO].m_v,
typename Operations::template scale_sum1<double>(dt*(m_coeff.g[eO])));
};
const coeff_type& coeff() const { return m_coeff; };
coeff_type & coeff() { return m_coeff; };
void reset() { m_coeff.reset(); };
const deriv_type & dxdt() const { return m_dxdt.m_v; };
private:
template< class StateType >
bool resize_dxdt_impl( const StateType &x )
{
return adjust_size_by_resizeability( m_dxdt, x, typename is_resizeable<deriv_type>::type() );
};
template< class StateType >
bool resize_xnew_impl( const StateType &x )
{
return adjust_size_by_resizeability( m_xnew, x, typename is_resizeable<state_type>::type() );
};
template< class StateType >
bool resize_xerr_impl( const StateType &x )
{
return adjust_size_by_resizeability( m_xerr, x, typename is_resizeable<state_type>::type() );
};
coeff_type m_coeff;
resizer_type m_dxdt_resizer;
resizer_type m_xnew_resizer;
resizer_type m_xerr_resizer;
wrapped_deriv_type m_dxdt;
wrapped_state_type m_xnew;
wrapped_state_type m_xerr;
};
} // odeint
} // numeric
} // boost
#endif

23
include/boost/numeric/odeint/stepper/bulirsch_stoer.hpp
View File

@ -98,6 +98,7 @@ public:
m_interval_sequence( m_k_max+1 ) ,
m_coeff( m_k_max+1 ) ,
m_cost( m_k_max+1 ) ,
m_facmin_table( m_k_max+1 ) ,
m_table( m_k_max ) ,
STEPFAC1( 0.65 ) , STEPFAC2( 0.94 ) , STEPFAC3( 0.02 ) , STEPFAC4( 4.0 ) , KFAC1( 0.8 ) , KFAC2( 0.9 )
{
@ -112,21 +113,14 @@ public:
else
m_cost[i] = m_cost[i-1] + m_interval_sequence[i];
m_coeff[i].resize(i);
m_facmin_table[i] = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , static_cast< value_type >(1) / static_cast< value_type >( 2*i+1 ) );
for( size_t k = 0 ; k < i ; ++k )
{
const value_type r = static_cast< value_type >( m_interval_sequence[i] ) / static_cast< value_type >( m_interval_sequence[k] );
m_coeff[i][k] = 1.0 / ( r*r - static_cast< value_type >( 1.0 ) ); // coefficients for extrapolation
}
// crude estimate of optimal order
m_current_k_opt = 4;
/* no calculation because log10 might not exist for value_type!
const value_type logfact( -log10( max BOOST_PREVENT_MACRO_SUBSTITUTION( eps_rel , static_cast< value_type >(1.0E-12) ) ) * 0.6 + 0.5 );
m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>( 1 ) , min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>( m_k_max-1 ) , logfact ));
*/
}
reset();
}
@ -282,7 +276,7 @@ public:
{
m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<int>(m_k_max-1) , static_cast<int>(m_current_k_opt)+1 );
new_h = h_opt[k];
new_h *= m_cost[m_current_k_opt]/m_cost[k];
new_h *= static_cast<value_type>(m_cost[m_current_k_opt])/static_cast<value_type>(m_cost[k]);
} else
new_h = h_opt[m_current_k_opt];
break;
@ -344,6 +338,12 @@ public:
{
m_first = true;
m_last_step_rejected = false;
// crude estimate of optimal order
m_current_k_opt = 4;
/* no calculation because log10 might not exist for value_type!
const value_type logfact( -log10( max BOOST_PREVENT_MACRO_SUBSTITUTION( eps_rel , static_cast< value_type >(1.0E-12) ) ) * 0.6 + 0.5 );
m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>( 1 ) , min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>( m_k_max-1 ) , logfact ));
*/
}
@ -417,7 +417,7 @@ private:
BOOST_USING_STD_MAX();
using std::pow;
value_type expo( 1.0/(2*k+1) );
value_type facmin = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , expo );
value_type facmin = m_facmin_table[k];
value_type fac;
if (error == 0.0)
fac=1.0/facmin;
@ -506,6 +506,7 @@ private:
int_vector m_interval_sequence; // stores the successive interval counts
value_matrix m_coeff;
int_vector m_cost; // costs for interval count
value_vector m_facmin_table; // for precomputed facmin to save pow calls
state_table_type m_table; // sequence of states for extrapolation

17
include/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp
View File

@ -110,6 +110,7 @@ public:
m_interval_sequence( m_k_max+1 ) ,
m_coeff( m_k_max+1 ) ,
m_cost( m_k_max+1 ) ,
m_facmin_table( m_k_max+1 ) ,
m_table( m_k_max ) ,
m_mp_states( m_k_max+1 ) ,
m_derivs( m_k_max+1 ) ,
@ -125,9 +126,13 @@ public:
m_interval_sequence[i] = 2 + 4*i; // 2 6 10 14 ...
m_derivs[i].resize( m_interval_sequence[i] );
if( i == 0 )
{
m_cost[i] = m_interval_sequence[i];
else
} else
{
m_cost[i] = m_cost[i-1] + m_interval_sequence[i];
}
m_facmin_table[i] = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , static_cast< value_type >(1) / static_cast< value_type >( 2*i+1 ) );
m_coeff[i].resize(i);
for( size_t k = 0 ; k < i ; ++k )
{
@ -386,6 +391,11 @@ public:
}
protected:
time_type m_max_dt;
private:
template< class StateInOut , class StateVector >
@ -429,7 +439,7 @@ private:
using std::pow;
value_type expo = static_cast<value_type>(1)/(m_interval_sequence[k-1]);
value_type facmin = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , expo );
value_type facmin = m_facmin_table[k];
value_type fac;
if (error == 0.0)
fac = static_cast<value_type>(1)/facmin;
@ -662,8 +672,6 @@ private:
default_error_checker< value_type, algebra_type , operations_type > m_error_checker;
modified_midpoint_dense_out< state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > m_midpoint;
time_type m_max_dt;
bool m_control_interpolation;
bool m_last_step_rejected;
@ -692,6 +700,7 @@ private:
int_vector m_interval_sequence; // stores the successive interval counts
value_matrix m_coeff;
int_vector m_cost; // costs for interval count
value_vector m_facmin_table; // for precomputed facmin to save pow calls
state_vector_type m_table; // sequence of states for extrapolation

322
include/boost/numeric/odeint/stepper/controlled_adams_bashforth_moulton.hpp Normal file
View File

@ -0,0 +1,322 @@
/*
boost/numeric/odeint/stepper/detail/controlled_adams_bashforth_moulton.hpp
[begin_description]
Implemetation of an controlled adams bashforth moulton stepper.
[end_description]
Copyright 2017 Valentin Noah Hartmann
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)
*/
#ifndef BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLED_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLED_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED
#include <boost/numeric/odeint/stepper/stepper_categories.hpp>
#include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
#include <boost/numeric/odeint/stepper/adaptive_adams_bashforth_moulton.hpp>
#include <boost/numeric/odeint/stepper/detail/pid_step_adjuster.hpp>
#include <boost/numeric/odeint/util/unwrap_reference.hpp>
#include <boost/numeric/odeint/util/is_resizeable.hpp>
#include <boost/numeric/odeint/util/resizer.hpp>
#include <boost/numeric/odeint/util/copy.hpp>
#include <boost/numeric/odeint/util/bind.hpp>
#include <iostream>
namespace boost {
namespace numeric {
namespace odeint {
template<
size_t MaxOrder,
class State,
class Value = double,
class Algebra = typename algebra_dispatcher< State >::algebra_type
>
class default_order_adjuster
{
public:
typedef State state_type;
typedef Value value_type;
typedef state_wrapper< state_type > wrapped_state_type;
typedef Algebra algebra_type;
default_order_adjuster( const algebra_type &algebra = algebra_type() )
: m_algebra( algebra )
{};
size_t adjust_order(size_t order, size_t init, boost::array<wrapped_state_type, 4> &xerr)
{
using std::abs;
value_type errc = abs(m_algebra.norm_inf(xerr[2].m_v));
value_type errm1 = 3*errc;
value_type errm2 = 3*errc;
if(order > 2)
{
errm2 = abs(m_algebra.norm_inf(xerr[0].m_v));
}
if(order >= 2)
{
errm1 = abs(m_algebra.norm_inf(xerr[1].m_v));
}
size_t o_new = order;
if(order == 2 && errm1 <= 0.5*errc)
{
o_new = order - 1;
}
else if(order > 2 && errm2 < errc && errm1 < errc)
{
o_new = order - 1;
}
if(init < order)
{
return order+1;
}
else if(o_new == order - 1)
{
return order-1;
}
else if(order <= MaxOrder)
{
value_type errp = abs(m_algebra.norm_inf(xerr[3].m_v));
if(order > 1 && errm1 < errc && errp)
{
return order-1;
}
else if(order < MaxOrder && errp < (0.5-0.25*order/MaxOrder) * errc)
{
return order+1;
}
}
return order;
};
private:
algebra_type m_algebra;
};
template<
class ErrorStepper,
class StepAdjuster = detail::pid_step_adjuster< typename ErrorStepper::state_type,
typename ErrorStepper::value_type,
typename ErrorStepper::deriv_type,
typename ErrorStepper::time_type,
typename ErrorStepper::algebra_type,
typename ErrorStepper::operations_type,
detail::H211PI
>,
class OrderAdjuster = default_order_adjuster< ErrorStepper::order_value,
typename ErrorStepper::state_type,
typename ErrorStepper::value_type,
typename ErrorStepper::algebra_type
>,
class Resizer = initially_resizer
>
class controlled_adams_bashforth_moulton
{
public:
typedef ErrorStepper stepper_type;
static const typename stepper_type::order_type order_value = stepper_type::order_value;
typedef typename stepper_type::state_type state_type;
typedef typename stepper_type::value_type value_type;
typedef typename stepper_type::deriv_type deriv_type;
typedef typename stepper_type::time_type time_type;
typedef typename stepper_type::algebra_type algebra_type;
typedef typename stepper_type::operations_type operations_type;
typedef Resizer resizer_type;
typedef StepAdjuster step_adjuster_type;
typedef OrderAdjuster order_adjuster_type;
typedef controlled_stepper_tag stepper_category;
typedef typename stepper_type::wrapped_state_type wrapped_state_type;
typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type;
typedef boost::array< wrapped_state_type , 4 > error_storage_type;
typedef typename stepper_type::coeff_type coeff_type;
typedef controlled_adams_bashforth_moulton< ErrorStepper , StepAdjuster , OrderAdjuster , Resizer > controlled_stepper_type;
controlled_adams_bashforth_moulton(step_adjuster_type step_adjuster = step_adjuster_type())
:m_stepper(),
m_dxdt_resizer(), m_xerr_resizer(), m_xnew_resizer(),
m_step_adjuster( step_adjuster ), m_order_adjuster()
{};
template< class ExplicitStepper, class System >
void initialize(ExplicitStepper stepper, System system, state_type &inOut, time_type &t, time_type dt)
{
m_stepper.initialize(stepper, system, inOut, t, dt);
};
template< class System >
void initialize(System system, state_type &inOut, time_type &t, time_type dt)
{
m_stepper.initialize(system, inOut, t, dt);
};
template< class ExplicitStepper, class System >
void initialize_controlled(ExplicitStepper stepper, System system, state_type &inOut, time_type &t, time_type &dt)
{
reset();
coeff_type &coeff = m_stepper.coeff();
m_dxdt_resizer.adjust_size( inOut , detail::bind( &controlled_stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
controlled_step_result res = fail;
for( size_t i=0 ; i<order_value; ++i )
{
do
{
res = stepper.try_step( system, inOut, t, dt );
}
while(res != success);
system( inOut , m_dxdt.m_v , t );
coeff.predict(t-dt, dt);
coeff.do_step(m_dxdt.m_v);
coeff.confirm();
if(coeff.m_eo < order_value)
{
++coeff.m_eo;
}
}
}
template< class System >
controlled_step_result try_step(System system, state_type & inOut, time_type &t, time_type &dt)
{
m_xnew_resizer.adjust_size( inOut , detail::bind( &controlled_stepper_type::template resize_xnew_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
controlled_step_result res = try_step(system, inOut, t, m_xnew.m_v, dt);
if(res == success)
{
boost::numeric::odeint::copy( m_xnew.m_v , inOut);
}
return res;
};
template< class System >
controlled_step_result try_step(System system, const state_type & in, time_type &t, state_type & out, time_type &dt)
{
m_xerr_resizer.adjust_size( in , detail::bind( &controlled_stepper_type::template resize_xerr_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
m_dxdt_resizer.adjust_size( in , detail::bind( &controlled_stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
m_stepper.do_step_impl(system, in, t, out, dt, m_xerr[2].m_v);
coeff_type &coeff = m_stepper.coeff();
time_type dtPrev = dt;
dt = m_step_adjuster.adjust_stepsize(coeff.m_eo, dt, m_xerr[2].m_v, out, m_stepper.dxdt() );
if(dt / dtPrev >= step_adjuster_type::threshold())
{
system(out, m_dxdt.m_v, t+dtPrev);
coeff.do_step(m_dxdt.m_v);
coeff.confirm();
t += dtPrev;
size_t eo = coeff.m_eo;
// estimate errors for next step
double factor = 1;
algebra_type m_algebra;
m_algebra.for_each2(m_xerr[2].m_v, coeff.phi[1][eo].m_v,
typename operations_type::template scale_sum1<double>(factor*dt*(coeff.gs[eo])));
if(eo > 1)
{
m_algebra.for_each2(m_xerr[1].m_v, coeff.phi[1][eo-1].m_v,
typename operations_type::template scale_sum1<double>(factor*dt*(coeff.gs[eo-1])));
}
if(eo > 2)
{
m_algebra.for_each2(m_xerr[0].m_v, coeff.phi[1][eo-2].m_v,
typename operations_type::template scale_sum1<double>(factor*dt*(coeff.gs[eo-2])));
}
if(eo < order_value && coeff.m_eo < coeff.m_steps_init-1)
{
m_algebra.for_each2(m_xerr[3].m_v, coeff.phi[1][eo+1].m_v,
typename operations_type::template scale_sum1<double>(factor*dt*(coeff.gs[eo+1])));
}
// adjust order
coeff.m_eo = m_order_adjuster.adjust_order(coeff.m_eo, coeff.m_steps_init-1, m_xerr);
return success;
}
else
{
return fail;
}
};
void reset() { m_stepper.reset(); };
private:
template< class StateType >
bool resize_dxdt_impl( const StateType &x )
{
return adjust_size_by_resizeability( m_dxdt, x, typename is_resizeable<deriv_type>::type() );
};
template< class StateType >
bool resize_xerr_impl( const StateType &x )
{
bool resized( false );
for(size_t i=0; i<m_xerr.size(); ++i)
{
resized |= adjust_size_by_resizeability( m_xerr[i], x, typename is_resizeable<state_type>::type() );
}
return resized;
};
template< class StateType >
bool resize_xnew_impl( const StateType &x )
{
return adjust_size_by_resizeability( m_xnew, x, typename is_resizeable<state_type>::type() );
};
stepper_type m_stepper;
wrapped_deriv_type m_dxdt;
error_storage_type m_xerr;
wrapped_state_type m_xnew;
resizer_type m_dxdt_resizer;
resizer_type m_xerr_resizer;
resizer_type m_xnew_resizer;
step_adjuster_type m_step_adjuster;
order_adjuster_type m_order_adjuster;
};
} // odeint
} // numeric
} // boost
#endif

26
include/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp
View File

@ -143,7 +143,7 @@ public:
error = max BOOST_PREVENT_MACRO_SUBSTITUTION (
static_cast<value_type>( pow( static_cast<value_type>(5.0) , -static_cast<value_type>(stepper_order) ) ) ,
error);
time_type dt_old = dt;
// time_type dt_old = dt; unused variable warning
//error too small - increase dt and keep the evolution and limit scaling factor to 5.0
dt *= static_cast<value_type>(9)/static_cast<value_type>(10) *
pow(error, static_cast<value_type>(-1) / stepper_order);
@ -163,7 +163,7 @@ public:
time_type get_max_dt() { return m_max_dt; }
private:
protected:
time_type m_max_dt;
};
@ -411,10 +411,11 @@ public:
template< class System , class StateIn , class DerivIn , class StateOut >
controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , time_type &dt )
{
if( !m_step_adjuster.check_step_size_limit(dt) )
unwrapped_step_adjuster &step_adjuster = m_step_adjuster;
if( !step_adjuster.check_step_size_limit(dt) )
{
// given dt was above step size limit - adjust and return fail;
dt = m_step_adjuster.get_max_dt();
dt = step_adjuster.get_max_dt();
return fail;
}
@ -428,13 +429,13 @@ public:
if( max_rel_err > 1.0 )
{
// error too big, decrease step size and reject this step
dt = m_step_adjuster.decrease_step(dt, max_rel_err, m_stepper.error_order());
dt = step_adjuster.decrease_step(dt, max_rel_err, m_stepper.error_order());
return fail;
} else
{
// otherwise, increase step size and accept
t += dt;
dt = m_step_adjuster.increase_step(dt, max_rel_err, m_stepper.stepper_order());
dt = step_adjuster.increase_step(dt, max_rel_err, m_stepper.stepper_order());
return success;
}
}
@ -505,6 +506,7 @@ private:
stepper_type m_stepper;
error_checker_type m_error_checker;
step_adjuster_type m_step_adjuster;
typedef typename unwrap_reference< step_adjuster_type >::type unwrapped_step_adjuster;
resizer_type m_dxdt_resizer;
resizer_type m_xerr_resizer;
@ -584,7 +586,7 @@ public:
const step_adjuster_type &step_adjuster = step_adjuster_type() ,
const stepper_type &stepper = stepper_type()
)
: m_stepper( stepper ) , m_error_checker( error_checker ) , m_step_adjuster(step_adjuster) ,
: m_stepper( stepper ) , m_error_checker( error_checker ) , m_step_adjuster(step_adjuster) ,
m_first_call( true )
{ }
@ -751,10 +753,11 @@ public:
controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt_in , time_type &t ,
StateOut &out , DerivOut &dxdt_out , time_type &dt )
{
if( !m_step_adjuster.check_step_size_limit(dt) )
unwrapped_step_adjuster &step_adjuster = m_step_adjuster;
if( !step_adjuster.check_step_size_limit(dt) )
{
// given dt was above step size limit - adjust and return fail;
dt = m_step_adjuster.get_max_dt();
dt = step_adjuster.get_max_dt();
return fail;
}
@ -770,12 +773,12 @@ public:
if( max_rel_err > 1.0 )
{
// error too big, decrease step size and reject this step
dt = m_step_adjuster.decrease_step(dt, max_rel_err, m_stepper.error_order());
dt = step_adjuster.decrease_step(dt, max_rel_err, m_stepper.error_order());
return fail;
}
// otherwise, increase step size and accept
t += dt;
dt = m_step_adjuster.increase_step(dt, max_rel_err, m_stepper.stepper_order());
dt = step_adjuster.increase_step(dt, max_rel_err, m_stepper.stepper_order());
return success;
}
@ -903,6 +906,7 @@ private:
stepper_type m_stepper;
error_checker_type m_error_checker;
step_adjuster_type m_step_adjuster;
typedef typename unwrap_reference< step_adjuster_type >::type unwrapped_step_adjuster;
resizer_type m_dxdt_resizer;
resizer_type m_xerr_resizer;

2
include/boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp
View File

@ -126,7 +126,7 @@ public:
m_t_old = m_t;
m_t += m_dt;
toggle_current_state();
return std::make_pair( m_t_old , m_dt );
return std::make_pair( m_t_old , m_t );
}
/*

207
include/boost/numeric/odeint/stepper/detail/adaptive_adams_coefficients.hpp Normal file
View File

@ -0,0 +1,207 @@
/*
boost/numeric/odeint/stepper/detail/adaptive_adams_coefficients.hpp
[begin_description]
Calculation of the coefficients for the adaptive adams stepper.
[end_description]
Copyright 2017 Valentin Noah Hartmann
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)
*/
#ifndef BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_ADAPTIVE_ADAMS_COEFFICIENTS_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_ADAPTIVE_ADAMS_COEFFICIENTS_HPP_INCLUDED
#include <boost/numeric/odeint/stepper/detail/rotating_buffer.hpp>
#include <boost/numeric/odeint/util/state_wrapper.hpp>
#include <boost/numeric/odeint/util/is_resizeable.hpp>
#include <boost/numeric/odeint/util/resizer.hpp>
#include <boost/numeric/odeint/util/unwrap_reference.hpp>
#include <boost/numeric/odeint/util/bind.hpp>
#include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp>
#include <boost/numeric/odeint/algebra/operations_dispatcher.hpp>
#include <boost/array.hpp>
namespace boost {
namespace numeric {
namespace odeint {
namespace detail {
template<
size_t Steps,
class Deriv,
class Value = double,
class Time = double,
class Algebra = typename algebra_dispatcher< Deriv >::algebra_type,
class Operations = typename operations_dispatcher< Deriv >::operations_type,
class Resizer = initially_resizer
>
class adaptive_adams_coefficients
{
public:
static const size_t steps = Steps;
typedef unsigned short order_type;
static const order_type order_value = steps;
typedef Value value_type;
typedef Deriv deriv_type;
typedef Time time_type;
typedef state_wrapper< deriv_type > wrapped_deriv_type;
typedef rotating_buffer< time_type , steps+1 > time_storage_type;
typedef Algebra algebra_type;
typedef Operations operations_type;
typedef Resizer resizer_type;
typedef adaptive_adams_coefficients< Steps , Deriv , Value , Time , Algebra , Operations , Resizer > aac_type;
adaptive_adams_coefficients( const algebra_type &algebra = algebra_type())
:m_eo(1), m_steps_init(1), beta(), phi(), m_ns(0), m_time_storage(),
m_algebra(algebra),
m_phi_resizer()
{
for (size_t i=0; i<order_value+2; ++i)
{
c[i] = 1.0/(i+1);
c[c_size+i] = 1.0/((i+1)*(i+2));
}
g[0] = c[0];
g[1] = c[c_size];
beta[0][0] = 1;
beta[1][0] = 1;
gs[0] = 1.0;
gs[1] = -1.0/2;
gs[2] = -1.0/12;
gs[3] = -1.0/24;
gs[4] = -19.0/720;
gs[5] = -3.0/160;
gs[6] = -863.0/60480;
gs[7] = -275.0/24192;
gs[8] = -33953.0/3628800;
gs[9] = 35.0/4436;
gs[10] = 40.0/5891;
gs[11] = 37.0/6250;
gs[12] = 25.0/4771;
gs[13] = 40.0/8547;
};
void predict(time_type t, time_type dt)
{
using std::abs;
m_time_storage[0] = t;
if (abs(m_time_storage[0] - m_time_storage[1] - dt) > 1e-16 || m_eo >= m_ns)
{
m_ns = 0;
}
else if (m_ns < order_value + 2)
{
m_ns++;
}
for(size_t i=1+m_ns; i<m_eo+1 && i<m_steps_init; ++i)
{
time_type diff = m_time_storage[0] - m_time_storage[i];
beta[0][i] = beta[0][i-1]*(m_time_storage[0] + dt - m_time_storage[i-1])/diff;
}
for(size_t i=2+m_ns; i<m_eo+2 && i<m_steps_init+1; ++i)
{
time_type diff = m_time_storage[0] + dt - m_time_storage[i-1];
for(size_t j=0; j<m_eo+1-i+1; ++j)
{
c[c_size*i+j] = c[c_size*(i-1)+j] - c[c_size*(i-1)+j+1]*dt/diff;
}
g[i] = c[c_size*i];
}
};
void do_step(const deriv_type &dxdt, const int o = 0)
{
m_phi_resizer.adjust_size( dxdt , detail::bind( &aac_type::template resize_phi_impl< deriv_type > , detail::ref( *this ) , detail::_1 ) );
phi[o][0].m_v = dxdt;
for(size_t i=1; i<m_eo+3 && i<m_steps_init+2 && i<order_value+2; ++i)
{
if (o == 0)
{
this->m_algebra.for_each3(phi[o][i].m_v, phi[o][i-1].m_v, phi[o+1][i-1].m_v,
typename Operations::template scale_sum2<value_type, value_type>(1.0, -beta[o][i-1]));
}
else
{
this->m_algebra.for_each2(phi[o][i].m_v, phi[o][i-1].m_v,
typename Operations::template scale_sum1<value_type>(1.0));
}
}
};
void confirm()
{
beta.rotate();
phi.rotate();
m_time_storage.rotate();
if(m_steps_init < order_value+1)
{
++m_steps_init;
}
};
void reset() { m_eo = 1; m_steps_init = 1; };
size_t m_eo;
size_t m_steps_init;
rotating_buffer<boost::array<value_type, order_value+1>, 2> beta; // beta[0] = beta(n)
rotating_buffer<boost::array<wrapped_deriv_type, order_value+2>, 3> phi; // phi[0] = phi(n+1)
boost::array<value_type, order_value + 2> g;
boost::array<value_type, 14> gs;
private:
template< class StateType >
bool resize_phi_impl( const StateType &x )
{
bool resized( false );
for(size_t i=0; i<(order_value + 2); ++i)
{
resized |= adjust_size_by_resizeability( phi[0][i], x, typename is_resizeable<deriv_type>::type() );
resized |= adjust_size_by_resizeability( phi[1][i], x, typename is_resizeable<deriv_type>::type() );
resized |= adjust_size_by_resizeability( phi[2][i], x, typename is_resizeable<deriv_type>::type() );
}
return resized;
};
size_t m_ns;
time_storage_type m_time_storage;
static const size_t c_size = order_value + 2;
boost::array<value_type, c_size*c_size> c;
algebra_type m_algebra;
resizer_type m_phi_resizer;
};
} // detail
} // odeint
} // numeric
} // boost
#endif

199
include/boost/numeric/odeint/stepper/detail/pid_step_adjuster.hpp Normal file
View File

@ -0,0 +1,199 @@
/*
boost/numeric/odeint/stepper/detail/pid_step_adjuster.hpp
[begin_description]
Implementation of the stepsize controller for the controlled adams bashforth moulton stepper.
[end_description]
Copyright 2017 Valentin Noah Hartmann
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)
*/
#ifndef BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_PID_STEP_ADJUSTER_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_PID_STEP_ADJUSTER_HPP_INCLUDED
#include <boost/numeric/odeint/stepper/detail/rotating_buffer.hpp>
#include <boost/numeric/odeint/stepper/detail/pid_step_adjuster_coefficients.hpp>
#include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp>
#include <boost/numeric/odeint/algebra/operations_dispatcher.hpp>
#include <math.h>
namespace boost {
namespace numeric {
namespace odeint {
namespace detail {
template<
class Value = double,
class Time = double
>
struct pid_op
{
public:
typedef Value value_type;
typedef Time time_type;
const double beta1;
const double beta2;
const double beta3;
const double alpha1;
const double alpha2;
const time_type dt1;
const time_type dt2;
const time_type dt3;
const size_t m_steps;
pid_op(const size_t steps, const double _dt1, const double _dt2, const double _dt3,
const double b1 = 1, const double b2 = 0, const double b3 = 0, const double a1 = 0, const double a2 = 0)
:beta1(b1), beta2(b2), beta3(b3), alpha1(a1), alpha2(a2),
dt1(_dt1), dt2(_dt2), dt3(_dt3),
m_steps(steps)
{};
template<class T1, class T2, class T3, class T4>
void operator()(T1 &t1, const T2 &t2, const T3 &t3, const T4 &t4)
{
using std::abs;
t1 = adapted_pow(abs(t2), -beta1/(m_steps + 1)) *
adapted_pow(abs(t3), -beta2/(m_steps + 1)) *
adapted_pow(abs(t4), -beta3/(m_steps + 1)) *
adapted_pow(abs(dt1/dt2), -alpha1/(m_steps + 1))*
adapted_pow(abs(dt2/dt3), -alpha2/(m_steps + 1));
t1 = 1/t1;
};
template<class T1, class T2>
void operator()(T1 &t1, const T2 &t2)
{
using std::abs;
t1 = adapted_pow(abs(t2), -beta1/(m_steps + 1));
t1 = 1/t1;
};
private:
template<class T>
inline value_type adapted_pow(T base, double exp)
{
if(exp == 0)
{
return 1;
}
else if (exp > 0)
{
return pow(base, exp);
}
else
{
return 1/pow(base, -exp);
}
};
};
template<
class State,
class Value = double,
class Deriv = State,
class Time = double,
class Algebra = typename algebra_dispatcher< State >::algebra_type,
class Operations = typename operations_dispatcher< Deriv >::operations_type,
size_t Type = BASIC
>
struct pid_step_adjuster
{
public:
static double threshold() { return 0.9; };
typedef State state_type;
typedef Value value_type;
typedef Deriv deriv_type;
typedef Time time_type;
typedef Algebra algebra_type;
typedef Operations operations_type;
typedef rotating_buffer<state_type, 3> error_storage_type;
typedef rotating_buffer<time_type, 3> time_storage_type;
typedef pid_step_adjuster_coefficients<Type> coeff_type;
pid_step_adjuster(double abs_tol = 1e-6, double rel_tol = 1e-6, time_type dtmax = 1.0)
:m_dtmax(dtmax), m_error_storage(), m_dt_storage(), m_init(0),
m_abs_tol(abs_tol), m_rel_tol(rel_tol)
{};
time_type adjust_stepsize(const size_t steps, time_type dt, state_type &err, const state_type &x, const deriv_type &dxdt)
{
using std::abs;
m_algebra.for_each3( err , x , dxdt ,
typename operations_type::template rel_error< value_type >( m_abs_tol , m_rel_tol , 1.0 , 1.0 * abs(get_unit_value( dt )) ) );
m_error_storage[0] = err;
m_dt_storage[0] = dt;
if(m_init >= 2)
{
m_algebra.for_each4(err, m_error_storage[0], m_error_storage[1], m_error_storage[2],
pid_op<>(steps, m_dt_storage[0], m_dt_storage[1], m_dt_storage[2],
m_coeff[0], m_coeff[1], m_coeff[2], m_coeff[3], m_coeff[4]));
}
else
{
m_algebra.for_each2(err, m_error_storage[0],
pid_op<>(steps, m_dt_storage[0], m_dt_storage[1], m_dt_storage[2], 0.7));
}
value_type ratio = 1 / m_algebra.norm_inf(err);
value_type kappa = 1.0;
ratio = 1.0 + kappa*atan((ratio - 1) / kappa);
if(ratio*dt >= m_dtmax)
{
ratio = m_dtmax / dt;
}
if(ratio >= threshold())
{
m_error_storage.rotate();
m_dt_storage.rotate();
++m_init;
}
else
{
m_init = 0;
}
return dt * static_cast<time_type>(ratio);
};
private:
algebra_type m_algebra;
time_type m_dtmax;
error_storage_type m_error_storage;
time_storage_type m_dt_storage;
size_t m_init;
double m_abs_tol;
double m_rel_tol;
coeff_type m_coeff;
};
} // detail
} // odeint
} // numeric
} // boost
#endif

180
include/boost/numeric/odeint/stepper/detail/pid_step_adjuster_coefficients.hpp Normal file
View File

@ -0,0 +1,180 @@
/*
boost/numeric/odeint/stepper/detail/pid_step_adjuster_coefficients.hpp
[begin_description]
Coefficients for the PID stepsize controller.
[end_description]
Copyright 2017 Valentin Noah Hartmann
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)
*/
#ifndef BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_PID_STEP_ADJUSTER_COEFFICIENTS_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_PID_STEP_ADJUSTER_COEFFICIENTS_HPP_INCLUDED
#include <boost/array.hpp>
namespace boost {
namespace numeric {
namespace odeint {
namespace detail {
enum adjuster_type{
BASIC,
H0211,
H211b,
H211PI,
H0312,
H312b,
H312PID,
H0321,
H321
};
template<int Type>
class pid_step_adjuster_coefficients;
template<>
class pid_step_adjuster_coefficients<BASIC> : public boost::array<double, 5>
{
public:
pid_step_adjuster_coefficients()
: boost::array<double, 5>()
{
(*this)[0] = 1.0;
(*this)[1] = 0.0;
(*this)[2] = 0.0;
(*this)[3] = 0.0;
(*this)[4] = 0.0;
}
};
template<>
class pid_step_adjuster_coefficients<H0211> : public boost::array<double, 5>
{
public:
pid_step_adjuster_coefficients()
: boost::array<double, 5>()
{
(*this)[0] = 1.0 / 2.0;
(*this)[1] = 1.0 / 2.0;
(*this)[2] = 0.0;
(*this)[3] = 1.0 / 2.0;
(*this)[4] = 0.0;
}
};
template<>
class pid_step_adjuster_coefficients<H211b> : public boost::array<double, 5>
{
public:
pid_step_adjuster_coefficients()
: boost::array<double, 5>()
{
(*this)[0] = 1.0 / 5.0;
(*this)[1] = 2.0 / 5.0;
(*this)[2] = 0.0;
(*this)[3] = 1.0 / 5.0;
(*this)[4] = 0.0;
}
};
template<>
class pid_step_adjuster_coefficients<H211PI> : public boost::array<double, 5>
{
public:
pid_step_adjuster_coefficients()
: boost::array<double, 5>()
{
(*this)[0] = 1.0 / 6.0;
(*this)[1] = 2.0 / 6.0;
(*this)[2] = 0.0;
(*this)[3] = 0.0;
(*this)[4] = 0.0;
}
};
template<>
class pid_step_adjuster_coefficients<H0312> : public boost::array<double, 5>
{
public:
pid_step_adjuster_coefficients()
: boost::array<double, 5>()
{
(*this)[0] = 1.0 / 4.0;
(*this)[1] = 2.0 / 2.0;
(*this)[2] = 1.0 / 4.0;
(*this)[3] = 3.0 / 4.0;
(*this)[4] = 1.0 / 4.0;
}
};
template<>
class pid_step_adjuster_coefficients<H312b> : public boost::array<double, 5>
{
public:
pid_step_adjuster_coefficients()
: boost::array<double, 5>()
{
(*this)[0] = 1.0 / 6.0;
(*this)[1] = 2.0 / 6.0;
(*this)[2] = 1.0 / 6.0;
(*this)[3] = 3.0 / 6.0;
(*this)[4] = 1.0 / 6.0;
}
};
template<>
class pid_step_adjuster_coefficients<H312PID> : public boost::array<double, 5>
{
public:
pid_step_adjuster_coefficients()
: boost::array<double, 5>()
{
(*this)[0] = 1.0 / 18.0;
(*this)[1] = 2.0 / 9.0;
(*this)[2] = 1.0 / 18.0;
(*this)[3] = 0.0;
(*this)[4] = 0.0;
}
};
template<>
class pid_step_adjuster_coefficients<H0321> : public boost::array<double, 5>
{
public:
pid_step_adjuster_coefficients()
: boost::array<double, 5>()
{
(*this)[0] = 5.0 / 4.0;
(*this)[1] = 1.0 / 2.0;
(*this)[2] = -3.0 / 4.0;
(*this)[3] = -1.0 / 4.0;
(*this)[4] = -3.0 / 4.0;
}
};
template<>
class pid_step_adjuster_coefficients<H321> : public boost::array<double, 5>
{
public:
pid_step_adjuster_coefficients()
: boost::array<double, 5>()
{
(*this)[0] = 1.0 / 3.0;
(*this)[1] = 1.0 / 18.0;
(*this)[2] = -5.0 / 18.0;
(*this)[3] = -5.0 / 16.0;
(*this)[4] = -1.0 / 6.0;
}
};
} // detail
} // odeint
} // numeric
} // boost
#endif

1
include/boost/numeric/odeint/stepper/generation.hpp
View File

@ -29,6 +29,7 @@
#include <boost/numeric/odeint/stepper/generation/generation_runge_kutta_dopri5.hpp>
#include <boost/numeric/odeint/stepper/generation/generation_runge_kutta_fehlberg78.hpp>
#include <boost/numeric/odeint/stepper/generation/generation_controlled_adams_bashforth_moulton.hpp>
#include <boost/numeric/odeint/stepper/generation/generation_rosenbrock4.hpp>

59
include/boost/numeric/odeint/stepper/generation/generation_controlled_adams_bashforth_moulton.hpp Normal file
View File

@ -0,0 +1,59 @@
/*
boost/numeric/odeint/stepper/detail/generation_controlled_adams_bashforth_moulton.hpp
[begin_description]
Spezialization of the generation functions for creation of the controlled adams bashforth moulton stepper.
[end_description]
Copyright 2017 Valentin Noah Hartmann
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)
*/
#ifndef GENERATION_CONTROLLED_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED
#define GENERATION_CONTROLLED_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED
#include <boost/numeric/odeint/stepper/adaptive_adams_bashforth_moulton.hpp>
#include <boost/numeric/odeint/stepper/controlled_adams_bashforth_moulton.hpp>
#include <boost/numeric/odeint/stepper/generation/make_controlled.hpp>
namespace boost {
namespace numeric {
namespace odeint {
template< size_t Steps, class State , class Value , class Deriv , class Time , class Algebra , class Operations , class Resizer >
struct get_controller< adaptive_adams_bashforth_moulton< Steps, State , Value , Deriv , Time , Algebra , Operations , Resizer > >
{
typedef adaptive_adams_bashforth_moulton<Steps, State, Value, Deriv, Time, Algebra, Operations, Resizer> stepper_type;
typedef controlled_adams_bashforth_moulton< stepper_type > type;
};
// controller factory for controlled_adams_bashforth_moulton
template< class Stepper >
struct controller_factory< Stepper , controlled_adams_bashforth_moulton< Stepper > >
{
typedef Stepper stepper_type;
typedef controlled_adams_bashforth_moulton< stepper_type > controller_type;
typedef typename controller_type::step_adjuster_type step_adjuster_type;
typedef typename stepper_type::value_type value_type;
typedef typename stepper_type::value_type time_type;
controller_type operator()( value_type abs_error , value_type rel_error , const stepper_type &stepper )
{
return controller_type(step_adjuster_type(abs_error, rel_error));
}
controller_type operator()( value_type abs_error , value_type rel_error ,
time_type max_dt, const stepper_type &stepper )
{
return controller_type( step_adjuster_type(abs_error, rel_error, max_dt));
}
};
}
}
}
#endif

4
include/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp
View File

@ -163,6 +163,8 @@ public:
if( m_max_dt != static_cast<time_type>(0) )
{
dt = detail::min_abs(m_max_dt, dt_new);
} else {
dt = dt_new;
}
m_last_rejected = false;
return success;
@ -198,7 +200,7 @@ public:
private:
protected:
template< class StateIn >
bool resize_m_xerr( const StateIn &x )

19
include/boost/numeric/odeint/stepper/rosenbrock4_dense_output.hpp
View File

@ -41,7 +41,8 @@ class rosenbrock4_dense_output
public:
typedef ControlledStepper controlled_stepper_type;
typedef typename controlled_stepper_type::stepper_type stepper_type;
typedef typename unwrap_reference< controlled_stepper_type >::type unwrapped_controlled_stepper_type;
typedef typename unwrapped_controlled_stepper_type::stepper_type stepper_type;
typedef typename stepper_type::value_type value_type;
typedef typename stepper_type::state_type state_type;
typedef typename stepper_type::wrapped_state_type wrapped_state_type;
@ -54,7 +55,7 @@ public:
typedef rosenbrock4_dense_output< ControlledStepper > dense_output_stepper_type;
rosenbrock4_dense_output( const controlled_stepper_type &stepper = controlled_stepper_type() )
: m_stepper( stepper ) ,
: m_stepper( stepper ) ,
m_x1() , m_x2() ,
m_current_state_x1( true ) ,
m_t() , m_t_old() , m_dt()
@ -75,16 +76,17 @@ public:
template< class System >
std::pair< time_type , time_type > do_step( System system )
{
unwrapped_controlled_stepper_type &stepper = m_stepper;
failed_step_checker fail_checker; // to throw a runtime_error if step size adjustment fails
controlled_step_result res = fail;
m_t_old = m_t;
do
{
res = m_stepper.try_step( system , get_current_state() , m_t , get_old_state() , m_dt );
res = stepper.try_step( system , get_current_state() , m_t , get_old_state() , m_dt );
fail_checker(); // check for overflow of failed steps
}
while( res == fail );
m_stepper.stepper().prepare_dense_output();
stepper.stepper().prepare_dense_output();
this->toggle_current_state();
return std::make_pair( m_t_old , m_t );
}
@ -96,20 +98,23 @@ public:
template< class StateOut >
void calc_state( time_type t , StateOut &x )
{
m_stepper.stepper().calc_state( t , x , get_old_state() , m_t_old , get_current_state() , m_t );
unwrapped_controlled_stepper_type &stepper = m_stepper;
stepper.stepper().calc_state( t , x , get_old_state() , m_t_old , get_current_state() , m_t );
}
template< class StateOut >
void calc_state( time_type t , const StateOut &x )
{
m_stepper.stepper().calc_state( t , x , get_old_state() , m_t_old , get_current_state() , m_t );
unwrapped_controlled_stepper_type &stepper = m_stepper;
stepper.stepper().calc_state( t , x , get_old_state() , m_t_old , get_current_state() , m_t );
}
template< class StateType >
void adjust_size( const StateType &x )
{
m_stepper.adjust_size( x );
unwrapped_controlled_stepper_type &stepper = m_stepper;
stepper.adjust_size( x );
resize_impl( x );
}

3
test/Jamfile.v2
View File

@ -48,6 +48,8 @@ test-suite "odeint"
[ run adams_bashforth.cpp ]
[ run adams_moulton.cpp ]
[ run adams_bashforth_moulton.cpp ]
[ run controlled_adams_bashforth_moulton.cpp ]
[ run adaptive_adams_coefficients.cpp ]
[ run generic_stepper.cpp ]
[ run generic_error_stepper.cpp ]
[ run bulirsch_stoer.cpp ]
@ -82,7 +84,6 @@ test-suite "odeint"
[ run integrate_overflow.cpp ]
[ compile unwrap_boost_reference.cpp ]
[ compile unwrap_reference.cpp : <cxxflags>-std=c++0x : unwrap_reference_C++11 ]
[ compile-fail unwrap_reference.cpp : <cxxflags>-std=c++98 : unwrap_reference_C++98 ]
[ compile std_array.cpp : <cxxflags>-std=c++0x ]
:
<testing.launcher>valgrind

7
test/adams_bashforth.cpp
View File

@ -208,8 +208,13 @@ BOOST_AUTO_TEST_CASE( test_auto_initialization )
adams.do_step( lorenz() , x , 0.0 , x , 0.1 );
BOOST_CHECK_EQUAL( adams.initializing_stepper().do_count , size_t( 2 ) );
adams.reset();
adams.do_step( lorenz() , x , 0.0 , x , 0.1 );
BOOST_CHECK_EQUAL( adams.initializing_stepper().do_count , size_t( 2 ) );
BOOST_CHECK_EQUAL( adams.initializing_stepper().do_count , size_t( 3 ) );
adams.do_step( lorenz() , x , 0.0 , x , 0.1 );
BOOST_CHECK_EQUAL( adams.initializing_stepper().do_count , size_t( 4 ) );
}
BOOST_AUTO_TEST_CASE( test_manual_initialization )

4
test/adams_bashforth_moulton.cpp
View File

@ -104,8 +104,8 @@ BOOST_AUTO_TEST_CASE( test_instantiation )
s4.do_step( lorenz() , x , t , dt );
s5.do_step( lorenz() , x , t , dt );
s6.do_step( lorenz() , x , t , dt );
// s7.do_step( lorenz() , x , t , dt );
// s8.do_step( lorenz() , x , t , dt );
// s7.do_step( lorenz() , x , t , dt );
// s8.do_step( lorenz() , x , t , dt );
}

120
test/adaptive_adams_coefficients.cpp Normal file
View File

@ -0,0 +1,120 @@
#include <boost/config.hpp>
#ifdef BOOST_MSVC
#pragma warning(disable:4996)
#endif
#define BOOST_TEST_MODULE odeint_adaptive_adams_coefficients
#include <boost/test/unit_test.hpp>
#include <boost/numeric/odeint/stepper/detail/adaptive_adams_coefficients.hpp>
#include <vector>
#include <boost/mpl/list.hpp>
#include <boost/mpl/size_t.hpp>
#include <boost/mpl/range_c.hpp>
using namespace boost::unit_test;
using namespace boost::numeric::odeint;
typedef double value_type;
BOOST_AUTO_TEST_SUITE( adaptive_adams_coefficients_test )
typedef boost::mpl::range_c< size_t , 2 , 10 > vector_of_steps;
BOOST_AUTO_TEST_CASE_TEMPLATE( test_step, step_type, vector_of_steps )
{
const static size_t steps = step_type::value;
typedef std::vector<double> deriv_type;
typedef double time_type;
typedef detail::adaptive_adams_coefficients<steps, deriv_type, time_type> aac_type;
std::vector<double> deriv;
deriv.push_back(-1);
time_type t = 0.0;
time_type dt = 0.1;
aac_type coeff;
for(size_t i=0; i<steps; ++i)
{
coeff.predict(t, dt);
coeff.do_step(deriv);
coeff.confirm();
t+= dt;
if(coeff.m_eo < steps)
coeff.m_eo ++;
}
std::vector<value_type> v(10);
v[0] = 1.0/1.0;
v[1] = 1.0/2.0;
v[2] = 5.0/12.0;
v[3] = 9.0/24.0;
v[4] = 251.0/720.0;
v[5] = 95.0/288.0;
v[6] = 19087.0/60480.0;
v[7] = 5257.0/17280.0;
v[8] = 5311869667636789.0/18014398509481984.0;
for(size_t i=0; i<steps; ++i)
{
BOOST_CHECK_SMALL(coeff.beta[1][i] - 1.0, 1e-15);
if(i == 0)
BOOST_CHECK_SMALL(coeff.phi[2][i].m_v[0] + 1, 1e-15);
else if (i == steps-1 && steps%2 == 1)
BOOST_CHECK_SMALL(coeff.phi[2][i].m_v[0] - 1, 1e-14);
else if (i == steps-1 && steps%2 == 0)
BOOST_CHECK_SMALL(coeff.phi[2][i].m_v[0] + 1, 1e-14);
else
BOOST_CHECK_SMALL(coeff.phi[2][i].m_v[0], 1e-15);
BOOST_CHECK_SMALL(coeff.g[i] - v[i], 1e-15);
}
}
BOOST_AUTO_TEST_CASE( test_copy )
{
typedef std::vector<double> deriv_type;
typedef double time_type;
typedef detail::adaptive_adams_coefficients<3, deriv_type, time_type> aac_type;
aac_type c1;
deriv_type deriv(1);
deriv[0] = 1.0;
time_type t = 0.0;
time_type dt = 0.01;
for(size_t i=0; i<3; ++i)
{
c1.predict(t, dt);
c1.do_step(deriv);
c1.confirm();
t+= dt;
if(c1.m_eo < 3)
c1.m_eo ++;
}
aac_type c2(c1);
BOOST_CHECK_EQUAL(c1.phi[0][0].m_v[0], c2.phi[0][0].m_v[0]);
BOOST_CHECK(&(c1.phi[0][0].m_v) != &(c2.phi[0][0].m_v));
aac_type c3;
deriv_type *p1 = &(c3.phi[0][0].m_v);
c3 = c1;
BOOST_CHECK(p1 == (&(c3.phi[0][0].m_v)));
BOOST_CHECK_EQUAL(c1.phi[0][0].m_v[0], c3.phi[0][0].m_v[0]);
}
BOOST_AUTO_TEST_SUITE_END()

56
test/controlled_adams_bashforth_moulton.cpp Normal file
View File

@ -0,0 +1,56 @@
#include <boost/config.hpp>
#ifdef BOOST_MSVC
#pragma warning(disable:4996)
#endif
#define BOOST_TEST_MODULE odeint_controlled_adams_bashforth_moulton
#include <boost/test/unit_test.hpp>
#include <boost/numeric/odeint/stepper/adaptive_adams_bashforth_moulton.hpp>
#include <boost/numeric/odeint/stepper/controlled_adams_bashforth_moulton.hpp>
using namespace boost::unit_test;
using namespace boost::numeric::odeint;
struct const_sys
{
template< class State , class Deriv , class Value >
void operator()( const State &x , Deriv &dxdt , const Value &dt ) const
{
dxdt[0] = 1;
}
};
typedef boost::array< double , 1 > state_type;
typedef double value_type;
BOOST_AUTO_TEST_SUITE( controlled_adams_bashforth_moulton_test )
BOOST_AUTO_TEST_CASE( test_instantiation )
{
controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<1, state_type> > s1;
controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<2, state_type> > s2;
controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<3, state_type> > s3;
controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<4, state_type> > s4;
controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<5, state_type> > s5;
controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<6, state_type> > s6;
controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<7, state_type> > s7;
controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<8, state_type> > s8;
controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<9, state_type> > s9;
state_type x = {{ 10.0 }};
value_type t = 0.0 , dt = 0.01;
s1.try_step(const_sys(), x, t, dt);
s2.try_step(const_sys(), x, t, dt);
s3.try_step(const_sys(), x, t, dt);
s4.try_step(const_sys(), x, t, dt);
s5.try_step(const_sys(), x, t, dt);
s6.try_step(const_sys(), x, t, dt);
s7.try_step(const_sys(), x, t, dt);
s8.try_step(const_sys(), x, t, dt);
s9.try_step(const_sys(), x, t, dt);
}
BOOST_AUTO_TEST_SUITE_END()

28
test/integrate.cpp
View File

@ -17,6 +17,8 @@
#define BOOST_TEST_MODULE odeint_integrate_functions
#include <boost/type_traits.hpp>
#include <vector>
#include <cmath>
#include <iostream>
@ -54,6 +56,9 @@
#include <boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp>
#include <boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp>
#include <boost/numeric/odeint/stepper/adaptive_adams_bashforth_moulton.hpp>
#include <boost/numeric/odeint/stepper/controlled_adams_bashforth_moulton.hpp>
#include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
using namespace boost::unit_test;
@ -237,7 +242,12 @@ class stepper_methods : public mpl::vector<
controlled_runge_kutta< runge_kutta_dopri5< state_type > > ,
controlled_runge_kutta< runge_kutta_fehlberg78< state_type > > ,
bulirsch_stoer< state_type > ,
dense_output_runge_kutta< controlled_runge_kutta< runge_kutta_dopri5< state_type > > >
dense_output_runge_kutta< controlled_runge_kutta< runge_kutta_dopri5< state_type > > >,
adaptive_adams_bashforth_moulton<3, state_type>,
adaptive_adams_bashforth_moulton<5, state_type>,
adaptive_adams_bashforth_moulton<7, state_type>,
controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<3, state_type> >,
controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<5, state_type> >
//bulirsch_stoer_dense_out< state_type >
> { };
@ -274,12 +284,16 @@ BOOST_AUTO_TEST_CASE_TEMPLATE( integrate_times_test_case , Stepper, stepper_meth
BOOST_AUTO_TEST_CASE_TEMPLATE( integrate_n_steps_test_case , Stepper, stepper_methods )
{
perform_integrate_n_steps_test< Stepper > tester;
tester();
tester( 200 , 0.01 );
tester( 200 , 0.01 );
tester( 200 , 0.01 );
tester( 200 , -0.01 );
if(!boost::is_same<Stepper, controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<3, state_type> > >::value &&
!boost::is_same<Stepper, controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<5, state_type> > >::value)
{
perform_integrate_n_steps_test< Stepper > tester;
tester();
tester( 200 , 0.01 );
tester( 200 , 0.01 );
tester( 200 , 0.01 );
tester( 200 , -0.01 );
}
}
BOOST_AUTO_TEST_SUITE_END()

1
test/numeric/Jamfile.v2
View File

@ -28,6 +28,7 @@ test-suite "odeint"
[ run rosenbrock.cpp ]
[ run adams_bashforth.cpp ]
[ run adams_bashforth_moulton.cpp ]
[ run adaptive_adams_bashforth_moulton.cpp ]
[ run abm_time_dependent.cpp ]
[ run order_quadrature_formula.cpp ]
[ run velocity_verlet.cpp ]

113
test/numeric/adaptive_adams_bashforth_moulton.cpp Normal file
View File

@ -0,0 +1,113 @@
/* Boost numeric test of the adams-bashforth steppers test file
Copyright 2013 Karsten Ahnert
Copyright 2013-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)
*/
// disable checked iterator warning for msvc
#include <boost/config.hpp>
#ifdef BOOST_MSVC
#pragma warning(disable:4996)
#endif
#define BOOST_TEST_MODULE numeric_adaptive_adams_bashforth_moulton
#include <iostream>
#include <cmath>
#include <boost/array.hpp>
#include <boost/test/unit_test.hpp>
#include <boost/mpl/vector.hpp>
#include <boost/numeric/odeint.hpp>
using namespace boost::unit_test;
using namespace boost::numeric::odeint;
namespace mpl = boost::mpl;
typedef double value_type;
typedef boost::array< double , 2 > state_type;
typedef runge_kutta_fehlberg78<state_type> initializing_stepper;
// harmonic oscillator, analytic solution x[0] = sin( t )
struct osc
{
void operator()( const state_type &x , state_type &dxdt , const double t ) const
{
dxdt[0] = x[1];
dxdt[1] = -x[0];
}
};
BOOST_AUTO_TEST_SUITE( numeric_adaptive_adams_bashforth_moulton_test )
/* generic test for all adams bashforth steppers */
template< class Stepper >
struct perform_adaptive_adams_bashforth_moulton_test
{
void operator()( void )
{
Stepper stepper;
initializing_stepper init_stepper;
const int o = stepper.order()+1; //order of the error is order of approximation + 1
const state_type x0 = {{ 0.0 , 1.0 }};
state_type x1 = x0;
double t = 0.0;
double dt = 0.25;
// initialization, does a number of steps to self-start the stepper with a small stepsize
stepper.initialize( init_stepper, osc() , x1 , t , dt);
double A = std::sqrt( x1[0]*x1[0] + x1[1]*x1[1] );
double phi = std::asin(x1[0]/A) - t;
// now we do the actual step
stepper.do_step( osc() , x1 , t , dt );
// only examine the error of the adams-bashforth step, not the initialization
const double f = 2.0 * std::abs( A*sin(t+dt+phi) - x1[0] ) / std::pow( dt , o ); // upper bound
std::cout << o << " , " << f << std::endl;
/* as long as we have errors above machine precision */
while( f*std::pow( dt , o ) > 1E-16 )
{
x1 = x0;
t = 0.0;
stepper.initialize( init_stepper, osc() , x1 , t , dt );
A = std::sqrt( x1[0]*x1[0] + x1[1]*x1[1] );
phi = std::asin(x1[0]/A) - t;
// now we do the actual step
stepper.do_step( osc() , x1 , t , dt );
stepper.reset();
// only examine the error of the adams-bashforth step, not the initialization
std::cout << "Testing dt=" << dt << " , " << std::abs( A*sin(t+dt+phi) - x1[0] ) << std::endl;
BOOST_CHECK_LT( std::abs( A*sin(t+dt+phi) - x1[0] ) , f*std::pow( dt , o ) );
dt *= 0.5;
}
}
};
typedef mpl::vector<
adaptive_adams_bashforth_moulton< 2 , state_type > ,
adaptive_adams_bashforth_moulton< 3 , state_type > ,
adaptive_adams_bashforth_moulton< 4 , state_type > ,
adaptive_adams_bashforth_moulton< 5 , state_type > ,
adaptive_adams_bashforth_moulton< 6 , state_type > ,
adaptive_adams_bashforth_moulton< 7 , state_type >
> adaptive_adams_bashforth_moulton_steppers;
BOOST_AUTO_TEST_CASE_TEMPLATE( adaptive_adams_bashforth_moulton_test , Stepper, adaptive_adams_bashforth_moulton_steppers )
{
perform_adaptive_adams_bashforth_moulton_test< Stepper > tester;
tester();
}
BOOST_AUTO_TEST_SUITE_END()

2
test/numeric/order_quadrature_formula.cpp
View File

@ -20,6 +20,8 @@
#include <iostream>
#include <cmath>
#include "boost/format.hpp"
#include <boost/test/unit_test.hpp>
#include <boost/mpl/vector.hpp>

1
test/regression/Jamfile.v2
View File

@ -26,5 +26,6 @@ test-suite "odeint"
[ run regression_147.cpp ]
[ compile regression_149.cpp : <cxxflags>-std=c++0x ]
[ run regression_168.cpp ]
[ run regression_189.cpp ]
: <testing.launcher>valgrind
;

72
test/regression/regression_189.cpp Normal file
View File

@ -0,0 +1,72 @@
/*
[begin_description]
Test case for issue 189:
Controlled Rosenbrock stepper fails to increase step size
[end_description]
Copyright 2016 Karsten Ahnert
Copyright 2016 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_regression_189
#include <boost/numeric/odeint.hpp>
#include <boost/test/unit_test.hpp>
#include <boost/phoenix/core.hpp>
#include <boost/phoenix/operator.hpp>
using namespace boost::numeric::odeint;
namespace phoenix = boost::phoenix;
typedef boost::numeric::ublas::vector< double > vector_type;
typedef boost::numeric::ublas::matrix< double > matrix_type;
struct stiff_system
{
void operator()( const vector_type &x , vector_type &dxdt , double /* t */ )
{
dxdt[ 0 ] = -101.0 * x[ 0 ] - 100.0 * x[ 1 ];
dxdt[ 1 ] = x[ 0 ];
}
};
struct stiff_system_jacobi
{
void operator()( const vector_type & /* x */ , matrix_type &J , const double & /* t */ , vector_type &dfdt )
{
J( 0 , 0 ) = -101.0;
J( 0 , 1 ) = -100.0;
J( 1 , 0 ) = 1.0;
J( 1 , 1 ) = 0.0;
dfdt[0] = 0.0;
dfdt[1] = 0.0;
}
};
BOOST_AUTO_TEST_CASE( regression_189 )
{
vector_type x( 2 , 1.0 );
size_t num_of_steps = integrate_const( make_dense_output< rosenbrock4< double > >( 1.0e-6 , 1.0e-6 ) ,
std::make_pair( stiff_system() , stiff_system_jacobi() ) ,
x , 0.0 , 50.0 , 0.01 ,
std::cout << phoenix::arg_names::arg2 << " " << phoenix::arg_names::arg1[0] << "\n" );
// regression: number of steps should be 74
size_t num_of_steps_expected = 74;
BOOST_CHECK_EQUAL( num_of_steps , num_of_steps_expected );
vector_type x2( 2 , 1.0 );
size_t num_of_steps2 = integrate_const( make_dense_output< runge_kutta_dopri5< vector_type > >( 1.0e-6 , 1.0e-6 ) ,
stiff_system() , x2 , 0.0 , 50.0 , 0.01 ,
std::cout << phoenix::arg_names::arg2 << " " << phoenix::arg_names::arg1[0] << "\n" );
num_of_steps_expected = 1531;
BOOST_CHECK_EQUAL( num_of_steps2 , num_of_steps_expected );
}

33
test/rosenbrock4.cpp
View File

@ -127,6 +127,39 @@ BOOST_AUTO_TEST_CASE( test_rosenbrock4_dense_output )
stepper.calc_state( 0.5 * ( tr.first + tr.second ) , x );
}
class rosenbrock4_controller_max_dt_adaptable : public rosenbrock4_controller< rosenbrock4< value_type > >
{
public:
void set_max_dt(value_type max_dt)
{
m_max_dt = max_dt;
}
};
BOOST_AUTO_TEST_CASE( test_rosenbrock4_dense_output_ref )
{
typedef rosenbrock4_dense_output< boost::reference_wrapper< rosenbrock4_controller_max_dt_adaptable > > stepper_type;
rosenbrock4_controller_max_dt_adaptable c_stepper;
stepper_type stepper( boost::ref( c_stepper ) );
typedef stepper_type::state_type state_type;
typedef stepper_type::value_type stepper_value_type;
typedef stepper_type::deriv_type deriv_type;
typedef stepper_type::time_type time_type;
state_type x( 2 );
x( 0 ) = 0.0 ; x(1) = 1.0;
stepper.initialize( x , 0.0 , 0.1 );
std::pair< value_type , value_type > tr = stepper.do_step( std::make_pair( sys() , jacobi() ) );
stepper.calc_state( 0.5 * ( tr.first + tr.second ) , x );
// adapt the maximal step size to a small value (smaller than the step size) and make a step
const double max_dt = 1e-8;
c_stepper.set_max_dt(max_dt);
stepper.do_step( std::make_pair( sys() , jacobi() ) );
// check if the step was done with the requested maximal step size
BOOST_CHECK_CLOSE(max_dt, stepper.current_time_step(), 1e-14);
}
BOOST_AUTO_TEST_CASE( test_rosenbrock4_copy_dense_output )
{
typedef rosenbrock4_controller< rosenbrock4< value_type > > controlled_stepper_type;

1
test/runge_kutta_concepts.cpp
View File

@ -88,6 +88,7 @@ struct perform_stepper_test
typedef T vector_space_type;
void operator()( void ) const
{
using std::abs;
vector_space_type x;
x = 2.0;
Stepper stepper;

2
test/runge_kutta_controlled_concepts.cpp
View File

@ -87,6 +87,7 @@ struct perform_controlled_stepper_test
typedef T vector_space_type;
void operator()( void ) const
{
using std::abs;
vector_space_type x;
x = 2.0;
ControlledStepper controlled_stepper;
@ -134,6 +135,7 @@ struct perform_controlled_stepper_test< ControlledStepper , vector_space_type >
{
void operator()( void ) const
{
using std::abs;
vector_space_type x;
x = 2.0;
ControlledStepper controlled_stepper;

1
test_external/eigen/Jamfile.v2
View File

@ -31,6 +31,5 @@ test-suite "odeint"
[ run runge_kutta4.cpp ]
[ run runge_kutta_dopri5.cpp ]
[ run integrate.cpp ]
[ compile-fail fail_integrate.cpp ]
: <testing.launcher>valgrind
;

47
test_external/eigen/fail_integrate.cpp
View File

@ -1,47 +0,0 @@
/*
[auto_generated]
fail_integrate.cpp
[begin_description]
tba.
[end_description]
Copyright 2009-2012 Karsten Ahnert
Copyright 2009-2012 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)
*/
#include <boost/config.hpp>
#ifdef BOOST_MSVC
#pragma warning(disable:4996)
#endif
#define BOOST_TEST_MODULE odeint_dummy
#include <boost/numeric/odeint/integrate/integrate.hpp>
#include <boost/numeric/odeint/external/eigen/eigen.hpp>
#include <boost/test/unit_test.hpp>
#include "dummy_odes.hpp"
using namespace boost::unit_test;
using namespace boost::numeric::odeint;
BOOST_AUTO_TEST_SUITE( eigen_fail_integrate )
BOOST_AUTO_TEST_CASE( test )
{
typedef Eigen::Matrix< double , 1 , 1 > state_type;
state_type x;
x[0] = 10.0;
double t_start = 0.0 , t_end = 1.0 , dt = 0.1;
integrate( constant_system_functor_standard() , x , t_start , t_end , dt );
}
BOOST_AUTO_TEST_SUITE_END()