Subject: [Boost-commit] svn:boost r62953 - sandbox/odeint/boost/numeric/odeint
From: mario.mulansky_at_[hidden]
Date: 2010-06-14 18:04:04
Author: mariomulansky
Date: 2010-06-14 18:04:03 EDT (Mon, 14 Jun 2010)
New Revision: 62953
URL: http://svn.boost.org/trac/boost/changeset/62953
Log:
test
Text files modified:
sandbox/odeint/boost/numeric/odeint/hamiltonian_stepper_rk.hpp | 194 ++++++++++++++++++++++++++++++---------
1 files changed, 147 insertions(+), 47 deletions(-)
Modified: sandbox/odeint/boost/numeric/odeint/hamiltonian_stepper_rk.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/hamiltonian_stepper_rk.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/hamiltonian_stepper_rk.hpp 2010-06-14 18:04:03 EDT (Mon, 14 Jun 2010)
@@ -1,5 +1,5 @@
/*
- boost header: numeric/odeint/hamiltonian_stepper_euler.hpp
+ boost header: numeric/odeint/hamiltonian_stepper_rk.hpp
Copyright 2009 Karsten Ahnert
Copyright 2009 Mario Mulansky
@@ -23,6 +23,8 @@
namespace numeric {
namespace odeint {
+
+
template<
class Container ,
class Time = double ,
@@ -44,7 +46,7 @@
private:
- container_type m_dpdt;
+ container_type m_dqdt , m_dpdt;
/*
rk_a[0]=0.40518861839525227722;
@@ -64,57 +66,155 @@
public:
- template< class CoordinateFunction >
- void do_step( CoordinateFunction qfunc ,
- state_type &x ,
- time_type t ,
- time_type dt )
+ template< class SystemFunction >
+ void do_step( SystemFunction func ,
+ state_type &state ,
+ time_type t ,
+ time_type dt )
{
- const size_t order = 6;
+ const size_t order = 6;
- const time_type rk_a[order] = {
- static_cast<time_type>( 0.40518861839525227722 ) ,
- static_cast<time_type>( -0.28714404081652408900 ) ,
- static_cast<time_type>( 0.3819554224212718118 ) ,
- static_cast<time_type>( 0.3819554224212718118 ) ,
- static_cast<time_type>( -0.28714404081652408900 ) ,
- static_cast<time_type>( 0.40518861839525227722 )
- };
- const time_type rk_b[order] = {
- static_cast<time_type>( -3.0/73.0 ) ,
- static_cast<time_type>( 17.0/59.0 ) ,
- static_cast<time_type>( 0.50592059438123984212 ) ,
- static_cast<time_type>( 17.0/59.0 ) ,
- static_cast<time_type>( -3.0/73.0 ) ,
- static_cast<time_type>( 0.0 )
- };
-
- container_type &q = x.first;
- container_type &p = x.second;
-
- if( !traits_type::same_size( q , p ) )
- {
- std::string msg( "hamiltonian_stepper_rk::do_step(): " );
- msg += "q and p have different sizes";
- throw std::invalid_argument( msg );
- }
+ const time_type rk_a[order] = {
+ static_cast<time_type>( 0.40518861839525227722 ) ,
+ static_cast<time_type>( -0.28714404081652408900 ) ,
+ static_cast<time_type>( 0.3819554224212718118 ) ,
+ static_cast<time_type>( 0.3819554224212718118 ) ,
+ static_cast<time_type>( -0.28714404081652408900 ) ,
+ static_cast<time_type>( 0.40518861839525227722 )
+ };
+ const time_type rk_b[order] = {
+ static_cast<time_type>( -3.0/73.0 ) ,
+ static_cast<time_type>( 17.0/59.0 ) ,
+ static_cast<time_type>( 0.50592059438123984212 ) ,
+ static_cast<time_type>( 17.0/59.0 ) ,
+ static_cast<time_type>( -3.0/73.0 ) ,
+ static_cast<time_type>( 0.0 )
+ };
+
+ container_type &q = state.first;
+ container_type &p = state.second;
+
+ if( !traits_type::same_size( q , p ) )
+ {
+ std::string msg( "hamiltonian_stepper_rk::do_step(): " );
+ msg += "q and p have different sizes";
+ throw std::invalid_argument( msg );
+ }
+ traits_type::adjust_size( p , m_dqdt );
traits_type::adjust_size( p , m_dpdt );
- for( size_t l=0 ; l<order ; ++l )
- {
- detail::it_algebra::increment( traits_type::begin(q) ,
- traits_type::end(q) ,
- traits_type::begin(p) ,
- rk_a[l]*dt );
- qfunc( q , m_dpdt );
- detail::it_algebra::increment( traits_type::begin(p) ,
- traits_type::end(p) ,
- traits_type::begin(m_dpdt) ,
- rk_b[l]*dt );
- }
- }
+ for( size_t l=0 ; l<order ; ++l )
+ {
+ func.first( p , m_dqdt );
+ detail::it_algebra::increment( traits_type::begin(q) ,
+ traits_type::end(q) ,
+ traits_type::begin(m_dqdt) ,
+ rk_a[l]*dt );
+ func.second( q , m_dpdt );
+ detail::it_algebra::increment( traits_type::begin(p) ,
+ traits_type::end(p) ,
+ traits_type::begin(m_dpdt) ,
+ rk_b[l]*dt );
+ }
+ }
+
+ };
+
+
+
+ template<
+ class Container ,
+ class Time = double ,
+ class Traits = container_traits< Container >
+ >
+ class hamiltonian_stepper_rk_qfunc
+ {
+ // provide basic typedefs
+ public:
+
+ typedef unsigned short order_type;
+ typedef Time time_type;
+ typedef Traits traits_type;
+ typedef typename traits_type::container_type container_type;
+ typedef std::pair< container_type , container_type > state_type;
+ typedef typename traits_type::value_type value_type;
+ typedef typename traits_type::iterator iterator;
+ typedef typename traits_type::const_iterator const_iterator;
+
+ private:
+
+ container_type m_dpdt;
+
+/*
+ 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;
+*/
+
+ public:
+
+ template< class CoordinateFunction >
+ void do_step( CoordinateFunction qfunc ,
+ state_type &state ,
+ time_type t ,
+ time_type dt )
+ {
+ const size_t order = 6;
+
+ const time_type rk_a[order] = {
+ static_cast<time_type>( 0.40518861839525227722 ) ,
+ static_cast<time_type>( -0.28714404081652408900 ) ,
+ static_cast<time_type>( 0.3819554224212718118 ) ,
+ static_cast<time_type>( 0.3819554224212718118 ) ,
+ static_cast<time_type>( -0.28714404081652408900 ) ,
+ static_cast<time_type>( 0.40518861839525227722 )
+ };
+ const time_type rk_b[order] = {
+ static_cast<time_type>( -3.0/73.0 ) ,
+ static_cast<time_type>( 17.0/59.0 ) ,
+ static_cast<time_type>( 0.50592059438123984212 ) ,
+ static_cast<time_type>( 17.0/59.0 ) ,
+ static_cast<time_type>( -3.0/73.0 ) ,
+ static_cast<time_type>( 0.0 )
+ };
+
+ container_type &q = state.first;
+ container_type &p = state.second;
+
+ if( !traits_type::same_size( q , p ) )
+ {
+ std::string msg( "hamiltonian_stepper_rk::do_step(): " );
+ msg += "q and p have different sizes";
+ throw std::invalid_argument( msg );
+ }
+
+ traits_type::adjust_size( p , m_dpdt );
+ for( size_t l=0 ; l<order ; ++l )
+ {
+ detail::it_algebra::increment( traits_type::begin(q) ,
+ traits_type::end(q) ,
+ traits_type::begin(p) ,
+ rk_a[l]*dt );
+ qfunc( q , m_dpdt );
+ detail::it_algebra::increment( traits_type::begin(p) ,
+ traits_type::end(p) ,
+ traits_type::begin(m_dpdt) ,
+ rk_b[l]*dt );
+ }
+ }
+
};