[Boost-commit] svn:boost r62953 - sandbox/odeint/boost/numeric/odeint

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 );
+            }
+        }
+        
     };