Subject: [Boost-commit] svn:boost r57132 - in sandbox/odeint: boost/numeric boost/numeric/odeint boost/numeric/odeint/concepts libs/numeric/odeint/examples
From: karsten.ahnert_at_[hidden]
Date: 2009-10-24 11:11:04
Author: karsten
Date: 2009-10-24 11:11:03 EDT (Sat, 24 Oct 2009)
New Revision: 57132
URL: http://svn.boost.org/trac/boost/changeset/57132
Log:
added runge kutta 4th order, and small structural changes
Added:
sandbox/odeint/boost/numeric/odeint/concepts/state_concept.hpp (contents, props changed)
- copied, changed from r57126, /sandbox/odeint/boost/numeric/odeint/concepts/concepts.hpp
sandbox/odeint/boost/numeric/odeint/resizer.hpp (contents, props changed)
sandbox/odeint/boost/numeric/odeint/runge_kutta_4.hpp (contents, props changed)
Removed:
sandbox/odeint/boost/numeric/odeint/concepts/concepts.hpp
Text files modified:
sandbox/odeint/boost/numeric/odeint.hpp | 1
sandbox/odeint/boost/numeric/odeint/concepts/state_concept.hpp | 7 +++--
sandbox/odeint/boost/numeric/odeint/euler.hpp | 44 ++-------------------------------------
sandbox/odeint/libs/numeric/odeint/examples/lorenz_array.cpp | 10 ++++++--
4 files changed, 15 insertions(+), 47 deletions(-)
Modified: sandbox/odeint/boost/numeric/odeint.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint.hpp 2009-10-24 11:11:03 EDT (Sat, 24 Oct 2009)
@@ -16,5 +16,6 @@
#define BOOST_NUMERIC_ODEINT_HPP
#include <boost/numeric/odeint/euler.hpp>
+#include <boost/numeric/odeint/runge_kutta_4.hpp>
#endif // BOOST_NUMERIC_ODEINT_HPP
Deleted: sandbox/odeint/boost/numeric/odeint/concepts/concepts.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/concepts/concepts.hpp 2009-10-24 11:11:03 EDT (Sat, 24 Oct 2009)
+++ (empty file)
@@ -1,51 +0,0 @@
-/* Boost odeint/euler.hpp header file
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
- Copyright 2009 Andre Bergner
-
- This file contains the concepts used in the odeint library
-
- 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_CONCEPTS_HPP
-#define BOOST_NUMERIC_ODEINT_CONCEPTS_HPP
-
-#include <boost/concept_check.hpp>
-
-namespace boost {
-namespace numeric {
-namespace odeint {
-
- template<class X>
- struct StateType {
-
- public:
- typedef typename X::iterator iterator; // requires iterator typedef
-
- // requires iterator being ForwardIterator
- BOOST_CONCEPT_ASSERT((ForwardIterator<iterator>));
-
- BOOST_CONCEPT_USAGE(StateType)
- {
- same_type(state.begin(), it); //requires begin() method
- same_type(state.end(), it); // requires end() method
- }
-
- private:
- X state;
- iterator it;
-
- template<class T>
- void same_type( T const&, T const& );
-
- };
-
-} // namespace odeint
-} // namespace numeric
-} // namespace boost
-
-#endif
Copied: sandbox/odeint/boost/numeric/odeint/concepts/state_concept.hpp (from r57126, /sandbox/odeint/boost/numeric/odeint/concepts/concepts.hpp)
==============================================================================
--- /sandbox/odeint/boost/numeric/odeint/concepts/concepts.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/concepts/state_concept.hpp 2009-10-24 11:11:03 EDT (Sat, 24 Oct 2009)
@@ -1,4 +1,4 @@
-/* Boost odeint/euler.hpp header file
+/* Boost odeint/concepts/state_concept.hpp header file
Copyright 2009 Karsten Ahnert
Copyright 2009 Mario Mulansky
@@ -11,8 +11,8 @@
copy at http://www.boost.org/LICENSE_1_0.txt)
*/
-#ifndef BOOST_NUMERIC_ODEINT_CONCEPTS_HPP
-#define BOOST_NUMERIC_ODEINT_CONCEPTS_HPP
+#ifndef BOOST_NUMERIC_ODEINT_CONCEPTS_STATE_CONCEPT_HPP
+#define BOOST_NUMERIC_ODEINT_CONCEPTS_STATE_CONCEPT_HPP
#include <boost/concept_check.hpp>
@@ -36,6 +36,7 @@
}
private:
+
X state;
iterator it;
Modified: sandbox/odeint/boost/numeric/odeint/euler.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/euler.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/euler.hpp 2009-10-24 11:11:03 EDT (Sat, 24 Oct 2009)
@@ -18,44 +18,15 @@
#define BOOST_NUMERIC_ODEINT_EULER_HPP
#include <boost/concept_check.hpp>
-#include <boost/numeric/odeint/concept/concepts.hpp>
-#include <tr1/array>
+
+#include <boost/numeric/odeint/concepts/state_concept.hpp>
+#include <boost/numeric/odeint/resizer.hpp>
namespace boost {
namespace numeric {
namespace odeint {
- template< class ContainerType >
- class resizer
- {
- public:
- void resize( const ContainerType &x , ContainerType &dxdt ) const
- {
- dxdt.resize( x.size() );
- }
-
- bool same_size( const ContainerType &x1 , ContainerType &x2 ) const
- {
- return (x1.size() == x2.size());
- }
- };
-
- template< class T , size_t N >
- class resizer< std::tr1::array< T , N > >
- {
- public:
- void resize( const std::tr1::array<T,N> &x , std::tr1::array<T,N> &dxdt ) const
- {
- throw; // should never be called
- }
-
- const bool same_size( const std::tr1::array<T,N> &x1 , std::tr1::array<T,N> &x2 ) const
- {
- return true; // if this was false, the code wouldn't compile
- }
- };
-
template<
class ContainerType ,
class ResizeType = resizer< ContainerType >
@@ -86,15 +57,6 @@
}
};
-
-/* ToDo:
- Write stepper for
- * fixed size systems
- * array<T>
- * system( T* , T* , T )
-*/
-
-
} // namespace odeint
} // namespace numeric
} // namespace boost
Added: sandbox/odeint/boost/numeric/odeint/resizer.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/boost/numeric/odeint/resizer.hpp 2009-10-24 11:11:03 EDT (Sat, 24 Oct 2009)
@@ -0,0 +1,60 @@
+/* Boost odeint/resizer.hpp header file
+
+ Copyright 2009 Karsten Ahnert
+ Copyright 2009 Mario Mulansky
+
+ This file includes resizer functionality for containers
+
+ 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_RESIZER_HPP
+#define BOOST_NUMERIC_ODEINT_RESIZER_HPP
+
+#include <tr1/array>
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+ template< class ContainerType >
+ class resizer
+ {
+ public:
+ void resize( const ContainerType &x , ContainerType &dxdt ) const
+ {
+ dxdt.resize( x.size() );
+ }
+
+ bool same_size( const ContainerType &x1 , ContainerType &x2 ) const
+ {
+ return (x1.size() == x2.size());
+ }
+ };
+
+ template< class T , size_t N >
+ class resizer< std::tr1::array< T , N > >
+ {
+ public:
+ void resize( const std::tr1::array<T,N> &x ,
+ std::tr1::array<T,N> &dxdt ) const
+ {
+ throw; // should never be called
+ }
+
+ const bool same_size( const std::tr1::array<T,N> &x1 ,
+ std::tr1::array<T,N> &x2 ) const
+ {
+ return true; // if this was false, the code wouldn't compile
+ }
+ };
+
+
+} // namespace odeint
+} // namespace numeric
+} // namespace boost
+
+
+#endif // BOOST_NUMERIC_ODEINT_RESIZER_HPP
Added: sandbox/odeint/boost/numeric/odeint/runge_kutta_4.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/boost/numeric/odeint/runge_kutta_4.hpp 2009-10-24 11:11:03 EDT (Sat, 24 Oct 2009)
@@ -0,0 +1,101 @@
+/* Boost odeint/runge_kutta_4.hpp header file
+
+ Copyright 2009 Karsten Ahnert
+ Copyright 2009 Mario Mulansky
+ Copyright 2009 Andre Bergner
+
+ This file includes the explicit runge kutta solver for
+ ordinary differential equations.
+
+ It solves any ODE dx/dt = f(x,t).
+
+ 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_RUNGE_KUTTA_4_HPP
+#define BOOST_NUMERIC_ODEINT_RUNGE_KUTTA_4_HPP
+
+#include <boost/concept_check.hpp>
+
+#include <boost/numeric/odeint/concepts/state_concept.hpp>
+#include <boost/numeric/odeint/resizer.hpp>
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+
+ template<
+ class ContainerType ,
+ class ResizeType = resizer< ContainerType >
+ >
+ class ode_step_runge_kutta_4
+ {
+ BOOST_CLASS_REQUIRE( ContainerType , boost::numeric::odeint, StateType );
+
+ ContainerType dxdt;
+ ContainerType dxt;
+ ContainerType dxm;
+ ContainerType xt;
+
+ ResizeType resizer;
+
+ typedef typename ContainerType::iterator iterator;
+ typedef typename ContainerType::value_type value_type;
+
+ public:
+
+ template< class DynamicalSystem , class TimeType>
+ void next_step( DynamicalSystem system ,
+ ContainerType &x ,
+ TimeType t ,
+ TimeType dt )
+ {
+ const value_type val2 = value_type( 2.0 );
+
+ if( ! resizer.same_size( x , dxdt ) ) resizer.resize( x , dxdt );
+ if( ! resizer.same_size( x , dxt ) ) resizer.resize( x , dxt );
+ if( ! resizer.same_size( x , dxm ) ) resizer.resize( x , dxm );
+ if( ! resizer.same_size( x , xt ) ) resizer.resize( x , xt );
+
+ value_type dt_val = value_type( dt );
+ value_type dh = dt_val * value_type( 0.5 );
+ value_type d6 = dt_val / value_type( 6.0 );
+ value_type th = dh + value_type( t );
+
+ iterator iter1 , iter2 ,iter3 , iter4;
+ iterator x_end = x.end() , xt_end = xt.end();
+
+ system( x , dxdt , t );
+ iter1 = xt.begin() ; iter2 = x.begin() ; iter3 = dxdt.begin();
+ while( iter1 != xt_end )
+ (*iter1++) = (*iter2++) + dh * (*iter3++);
+
+ system( xt , dxt , th );
+ iter1 = xt.begin() ; iter2 = x.begin() ; iter3 = dxt.begin();
+ while( iter1 != xt_end )
+ (*iter1++) = (*iter2++) + dh * (*iter3++);
+
+ system( xt , dxm , th );
+ iter1 = xt.begin() ; iter2 = x.begin() ; iter3 = dxm.begin() ; iter4 = dxt.begin();
+ while( iter1 != xt_end )
+ {
+ (*iter1++) = (*iter2++) + dt_val * (*iter3);
+ (*iter3++) += (*iter4++);
+ }
+
+ system( xt , dxt , value_type( t + dt ) );
+ iter1 = x.begin() ; iter2 = dxdt.begin() ; iter3 = dxt.begin() ; iter4 = dxm.begin();
+ while( iter1 != x_end )
+ (*iter1++) += d6 * ( (*iter2++) + (*iter3++) + val2 * (*iter4++) );
+ }
+ };
+
+} // namespace odeint
+} // namespace numeric
+} // namespace boost
+
+
+#endif // BOOST_NUMERIC_ODEINT_RUNGE_KUTTA_4_HPP
Modified: sandbox/odeint/libs/numeric/odeint/examples/lorenz_array.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/lorenz_array.cpp (original)
+++ sandbox/odeint/libs/numeric/odeint/examples/lorenz_array.cpp 2009-10-24 11:11:03 EDT (Sat, 24 Oct 2009)
@@ -47,18 +47,22 @@
int main( int argc , char **argv )
{
- state_type x;
+ state_type x , x2;
x[0] = 1.0;
x[1] = 0.0;
- x[2] = 0.0;
+ x[2] = 20.0;
+ x2 = x;
ode_step_euler< state_type > euler;
+ ode_step_runge_kutta_4< state_type > rk4;
double t = 0.0;
for( size_t oi=0 ; oi<olen ; ++oi,t+=dt )
{
- cout << t << tab << x[0] << tab << x[1] << tab << x[2] << endl;
+ cout << t << tab << x[0] << tab << x[1] << tab << x[2] << tab;
+ cout << x2[0] << tab << x2[1] << tab << x2[2] << endl;
euler.next_step( lorenz , x , t , dt );
+ rk4.next_step( lorenz , x2 , t , dt );
}
return 0;