|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r58269 - in sandbox/odeint/branches/karsten: boost/numeric boost/numeric/odeint libs/numeric/odeint/examples
From: karsten.ahnert_at_[hidden]
Date: 2009-12-10 12:06:59
Author: karsten
Date: 2009-12-10 12:06:58 EST (Thu, 10 Dec 2009)
New Revision: 58269
URL: http://svn.boost.org/trac/boost/changeset/58269
Log:
added two hamiltonian solvers
Added:
sandbox/odeint/branches/karsten/boost/numeric/odeint/hamiltonian_stepper_euler.hpp (contents, props changed)
sandbox/odeint/branches/karsten/boost/numeric/odeint/hamiltonian_stepper_rk.hpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/hamiltonian_test.cpp (contents, props changed)
Text files modified:
sandbox/odeint/branches/karsten/boost/numeric/odeint.hpp | 3 +++
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_euler.hpp | 8 ++++----
sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/Jamfile | 4 +++-
3 files changed, 10 insertions(+), 5 deletions(-)
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint.hpp 2009-12-10 12:06:58 EST (Thu, 10 Dec 2009)
@@ -29,6 +29,9 @@
#include <boost/numeric/odeint/stepper_midpoint.hpp>
#include <boost/numeric/odeint/stepper_rk78_fehlberg.hpp>
+#include <boost/numeric/odeint/hamiltonian_stepper_euler.hpp>
+#include <boost/numeric/odeint/hamiltonian_stepper_rk.hpp>
+
#include <boost/numeric/odeint/controlled_stepper_standard.hpp>
#include <boost/numeric/odeint/controlled_stepper_bs.hpp>
Added: sandbox/odeint/branches/karsten/boost/numeric/odeint/hamiltonian_stepper_euler.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/hamiltonian_stepper_euler.hpp 2009-12-10 12:06:58 EST (Thu, 10 Dec 2009)
@@ -0,0 +1,88 @@
+/*
+ boost header: numeric/odeint/hamiltonian_stepper_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_EULER_HPP_INCLUDED
+#define BOOST_NUMERIC_ODEINT_HAMILTONIAN_STEPPER_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_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_EULER_HPP_INCLUDED
Added: sandbox/odeint/branches/karsten/boost/numeric/odeint/hamiltonian_stepper_rk.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/hamiltonian_stepper_rk.hpp 2009-12-10 12:06:58 EST (Thu, 10 Dec 2009)
@@ -0,0 +1,126 @@
+/*
+ boost header: numeric/odeint/hamiltonian_stepper_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_RK_HPP_INCLUDED
+#define BOOST_NUMERIC_ODEINT_HAMILTONIAN_STEPPER_RK_HPP_INCLUDED
+
+#include <stdexcept>
+#include <tr1/array>
+
+#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_rk
+ {
+ // 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;
+
+/*
+ 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 ,
+ container_type &q ,
+ container_type &p ,
+ time_type dt )
+ {
+ const size_t order = 6;
+ const std::tr1::array< time_type , order > rk_a = {{
+ 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 std::tr1::array< time_type , order > rk_b = {{
+ 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 )
+ }};
+
+
+
+ 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 );
+
+ 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 );
+ }
+ }
+
+ };
+
+
+} // namespace odeint
+} // namespace numeric
+} // namespace boost
+
+
+
+
+
+#endif //BOOST_NUMERIC_ODEINT_HAMILTONIAN_STEPPER_RK_HPP_INCLUDED
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_euler.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_euler.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_euler.hpp 2009-12-10 12:06:58 EST (Thu, 10 Dec 2009)
@@ -62,10 +62,10 @@
template< class DynamicalSystem >
void do_step( DynamicalSystem &system ,
- container_type &x ,
- const container_type &dxdt ,
- time_type t ,
- time_type dt )
+ container_type &x ,
+ const container_type &dxdt ,
+ time_type t ,
+ time_type dt )
{
//x = x + dt*dxdt
detail::it_algebra::increment( traits_type::begin(x) ,
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/Jamfile
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/Jamfile (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/Jamfile 2009-12-10 12:06:58 EST (Thu, 10 Dec 2009)
@@ -15,7 +15,9 @@
;
exe harmonic_oscillator : harmonic_oscillator.cpp ;
-exe test_container_and_stepper : test_container_and_stepper.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 ;
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/hamiltonian_test.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/hamiltonian_test.cpp 2009-12-10 12:06:58 EST (Thu, 10 Dec 2009)
@@ -0,0 +1,67 @@
+/* Boost numeric/odeint/examples/lorenz_stepper.cpp
+
+ Copyright 2009 Karsten Ahnert
+ Copyright 2009 Mario Mulansky
+
+ Testing of the Hamiltonian solvers
+
+ Furthermore, the usage of std::tr1::array and std::vector in odeint is
+ shown and the performance of both containers is compared.
+
+ 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)
+*/
+
+#include <iostream>
+#include <vector>
+#include <list>
+#include <tr1/array>
+
+#include <boost/numeric/odeint.hpp>
+
+#define tab "\t"
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+
+typedef std::tr1::array< double , 1 > state_type;
+
+void harm_osc( state_type& q , state_type &dpdt )
+{
+ dpdt[0] = - q[0];
+}
+
+
+
+
+int main( int argc , char **argv )
+{
+ const double dt = 0.25;
+ const size_t olen = 10000;
+
+ state_type q1 = {{ 1.0 }} , p1 = {{ 0.0 }};
+ state_type q2 = {{ 1.0 }} , p2 = {{ 0.0 }};
+
+ hamiltonian_stepper_euler< state_type > euler;
+ hamiltonian_stepper_rk< state_type > rk;
+
+ double t = 0.0;
+ for( size_t oi=0 ; oi<olen ; ++oi,t+=dt )
+ {
+ euler.do_step( harm_osc , q1 , p1 , dt );
+ rk.do_step( harm_osc , q2 , p2 , dt );
+ cout << t << tab << q1[0] << tab << p1[0] << tab;
+ cout << q2[0] << tab << p2[0] << "\n";
+ }
+
+ return 0;
+}
+
+
+
+/*
+ Compile with
+ g++ -Wall -O3 -I$BOOST_ROOT -I../../../../ lorenz_stepper.cpp
+*/
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