[Boost-commit] svn:boost r62905 - in sandbox/odeint: boost/numeric/odeint libs/numeric/odeint/examples

Subject: [Boost-commit] svn:boost r62905 - in sandbox/odeint: boost/numeric/odeint libs/numeric/odeint/examples
From: mario.mulansky_at_[hidden]
Date: 2010-06-13 10:01:08

Author: mariomulansky
Date: 2010-06-13 10:01:06 EDT (Sun, 13 Jun 2010)
New Revision: 62905
URL: http://svn.boost.org/trac/boost/changeset/62905

Log:
solar system
Added:
   sandbox/odeint/boost/numeric/odeint/hamiltonian_stepper_qfunc.hpp   (contents, props changed)
Text files modified: 
   sandbox/odeint/boost/numeric/odeint/hamiltonian_stepper_euler.hpp |    78 +++++++++++++++--------                 
   sandbox/odeint/libs/numeric/odeint/examples/solar_system.cpp      |   133 +++++++++++++++++++++++---------------- 
   2 files changed, 127 insertions(+), 84 deletions(-)

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-13 10:01:06 EDT (Sun, 13 Jun 2010)
@@ -44,42 +44,64 @@
 
     private:
 
-	container_type m_dpdt;
-
-
+        state_type m_dxdt;
+        container_type m_dqdt , m_dpdt ;
 
     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 )
         {
+            if( !traits_type::same_size( state.first , state.second ) )
+            {
+                std::string msg( "hamiltonian_stepper_euler::do_step(): " );
+                msg += "q and p have different sizes";
+                throw std::invalid_argument( msg );
+            }
             
-            container_type &q = x.first;
-            container_type &p = x.second;
+            container_type &q = state.first , &p = state.second;
+            container_type &dqdt = m_dxdt.first , &dpdt = m_dxdt.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( q , dqdt );
+            traits_type::adjust_size( p , dpdt );
 
+            func.first( p , dqdt );
+            detail::it_algebra::increment( traits_type::begin(q) ,
+                                           traits_type::end(q) ,
+                                           traits_type::begin(dqdt) ,
+                                           dt );
+            func.second( q , dpdt );
+            detail::it_algebra::increment( traits_type::begin(p) ,
+                                           traits_type::end(p) ,
+                                           traits_type::begin(dpdt) ,
+                                           dt );
+
+        }
+
+        template< class CoordinateFunction , class MomentumFunction >
+        void do_step( CoordinateFunction dqdt_func ,
+                      MomentumFunction dpdt_func ,
+                      container_type &q ,
+                      container_type &p ,
+                      time_type dt )
+        {
+            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_dqdt );
             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 );
-	}
-
+            
+            dqdt_func( p , m_dqdt );
+            detail::it_algebra::increment( traits_type::begin(q) ,
+                                           traits_type::end(q) ,
+                                           traits_type::begin(m_dqdt) ,
+                                           dt );
+        }
     };
 
 

