|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r57251 - in sandbox/odeint: boost/numeric boost/numeric/odeint boost/numeric/odeint/concepts boost/numeric/odeint/detail libs/numeric/odeint/examples
From: mario.mulansky_at_[hidden]
Date: 2009-10-30 13:56:40
Author: mariomulansky
Date: 2009-10-30 13:56:39 EDT (Fri, 30 Oct 2009)
New Revision: 57251
URL: http://svn.boost.org/trac/boost/changeset/57251
Log:
+stepsize_controller_standard
+lorenz_controlled example
Added:
sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp
- copied, changed from r57244, /sandbox/odeint/boost/numeric/odeint/detail/accumulators.hpp
sandbox/odeint/libs/numeric/odeint/examples/lorenz_controlled.cpp (contents, props changed)
Removed:
sandbox/odeint/boost/numeric/odeint/detail/accumulators.hpp
Text files modified:
sandbox/odeint/boost/numeric/odeint.hpp | 1
sandbox/odeint/boost/numeric/odeint/concepts/state_concept.hpp | 2
sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp | 41 +++++++++++++-------
sandbox/odeint/boost/numeric/odeint/euler.hpp | 35 ++++++++++++++---
sandbox/odeint/boost/numeric/odeint/resizer.hpp | 4 +-
sandbox/odeint/boost/numeric/odeint/stepsize_controller_standard.hpp | 79 +++++++++++++++++++++++++++++++++++++--
6 files changed, 134 insertions(+), 28 deletions(-)
Modified: sandbox/odeint/boost/numeric/odeint.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint.hpp 2009-10-30 13:56:39 EDT (Fri, 30 Oct 2009)
@@ -17,5 +17,6 @@
#include <boost/numeric/odeint/euler.hpp>
#include <boost/numeric/odeint/runge_kutta_4.hpp>
+#include <boost/numeric/odeint/stepsize_controller_standard.hpp>
#endif // BOOST_NUMERIC_ODEINT_HPP
Modified: sandbox/odeint/boost/numeric/odeint/concepts/state_concept.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/concepts/state_concept.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/concepts/state_concept.hpp 2009-10-30 13:56:39 EDT (Fri, 30 Oct 2009)
@@ -53,7 +53,7 @@
public:
- //BOOST_CONCEPT_ASSERT((StateType<X>));
+ BOOST_CONCEPT_ASSERT((StateType<X>));
BOOST_CONCEPT_USAGE(Resizable)
{
Deleted: sandbox/odeint/boost/numeric/odeint/detail/accumulators.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/detail/accumulators.hpp 2009-10-30 13:56:39 EDT (Fri, 30 Oct 2009)
+++ (empty file)
@@ -1,54 +0,0 @@
-/* Boost odeint/detail/accumulators.hpp header file
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
- Copyright 2009 Andre Bergner
-
- Some accumulators for odeint
-
- 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_DETAIL_ACCUMULATORS_HPP
-#define BOOST_NUMERIC_ODEINT_DETAIL_ACCUMULATORS_HPP
-
-
-
-namespace boost {
-namespace numeric {
-namespace odeint {
-namespace detail {
-
- template< class InputIterator ,
- class OutputIterator ,
- class T >
- void multiply_and_add( InputIterator first1 ,
- InputIterator last1 ,
- OutputIterator first2 ,
- T dt )
- {
- while( first1 != last1 )
- (*first1++) += dt * (*first2++);
- }
-
- template< class InputIterator1 ,
- class InputIterator2 ,
- class OutputIterator >
- void substract_and_assign( InputIterator1 first1 ,
- InputIterator1 last1 ,
- InputIterator2 first2 ,
- OutputIterator first3 )
- {
- while( first1 != last1 )
- ( *first3++ ) = ( *first1++ ) - ( *first2++ );
- }
-
-}
-}
-}
-}
-
-
-#endif //BOOST_NUMERIC_ODEINT_DETAIL_ACCUMULATORS_HPP
Copied: sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp (from r57244, /sandbox/odeint/boost/numeric/odeint/detail/accumulators.hpp)
==============================================================================
--- /sandbox/odeint/boost/numeric/odeint/detail/accumulators.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp 2009-10-30 13:56:39 EDT (Fri, 30 Oct 2009)
@@ -4,7 +4,7 @@
Copyright 2009 Mario Mulansky
Copyright 2009 Andre Bergner
- Some accumulators for odeint
+ Some algebraic operations for iterators
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or
@@ -20,35 +20,46 @@
namespace numeric {
namespace odeint {
namespace detail {
+namespace it_algebra { // iterator algebra
- template< class InputIterator ,
- class OutputIterator ,
+ /* computes x_n += alpha * y_n for all n
+ where x and y are given by the iterators x_start, x_end, y_start
+ make sure that x and y have the same range
+ */
+ template< class Iterator1 ,
+ class Iterator2 ,
class T >
- void multiply_and_add( InputIterator first1 ,
- InputIterator last1 ,
- OutputIterator first2 ,
- T dt )
+ void scale_and_add( Iterator1 x_start ,
+ Iterator1 x_end ,
+ Iterator2 y_start ,
+ T alpha )
{
- while( first1 != last1 )
- (*first1++) += dt * (*first2++);
+ while( x_start != x_end )
+ (*x_start++) += alpha * (*y_start++);
}
+
+ /* computes out_n = x_n - 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 >
- void substract_and_assign( InputIterator1 first1 ,
- InputIterator1 last1 ,
- InputIterator2 first2 ,
- OutputIterator first3 )
+ void substract_and_assign( InputIterator1 x_start ,
+ InputIterator1 x_end ,
+ InputIterator2 y_start ,
+ OutputIterator out_start )
{
- while( first1 != last1 )
- ( *first3++ ) = ( *first1++ ) - ( *first2++ );
+ while( x_start != x_end )
+ ( *out_start++ ) = ( *x_start++ ) - ( *y_start++ );
}
}
}
}
}
+}
#endif //BOOST_NUMERIC_ODEINT_DETAIL_ACCUMULATORS_HPP
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-30 13:56:39 EDT (Fri, 30 Oct 2009)
@@ -19,7 +19,7 @@
#include <boost/concept_check.hpp>
-#include <boost/numeric/odeint/detail/accumulators.hpp>
+#include <boost/numeric/odeint/detail/iterator_algebra.hpp>
#include <boost/numeric/odeint/concepts/state_concept.hpp>
#include <boost/numeric/odeint/resizer.hpp>
@@ -47,6 +47,8 @@
public:
+ const unsigned int order() { return 1; }
+
template< class DynamicalSystem , class TimeType >
void next_step( DynamicalSystem system ,
ContainerType &x ,
@@ -54,7 +56,7 @@
TimeType t ,
TimeType dt )
{
- detail::multiply_and_add( x.begin() , x.end() , dxdt.begin() , dt );
+ detail::it_algebra::scale_and_add( x.begin() , x.end() , dxdt.begin() , dt );
}
template< class DynamicalSystem , class TimeType >
@@ -63,7 +65,7 @@
TimeType t ,
TimeType dt )
{
- resizer.check_size_and_resize( x , dxdt );
+ resizer.adjust_size( x , dxdt );
system( x , dxdt , t );
next_step( system , x , dxdt , t , dt );
}
@@ -76,8 +78,8 @@
TimeType dt ,
ContainerType &xerr )
{
- resizer.check_size_and_resize( x , dxdt );
- resizer.check_size_and_resize( x , xerr );
+ resizer.adjust_size( x , dxdt );
+ resizer.adjust_size( x , xerr );
xtemp = x;
TimeType dt2 = 0.5*dt;
@@ -88,9 +90,30 @@
next_step( system , xtemp , dxdt , t , dt2 );
next_step( system , xtemp , t+dt2 , dt2 );
- detail::substract_and_assign( x.begin() , x.end() , xtemp.begin() , xerr.begin() );
+ detail::it_algebra::substract_and_assign( x.begin() , x.end() , xtemp.begin() , xerr.begin() );
+ }
+
+ template< class DynamicalSystem , class TimeType >
+ void next_step( DynamicalSystem system ,
+ ContainerType &x ,
+ ContainerType &dxdt ,
+ TimeType t ,
+ TimeType dt ,
+ ContainerType &xerr )
+ {
+ resizer.adjust_size( x , xerr );
+
+ xtemp = x;
+ TimeType dt2 = 0.5*dt;
+
+ next_step( system , x , dxdt , t , dt );
+ next_step( system , xtemp , dxdt , t , dt2 );
+ next_step( system , xtemp , t+dt2 , dt2 );
+
+ detail::it_algebra::substract_and_assign( x.begin() , x.end() , xtemp.begin() , xerr.begin() );
}
+
};
} // namespace odeint
Modified: sandbox/odeint/boost/numeric/odeint/resizer.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/resizer.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/resizer.hpp 2009-10-30 13:56:39 EDT (Fri, 30 Oct 2009)
@@ -36,7 +36,7 @@
return (x1.size() == x2.size());
}
- void check_size_and_resize( const ContainerType &x1 , ContainerType &x2 ) const
+ void adjust_size( const ContainerType &x1 , ContainerType &x2 ) const
{
if( same_size( x1 , x2 ) ) resize( x1 , x2 );
}
@@ -59,7 +59,7 @@
return true; // if this was false, the code wouldn't compile
}
- void check_size_and_resize( const std::tr1::array<T,N> &x1 , std::tr1::array<T,N> &x2 ) const
+ void adjust_size( const std::tr1::array<T,N> &x1 , std::tr1::array<T,N> &x2 ) const
{
if( !same_size( x1 , x2 ) ) throw;
}
Modified: sandbox/odeint/boost/numeric/odeint/stepsize_controller_standard.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepsize_controller_standard.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/stepsize_controller_standard.hpp 2009-10-30 13:56:39 EDT (Fri, 30 Oct 2009)
@@ -1,13 +1,10 @@
-/* Boost odeint/stepsite_controller_standard.hpp header file
+/* Boost odeint/stepsize_controller_standard.hpp header file
Copyright 2009 Karsten Ahnert
Copyright 2009 Mario Mulansky
This file includes the standard step size controller
- It solves any ODE dx/dt = f(x,t) via
- x(t+dt) = x(t) + 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)
@@ -16,16 +13,90 @@
#ifndef BOOST_NUMERIC_ODEINT_STEPSIZE_CONTROLLER_STANDARD_HPP
#define BOOST_NUMERIC_ODEINT_STEPSIZE_CONTROLLER_STANDARD_HPP
+#include <cmath> // for pow( ) and abs()
+#include <complex>
+
#include <boost/concept_check.hpp>
+#include <boost/numeric/odeint/detail/iterator_algebra.hpp>
#include <boost/numeric/odeint/concepts/state_concept.hpp>
#include <boost/numeric/odeint/resizer.hpp>
+#include <iostream>
+
namespace boost {
namespace numeric {
namespace odeint {
+typedef enum{SUCCESS, STEP_SIZE_DECREASED, STEP_SIZE_INCREASED} controlled_step_result;
+ template<
+ class ContainerType,
+ class T,
+ class ResizeType = resizer< ContainerType > >
+ class step_controller_standard {
+
+ typedef typename ContainerType::iterator iterator;
+
+ T eps_abs;
+ T eps_rel;
+ T a_x;
+ T a_dxdt;
+ ContainerType dxdt;
+ ContainerType x_tmp;
+ ContainerType x_err;
+ ResizeType resizer;
+
+ public:
+ step_controller_standard( T abs_err, T rel_err, T factor_x, T factor_dxdt )
+ : eps_abs(abs_err), eps_rel(rel_err), a_x(factor_x), a_dxdt(factor_dxdt)
+ { }
+
+ template< class Step, class DynamicalSystem >
+ controlled_step_result controlled_step( Step &stepper,
+ DynamicalSystem &system,
+ ContainerType &x,
+ T &t,
+ T &dt )
+ {
+ resizer.adjust_size(x, x_err);
+
+ x_tmp = x; // copy current state
+ system( x, dxdt, t ); // compute dxdt
+ stepper.next_step(system, x, dxdt, t, dt, x_err ); // do step forward with error
+
+ iterator x_start = x_tmp.begin();
+ iterator dxdt_start = dxdt.begin();
+ iterator x_err_start = x_err.begin();
+ T max_rel_err = 0.0;
+
+ while( x_start != x_tmp.end() ) {
+ // get the maximal value of x_err/D where
+ // D = eps_abs + eps_rel * (a_x*|x| + a_dxdt*|dxdt|);
+ T err = eps_abs + eps_rel * (a_x * std::abs(*x_start++) +
+ a_dxdt * dt * std::abs(*dxdt_start++));
+ max_rel_err = max( std::abs(*x_err_start++)/err , max_rel_err );
+ }
+
+ //std::cout<<max_rel_err<<std::endl;
+
+ if( max_rel_err > 1.1 ) { // error too large - decrease dt
+ // limit scaling factor to 0.2
+ dt *= max( 0.9*pow(max_rel_err , -1.0/stepper.order()) , 0.2 );
+ // reset state
+ x = x_tmp;
+ return STEP_SIZE_DECREASED;
+ } else if( max_rel_err < 0.5 ) { //error too small - increase dt
+ t += dt; // we keep the evolution -> increase time
+ // limit scaling factor to 5.0
+ dt *= min( 0.9*pow(max_rel_err , -1.0/(stepper.order()+1.0)), 5.0 );
+ return STEP_SIZE_INCREASED;
+ } else {
+ t += dt;
+ return SUCCESS;
+ }
+ }
+ };
} // namespace odeint
} // namespace numeric
Added: sandbox/odeint/libs/numeric/odeint/examples/lorenz_controlled.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/libs/numeric/odeint/examples/lorenz_controlled.cpp 2009-10-30 13:56:39 EDT (Fri, 30 Oct 2009)
@@ -0,0 +1,89 @@
+/* Boost numeric/odeint/examples/lorenz_array.cpp
+
+ Copyright 2009 Karsten Ahnert
+
+ Shows, the usage of odeint, and integrates the Lorenz equations,
+ a deterministic chaotic system
+
+ dx/dt = sigma * ( x - y)
+ dy/dt = R*x - y - x*z
+ dz/dt = x*y - b*z
+
+ mit sigma = 10, r=28, b = 8/3
+
+ 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 std::tr1;
+using namespace boost::numeric::odeint;
+
+const double sigma = 10.0;
+const double R = 28.0;
+const double b = 8.0 / 3.0;
+
+const size_t olen = 10000;
+
+const double eps_abs = 1E-3;
+const double eps_rel = 1E-3;
+
+const double min_dt = 1E-10;
+
+typedef array<double, 3> state_type;
+
+void lorenz( state_type &x , state_type &dxdt , double t )
+{
+ dxdt[0] = sigma * ( x[1] - x[0] );
+ dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
+ dxdt[2] = x[0]*x[1] - b * x[2];
+}
+
+int main( int argc , char **argv )
+{
+ state_type x;
+ x[0] = 1.0;
+ x[1] = 0.0;
+ x[2] = 20.0;
+
+ ode_step_euler< state_type > euler;
+ step_controller_standard< state_type, double > controlled_euler(eps_abs,
+ eps_rel,
+ 1.0, 1.0);
+
+ double t = 0.0;
+ double dt = 0.01;
+ controlled_step_result result;
+
+ cout.precision(5);
+ cout.setf(ios::fixed,ios::floatfield);
+
+
+ while( (t<10.0) && (dt > min_dt) ) {
+ //cout << "current state: " ;
+ cout << t << tab << dt << tab << x[0] << tab << x[1] << tab << x[2] << endl;
+ result = controlled_euler.controlled_step( euler , lorenz , x , t , dt );
+ while( result != SUCCESS ) {
+ result = controlled_euler.controlled_step( euler , lorenz , x , t , dt );
+ if( dt < min_dt ) break;
+ }
+ //cout<<"SUCCESS with dt="<<dt<<endl;
+ }
+
+ return 0;
+}
+
+/*
+ Compile with
+ g++ -Wall -I$BOOST_ROOT -I../../../../ lorenz_array.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