|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r57410 - in sandbox/odeint: boost/numeric/odeint boost/numeric/odeint/detail libs/numeric/odeint/examples
From: karsten.ahnert_at_[hidden]
Date: 2009-11-05 11:56:07
Author: karsten
Date: 2009-11-05 11:56:06 EST (Thu, 05 Nov 2009)
New Revision: 57410
URL: http://svn.boost.org/trac/boost/changeset/57410
Log:
cosmetics in rk4
Text files modified:
sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp | 43 ++++++++++++++++
sandbox/odeint/boost/numeric/odeint/euler.hpp | 3
sandbox/odeint/boost/numeric/odeint/runge_kutta_4.hpp | 104 +++++++++++++++++++++++++---------------
sandbox/odeint/libs/numeric/odeint/examples/lorenz_integrate_constant_step.cpp | 10 +-
4 files changed, 115 insertions(+), 45 deletions(-)
Modified: sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp 2009-11-05 11:56:06 EST (Thu, 05 Nov 2009)
@@ -55,6 +55,49 @@
( *out_start++ ) = ( *x_start++ ) - ( *y_start++ );
}
+
+ /* computes out_n = x_n + alpha * y_n for all n
+ where x,y and out are given as iterators.
+ make sure, that x,y and out have the same range
+ */
+ template< class InputIterator1 ,
+ class InputIterator2 ,
+ class OutputIterator ,
+ class T >
+ void scale_and_add_and_assign( InputIterator1 x_start ,
+ InputIterator1 x_end ,
+ InputIterator2 y_start ,
+ OutputIterator out_start ,
+ T alpha )
+ {
+ while( x_start != x_end )
+ ( *out_start++) = ( *x_start++ ) + alpha * ( *y_start++ );
+ }
+
+
+
+
+ /* computes out_n = alpha2* ( x_n + y_n + alpha1*z_n )for all n
+ where x,y,z and out are given as iterators.
+ make sure, that x,y,z and out have the same range
+ */
+ template< class InputIterator1 ,
+ class InputIterator2 ,
+ class InputIterator3 ,
+ class OutputIterator ,
+ class T >
+ void scale_and_add_and_add_and_assign( InputIterator1 x_start ,
+ InputIterator1 x_end ,
+ InputIterator2 y_start ,
+ InputIterator3 z_start ,
+ OutputIterator out_start ,
+ T alpha1 ,
+ T alpha2 )
+ {
+ while( x_start != x_end )
+ (*out_start++) += alpha2 * ( (*x_start++) + (*y_start++) + alpha1*(*z_start++) );
+ }
+
}
}
}
Modified: sandbox/odeint/boost/numeric/odeint/euler.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/euler.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/euler.hpp 2009-11-05 11:56:06 EST (Thu, 05 Nov 2009)
@@ -53,7 +53,8 @@
// check the concept of the ContainerType
private:
- BOOST_CLASS_REQUIRE( container_type , boost::numeric::odeint, StateType );
+ BOOST_CLASS_REQUIRE( container_type ,
+ boost::numeric::odeint, StateType );
Modified: sandbox/odeint/boost/numeric/odeint/runge_kutta_4.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/runge_kutta_4.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/runge_kutta_4.hpp 2009-11-05 11:56:06 EST (Thu, 05 Nov 2009)
@@ -21,6 +21,7 @@
#include <boost/numeric/odeint/concepts/state_concept.hpp>
#include <boost/numeric/odeint/resizer.hpp>
+#include <iostream>
namespace boost {
namespace numeric {
@@ -28,73 +29,98 @@
template<
- class ContainerType ,
- class ResizeType = resizer< ContainerType >
+ class Container ,
+ class Time = double ,
+ class Resizer = resizer< Container >
>
class ode_step_runge_kutta_4
{
- BOOST_CLASS_REQUIRE( ContainerType , boost::numeric::odeint, StateType );
- ContainerType dxdt;
- ContainerType dxt;
- ContainerType dxm;
- ContainerType xt;
+ // provide basic typedefs
+ public:
+
+ typedef Container container_type;
+ typedef Resizer resizer_type;
+ typedef Time time_type;
+ typedef typename container_type::value_type value_type;
+ typedef typename container_type::iterator iterator;
+
+
+
+
+
+ // check the concept of the ContainerType
+ private:
+
+ BOOST_CLASS_REQUIRE( container_type ,
+ boost::numeric::odeint, StateType );
+
+
- ResizeType resizer;
- typedef typename ContainerType::iterator iterator;
- typedef typename ContainerType::value_type value_type;
+ // private members
+ private:
+
+ container_type m_dxdt;
+ container_type m_dxt;
+ container_type m_dxm;
+ container_type m_xt;
+ resizer_type m_resizer;
+
+
+
+
+ // public interface
public:
- template< class DynamicalSystem , class TimeType>
+ template< class DynamicalSystem >
void next_step( DynamicalSystem system ,
- ContainerType &x ,
- TimeType t ,
- TimeType dt )
+ container_type &x ,
+ time_type t ,
+ time_type dt )
{
- const TimeType val2 = TimeType( 2.0 );
+ using namespace detail::it_algebra;
- 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 );
-
- TimeType dh = TimeType( 0.5 ) * dt;
- TimeType d6 = dt / TimeType( 6.0 );
- TimeType th = t + dh;
+ const time_type val2 = time_type( 2.0 );
+
+ m_resizer.adjust_size( x , m_dxdt );
+ m_resizer.adjust_size( x , m_dxt );
+ m_resizer.adjust_size( x , m_dxm );
+ m_resizer.adjust_size( x , m_xt );
+
+ time_type dh = time_type( 0.5 ) * dt;
+ time_type d6 = dt / time_type( 6.0 );
+ time_type th = t + dh;
iterator iter1 , iter2 ,iter3 , iter4;
- iterator x_end = x.end() , xt_end = xt.end();
+ iterator x_end = x.end() , xt_end = m_xt.end();
- system( x , dxdt , t );
- iter1 = xt.begin() ; iter2 = x.begin() ; iter3 = dxdt.begin();
- while( iter1 != xt_end )
- (*iter1++) = (*iter2++) + dh * (*iter3++);
+ system( x , m_dxdt , t );
+ scale_and_add_and_assign( x.begin() , x.end() , m_dxdt.begin() , m_xt.begin() , dh );
- system( xt , dxt , th );
- iter1 = xt.begin() ; iter2 = x.begin() ; iter3 = dxt.begin();
- while( iter1 != xt_end )
- (*iter1++) = (*iter2++) + dh * (*iter3++);
+ system( m_xt , m_dxt , th );
+ scale_and_add_and_assign( x.begin() , x.end() , m_dxt.begin() , m_xt.begin() , dh );
- system( xt , dxm , th );
- iter1 = xt.begin() ; iter2 = x.begin() ; iter3 = dxm.begin() ; iter4 = dxt.begin();
+ system( m_xt , m_dxm , th );
+ iter1 = m_xt.begin() ; iter2 = x.begin() ; iter3 = m_dxm.begin() ; iter4 = m_dxt.begin();
while( iter1 != xt_end )
{
(*iter1++) = (*iter2++) + dt * (*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++) );
+ system( m_xt , m_dxt , value_type( t + dt ) );
+ scale_and_add_and_add_and_assign( m_dxdt.begin() , m_dxdt.end() , m_dxt.begin() , m_dxm.begin() , x.begin() , val2 , d6 );
+
}
+
+
};
} // namespace odeint
} // namespace numeric
- } // namespace boost
+} // namespace boost
#endif // BOOST_NUMERIC_ODEINT_RUNGE_KUTTA_4_HPP
Modified: sandbox/odeint/libs/numeric/odeint/examples/lorenz_integrate_constant_step.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/lorenz_integrate_constant_step.cpp (original)
+++ sandbox/odeint/libs/numeric/odeint/examples/lorenz_integrate_constant_step.cpp 2009-11-05 11:56:06 EST (Thu, 05 Nov 2009)
@@ -58,20 +58,20 @@
int main( int argc , char **argv )
{
-
state_type x(3);
x[0] = 1.0;
x[1] = 0.0;
x[2] = 0.0;
+ ode_step_runge_kutta_4< state_type , double > rk4;
ode_step_euler< state_type , double > euler;
-// integrate( euler , lorenz , 0.0 , 0.01 , x , 10.0 , my_observer );
- integrate( euler , lorenz , 0.0 , 0.01 , x , 1.0 , cout << _1 << tab << _2[0] << "\n" );
+ integrate( rk4 , lorenz , 0.0 , 0.01 , x , 100.0 , cout << _1 << tab << _2[0] << "\n" );
+// integrate( euler , lorenz , 0.0 , 0.01 , x , 100.0 , cout << _1 << tab << _2[0] << "\n" );
- vector<double> traj;
+/* vector<double> traj;
back_insert_iterator< vector<double> > iter(traj);
integrate( euler , lorenz , 0.0 , 0.01 , x , 1.0 , var(*iter++) = _2[1] );
- copy( traj.begin() , traj.end() , ostream_iterator<double>( cout , "\n" ) );
+ copy( traj.begin() , traj.end() , ostream_iterator<double>( cout , "\n" ) );*/
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