Added: sandbox/odeint/boost/numeric/odeint/hamiltonian_stepper_qfunc.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/boost/numeric/odeint/hamiltonian_stepper_qfunc.hpp	2010-06-13 10:01:06 EDT (Sun, 13 Jun 2010)
@@ -0,0 +1,88 @@
+/*
+ boost header: numeric/odeint/hamiltonian_stepper_qfunc_euler.hpp
+
+ Copyright 2009 Karsten Ahnert
+ Copyright 2009 Mario Mulansky
+ Copyright 2009 Andre Bergner
+
+ 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_HAMILTONIAN_STEPPER_QFUNC_EULER_HPP_INCLUDED
+#define BOOST_NUMERIC_ODEINT_HAMILTONIAN_STEPPER_QFUNC_EULER_HPP_INCLUDED
+
+#include <stdexcept>
+
+#include <boost/numeric/odeint/detail/iterator_algebra.hpp>
+#include <boost/numeric/odeint/container_traits.hpp>
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+    template<
+        class Container ,
+        class Time = double ,
+        class Traits = container_traits< Container >
+        >
+    class hamiltonian_stepper_qfunc_euler
+    {
+        // 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 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 ,
+                  container_type &q ,
+                  container_type &p ,
+                  time_type dt )
+        {
+            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
+
+
+
+
+
+#endif //BOOST_NUMERIC_ODEINT_HAMILTONIAN_STEPPER_QFUNC_EULER_HPP_INCLUDED

Modified: sandbox/odeint/libs/numeric/odeint/examples/solar_system.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/solar_system.cpp	(original)
+++ sandbox/odeint/libs/numeric/odeint/examples/solar_system.cpp	2010-06-13 10:01:06 EDT (Sun, 13 Jun 2010)
@@ -16,10 +16,9 @@
 
 #include <boost/circular_buffer.hpp>
 #include <boost/ref.hpp>
-#include <boost/numeric/odeint.hpp>
-#include <boost/numeric/odeint/container_traits_tr1_array.hpp>
-#include <utility>
 
+#include <boost/numeric/odeint/container_traits_tr1_array.hpp>
+#include <boost/numeric/odeint/hamiltonian_stepper_euler.hpp>
 #include "point_type.hpp"
 
 #define tab "\t" 
@@ -32,9 +31,6 @@
 
 typedef point< double , 3 > point_type;
 typedef std::tr1::array< point_type , n > state_type;
-
-typedef std::pair< state_type , state_type > state_pair;
-
 typedef std::tr1::array< double , n > mass_type;
 typedef hamiltonian_stepper_euler< state_type > stepper_type;
 
@@ -57,23 +53,30 @@
     point_type mean( 0.0 );
     for( size_t i=0 ; i<x.size() ; ++i )
     {
-	overall_mass += m[i];
-	mean += m[i] * x[i];
+        overall_mass += m[i];
+        mean += m[i] * x[i];
     }
     if( x.size() != 0 ) mean /= overall_mass;
     return mean;
 }
 
+point_type center( const state_type &x )
+{
+    point_type mean( 0.0 );
+    for( size_t i=0 ; i<x.size() ; ++i ) mean += x[i];
+    if( !x.empty() ) mean /= x.size();
+    return mean;
+}
 
-double energy( const state_type &q , const state_type &p ,
-               const mass_type &masses )
+
+double energy( const state_type &q , const state_type &p , const mass_type &masses )
 {
     const size_t n = q.size();
     double en = 0.0;
     for( size_t i=0 ; i<n ; ++i )
     {
         en += 0.5 * norm( p[i] ) / masses[i];
-	for( size_t j=0 ; j<i ; ++j )
+        for( size_t j=0 ; j<i ; ++j )
         {
             double diff = abs( q[i] - q[j] );
             en -= gravitational_constant * masses[j] * masses[i] / diff;
@@ -84,32 +87,46 @@
 
 struct solar_system
 {
-    mass_type &m_masses;
+    const mass_type &m_masses;
 
-    solar_system( mass_type &masses ) : m_masses( masses ) { }
+    solar_system( const mass_type &masses ) : m_masses( masses ) { }
 
-    void operator()( const state_type &q , state_type &dpdt )
+    void operator()( const state_type &q , state_type &dpdt ) const
     {
         const size_t n = q.size();
-        fill( dpdt.begin() , dpdt.end() , 0.0 );
+        fill( dpdt.begin() , dpdt.end() , point_type( 0.0 , 0.0 , 0.0 ) );
         for( size_t i=0 ; i<n ; ++i )
         {
             for( size_t j=i+1 ; j<n ; ++j )
             {
                 point_type diff = q[j] - q[i];
-		double d = abs( diff );
-                diff = gravitational_constant * diff / d / d / d;
-                dpdt[i] += diff * m_masses[j];
-                dpdt[j] -= diff * m_masses[i];
+                double d = abs( diff );
+                diff *= ( gravitational_constant / d / d / d * m_masses[i] * m_masses[j] ) ;
+                dpdt[i] += diff ;
+                dpdt[j] -= diff ;
             }
         }
     }
 };
 
+struct solar_system_momentum
+{
+    const mass_type &m_masses;
+
+    solar_system_momentum( const mass_type &masses ) : m_masses( masses ) { }
+
+    void operator()( const state_type &p , state_type &dqdt ) const
+    {
+        for( size_t i=0 ; i<p.size() ; ++i ) dqdt[i] = p[i] / m_masses[i];
+    }
+};
+
+
+
 
 int main( int argc , char **argv )
 {
-    mass_type masses = {{
+    const mass_type masses = {{
             1.00000597682 ,      // sun
             0.000954786104043 ,  // jupiter
             0.000285583733151 ,  // saturn
@@ -118,46 +135,50 @@
             1.0 / ( 1.3e8 )      // pluto
         }};
 
-    state_pair state;
-    state_type &q = state.first;
-    state_type &p = state.second;
-
-    q[0] = point_type( 0.0 , 0.0 , 0.0 );                         // sun
-    q[1] = point_type( -3.5023653 , -3.8169847 , -1.5507963 );    // jupiter
-    q[2] = point_type( 9.0755314 , -3.0458353 , -1.6483708 );     // saturn
-    q[3] = point_type( 8.3101420 , -16.2901086 , -7.2521278 );    // uranus
-    q[4] = point_type( 11.4707666 , -25.7294829 , -10.8169456 );  // neptune
-    q[5] = point_type( -15.5387357 , -25.2225594 , -3.1902382 );  // pluto
-    
-
-    p[0] = point_type( 0.0 , 0.0 , 0.0 );                        // sun
-    p[1] = point_type( 0.00565429 , -0.00412490 , -0.00190589 ); // jupiter
-    p[2] = point_type( 0.00168318 , 0.00483525 , 0.00192462 );   // saturn
-    p[3] = point_type( 0.00354178 , 0.00137102 , 0.00055029 );   // uranus
-    p[4] = point_type( 0.00288930 , 0.00114527 , 0.00039677 );   // neptune
-    p[5] = point_type( 0.00276725 , -0.00170702 , -0.00136504 ); // pluto
-
-    point_type qmean = center_of_mass( q , masses );
-    point_type pmean = center_of_mass( p , masses );
-    for( size_t i=0 ; i<n ; ++i ) { q[i] -= qmean ; p[i] -= pmean; }
+    state_type q = {{
+            point_type( 0.0 , 0.0 , 0.0 ) ,                        // sun
+            point_type( -3.5023653 , -3.8169847 , -1.5507963 ) ,   // jupiter
+            point_type( 9.0755314 , -3.0458353 , -1.6483708 ) ,    // saturn
+            point_type( 8.3101420 , -16.2901086 , -7.2521278 ) ,   // uranus
+            point_type( 11.4707666 , -25.7294829 , -10.8169456 ) , // neptune
+            point_type( -15.5387357 , -25.2225594 , -3.1902382 )   // pluto
+        }};
+
+    state_type p = {{
+            point_type( 0.0 , 0.0 , 0.0 ) ,                        // sun
+            point_type( 0.00565429 , -0.00412490 , -0.00190589 ) , // jupiter
+            point_type( 0.00168318 , 0.00483525 , 0.00192462 ) ,   // saturn
+            point_type( 0.00354178 , 0.00137102 , 0.00055029 ) ,   // uranus
+            point_type( 0.00288930 , 0.00114527 , 0.00039677 ) ,   // neptune
+            point_type( 0.00276725 , -0.00170702 , -0.00136504 )   // pluto
+        }};
+    for( size_t i=0 ; i<n ; ++i ) p[i] *= masses[i];
 
-    stepper_type stepper;
-    solar_system sol( masses );
 
-    const double dt = 100.0;
+/*
+ * ToDo : remove center of mass velocity
+*/
+
+    stepper_type stepper;
+    const double dt = 1.0;
     double t = 0.0;
-    while( t < 10000000.0 )
-    {
-	clog << t << tab << energy( q , p , masses ) << tab;
-        clog << center_of_mass( q , masses ) << tab;
-        clog << center_of_mass( p , masses ) << endl;
-
-	cout << t;
-        for( size_t i=0 ; i<n ; ++i ) cout << tab << q[i];
-	cout << endl;
 
-        for( size_t i=0 ; i<1 ; ++i,t+=dt )
-	    stepper.do_step( sol , state , t , dt );
+
+    pair< state_type , state_type > state = make_pair( q , p );
+    for( size_t c = 0 ; c<2000 ; ++c )
+    {
+        clog << t << tab << energy( state.first , state.second , masses ) << tab;
+        clog << center_of_mass( state.first , masses ) << tab;
+        clog << center_of_mass( state.second , masses ) << endl;
+
+        cout << t;
+        for( size_t i=0 ; i<n ; ++i ) cout << tab << state.first[i];
+        cout << endl;
+
+        for( size_t i=0 ; i<100 ; ++i,t+=dt )
+            stepper.do_step( make_pair( solar_system_momentum( masses ) ,
+                                        solar_system( masses ) ) ,
+                             state , t , dt );
         t += dt;
     }