[Boost-commit] svn:boost r62954 - in sandbox/odeint: . boost/numeric/odeint libs/numeric/odeint/stuff

Subject: [Boost-commit] svn:boost r62954 - in sandbox/odeint: . boost/numeric/odeint libs/numeric/odeint/stuff
From: mario.mulansky_at_[hidden]
Date: 2010-06-14 18:05:13

Author: mariomulansky
Date: 2010-06-14 18:05:12 EDT (Mon, 14 Jun 2010)
New Revision: 62954
URL: http://svn.boost.org/trac/boost/changeset/62954

Log:
reorganized hamiltonian steppers
Added:
   sandbox/odeint/libs/numeric/odeint/stuff/Jamfile
      - copied, changed from r62936, /sandbox/odeint/libs/numeric/odeint/examples/Jamfile
   sandbox/odeint/libs/numeric/odeint/stuff/test_hamiltonian_steppers.cpp   (contents, props changed)
Text files modified: 
   sandbox/odeint/Jamroot                                                  |     1                                         
   sandbox/odeint/boost/numeric/odeint/hamiltonian_stepper_euler.hpp       |    65 ++++++++++++++++++++++++++++++++++++++- 
   sandbox/odeint/boost/numeric/odeint/hamiltonian_stepper_euler_qfunc.hpp |     2                                         
   sandbox/odeint/libs/numeric/odeint/stuff/Jamfile                        |    17 ---------                               
   4 files changed, 66 insertions(+), 19 deletions(-)

Modified: sandbox/odeint/Jamroot
==============================================================================
--- sandbox/odeint/Jamroot	(original)
+++ sandbox/odeint/Jamroot	2010-06-14 18:05:12 EDT (Mon, 14 Jun 2010)
@@ -13,6 +13,7 @@
 
 build-project libs/numeric/odeint/examples ;
 build-project libs/numeric/odeint/test ;
+build-project libs/numeric/odeint/stuff ;
 # build-project libs/numeric/odeint/doc ;
 
 

Modified: sandbox/odeint/boost/numeric/odeint/hamiltonian_stepper_euler.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/hamiltonian_stepper_euler.hpp	(original)
+++ sandbox/odeint/boost/numeric/odeint/hamiltonian_stepper_euler.hpp	2010-06-14 18:05:12 EDT (Mon, 14 Jun 2010)
@@ -50,8 +50,10 @@
     public:
 
         template< class SystemFunction >
-        void do_step( SystemFunction func , state_type &state ,
-                      time_type t , time_type dt )
+        void do_step( SystemFunction func , 
+                      state_type &state ,
+                      time_type t , 
+                      time_type dt )
         {
             if( !traits_type::same_size( state.first , state.second ) )
             {
@@ -81,6 +83,65 @@
     };
 
 
+    template<
+        class Container ,
+        class Time = double ,
+        class Traits = container_traits< Container >
+        >
+    class hamiltonian_stepper_euler_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;
+
+    public:
+
+        template< class CoordinateFunction >
+        void do_step( CoordinateFunction qfunc ,
+                      state_type &state ,
+                      time_type t ,
+                      time_type dt )
+        {
+
+            container_type &q = state.first;
+            container_type &p = state.second;
+
+            if( !traits_type::same_size( q , p ) )
+            {
+                std::string msg( "hamiltonian_stepper_euler::do_step(): " );
+                msg += "q and p have different sizes";
+                throw std::invalid_argument( msg );
+            }
+            
+            traits_type::adjust_size( p , m_dpdt );
+            
+            detail::it_algebra::increment( traits_type::begin(q) ,
+                                           traits_type::end(q) ,
+                                           traits_type::begin(p) ,
+                                           dt );
+            qfunc( q , m_dpdt );
+            detail::it_algebra::increment( traits_type::begin(p) ,
+                                           traits_type::end(p) ,
+                                           traits_type::begin(m_dpdt) ,
+                                           dt );
+        }
+
+    };
+
+
+
 } // namespace odeint
 } // namespace numeric
 } // namespace boost

Modified: sandbox/odeint/boost/numeric/odeint/hamiltonian_stepper_euler_qfunc.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/hamiltonian_stepper_euler_qfunc.hpp	(original)
+++ sandbox/odeint/boost/numeric/odeint/hamiltonian_stepper_euler_qfunc.hpp	2010-06-14 18:05:12 EDT (Mon, 14 Jun 2010)
@@ -27,7 +27,7 @@
         class Time = double ,
         class Traits = container_traits< Container >
         >
