summaryrefslogtreecommitdiff
path: root/boost/numeric/odeint/stepper
diff options
context:
space:
mode:
Diffstat (limited to 'boost/numeric/odeint/stepper')
-rw-r--r--boost/numeric/odeint/stepper/adams_bashforth.hpp28
-rw-r--r--boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp10
-rw-r--r--boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp35
-rw-r--r--boost/numeric/odeint/stepper/base/explicit_error_stepper_fsal_base.hpp30
-rw-r--r--boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp31
-rw-r--r--boost/numeric/odeint/stepper/base/symplectic_rkn_stepper_base.hpp2
-rw-r--r--boost/numeric/odeint/stepper/bulirsch_stoer.hpp7
-rw-r--r--boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp5
-rw-r--r--boost/numeric/odeint/stepper/detail/adams_bashforth_call_algebra.hpp4
-rw-r--r--boost/numeric/odeint/stepper/detail/adams_moulton_call_algebra.hpp2
-rw-r--r--boost/numeric/odeint/stepper/euler.hpp4
-rw-r--r--boost/numeric/odeint/stepper/extrapolation_stepper.hpp293
-rw-r--r--boost/numeric/odeint/stepper/implicit_euler.hpp2
-rw-r--r--boost/numeric/odeint/stepper/rosenbrock4_controller.hpp1
-rw-r--r--boost/numeric/odeint/stepper/runge_kutta4_classic.hpp11
-rw-r--r--boost/numeric/odeint/stepper/velocity_verlet.hpp2
16 files changed, 434 insertions, 33 deletions
diff --git a/boost/numeric/odeint/stepper/adams_bashforth.hpp b/boost/numeric/odeint/stepper/adams_bashforth.hpp
index 59ed1c9..5ff1e83 100644
--- a/boost/numeric/odeint/stepper/adams_bashforth.hpp
+++ b/boost/numeric/odeint/stepper/adams_bashforth.hpp
@@ -37,6 +37,7 @@
#include <boost/numeric/odeint/stepper/stepper_categories.hpp>
#include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
+#include <boost/numeric/odeint/stepper/extrapolation_stepper.hpp>
#include <boost/numeric/odeint/stepper/base/algebra_stepper_base.hpp>
@@ -44,11 +45,28 @@
#include <boost/numeric/odeint/stepper/detail/adams_bashforth_call_algebra.hpp>
#include <boost/numeric/odeint/stepper/detail/rotating_buffer.hpp>
+#include <boost/mpl/arithmetic.hpp>
+#include <boost/mpl/min_max.hpp>
+#include <boost/mpl/equal_to.hpp>
+
+namespace mpl = boost::mpl;
+
namespace boost {
namespace numeric {
namespace odeint {
+ using mpl::int_;
+
+ /* if N >= 4, returns the smallest even number > N, otherwise returns 4 */
+ template < int N >
+ struct order_helper
+ : mpl::max< typename mpl::eval_if<
+ mpl::equal_to< mpl::modulus< int_< N >, int_< 2 > >,
+ int_< 0 > >,
+ int_< N >, int_< N + 1 > >::type,
+ int_< 4 > >::type
+ { };
template<
size_t Steps ,
@@ -59,7 +77,9 @@ class Time = Value ,
class Algebra = typename algebra_dispatcher< State >::algebra_type ,
class Operations = typename operations_dispatcher< State >::operations_type ,
class Resizer = initially_resizer ,
-class InitializingStepper = runge_kutta4< State , Value , Deriv , Time , Algebra , Operations, Resizer >
+class InitializingStepper = extrapolation_stepper< order_helper<Steps>::value,
+ State, Value, Deriv, Time,
+ Algebra, Operations, Resizer >
>
class adams_bashforth : public algebra_stepper_base< Algebra , Operations >
{
@@ -181,7 +201,8 @@ public :
{
if( i != 0 ) m_step_storage.rotate();
sys( x , m_step_storage[0].m_v , t );
- stepper.do_step( system , x , m_step_storage[0].m_v , t , dt );
+ stepper.do_step_dxdt_impl( system, x, m_step_storage[0].m_v, t,
+ dt );
t += dt;
}
m_steps_initialized = steps;
@@ -222,7 +243,8 @@ private:
{
if( m_steps_initialized != 0 ) m_step_storage.rotate();
sys( in , m_step_storage[0].m_v , t );
- m_initializing_stepper.do_step( system , in , m_step_storage[0].m_v , t , out , dt );
+ m_initializing_stepper.do_step_dxdt_impl(
+ system, in, m_step_storage[0].m_v, t, out, dt );
++m_steps_initialized;
}
else
diff --git a/boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp b/boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp
index 19b0239..2f7cc4c 100644
--- a/boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp
+++ b/boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp
@@ -49,7 +49,8 @@ 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 Resizer = initially_resizer,
+class InitializingStepper = runge_kutta4< State , Value , Deriv , Time , Algebra , Operations, Resizer >
>
class adams_bashforth_moulton
{
@@ -71,12 +72,13 @@ public :
typedef Operations operations_type;
typedef Resizer resizer_type;
typedef stepper_tag stepper_category;
+ typedef InitializingStepper initializing_stepper_type;
static const size_t steps = Steps;
#ifndef DOXYGEN_SKIP
- typedef adams_bashforth< steps , state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > adams_bashforth_type;
+ typedef adams_bashforth< steps , state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type, initializing_stepper_type > adams_bashforth_type;
typedef adams_moulton< steps , state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > adams_moulton_type;
- typedef adams_bashforth_moulton< steps , state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > stepper_type;
+ typedef adams_bashforth_moulton< steps , state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type , initializing_stepper_type> stepper_type;
#endif //DOXYGEN_SKIP
typedef unsigned short order_type;
static const order_type order_value = steps;
@@ -158,7 +160,7 @@ private:
{
m_resizer.adjust_size( x , detail::bind( &stepper_type::template resize_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) );
m_adams_bashforth.do_step( system , x , t , m_x.m_v , dt );
- m_adams_moulton.do_step( system , x , m_x.m_v , t , x , dt , m_adams_bashforth.step_storage() );
+ m_adams_moulton.do_step( system , x , m_x.m_v , t+dt , x , dt , m_adams_bashforth.step_storage() );
}
else
{
diff --git a/boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp b/boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp
index 3ba11a4..08009dc 100644
--- a/boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp
+++ b/boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp
@@ -152,6 +152,22 @@ public:
/*
+ * named Version 2: do_step_dxdt_impl( sys , in , dxdt , t , dt )
+ *
+ * this version is needed when this stepper is used for initializing
+ * multistep stepper like adams-bashforth. Hence we provide an explicitely
+ * named version that is not disabled. Meant for internal use only.
+ */
+ template < class System, class StateInOut, class DerivIn >
+ void do_step_dxdt_impl( System system, StateInOut &x, const DerivIn &dxdt,
+ time_type t, time_type dt )
+ {
+ this->stepper().do_step_impl( system , x , dxdt , t , x , dt );
+ }
+
+
+
+ /*
* Version 3 : do_step( sys , in , t , out , dt )
*
* this version does not solve the forwarding problem, boost.range can not be used
@@ -181,10 +197,21 @@ public:
{
this->stepper().do_step_impl( system , in , dxdt , t , out , dt );
}
-
-
-
-
+
+ /*
+ * named Version 4: do_step_dxdt_impl( sys , in , dxdt , t , out, dt )
+ *
+ * this version is needed when this stepper is used for initializing
+ * multistep stepper like adams-bashforth. Hence we provide an explicitely
+ * named version that is not disabled. Meant for internal use only.
+ */
+ template < class System, class StateIn, class DerivIn, class StateOut >
+ void do_step_dxdt_impl( System system, const StateIn &in,
+ const DerivIn &dxdt, time_type t, StateOut &out,
+ time_type dt )
+ {
+ this->stepper().do_step_impl( system , in , dxdt , t , out , dt );
+ }
/*
* Version 5 :do_step( sys , x , t , dt , xerr )
diff --git a/boost/numeric/odeint/stepper/base/explicit_error_stepper_fsal_base.hpp b/boost/numeric/odeint/stepper/base/explicit_error_stepper_fsal_base.hpp
index a055c7f..b1d751a 100644
--- a/boost/numeric/odeint/stepper/base/explicit_error_stepper_fsal_base.hpp
+++ b/boost/numeric/odeint/stepper/base/explicit_error_stepper_fsal_base.hpp
@@ -151,11 +151,27 @@ public:
/*
+ * named Version 2: do_step_dxdt_impl( sys , in , dxdt , t , dt )
+ *
+ * this version is needed when this stepper is used for initializing
+ * multistep stepper like adams-bashforth. Hence we provide an explicitely
+ * named version that is not disabled. Meant for internal use only.
+ */
+ template< class System , class StateInOut , class DerivInOut >
+ void do_step_dxdt_impl( System system , StateInOut &x , DerivInOut &dxdt , time_type t , time_type dt )
+ {
+ m_first_call = true;
+ this->stepper().do_step_impl( system , x , dxdt , t , x , dxdt , dt );
+ }
+
+ /*
* version 3 : do_step( sys , in , t , out , dt )
*
- * this version does not solve the forwarding problem, boost.range can not be used
+ * this version does not solve the forwarding problem, boost.range can not
+ * be used.
*
- * the disable is needed to avoid ambiguous overloads if state_type = time_type
+ * the disable is needed to avoid ambiguous overloads if
+ * state_type = time_type
*/
template< class System , class StateIn , class StateOut >
typename boost::disable_if< boost::is_same< StateIn , time_type > , void >::type
@@ -174,12 +190,14 @@ public:
*
* this version does not solve the forwarding problem, boost.range can not be used
*/
- template< class System , class StateIn , class DerivIn , class StateOut , class DerivOut >
- void do_step( System system , const StateIn &in , const DerivIn &dxdt_in , time_type t ,
- StateOut &out , DerivOut &dxdt_out , time_type dt )
+ template< class System, class StateIn, class DerivIn, class StateOut,
+ class DerivOut >
+ void do_step( System system, const StateIn &in, const DerivIn &dxdt_in,
+ time_type t, StateOut &out, DerivOut &dxdt_out, time_type dt )
{
m_first_call = true;
- this->stepper().do_step_impl( system , in , dxdt_in , t , out , dxdt_out , dt );
+ this->stepper().do_step_impl( system, in, dxdt_in, t, out, dxdt_out,
+ dt );
}
diff --git a/boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp b/boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp
index 40aab80..d81c8c7 100644
--- a/boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp
+++ b/boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp
@@ -139,6 +139,21 @@ public:
/*
+ * named Version 2: do_step_dxdt_impl( sys , in , dxdt , t , dt )
+ *
+ * this version is needed when this stepper is used for initializing
+ * multistep stepper like adams-bashforth. Hence we provide an explicitely
+ * named version that is not disabled. Meant for internal use only.
+ */
+ template < class System, class StateInOut, class DerivIn >
+ void do_step_dxdt_impl( System system, StateInOut &x, const DerivIn &dxdt,
+ time_type t, time_type dt )
+ {
+ this->stepper().do_step_impl( system , x , dxdt , t , x , dt );
+ }
+
+
+ /*
* Version 3 : do_step( sys , in , t , out , dt )
*
* this version does not solve the forwarding problem, boost.range can not be used
@@ -164,6 +179,22 @@ public:
this->stepper().do_step_impl( system , in , dxdt , t , out , dt );
}
+
+ /*
+ * named Version 4: do_step_dxdt_impl( sys , in , dxdt , t , out, dt )
+ *
+ * this version is needed when this stepper is used for initializing
+ * multistep stepper like adams-bashforth. Hence we provide an explicitely
+ * named version. Meant for internal use only.
+ */
+ template < class System, class StateIn, class DerivIn, class StateOut >
+ void do_step_dxdt_impl( System system, const StateIn &in,
+ const DerivIn &dxdt, time_type t, StateOut &out,
+ time_type dt )
+ {
+ this->stepper().do_step_impl( system , in , dxdt , t , out , dt );
+ }
+
template< class StateIn >
void adjust_size( const StateIn &x )
{
diff --git a/boost/numeric/odeint/stepper/base/symplectic_rkn_stepper_base.hpp b/boost/numeric/odeint/stepper/base/symplectic_rkn_stepper_base.hpp
index ea3523e..eb09aef 100644
--- a/boost/numeric/odeint/stepper/base/symplectic_rkn_stepper_base.hpp
+++ b/boost/numeric/odeint/stepper/base/symplectic_rkn_stepper_base.hpp
@@ -182,7 +182,7 @@ private:
// stepper for systems with function for dq/dt = f(p) and dp/dt = -f(q)
template< class System , class StateIn , class StateOut >
- void do_step_impl( System system , const StateIn &in , time_type t , StateOut &out , time_type dt , boost::mpl::true_ )
+ void do_step_impl( System system , const StateIn &in , time_type /* t */ , StateOut &out , time_type dt , boost::mpl::true_ )
{
typedef typename odeint::unwrap_reference< System >::type system_type;
typedef typename odeint::unwrap_reference< typename system_type::first_type >::type coor_deriv_func_type;
diff --git a/boost/numeric/odeint/stepper/bulirsch_stoer.hpp b/boost/numeric/odeint/stepper/bulirsch_stoer.hpp
index b604f86..e71008a 100644
--- a/boost/numeric/odeint/stepper/bulirsch_stoer.hpp
+++ b/boost/numeric/odeint/stepper/bulirsch_stoer.hpp
@@ -194,7 +194,6 @@ public:
static const value_type val1( 1.0 );
- typename odeint::unwrap_reference< System >::type &sys = system;
if( m_resizer.adjust_size( in , detail::bind( &controlled_error_bs_type::template resize_impl< StateIn > , detail::ref( *this ) , detail::_1 ) ) )
{
reset(); // system resized -> reset
@@ -219,12 +218,12 @@ public:
m_midpoint.set_steps( m_interval_sequence[k] );
if( k == 0 )
{
- m_midpoint.do_step( sys , in , dxdt , t , out , dt );
+ m_midpoint.do_step( system , in , dxdt , t , out , dt );
/* the first step, nothing more to do */
}
else
{
- m_midpoint.do_step( sys , in , dxdt , t , m_table[k-1].m_v , dt );
+ m_midpoint.do_step( system , in , dxdt , t , m_table[k-1].m_v , dt );
extrapolate( k , m_table , m_coeff , out );
// get error estimate
m_algebra.for_each3( m_err.m_v , out , m_table[0].m_v ,
@@ -341,7 +340,7 @@ public:
resize_m_dxdt( x );
resize_m_xnew( x );
resize_impl( x );
- m_midpoint.adjust_size();
+ m_midpoint.adjust_size( x );
}
diff --git a/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp b/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp
index c352c96..0036d1d 100644
--- a/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp
+++ b/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp
@@ -155,8 +155,6 @@ public:
static const value_type val1( 1.0 );
- typename odeint::unwrap_reference< System >::type &sys = system;
-
bool reject( true );
time_vector h_opt( m_k_max+1 );
@@ -261,6 +259,7 @@ public:
{
//calculate dxdt for next step and dense output
+ typename odeint::unwrap_reference< System >::type &sys = system;
sys( out , dxdt_new , t+dt );
//prepare dense output
@@ -367,7 +366,7 @@ public:
void adjust_size( const StateIn &x )
{
resize_impl( x );
- m_midpoint.adjust_size();
+ m_midpoint.adjust_size( x );
}
diff --git a/boost/numeric/odeint/stepper/detail/adams_bashforth_call_algebra.hpp b/boost/numeric/odeint/stepper/detail/adams_bashforth_call_algebra.hpp
index 168bc69..cc279ab 100644
--- a/boost/numeric/odeint/stepper/detail/adams_bashforth_call_algebra.hpp
+++ b/boost/numeric/odeint/stepper/detail/adams_bashforth_call_algebra.hpp
@@ -117,7 +117,7 @@ struct adams_bashforth_call_algebra< 7 , Algebra , Operations >
{
//BOOST_ASSERT( false ); // not implemented
typedef typename Coefficients::value_type value_type;
- Algebra::for_each9( out , in , steps[0].m_v , steps[1].m_v , steps[2].m_v , steps[3].m_v , steps[4].m_v , steps[5].m_v , steps[6].m_v ,
+ algebra.for_each9( out , in , steps[0].m_v , steps[1].m_v , steps[2].m_v , steps[3].m_v , steps[4].m_v , steps[5].m_v , steps[6].m_v ,
typename Operations::template scale_sum8< value_type , Time , Time , Time , Time , Time , Time >(
1.0 , dt * coef[0] , dt * coef[1] , dt * coef[2] , dt * coef[3] , dt * coef[4] , dt * coef[5] , dt * coef[6] ) );
}
@@ -132,7 +132,7 @@ struct adams_bashforth_call_algebra< 8 , Algebra , Operations >
{
//BOOST_ASSERT( false ); // not implemented
typedef typename Coefficients::value_type value_type;
- Algebra::for_each10( out , in , steps[0].m_v , steps[1].m_v , steps[2].m_v , steps[3].m_v , steps[4].m_v , steps[5].m_v , steps[6].m_v , steps[7].m_v ,
+ algebra.for_each10( out , in , steps[0].m_v , steps[1].m_v , steps[2].m_v , steps[3].m_v , steps[4].m_v , steps[5].m_v , steps[6].m_v , steps[7].m_v ,
typename Operations::template scale_sum9< value_type , Time , Time , Time , Time , Time , Time , Time >(
1.0 , dt * coef[0] , dt * coef[1] , dt * coef[2] , dt * coef[3] , dt * coef[4] , dt * coef[5] , dt * coef[6] , dt * coef[7] ) );
}
diff --git a/boost/numeric/odeint/stepper/detail/adams_moulton_call_algebra.hpp b/boost/numeric/odeint/stepper/detail/adams_moulton_call_algebra.hpp
index f8a1505..b6f5f2a 100644
--- a/boost/numeric/odeint/stepper/detail/adams_moulton_call_algebra.hpp
+++ b/boost/numeric/odeint/stepper/detail/adams_moulton_call_algebra.hpp
@@ -32,7 +32,7 @@ template< class Algebra , class Operations >
struct adams_moulton_call_algebra< 1 , Algebra , Operations >
{
template< class StateIn , class StateOut , class DerivIn , class StepStorage , class Coefficients , class Time >
- void operator()( Algebra &algebra , const StateIn &in , StateOut &out , const DerivIn &dxdt , const StepStorage &steps , const Coefficients &coef , Time dt ) const
+ void operator()( Algebra &algebra , const StateIn &in , StateOut &out , const DerivIn &dxdt , const StepStorage& /* steps */ , const Coefficients &coef , Time dt ) const
{
typedef typename Coefficients::value_type value_type;
algebra.for_each3( out , in , dxdt , typename Operations::template scale_sum2< value_type , Time >( 1.0 , dt * coef[0] ) );
diff --git a/boost/numeric/odeint/stepper/euler.hpp b/boost/numeric/odeint/stepper/euler.hpp
index 443f942..1c7c126 100644
--- a/boost/numeric/odeint/stepper/euler.hpp
+++ b/boost/numeric/odeint/stepper/euler.hpp
@@ -76,7 +76,7 @@ public :
{ }
template< class System , class StateIn , class DerivIn , class StateOut >
- void do_step_impl( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt )
+ void do_step_impl( System /* system */ , const StateIn &in , const DerivIn &dxdt , time_type /* t */ , StateOut &out , time_type dt )
{
stepper_base_type::m_algebra.for_each3( out , in , dxdt ,
typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , dt ) );
@@ -84,7 +84,7 @@ public :
}
template< class StateOut , class StateIn1 , class StateIn2 >
- void calc_state( StateOut &x , time_type t , const StateIn1 &old_state , time_type t_old , const StateIn2 &current_state , time_type t_new ) const
+ void calc_state( StateOut &x , time_type t , const StateIn1 &old_state , time_type t_old , const StateIn2 & /*current_state*/ , time_type /* t_new */ ) const
{
const time_type delta = t - t_old;
stepper_base_type::m_algebra.for_each3( x , old_state , stepper_base_type::m_dxdt.m_v ,
diff --git a/boost/numeric/odeint/stepper/extrapolation_stepper.hpp b/boost/numeric/odeint/stepper/extrapolation_stepper.hpp
new file mode 100644
index 0000000..3b688e9
--- /dev/null
+++ b/boost/numeric/odeint/stepper/extrapolation_stepper.hpp
@@ -0,0 +1,293 @@
+/*
+ [auto_generated]
+ boost/numeric/odeint/stepper/extrapolation_stepper.hpp
+
+ [begin_description]
+ extrapolation stepper
+ [end_description]
+
+ Copyright 2009-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)
+*/
+
+#ifndef BOOST_NUMERIC_ODEINT_STEPPER_EXTRAPOLATION_STEPPER_HPP_INCLUDED
+#define BOOST_NUMERIC_ODEINT_STEPPER_EXTRAPOLATION_STEPPER_HPP_INCLUDED
+
+#include <iostream>
+
+#include <algorithm>
+
+#include <boost/config.hpp> // for min/max guidelines
+#include <boost/static_assert.hpp>
+
+#include <boost/numeric/odeint/util/bind.hpp>
+#include <boost/numeric/odeint/util/unwrap_reference.hpp>
+
+#include <boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp>
+#include <boost/numeric/odeint/stepper/modified_midpoint.hpp>
+#include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
+#include <boost/numeric/odeint/algebra/range_algebra.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/util/unit_helper.hpp>
+#include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
+
+namespace boost
+{
+namespace numeric
+{
+namespace odeint
+{
+
+template < unsigned short Order, 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 >
+#ifndef DOXYGEN_SKIP
+class extrapolation_stepper
+ : public explicit_error_stepper_base<
+ extrapolation_stepper< Order, State, Value, Deriv, Time, Algebra,
+ Operations, Resizer >,
+ Order, Order, Order - 2, State, Value, Deriv, Time, Algebra,
+ Operations, Resizer >
+#else
+class extrapolation_stepper : public explicit_error_stepper_base
+#endif
+{
+
+ private:
+ // check for Order being odd
+ BOOST_STATIC_ASSERT_MSG(
+ ( ( Order % 2 ) == 0 ) && ( Order > 2 ),
+ "extrapolation_stepper requires even Order larger than 2" );
+
+ public:
+#ifndef DOXYGEN_SKIP
+ typedef explicit_error_stepper_base<
+ extrapolation_stepper< Order, State, Value, Deriv, Time, Algebra,
+ Operations, Resizer >,
+ Order, Order, Order - 2, State, Value, Deriv, Time, Algebra, Operations,
+ Resizer > stepper_base_type;
+#else
+ typedef explicit_error_stepper_base< extrapolation_stepper< ... >, ... >
+ stepper_base_type;
+#endif
+
+ typedef typename stepper_base_type::state_type state_type;
+ typedef typename stepper_base_type::value_type value_type;
+ typedef typename stepper_base_type::deriv_type deriv_type;
+ typedef typename stepper_base_type::time_type time_type;
+ typedef typename stepper_base_type::algebra_type algebra_type;
+ typedef typename stepper_base_type::operations_type operations_type;
+ typedef typename stepper_base_type::resizer_type resizer_type;
+
+#ifndef DOXYGEN_SKIP
+ typedef typename stepper_base_type::stepper_type stepper_type;
+ typedef typename stepper_base_type::wrapped_state_type wrapped_state_type;
+ typedef typename stepper_base_type::wrapped_deriv_type wrapped_deriv_type;
+
+ typedef std::vector< value_type > value_vector;
+ typedef std::vector< value_vector > value_matrix;
+ typedef std::vector< size_t > int_vector;
+ typedef std::vector< wrapped_state_type > state_table_type;
+ typedef modified_midpoint< state_type, value_type, deriv_type, time_type,
+ algebra_type, operations_type,
+ resizer_type > midpoint_stepper_type;
+
+#endif // DOXYGEN_SKIP
+
+ typedef unsigned short order_type;
+ static const order_type order_value = stepper_base_type::order_value;
+ static const order_type stepper_order_value =
+ stepper_base_type::stepper_order_value;
+ static const order_type error_order_value =
+ stepper_base_type::error_order_value;
+
+ const static size_t m_k_max = ( order_value - 2 ) / 2;
+
+ extrapolation_stepper( const algebra_type &algebra = algebra_type() )
+ : stepper_base_type( algebra ), m_interval_sequence( m_k_max + 1 ),
+ m_coeff( m_k_max + 1 ), m_table( m_k_max )
+ {
+ for ( unsigned short i = 0; i < m_k_max + 1; i++ )
+ {
+ m_interval_sequence[i] = 2 * ( i + 1 );
+ m_coeff[i].resize( i );
+ 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] =
+ static_cast< value_type >( 1 ) /
+ ( r * r - static_cast< value_type >(
+ 1 ) ); // coefficients for extrapolation
+ }
+ }
+ }
+
+ template < class System, class StateIn, class DerivIn, class StateOut,
+ class Err >
+ void do_step_impl( System system, const StateIn &in, const DerivIn &dxdt,
+ time_type t, StateOut &out, time_type dt, Err &xerr )
+ {
+ // std::cout << "dt: " << dt << std::endl;
+ // normal step
+ do_step_impl( system, in, dxdt, t, out, dt );
+
+ static const value_type val1( 1.0 );
+ // additionally, perform the error calculation
+ stepper_base_type::m_algebra.for_each3(
+ xerr, out, m_table[0].m_v,
+ typename operations_type::template scale_sum2<
+ value_type, value_type >( val1, -val1 ) );
+ }
+
+ template < class System, class StateInOut, class DerivIn, class Err >
+ void do_step_impl_io( System system, StateInOut &inout, const DerivIn &dxdt,
+ time_type t, time_type dt, Err &xerr )
+ {
+ // normal step
+ do_step_impl_io( system, inout, dxdt, t, dt );
+
+ static const value_type val1( 1.0 );
+ // additionally, perform the error calculation
+ stepper_base_type::m_algebra.for_each3(
+ xerr, inout, m_table[0].m_v,
+ typename operations_type::template scale_sum2<
+ value_type, value_type >( val1, -val1 ) );
+ }
+
+ template < class System, class StateIn, class DerivIn, class StateOut >
+ void do_step_impl( System system, const StateIn &in, const DerivIn &dxdt,
+ time_type t, StateOut &out, time_type dt )
+ {
+ m_resizer.adjust_size(
+ in, detail::bind( &stepper_type::template resize_impl< StateIn >,
+ detail::ref( *this ), detail::_1 ) );
+ size_t k = 0;
+ m_midpoint.set_steps( m_interval_sequence[k] );
+ m_midpoint.do_step( system, in, dxdt, t, out, dt );
+ for ( k = 1; k <= m_k_max; ++k )
+ {
+ m_midpoint.set_steps( m_interval_sequence[k] );
+ m_midpoint.do_step( system, in, dxdt, t, m_table[k - 1].m_v, dt );
+ extrapolate( k, m_table, m_coeff, out );
+ }
+ }
+
+ template < class System, class StateInOut, class DerivIn >
+ void do_step_impl_io( System system, StateInOut &inout, const DerivIn &dxdt,
+ time_type t, time_type dt )
+ {
+ // special care for inout
+ m_xout_resizer.adjust_size(
+ inout,
+ detail::bind( &stepper_type::template resize_m_xout< StateInOut >,
+ detail::ref( *this ), detail::_1 ) );
+ do_step_impl( system, inout, dxdt, t, m_xout.m_v, dt );
+ boost::numeric::odeint::copy( m_xout.m_v, inout );
+ }
+
+ template < class System, class StateInOut, class DerivIn >
+ void do_step_dxdt_impl( System system, StateInOut &x, const DerivIn &dxdt,
+ time_type t, time_type dt )
+ {
+ do_step_impl_io( system , x , dxdt , t , dt );
+ }
+
+ template < class System, class StateIn, class DerivIn, class StateOut >
+ void do_step_dxdt_impl( System system, const StateIn &in,
+ const DerivIn &dxdt, time_type t, StateOut &out,
+ time_type dt )
+ {
+ do_step_impl( system , in , dxdt , t , out , dt );
+ }
+
+
+ template < class StateIn > void adjust_size( const StateIn &x )
+ {
+ resize_impl( x );
+ m_midpoint.adjust_size( x );
+ }
+
+ private:
+ template < class StateIn > bool resize_impl( const StateIn &x )
+ {
+ bool resized( false );
+ for ( size_t i = 0; i < m_k_max; ++i )
+ resized |= adjust_size_by_resizeability(
+ m_table[i], x, typename is_resizeable< state_type >::type() );
+ return resized;
+ }
+
+ template < class StateIn > bool resize_m_xout( const StateIn &x )
+ {
+ return adjust_size_by_resizeability(
+ m_xout, x, typename is_resizeable< state_type >::type() );
+ }
+
+ template < class StateInOut >
+ void extrapolate( size_t k, state_table_type &table,
+ const value_matrix &coeff, StateInOut &xest )
+ /* polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf
+ uses the obtained intermediate results to extrapolate to dt->0
+ */
+ {
+ static const value_type val1 = static_cast< value_type >( 1.0 );
+
+ for ( int j = k - 1; j > 0; --j )
+ {
+ stepper_base_type::m_algebra.for_each3(
+ table[j - 1].m_v, table[j].m_v, table[j - 1].m_v,
+ typename operations_type::template scale_sum2<
+ value_type, value_type >( val1 + coeff[k][j],
+ -coeff[k][j] ) );
+ }
+ stepper_base_type::m_algebra.for_each3(
+ xest, table[0].m_v, xest,
+ typename operations_type::template scale_sum2<
+ value_type, value_type >( val1 + coeff[k][0], -coeff[k][0] ) );
+ }
+
+ private:
+ midpoint_stepper_type m_midpoint;
+
+ resizer_type m_resizer;
+ resizer_type m_xout_resizer;
+
+ int_vector m_interval_sequence; // stores the successive interval counts
+ value_matrix m_coeff;
+
+ wrapped_state_type m_xout;
+ state_table_type m_table; // sequence of states for extrapolation
+};
+
+/******** DOXYGEN *******/
+
+/**
+ * \class extrapolation_stepper
+ * \brief Extrapolation stepper with configurable order, and error estimation.
+ *
+ * The extrapolation stepper is a stepper with error estimation and configurable
+ * order. The order is given as template parameter and needs to be an _odd_
+ * number. The stepper is based on several executions of the modified midpoint
+ * method and a Richardson extrapolation. This is essentially the same technique
+ * as for bulirsch_stoer, but without the variable order.
+ *
+ * \note The Order parameter has to be an even number greater 2.
+ */
+}
+}
+}
+#endif
diff --git a/boost/numeric/odeint/stepper/implicit_euler.hpp b/boost/numeric/odeint/stepper/implicit_euler.hpp
index f301988..e1c6416 100644
--- a/boost/numeric/odeint/stepper/implicit_euler.hpp
+++ b/boost/numeric/odeint/stepper/implicit_euler.hpp
@@ -144,7 +144,7 @@ private:
void solve( state_type &x , matrix_type &m )
{
int res = boost::numeric::ublas::lu_factorize( m , m_pm.m_v );
- if( res != 0 ) exit(0);
+ if( res != 0 ) std::exit(0);
boost::numeric::ublas::lu_substitute( m , m_pm.m_v , x );
}
diff --git a/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp b/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp
index 4b3824f..df4e6c4 100644
--- a/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp
+++ b/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp
@@ -66,6 +66,7 @@ public:
{
BOOST_USING_STD_MAX();
using std::abs;
+ using std::sqrt;
const size_t n = x.size();
value_type err = 0.0 , sk = 0.0;
diff --git a/boost/numeric/odeint/stepper/runge_kutta4_classic.hpp b/boost/numeric/odeint/stepper/runge_kutta4_classic.hpp
index 2b8e70b..32bda0b 100644
--- a/boost/numeric/odeint/stepper/runge_kutta4_classic.hpp
+++ b/boost/numeric/odeint/stepper/runge_kutta4_classic.hpp
@@ -121,11 +121,20 @@ public :
// dt * m_dxh = k4
sys( m_x_tmp.m_v , m_dxh.m_v , t + dt );
+
//x += dt/6 * ( m_dxdt + m_dxt + val2*m_dxm )
time_type dt6 = dt / static_cast< value_type >( 6 );
time_type dt3 = dt / static_cast< value_type >( 3 );
stepper_base_type::m_algebra.for_each6( out , in , dxdt , m_dxt.m_v , m_dxm.m_v , m_dxh.m_v ,
- typename operations_type::template scale_sum5< value_type , time_type , time_type , time_type , time_type >( 1.0 , dt6 , dt3 , dt3 , dt6 ) );
+ typename operations_type::template scale_sum5< value_type , time_type , time_type , time_type , time_type >( 1.0 , dt6 , dt3 , dt3 , dt6 ) );
+
+ // x += dt/6 * m_dxdt + dt/3 * m_dxt )
+ // stepper_base_type::m_algebra.for_each4( out , in , dxdt , m_dxt.m_v ,
+ // typename operations_type::template scale_sum3< value_type , time_type , time_type >( 1.0 , dt6 , dt3 ) );
+ // // x += dt/3 * m_dxm + dt/6 * m_dxh )
+ // stepper_base_type::m_algebra.for_each4( out , out , m_dxm.m_v , m_dxh.m_v ,
+ // typename operations_type::template scale_sum3< value_type , time_type , time_type >( 1.0 , dt3 , dt6 ) );
+
}
template< class StateType >
diff --git a/boost/numeric/odeint/stepper/velocity_verlet.hpp b/boost/numeric/odeint/stepper/velocity_verlet.hpp
index 24f3c0e..3a20fc2 100644
--- a/boost/numeric/odeint/stepper/velocity_verlet.hpp
+++ b/boost/numeric/odeint/stepper/velocity_verlet.hpp
@@ -114,7 +114,7 @@ public:
algebra_stepper_base_type::m_algebra.for_each4(
qout , qin , pin , ain ,
- typename operations_type::template scale_sum3< value_type , time_type , time_square_type >( one , one_half * dt , one * dt * dt ) );
+ typename operations_type::template scale_sum3< value_type , time_type , time_square_type >( one , one * dt , one_half * dt * dt ) );
typename odeint::unwrap_reference< System >::type & sys = system;