Subject: [Boost-commit] svn:boost r68831 - in sandbox/odeint: boost/numeric/odeint branches/karsten branches/karsten/boost/numeric/odeint/stepper branches/karsten/boost/numeric/odeint/stepper/base branches/karsten/boost/numeric/odeint/stepper/detail branches/karsten/libs/numeric/odeint/regression_test branches/karsten/libs/numeric/odeint/test
From: karsten.ahnert_at_[hidden]
Date: 2011-02-13 08:00:00
Author: karsten
Date: 2011-02-13 07:59:58 EST (Sun, 13 Feb 2011)
New Revision: 68831
URL: http://svn.boost.org/trac/boost/changeset/68831
Log:
symplectic runge kutta nystroem methods created
Added:
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/symplectic_nystroem_stepper_base.hpp (contents, props changed)
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/detail/is_pair.hpp (contents, props changed)
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/symplectic_euler.hpp (contents, props changed)
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/symplectic_rkn_sb3a_mclachlan.hpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/regression_test/hamiltonian_steppers.cpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/test/is_pair.cpp (contents, props changed)
Properties modified:
sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp (contents, props changed)
Text files modified:
sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp | 1 -
sandbox/odeint/branches/karsten/Jamroot | 1 +
sandbox/odeint/branches/karsten/TODO | 3 +++
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/detail/macros.hpp | 20 ++++++++++++++++++++
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/stepper_categories.hpp | 1 +
sandbox/odeint/branches/karsten/libs/numeric/odeint/regression_test/Jamfile | 5 +++++
sandbox/odeint/branches/karsten/libs/numeric/odeint/test/Jamfile | 1 +
7 files changed, 31 insertions(+), 1 deletions(-)
Modified: sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp 2011-02-13 07:59:58 EST (Sun, 13 Feb 2011)
@@ -59,7 +59,6 @@
stepper_type m_stepper;
const time_type m_fac;
-
//
// public interface
//
Modified: sandbox/odeint/branches/karsten/Jamroot
==============================================================================
--- sandbox/odeint/branches/karsten/Jamroot (original)
+++ sandbox/odeint/branches/karsten/Jamroot 2011-02-13 07:59:58 EST (Sun, 13 Feb 2011)
@@ -31,6 +31,7 @@
build-project libs/numeric/odeint/ideas/generic_stepper ;
build-project libs/numeric/odeint/ideas/rosenbrock4 ;
build-project libs/numeric/odeint/ideas/units ;
+build-project libs/numeric/odeint/ideas/algebra ;
# docs:
Modified: sandbox/odeint/branches/karsten/TODO
==============================================================================
--- sandbox/odeint/branches/karsten/TODO (original)
+++ sandbox/odeint/branches/karsten/TODO 2011-02-13 07:59:58 EST (Sun, 13 Feb 2011)
@@ -25,6 +25,9 @@
* check copyright note
* documente every file in the preamble
* check once more, if all contructor, destructors and assign-operators are present
+* symplecit_stepper
+ * find an appropriate name, (symplectic stroemer nystroem)
+ * check is the coefficients are named good
* Integrate functions
* skript for setting the include defines according to the position in file system an writing a general copyright comment at the beginning
Added: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/symplectic_nystroem_stepper_base.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/symplectic_nystroem_stepper_base.hpp 2011-02-13 07:59:58 EST (Sun, 13 Feb 2011)
@@ -0,0 +1,239 @@
+/*
+ * symplectic_nystroem.hpp
+ *
+ * Created on: Feb 12, 2011
+ * Author: karsten
+ */
+
+#ifndef BOOST_NUMERIC_ODEINT_STEPPER_SYMPLECTIC_NYSTROEM_STEPPER_BASE_HPP_INCLUDED
+#define BOOST_NUMERIC_ODEINT_STEPPER_SYMPLECTIC_NYSTROEM_STEPPER_BASE_HPP_INCLUDED
+
+#include <boost/ref.hpp>
+#include <boost/array.hpp>
+
+#include <boost/numeric/odeint/util/size_adjuster.hpp>
+#include <boost/numeric/odeint/util/construct.hpp>
+#include <boost/numeric/odeint/util/destruct.hpp>
+#include <boost/numeric/odeint/util/copy.hpp>
+
+#include <boost/numeric/odeint/stepper/stepper_categories.hpp>
+
+#include <boost/numeric/odeint/stepper/detail/is_pair.hpp>
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+template<
+ size_t NumOfStages ,
+ class Stepper ,
+ class State ,
+ class Value ,
+ class Deriv ,
+ class Time ,
+ class Algebra ,
+ class Operations ,
+ class AdjustSizePolicy
+ >
+class symplectic_nystroem_stepper_base
+{
+
+ void initialize( void )
+ {
+ boost::numeric::odeint::construct( m_dqdt );
+ boost::numeric::odeint::construct( m_dpdt );
+ m_coor_deriv_adjuster.register_state( 0 , m_dqdt );
+ m_momentum_deriv_adjuster.register_state( 0 , m_dpdt );
+ }
+
+ void copy( const symplectic_nystroem_stepper_base &b )
+ {
+ boost::numeric::odeint::copy( b.m_dqdt , m_dqdt );
+ boost::numeric::odeint::copy( b.m_dpdt , m_dpdt );
+ }
+
+public:
+
+ const static size_t num_of_stages = NumOfStages;
+ 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 AdjustSizePolicy adjust_size_policy;
+ typedef Stepper stepper_type;
+ typedef stepper_tag stepper_category;
+
+ typedef typename state_type::first_type coor_type;
+ typedef typename state_type::second_type momentum_type;
+
+ typedef typename deriv_type::first_type coor_deriv_type;
+ typedef typename deriv_type::second_type momentum_deriv_type;
+
+ typedef boost::array< value_type , num_of_stages > coef_type;
+
+ symplectic_nystroem_stepper_base( const coef_type &coef_a , const coef_type &coef_b )
+ : m_coef_a( coef_a ) , m_coef_b( coef_b ) ,
+ m_coor_deriv_adjuster() , m_momentum_deriv_adjuster() ,
+ m_dqdt() , m_dpdt()
+ {
+ initialize();
+ }
+
+ symplectic_nystroem_stepper_base( const symplectic_nystroem_stepper_base &b )
+ : m_coef_a( b.m_coef_a ) , m_coef_b( b.m_coef_b ) ,
+ m_coor_deriv_adjuster() , m_momentum_deriv_adjuster(),
+ m_dqdt() , m_dpdt()
+ {
+ initialize();
+ copy( b );
+ }
+
+ ~symplectic_nystroem_stepper_base( void )
+ {
+ boost::numeric::odeint::destruct( m_dqdt );
+ boost::numeric::odeint::destruct( m_dpdt );
+ }
+
+ symplectic_nystroem_stepper_base& operator=( const symplectic_nystroem_stepper_base &b )
+ {
+ copy( b );
+ return *this;
+ }
+
+
+ template< class System , class StateIn , class StateOut >
+ void do_step( System system , const StateIn &in , const time_type &t , const StateOut &out , const time_type &dt )
+ {
+ do_step( system , in , t , out , dt );
+ }
+
+
+ template< class System , class StateIn , class StateOut >
+ void do_step( System system , const StateIn &in , const time_type &t , StateOut &out , const time_type &dt )
+ {
+ typedef typename boost::unwrap_reference< System >::type system_type;
+ do_step_impl( system , in , t , out , dt , typename detail::is_pair< system_type >::type() );
+ }
+
+
+
+
+ template< class System , class StateInOut >
+ void do_step( System system , const StateInOut &state , const time_type &t , const time_type &dt )
+ {
+ do_step( system , state , t , dt );
+ }
+
+ template< class System , class StateInOut >
+ void do_step( System system , StateInOut &state , const time_type &t , const time_type &dt )
+ {
+ typedef typename boost::unwrap_reference< System >::type system_type;
+ do_step_impl( system , state , t , state , dt , typename detail::is_pair< system_type >::type() );
+ }
+
+ template< class StateType >
+ void adjust_size( const StateType &x )
+ {
+ m_coor_deriv_adjuster.adjust_size( x );
+ m_momentum_deriv_adjuster.adjust_size( x );
+ }
+
+ const coef_type& coef_a( void ) const { return m_coef_a; }
+ const coef_type& coef_b( void ) const { return m_coef_b; }
+
+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 , const time_type &t , StateOut &out , const time_type &dt , boost::mpl::true_ )
+ {
+ typedef typename boost::unwrap_reference< System >::type system_type;
+ typedef typename boost::unwrap_reference< typename system_type::first_type >::type coor_deriv_func_type;
+ typedef typename boost::unwrap_reference< typename system_type::second_type >::type momentum_deriv_func_type;
+ system_type &sys = system;
+ coor_deriv_func_type &coor_func = sys.first;
+ momentum_deriv_func_type &momentum_func = sys.second;
+
+
+ m_coor_deriv_adjuster.adjust_size_by_policy( in , adjust_size_policy() );
+ m_momentum_deriv_adjuster.adjust_size_by_policy( in , adjust_size_policy() );
+
+ // ToDo: check sizes?
+
+ for( size_t l=0 ; l<num_of_stages ; ++l )
+ {
+ if( l == 0 )
+ {
+ coor_func( in.second , m_dqdt );
+ typename algebra_type::for_each2()( out.first , in.first , m_dqdt ,
+ typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , coef_a[l] * dt ) );
+ momentum_func( out.first , m_dqdt );
+ typename algebra_type::for_each2()( out.second , in.second , m_dqdt ,
+ typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , coef_b[l] * dt ) );
+ }
+ else
+ {
+ coor_func( out.second , m_dqdt );
+ typename algebra_type::for_each2()( out.first , out.first , m_dqdt ,
+ typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , coef_a[l] * dt ) );
+ momentum_func( out.first , m_dqdt );
+ typename algebra_type::for_each2()( out.second , out.second , m_dqdt ,
+ typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , coef_b[l] * dt ) );
+ }
+ }
+ }
+
+
+ // stepper for systems with only function dp /dt = -f(q), dq/dt = p
+ template< class System , class StateIn , class StateOut >
+ void do_step_impl( System system , const StateIn &in , const time_type &t , StateOut &out , const time_type &dt , boost::mpl::false_ )
+ {
+ typedef typename boost::unwrap_reference< System >::type momentum_deriv_func_type;
+ momentum_deriv_func_type &momentum_func = system;
+
+
+ m_coor_deriv_adjuster.adjust_size_by_policy( in , adjust_size_policy() );
+ m_momentum_deriv_adjuster.adjust_size_by_policy( in , adjust_size_policy() );
+
+ // ToDo: check sizes?
+
+
+ for( size_t l=0 ; l<num_of_stages ; ++l )
+ {
+ if( l == 0 )
+ {
+ typename algebra_type::for_each3()( out.first , in.first , in.second ,
+ typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_a[l] * dt ) );
+ momentum_func( out.first , m_dqdt );
+ typename algebra_type::for_each3()( out.second , in.second , m_dqdt ,
+ typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_b[l] * dt ) );
+ }
+ else
+ {
+ typename algebra_type::for_each3()( out.first , out.first , out.second ,
+ typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_a[l] * dt ) );
+ momentum_func( out.first , m_dqdt );
+ typename algebra_type::for_each3()( out.second , out.second , m_dqdt ,
+ typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_b[l] * dt ) );
+ }
+ }
+ }
+
+ // ToDo : check if a and b are appropriate names for the coefficients
+ const coef_type m_coef_a;
+ const coef_type m_coef_b;
+
+ size_adjuster< coor_deriv_type , 1 > m_coor_deriv_adjuster;
+ size_adjuster< momentum_deriv_type , 1 > m_momentum_deriv_adjuster;
+ coor_deriv_type m_dqdt;
+ momentum_deriv_type m_dpdt;
+};
+
+} // namespace odeint
+} // namespace numeric
+} // namespace boost
+
+
+#endif // BOOST_NUMERIC_ODEINT_STEPPER_SYMPLECTIC_NYSTROEM_STEPPER_BASE_HPP_INCLUDED
Added: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/detail/is_pair.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/detail/is_pair.hpp 2011-02-13 07:59:58 EST (Sun, 13 Feb 2011)
@@ -0,0 +1,35 @@
+/*
+ * is_pair.hpp
+ *
+ * Created on: Feb 12, 2011
+ * Author: karsten
+ */
+
+#ifndef BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_IS_PAIR_HPP_
+#define BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_IS_PAIR_HPP_
+
+#include <boost/mpl/bool.hpp>
+#include <utility>
+
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+namespace detail {
+
+template< class T >
+struct is_pair : public boost::mpl::false_
+{
+};
+
+template< class T1 , class T2 >
+struct is_pair< std::pair< T1 , T2 > > : public boost::mpl::true_
+{
+};
+
+} // namespace detail
+} // namespace odeint
+} // namespace numeric
+} // namespace boost
+
+#endif /* BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_IS_PAIR_HPP_ */
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/detail/macros.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/detail/macros.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/detail/macros.hpp 2011-02-13 07:59:58 EST (Sun, 13 Feb 2011)
@@ -69,4 +69,24 @@
typedef typename stepper_base_type::stepper_type stepper_type
+#define BOOST_ODEINT_SYMPLECTIC_NYSTROEM_STEPPER_TYPEDEFS( STEPPER , STAGES ) \
+ typedef symplectic_nystroem_stepper_base< \
+ STAGES , \
+ STEPPER< State , Value , Deriv , Time , Algebra , Operations , AdjustSizePolicy > , \
+ State , Value , Deriv , Time , Algebra , Operations , AdjustSizePolicy > stepper_base_type; \
+ 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::adjust_size_policy adjust_size_policy; \
+ typedef typename stepper_base_type::stepper_type stepper_type; \
+ typedef typename stepper_base_type::coor_type coor_type; \
+ typedef typename stepper_base_type::momentum_type momentum_type; \
+ typedef typename stepper_base_type::coor_deriv_type coor_deriv_type; \
+ typedef typename stepper_base_type::momentum_deriv_type momentum_deriv_type; \
+ typedef typename stepper_base_type::coef_type coef_type
+
+
#endif //BOOST_BOOST_NUMERIC_ODEINT_DETAIL_MACROS_HPP_INCLUDED
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/stepper_categories.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/stepper_categories.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/stepper_categories.hpp 2011-02-13 07:59:58 EST (Sun, 13 Feb 2011)
@@ -27,6 +27,7 @@
//struct explicit_stepper_tag : stepper_tag {};
//struct implicit_stepper_tag : stepper_tag {};
+
struct error_stepper_tag {};
struct explicit_error_stepper_tag : error_stepper_tag {};
struct explicit_error_stepper_fsal_tag : error_stepper_tag {};
Added: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/symplectic_euler.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/symplectic_euler.hpp 2011-02-13 07:59:58 EST (Sun, 13 Feb 2011)
@@ -0,0 +1,84 @@
+/*
+ * symplectic_euler.hpp
+ *
+ * Created on: Feb 12, 2011
+ * Author: karsten
+ */
+
+#ifndef BOOST_NUMERIC_ODEINT_STEPPER_SYMPLECTIC_EULER_HPP_
+#define BOOST_NUMERIC_ODEINT_STEPPER_SYMPLECTIC_EULER_HPP_
+
+#include <boost/numeric/odeint/stepper/base/symplectic_nystroem_stepper_base.hpp>
+
+#include <boost/numeric/odeint/algebra/range_algebra.hpp>
+#include <boost/numeric/odeint/algebra/default_operations.hpp>
+
+#include <boost/numeric/odeint/stepper/detail/macros.hpp>
+
+#include <boost/array.hpp>
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+namespace detail {
+namespace symplectic_euler_coef {
+
+ template< class Value >
+ struct coef_a_type : public boost::array< Value , 1 >
+ {
+ coef_a_type( void )
+ {
+ (*this)[0] = 1.0;
+ }
+ };
+
+ template< class Value >
+ struct coef_b_type : public boost::array< Value , 1 >
+ {
+ coef_b_type( void )
+ {
+ (*this)[0] = 1.0;
+ }
+ };
+
+} // namespace symplectic_euler_coef
+} // namespace detail
+
+
+
+template<
+ class State ,
+ class Value = double ,
+ class Deriv = State ,
+ class Time = Value ,
+ class Algebra = range_algebra ,
+ class Operations = default_operations ,
+ class AdjustSizePolicy = adjust_size_initially_tag
+ >
+class symplectic_euler :
+ public symplectic_nystroem_stepper_base
+ <
+ 1 ,
+ symplectic_euler< State , Value , Deriv , Time , Algebra , Operations , AdjustSizePolicy > ,
+ State , Value , Deriv , Time , Algebra , Operations , AdjustSizePolicy
+ >
+{
+public:
+
+ BOOST_ODEINT_SYMPLECTIC_NYSTROEM_STEPPER_TYPEDEFS( symplectic_euler , 1 );
+
+ symplectic_euler( void )
+ : stepper_base_type( detail::symplectic_euler_coef::coef_a_type< value_type >() , detail::symplectic_euler_coef::coef_b_type< value_type >() )
+ {
+ }
+};
+
+
+
+} // namespace odeint
+} // namespace numeric
+} // namespace boost
+
+
+#endif /* BOOST_NUMERIC_ODEINT_STEPPER_SYMPLECTIC_EULER_HPP_ */
Added: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/symplectic_rkn_sb3a_mclachlan.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/symplectic_rkn_sb3a_mclachlan.hpp 2011-02-13 07:59:58 EST (Sun, 13 Feb 2011)
@@ -0,0 +1,114 @@
+/*
+ * symplectic_rkn_sb3a_mclachlan.hpp
+ *
+ * Created on: Feb 13, 2011
+ * Author: karsten
+ */
+
+#ifndef SYMPLECTIC_RKN_SB3A_MCLACHLAN_HPP_
+#define SYMPLECTIC_RKN_SB3A_MCLACHLAN_HPP_
+
+#include <boost/numeric/odeint/stepper/base/symplectic_nystroem_stepper_base.hpp>
+
+#include <boost/numeric/odeint/algebra/range_algebra.hpp>
+#include <boost/numeric/odeint/algebra/default_operations.hpp>
+
+#include <boost/numeric/odeint/stepper/detail/macros.hpp>
+
+#include <boost/array.hpp>
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+namespace detail {
+namespace symplectic_rkn_sb3a_mclachlan {
+
+/*
+ rk_a[0]=0.40518861839525227722;
+ rk_a[1]=-0.28714404081652408900;
+ rk_a[2]=0.5-(rk_a[0]+rk_a[1]);
+ rk_a[3]=rk_a[2];
+ rk_a[4]=rk_a[1];
+ rk_a[5]=rk_a[0];
+
+ rk_b[0]=-3.0/73.0;
+ rk_b[1]=17.0/59.0;
+ rk_b[2]=1.0-2.0*(rk_b[0]+rk_b[1]);
+ rk_b[3]=rk_b[1];
+ rk_b[4]=rk_b[0];
+ rk_b[5]=0.0;
+*/
+
+
+ template< class Value >
+ struct coef_a_type : public boost::array< Value , 6 >
+ {
+ coef_a_type( void )
+ {
+ (*this)[0] = 0.40518861839525227722;
+ (*this)[1] = -0.28714404081652408900;
+ (*this)[2] = 0.5 - ( (*this)[0] + (*this)[1] );
+ (*this)[3] = (*this)[2];
+ (*this)[4] = (*this)[1];
+ (*this)[5] = (*this)[0];
+
+ }
+ };
+
+ template< class Value >
+ struct coef_b_type : public boost::array< Value , 6 >
+ {
+ coef_b_type( void )
+ {
+ (*this)[0] = -3.0 / 73.0;
+ (*this)[1] = 17.0 / 59.0;
+ (*this)[2] = 1.0 - 2.0 * ( (*this)[0] + (*this)[1] );
+ (*this)[3] = (*this)[1];
+ (*this)[4] = (*this)[0];
+ (*this)[5] = 0.0;
+ }
+ };
+
+} // namespace symplectic_euler_coef
+} // namespace detail
+
+
+
+template<
+ class State ,
+ class Value = double ,
+ class Deriv = State ,
+ class Time = Value ,
+ class Algebra = range_algebra ,
+ class Operations = default_operations ,
+ class AdjustSizePolicy = adjust_size_initially_tag
+ >
+class symplectic_rkn_sb3a_mclachlan :
+ public symplectic_nystroem_stepper_base
+ <
+ 6 ,
+ symplectic_euler< State , Value , Deriv , Time , Algebra , Operations , AdjustSizePolicy > ,
+ State , Value , Deriv , Time , Algebra , Operations , AdjustSizePolicy
+ >
+{
+public:
+
+ BOOST_ODEINT_SYMPLECTIC_NYSTROEM_STEPPER_TYPEDEFS( symplectic_euler , 6 );
+
+ symplectic_rkn_sb3a_mclachlan( void )
+ : stepper_base_type(
+ detail::symplectic_rkn_sb3a_mclachlan::coef_a_type< value_type >() ,
+ detail::symplectic_rkn_sb3a_mclachlan::coef_b_type< value_type >()
+ )
+ {
+ }
+};
+
+
+
+} // namespace odeint
+} // namespace numeric
+} // namespace boost
+
+#endif /* SYMPLECTIC_RKN_SB3A_MCLACHLAN_HPP_ */
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/regression_test/Jamfile
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/regression_test/Jamfile (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/regression_test/Jamfile 2011-02-13 07:59:58 EST (Sun, 13 Feb 2011)
@@ -7,6 +7,7 @@
project
: requirements
<include>../../../..
+ <include>../../../../../..
;
exe controlled_stepper_evolution
@@ -31,4 +32,8 @@
exe integrate_functions
: integrate_functions.cpp
+ ;
+
+exe hamiltonian_steppers
+ : hamiltonian_steppers.cpp
;
\ No newline at end of file
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/regression_test/hamiltonian_steppers.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/regression_test/hamiltonian_steppers.cpp 2011-02-13 07:59:58 EST (Sun, 13 Feb 2011)
@@ -0,0 +1,184 @@
+/*
+ * hamiltonian_steppers.cpp
+ *
+ * Created on: Feb 12, 2011
+ * Author: karsten
+ */
+
+#include <iostream>
+#include <tr1/array>
+
+#include <boost/accumulators/accumulators.hpp>
+#include <boost/accumulators/statistics.hpp>
+#include <boost/iterator/zip_iterator.hpp>
+
+#include <boost/timer.hpp>
+
+#include <boost/numeric/odeint/stepper/symplectic_euler.hpp>
+#include <boost/numeric/odeint/stepper/symplectic_rkn_sb3a_mclachlan.hpp>
+
+#include <boost/numeric/odeint/hamiltonian_stepper_rk.hpp>
+#include <boost/numeric/odeint/container_traits_tr1_array.hpp>
+
+
+using namespace std;
+using namespace boost::accumulators;
+
+
+typedef accumulator_set< double , stats< tag::mean > > accumulator_type;
+
+ostream& operator<<( ostream& out , accumulator_type &acc )
+{
+ out << boost::accumulators::mean( acc ) << "\t";
+ return out;
+}
+
+typedef boost::timer timer_type;
+
+
+
+
+
+
+typedef std::tr1::array< double , 1 > container_type;
+typedef std::pair< container_type , container_type > state_type;
+
+
+
+struct harm_osc
+{
+ const double m_omega_sq;
+
+ harm_osc( double omega_sq ) : m_omega_sq( omega_sq ) { }
+
+ template< class Coor , class MomentumDeriv >
+ void operator()( const Coor &coor , MomentumDeriv &deriv )
+ {
+ deriv[0] = -m_omega_sq * coor[0];
+ }
+};
+
+
+void test_euler( void )
+{
+ typedef boost::numeric::odeint::symplectic_euler< state_type > stepper_type;
+ stepper_type stepper;
+ stepper_type stepper2( stepper );
+ stepper_type stepper3;
+ stepper3 = stepper;
+
+ state_type state;
+ state.first[0] = 1.0;
+ state.second[0] = 0.0;
+
+ double t = 0.0;
+ const double dt = 0.1;
+ const double omega_sq = 4.0;
+
+ stepper.do_step( harm_osc( omega_sq ) , state , t , dt );
+}
+
+void test_rkn_sb3a_mclachlan( void )
+{
+ typedef boost::numeric::odeint::symplectic_rkn_sb3a_mclachlan< state_type > stepper_type;
+ stepper_type stepper;
+ stepper_type stepper2( stepper );
+ stepper_type stepper3;
+ stepper3 = stepper;
+
+ state_type state;
+ state.first[0] = 1.0;
+ state.second[0] = 0.0;
+
+ double t = 0.0;
+ const double dt = 0.1;
+ const double omega_sq = 4.0;
+
+ stepper.do_step( harm_osc( omega_sq ) , state , t , dt );
+}
+
+void compare_euler_rkn( void )
+{
+ typedef boost::numeric::odeint::symplectic_euler< state_type > stepper_type1;
+ typedef boost::numeric::odeint::symplectic_rkn_sb3a_mclachlan< state_type > stepper_type2;
+
+ stepper_type1 stepper1;
+ stepper_type2 stepper2;
+
+ state_type state1;
+ state1.first[0] = 1.0;
+ state1.second[0] = 0.0;
+ state_type state2 = state1;
+
+ double t = 0.0;
+ const double dt = 0.1;
+ const double omega_sq = 4.0;
+
+ const size_t ilen = 100;
+ for( size_t i=0 ; i<10000 ; ++i )
+ {
+ double q1 = state1.first[0] , p1 = state1.second[0] , q2 = state2.first[0] , p2 = state2.second[0] ;
+ double energy1 = 0.5 * p1 * p1 + 0.5 * omega_sq * q1 * q1;
+ double energy2 = 0.5 * p2 * p2 + 0.5 * omega_sq * q2 * q2;
+ cout << t << "\t" << q1 << "\t" << p1 << "\t" << q2 << "\t" << p2 << "\t" << energy1 << "\t" << energy2 << "\n";
+ for( size_t ii=0 ; ii<ilen ; ++ii , t+= dt )
+ {
+ stepper1.do_step( harm_osc( omega_sq ) , state1 , t , dt );
+ stepper2.do_step( harm_osc( omega_sq ) , state2 , t , dt );
+ }
+ }
+}
+
+void performance_compare( void )
+{
+ typedef boost::numeric::odeint::hamiltonian_stepper_rk_qfunc< container_type > stepper_type1;
+ typedef boost::numeric::odeint::symplectic_rkn_sb3a_mclachlan< state_type > stepper_type2;
+
+ state_type state1;
+ state1.first[0] = 1.0;
+ state1.second[0] = 1.0;
+ state_type state2 = state1;
+
+ const size_t steps = 100000000;
+
+ timer_type timer;
+ accumulator_type acc1 , acc2;
+ harm_osc osc( 4.0 );
+ const double dt = 0.01;
+
+ for( size_t count=0 ; count<10000 ; ++count )
+ {
+ timer.restart();
+ stepper_type1 stepper1;
+ for( size_t i=0 ; i<steps ; ++i )
+ stepper1.do_step( osc , state1 , 0.0 , dt );
+ acc1( timer.elapsed() );
+
+ timer.restart();
+ stepper_type2 stepper2;
+ for( size_t i=0 ; i<steps ; ++i )
+ stepper2.do_step( osc , state2 , 0.0 , dt );
+ acc2( timer.elapsed() );
+
+
+
+ clog << count << "\t";
+ clog << acc1 << "\t" << acc2 << "\t";
+ clog << state1.first[0] << "\t" << state1.second[0] << "\t";
+ clog << state2.first[0] << "\t" << state2.second[0] << endl;
+ }
+
+
+}
+
+
+int main( int argc , char **argv )
+{
+ test_euler();
+ test_rkn_sb3a_mclachlan();
+ compare_euler_rkn();
+ performance_compare();
+
+
+ return 0;
+}
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/test/Jamfile
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/test/Jamfile (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/test/Jamfile 2011-02-13 07:59:58 EST (Sun, 13 Feb 2011)
@@ -28,5 +28,6 @@
[ run stepper_copying.cpp ]
[ run stepper_with_ranges.cpp ]
[ run rosenbrock4.cpp ]
+ [ run is_pair.cpp ]
: <testing.launcher>valgrind
;
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/test/is_pair.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/test/is_pair.cpp 2011-02-13 07:59:58 EST (Sun, 13 Feb 2011)
@@ -0,0 +1,40 @@
+/*
+ * is_pair.cpp
+ *
+ * Created on: Feb 12, 2011
+ * Author: karsten
+ */
+
+#define BOOST_TEST_MODULE odeint_standard_algebra
+
+#include <utility>
+
+#include <boost/test/unit_test.hpp>
+#include <boost/static_assert.hpp>
+
+#include <boost/numeric/odeint/stepper/detail/is_pair.hpp>
+
+using namespace boost::numeric::odeint;
+
+
+
+BOOST_AUTO_TEST_SUITE( is_pair_test )
+
+BOOST_AUTO_TEST_CASE( is_pair )
+{
+ typedef std::pair< int , int > type1;
+ typedef std::pair< int& , int > type2;
+ typedef std::pair< int , int& > type3;
+ typedef std::pair< int& , int& > type4;
+ typedef std::pair< const int , int > type5;
+ typedef std::pair< const int& , int > type6;
+
+ BOOST_STATIC_ASSERT(( detail::is_pair< type1 >::value ));
+ BOOST_STATIC_ASSERT(( detail::is_pair< type2 >::value ));
+ BOOST_STATIC_ASSERT(( detail::is_pair< type3 >::value ));
+ BOOST_STATIC_ASSERT(( detail::is_pair< type4 >::value ));
+ BOOST_STATIC_ASSERT(( detail::is_pair< type5 >::value ));
+ BOOST_STATIC_ASSERT(( detail::is_pair< type6 >::value ));
+}
+
+BOOST_AUTO_TEST_SUITE_END()