-    class hamiltonian_stepper_qfunc_euler
+    class hamiltonian_stepper_euler_qfunc
     {
         // provide basic typedefs
     public:

Copied: sandbox/odeint/libs/numeric/odeint/stuff/Jamfile (from r62936, /sandbox/odeint/libs/numeric/odeint/examples/Jamfile)
==============================================================================
--- /sandbox/odeint/libs/numeric/odeint/examples/Jamfile	(original)
+++ sandbox/odeint/libs/numeric/odeint/stuff/Jamfile	2010-06-14 18:05:12 EDT (Mon, 14 Jun 2010)
@@ -14,19 +14,4 @@
     : build-dir .
     ;
 
-exe harmonic_oscillator : harmonic_oscillator.cpp ;
-#exe test_container_and_stepper : test_container_and_stepper.cpp ;
-#exe hamiltonian_test : hamiltonian_test.cpp ;
-
-exe lorenz_cmp_rk4_rk_generic : lorenz_cmp_rk4_rk_generic.cpp ;
-exe lorenz_controlled : lorenz_controlled.cpp ;
-exe lorenz_integrate_constant_step : lorenz_integrate_constant_step.cpp ;
-exe lorenz_integrator : lorenz_integrator.cpp ;
-exe lorenz_stepper : lorenz_stepper.cpp ;
-exe pendulum_vibrating_pivot : pendulum_vibrating_pivot.cpp ;
-exe dnls : dnls.cpp ;
-
-exe solar_system : solar_system.cpp ;
-
-exe doc_harm_osc : doc_harm_osc.cpp ;
-exe doc_integrate : doc_integrate.cpp ;
+exe test_hamiltonian_steppers : test_hamiltonian_steppers.cpp ;

Added: sandbox/odeint/libs/numeric/odeint/stuff/test_hamiltonian_steppers.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/libs/numeric/odeint/stuff/test_hamiltonian_steppers.cpp	2010-06-14 18:05:12 EDT (Mon, 14 Jun 2010)
@@ -0,0 +1,84 @@
+#include <iostream>
+#include <vector>
+#include <utility>
+#include <boost/numeric/odeint.hpp>
+
+#define tab "\t"
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+typedef vector< double > coord_type;
+typedef pair< coord_type , coord_type > state_type;
+typedef hamiltonian_stepper_euler< coord_type > euler_stepper_type;
+typedef hamiltonian_stepper_euler_qfunc< coord_type > euler_stepper_type_qfunc;
+typedef hamiltonian_stepper_rk< coord_type > rk_stepper_type;
+typedef hamiltonian_stepper_rk_qfunc< coord_type > rk_stepper_type_qfunc;
+
+static const double omega0 = 1.0;
+static const double omega1 = 2.0;
+
+void qfunc(const coord_type &q, coord_type &dpdt )
+{
+    dpdt[0] = -omega0*omega0*q[0];// + q[1];
+    dpdt[1] = -omega1*omega1*q[1];// + q[0];
+}
+
+void pfunc(const coord_type &p, coord_type &dqdt )
+{
+    dqdt[0] = p[0];
+    dqdt[1] = p[1];
+}
+
+double energy( state_type &state )
+{
+    coord_type &q = state.first;
+    coord_type &p = state.second;
+    return p[0]*p[0]/2.0 + p[1]*p[1]/2.0 
+        + omega0*omega0*q[0]*q[0]/2.0 
+        + omega1*omega1*q[1]*q[1]/2.0;
+    //- q[0]*q[1];
+}
+
+
+int main( int argc , char **argv )
+{
+    coord_type x1(2) , p1(2 , 0.0);
+    x1[0] = 1.0;
+    x1[1] = -1.0;
+    coord_type x2(x1) , x3(x1) , x4(x1);
+    coord_type p2(p1) , p3(p1) , p4(p1);
+    state_type state1 = make_pair( x1 , p1 );
+    state_type state2 = make_pair( x2 , p2 );
+    state_type state3 = make_pair( x3 , p3 );
+    state_type state4 = make_pair( x4 , p4 );
+    double initial_energy = energy( state1 );
+    clog << "initial energy: " << initial_energy << endl;
+
+    euler_stepper_type euler;
+    euler_stepper_type_qfunc euler_qfunc;
+    rk_stepper_type rk;
+    rk_stepper_type_qfunc rk_qfunc;
+
+    const double dt = 0.1;
+    double t = 0.0;
+    cout.precision(16);
+    for( int i=0 ; i<100 ; ++i , t += dt)
+    {
+        euler.do_step( make_pair( &pfunc , &qfunc ) , state1 , t , dt );
+        euler_qfunc.do_step( &qfunc , state2 , t , dt );
+
+        rk.do_step( make_pair( &pfunc , &qfunc ) , state3 , t , dt );
+        rk_qfunc.do_step( &qfunc , state4 , t , dt );
+
+        cout << t << tab << energy( state1 ) - initial_energy << tab;
+        cout << energy( state2 ) - initial_energy << tab;
+        cout << energy( state3 ) - initial_energy << tab;
+        cout << energy( state4 ) - initial_energy << endl;
+    }
+    clog << "final energy : " << energy(state1) << tab;
+    clog << energy( state2 ) << tab;
+    clog << energy( state3 ) << tab;
+    clog << energy( state4 ) << endl;
+
+}