|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r61746 - in sandbox/odeint: boost/numeric boost/numeric/odeint libs/numeric/odeint/examples libs/numeric/odeint/performance
From: karsten.ahnert_at_[hidden]
Date: 2010-05-03 15:05:04
Author: karsten
Date: 2010-05-03 15:05:03 EDT (Mon, 03 May 2010)
New Revision: 61746
URL: http://svn.boost.org/trac/boost/changeset/61746
Log:
solar system example
Binary files modified:
sandbox/odeint/boost/numeric/odeint/stepper_rk5_ck.hpp
Text files modified:
sandbox/odeint/boost/numeric/odeint.hpp | 2
sandbox/odeint/boost/numeric/odeint/hamiltonian_stepper_euler.hpp | 2
sandbox/odeint/boost/numeric/odeint/hamiltonian_stepper_rk.hpp | 2
sandbox/odeint/libs/numeric/odeint/examples/Jamfile | 2
sandbox/odeint/libs/numeric/odeint/examples/point_type.hpp | 11 ---------
sandbox/odeint/libs/numeric/odeint/examples/solar_system.cpp | 45 ++++++++++++++++++++++++----------------
sandbox/odeint/libs/numeric/odeint/performance/gsl_compare_lorenz.cpp | 4 +-
7 files changed, 33 insertions(+), 35 deletions(-)
Modified: sandbox/odeint/boost/numeric/odeint.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint.hpp 2010-05-03 15:05:03 EDT (Mon, 03 May 2010)
@@ -15,7 +15,7 @@
#ifndef BOOST_NUMERIC_ODEINT_HPP
#define BOOST_NUMERIC_ODEINT_HPP
-#include <boost/config.hpp>
+// #include <boost/config.hpp>
#include <boost/numeric/odeint/container_traits.hpp>
#include <boost/numeric/odeint/container_traits_tr1_array.hpp>
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-05-03 15:05:03 EDT (Mon, 03 May 2010)
@@ -49,7 +49,7 @@
public:
template< class CoordinateFunction >
- void do_step( CoordinateFunction &qfunc ,
+ void do_step( CoordinateFunction qfunc ,
container_type &q ,
container_type &p ,
time_type dt )
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-05-03 15:05:03 EDT (Mon, 03 May 2010)
@@ -63,7 +63,7 @@
public:
template< class CoordinateFunction >
- void do_step( CoordinateFunction &qfunc ,
+ void do_step( CoordinateFunction qfunc ,
container_type &q ,
container_type &p ,
time_type dt )
Modified: sandbox/odeint/boost/numeric/odeint/stepper_rk5_ck.hpp
==============================================================================
Binary files. No diff available.
Modified: sandbox/odeint/libs/numeric/odeint/examples/Jamfile
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/Jamfile (original)
+++ sandbox/odeint/libs/numeric/odeint/examples/Jamfile 2010-05-03 15:05:03 EDT (Mon, 03 May 2010)
@@ -6,7 +6,7 @@
project
: requirements
- <include>../../../..
+ <include>../../../..
<include>$BOOST_ROOT
<include>$BLITZ_ROOT
<include>$MTL4_ROOT
Modified: sandbox/odeint/libs/numeric/odeint/examples/point_type.hpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/point_type.hpp (original)
+++ sandbox/odeint/libs/numeric/odeint/examples/point_type.hpp 2010-05-03 15:05:03 EDT (Mon, 03 May 2010)
@@ -49,17 +49,6 @@
if( dim > 2 ) m_val[2] = z;
}
-/* template< size_t i > T get( void ) const
- {
- BOOST_STATIC_ASSERT( i < dim );
- return m_val[i];
- }
- template< size_t i > T& get( void )
- {
- BOOST_STATIC_ASSERT( i < dim );
- return m_val[i];
- }*/
-
T operator[]( size_t i ) const { return m_val[i]; }
T& operator[]( size_t i ) { return m_val[i]; }
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-05-03 15:05:03 EDT (Mon, 03 May 2010)
@@ -15,6 +15,7 @@
#include <bits/stdtr1c++.h>
#include <boost/circular_buffer.hpp>
+#include <boost/ref.hpp>
#include <boost/numeric/odeint.hpp>
#include "point_type.hpp"
@@ -30,7 +31,7 @@
typedef point< double , 3 > point_type;
typedef std::tr1::array< point_type , n > state_type;
typedef std::tr1::array< double , n > mass_type;
-typedef hamiltonian_stepper_rk< state_type > stepper_type;
+typedef hamiltonian_stepper_euler< state_type > stepper_type;
const double gravitational_constant = 2.95912208286e-4;
@@ -63,26 +64,26 @@
const mass_type &masses )
{
const size_t n = q.size();
- double e = 0.0;
+ double en = 0.0;
for( size_t i=0 ; i<n ; ++i )
{
- e += 0.5 * norm( p[i] );
- for( size_t j=i+1 ; j<n ; ++j )
+ en += 0.5 * norm( p[i] ) / masses[i];
+ for( size_t j=0 ; j<i ; ++j )
{
- double diff = abs( q[j] - q[i] );
- e += gravitational_constant * masses[j] / diff;
+ double diff = abs( q[i] - q[j] );
+ en -= gravitational_constant * masses[j] * masses[i] / diff;
}
}
- return e;
+ return en;
}
struct solar_system
{
- const mass_type &m_masses;
+ mass_type &m_masses;
- solar_system( const mass_type &masses ) : m_masses( masses ) { }
+ solar_system( mass_type &masses ) : m_masses( masses ) { }
- void operator()( state_type &q , state_type &dpdt )
+ void operator()( const state_type &q , state_type &dpdt )
{
const size_t n = q.size();
fill( dpdt.begin() , dpdt.end() , 0.0 );
@@ -91,7 +92,8 @@
for( size_t j=i+1 ; j<n ; ++j )
{
point_type diff = q[j] - q[i];
- diff = gravitational_constant * diff / pow( abs( diff ) , 3.0 );
+ double d = abs( diff );
+ diff = gravitational_constant * diff / d / d / d;
dpdt[i] += diff * m_masses[j];
dpdt[j] -= diff * m_masses[i];
}
@@ -129,28 +131,35 @@
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; }
stepper_type stepper;
- solar_system sol( masses );
-
- const double dt = 0.001;
+ const double dt = 100.0;
double t = 0.0;
- while( t < 100000.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;
- for( size_t i=0 ; i<n ; ++i )
- cout << t << tab << q[i] << tab << p[i] << endl;
+ cout << t;
+ for( size_t i=0 ; i<n ; ++i ) cout << tab << q[i];
+ cout << endl;
- for( size_t i=0 ; i<1000 ; ++i ) stepper.do_step( sol , q , p , dt );
+ for( size_t i=0 ; i<1 ; ++i,t+=dt )
+ stepper.do_step( solar_system( masses ) , q , p , dt );
t += dt;
}
return 0;
}
+
+
+/*
+Plot with gnuplot:
+p "solar_system.dat" u 2:4 w l,"solar_system.dat" u 5:7 w l,"solar_system.dat" u 8:10 w l,"solar_system.dat" u 11:13 w l,"solar_system.dat" u 14:16 w l,"solar_system.dat" u 17:19 w l
+*/
Modified: sandbox/odeint/libs/numeric/odeint/performance/gsl_compare_lorenz.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/performance/gsl_compare_lorenz.cpp (original)
+++ sandbox/odeint/libs/numeric/odeint/performance/gsl_compare_lorenz.cpp 2010-05-03 15:05:03 EDT (Mon, 03 May 2010)
@@ -95,7 +95,7 @@
int main( int argc , char **argv )
{
const double dt = 0.01;
- const size_t olen = 10000000;
+ const size_t olen = 100000000;
state_type x_start = {{ -13.90595 , -15.54127 , 32.48918 }};
@@ -156,7 +156,7 @@
// just check, if rk4 is working correctly
- const size_t tslen = 10000;
+ const size_t tslen = 1000;
x1 = x_start;
copy( x_start.begin() , x_start.end() , x2 );
double x3[3] = { x_start[0] , x_start[1] , x_start[2] };
Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk