Subject: [Boost-commit] svn:boost r57915 - in sandbox/odeint: boost/numeric boost/numeric/odeint libs/numeric/odeint/examples
From: karsten.ahnert_at_[hidden]
Date: 2009-11-24 17:58:22
Author: karsten
Date: 2009-11-24 17:58:20 EST (Tue, 24 Nov 2009)
New Revision: 57915
URL: http://svn.boost.org/trac/boost/changeset/57915
Log:
starting implementation of rk78 fehlberg
Added:
sandbox/odeint/boost/numeric/odeint/stepper_rk78_fehlberg.hpp (contents, props changed)
Text files modified:
sandbox/odeint/boost/numeric/odeint.hpp | 1
sandbox/odeint/boost/numeric/odeint/controlled_stepper_standard.hpp | 120 ++++++++++++++++++++++++++-------------
sandbox/odeint/boost/numeric/odeint/stepper_euler.hpp | 1
sandbox/odeint/boost/numeric/odeint/stepper_midpoint.hpp | 22 +++----
sandbox/odeint/libs/numeric/odeint/examples/lorenz_cmp_rk4_rk_generic.cpp | 7 --
sandbox/odeint/libs/numeric/odeint/examples/lorenz_controlled.cpp | 6 +-
sandbox/odeint/libs/numeric/odeint/examples/lorenz_stepper.cpp | 3
7 files changed, 94 insertions(+), 66 deletions(-)
Modified: sandbox/odeint/boost/numeric/odeint.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint.hpp 2009-11-24 17:58:20 EST (Tue, 24 Nov 2009)
@@ -22,6 +22,7 @@
#include <boost/numeric/odeint/stepper_rk_generic.hpp>
#include <boost/numeric/odeint/stepper_half_step.hpp>
#include <boost/numeric/odeint/stepper_midpoint.hpp>
+#include <boost/numeric/odeint/stepper_rk78_fehlberg.hpp>
#include <boost/numeric/odeint/controlled_stepper_standard.hpp>
Modified: sandbox/odeint/boost/numeric/odeint/controlled_stepper_standard.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/controlled_stepper_standard.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/controlled_stepper_standard.hpp 2009-11-24 17:58:20 EST (Tue, 24 Nov 2009)
@@ -27,7 +27,11 @@
namespace numeric {
namespace odeint {
- typedef enum{success, step_size_decreased, step_size_increased} controlled_step_result;
+ typedef enum{
+ success ,
+ step_size_decreased ,
+ step_size_increased
+ } controlled_step_result;
/*
The initial state is given in x.
@@ -54,7 +58,8 @@
template<
class ErrorStepper,
class ResizeType = resizer< typename ErrorStepper::container_type > >
- class controlled_stepper_standard {
+ class controlled_stepper_standard
+ {
public:
@@ -65,20 +70,47 @@
typedef typename container_type::value_type value_type;
typedef typename container_type::iterator iterator;
- typedef const unsigned short order_type;
-
// private members
private:
+
ErrorStepper &m_stepper;
- time_type eps_abs;
- time_type eps_rel;
- time_type a_x;
- time_type a_dxdt;
- container_type dxdt;
- container_type x_tmp;
- container_type x_err;
- resizer_type resizer;
+ time_type m_eps_abs;
+ time_type m_eps_rel;
+ time_type m_a_x;
+ time_type m_a_dxdt;
+ container_type m_dxdt;
+ container_type m_x_tmp;
+ container_type m_x_err;
+ resizer_type m_resizer;
+
+
+
+ // private methods
+
+ time_type calc_max_rel_err(
+ iterator x_start ,
+ iterator x_end ,
+ iterator dxdt_start ,
+ iterator x_err_start ,
+ time_type dt
+ )
+ {
+ time_type max_rel_err = 0.0;
+ time_type err;
+
+ while( x_start != x_end )
+ {
+ // get the maximal value of x_err/D where
+ // D = eps_abs + eps_rel * (a_x*|x| + a_dxdt*|dxdt|);
+ err = m_eps_abs + m_eps_rel * (
+ m_a_x * std::abs(*x_start++) +
+ m_a_dxdt * dt * std::abs(*dxdt_start++) );
+ max_rel_err = max( std::abs(*x_err_start++)/err , max_rel_err );
+ }
+ return max_rel_err;
+ }
+
// public functions
@@ -89,7 +121,10 @@
time_type abs_err, time_type rel_err,
time_type factor_x, time_type factor_dxdt )
: m_stepper(stepper),
- eps_abs(abs_err), eps_rel(rel_err), a_x(factor_x), a_dxdt(factor_dxdt)
+ m_eps_abs(abs_err),
+ m_eps_rel(rel_err),
+ m_a_x(factor_x),
+ m_a_dxdt(factor_dxdt)
{ }
@@ -107,42 +142,45 @@
time_type &t,
time_type &dt )
{
- resizer.adjust_size(x, x_err);
+ m_resizer.adjust_size( x , m_x_err );
+ m_resizer.adjust_size( x , m_dxdt );
- x_tmp = x; // copy current state
- system( x, dxdt, t ); // compute dxdt
- m_stepper.do_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();
- time_type max_rel_err = 0.0;
+ m_x_tmp = x;
+ system( x , m_dxdt , t );
+ m_stepper.do_step( system , x , m_dxdt, t , dt , m_x_err );
- 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|);
- time_type 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 );
- }
+ time_type max_rel_err = calc_max_rel_err(
+ m_x_tmp.begin() , m_x_tmp.end() ,
+ m_dxdt.begin() , m_x_err.begin() , dt );
- //std::cout<<max_rel_err<<std::endl;
- if( max_rel_err > 1.1 ) { // error too large - decrease dt
+ 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/(m_stepper.order_error()-1.0)) , 0.2 );
+ dt *= max( 0.9*pow(max_rel_err , -1.0/(m_stepper.order_error()-1.0)),
+ 0.2 );
+
// reset state
- x = x_tmp;
+ x = m_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/m_stepper.order()), 5.0 );
- return step_size_increased;
- } else {
- t += dt;
- return success;
}
+ else
+ {
+ if( max_rel_err < 0.5 )
+ {
+ //error too small - increase dt and keep the evolution
+ t += dt;
+ // limit scaling factor to 5.0
+ dt *= min( 0.9*pow(max_rel_err , -1.0/m_stepper.order()), 5.0 );
+ return step_size_increased;
+ }
+ else
+ {
+ t += dt;
+ return success;
+ }
+ }
}
};
Modified: sandbox/odeint/boost/numeric/odeint/stepper_euler.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepper_euler.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/stepper_euler.hpp 2009-11-24 17:58:20 EST (Tue, 24 Nov 2009)
@@ -63,7 +63,6 @@
private:
container_type m_dxdt;
- container_type m_xtemp;
resizer_type m_resizer;
Modified: sandbox/odeint/boost/numeric/odeint/stepper_midpoint.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepper_midpoint.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/stepper_midpoint.hpp 2009-11-24 17:58:20 EST (Tue, 24 Nov 2009)
@@ -58,6 +58,7 @@
// private memebers
private:
+
resizer_type m_resizer;
unsigned short m_stepcount;
@@ -68,8 +69,7 @@
public:
- stepper_midpoint( unsigned short stepcount = 2 )
- { }
+ stepper_midpoint( unsigned short stepcount = 2 ) { }
order_type order() const { return 2; }
@@ -79,10 +79,7 @@
m_stepcount = stepcount;
}
- unsigned short get_step_count()
- {
- return m_stepcount;
- }
+ unsigned short get_step_count() const { return m_stepcount; }
template< class DynamicalSystem >
void do_step(
@@ -90,13 +87,13 @@
container_type &x ,
container_type &dxdt ,
time_type t ,
- time_type dt
- )
+ time_type dt )
{
- const time_type h = dt/static_cast<time_type>( m_stepcount );
- const time_type h2 = static_cast<time_type>( 2.0 )*h;
const time_type t_1 = static_cast<time_type>( 1.0 );
const time_type t_05 = static_cast<time_type>( 0.5 );
+
+ const time_type h = dt/static_cast<time_type>( m_stepcount );
+ const time_type h2 = static_cast<time_type>( 2.0 )*h;
time_type th = t + h;
m_resizer.adjust_size(x, m_x0);
@@ -139,12 +136,11 @@
DynamicalSystem &system ,
container_type &x ,
time_type t ,
- time_type dt ,
- unsigned short n = 2 )
+ time_type dt )
{
m_resizer.adjust_size(x, m_dxdt);
system( x, m_dxdt, t );
- do_step( system , x, m_dxdt, t, dt, n );
+ do_step( system , x, m_dxdt, t, dt );
}
Added: sandbox/odeint/boost/numeric/odeint/stepper_rk78_fehlberg.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/boost/numeric/odeint/stepper_rk78_fehlberg.hpp 2009-11-24 17:58:20 EST (Tue, 24 Nov 2009)
@@ -0,0 +1,230 @@
+/*
+ boost header: numeric/odeint/stepper_rk78_fehlberg.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_STEPPER_RK78_FEHLBERG_HPP_INCLUDED
+#define BOOST_NUMERIC_ODEINT_STEPPER_RK78_FEHLBERG_HPP_INCLUDED
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+ template<
+ class Container ,
+ class Time = double ,
+ class Resizer = resizer< Container >
+ >
+ class stepper_rk78_fehlberg
+ {
+ // provide basic typedefs
+ public:
+
+ typedef Container container_type;
+ typedef Resizer resizer_type;
+ typedef Time time_type;
+ typedef const unsigned short order_type;
+ typedef typename container_type::value_type value_type;
+ typedef typename container_type::iterator iterator;
+
+
+
+
+
+ // private members
+ private:
+
+ container_type m_dxdt;
+ container_type m_xt;
+ container_type m_k02 , m_k03 , m_k04 , m_k05 , m_k06 , m_k07 ,
+ m_k08 , m_k09 , m_k10 , m_k11 , m_k12;
+
+ resizer_type m_resizer;
+
+ // the times at which system is called
+ time_type m_t02 , m_t03 , m_t04 , m_t05 , m_t06 , m_t07 , m_t08 ,
+ m_t09 , m_t10 , m_t11 , m_t12 , m_t13;
+
+
+
+
+ // public interface
+ public:
+
+ order_type order() const { return 8; }
+
+
+ template< class DynamicalSystem >
+ void do_step( DynamicalSystem &system ,
+ container_type &x ,
+ const container_type &dxdt ,
+ time_type t ,
+ time_type dt )
+ {
+ // the time constant
+ const time_type a02 = static_cast<time_type>( 2.0 / 27.0 );
+ const time_type a03 = static_cast<time_type>( 1.0 / 9.0 );
+ const time_type a04 = static_cast<time_type>( 1.0 / 6.0 );
+ const time_type a05 = static_cast<time_type>( 5.0 / 12.0 );
+ const time_type a06 = static_cast<time_type>( 0.50 );
+ const time_type a07 = static_cast<time_type>( 5.0 / 6.0 );
+ const time_type a08 = static_cast<time_type>( 1.0 / 6.0 );
+ const time_type a09 = static_cast<time_type>( 2.0 / 3.0 );
+ const time_type a10 = static_cast<time_type>( 1.0 / 3.0 );
+
+ // the weights for each step
+ const time_type c06 = static_cast<time_type>( 34.0 / 105.0 );
+ const time_type c07 = static_cast<time_type>( 9.0 / 35.0 );
+ const time_type c08 = static_cast<time_type>( c07 );
+ const time_type c09 = static_cast<time_type>( 9.0 / 280.0 );
+ const time_type c10 = static_cast<time_type>( c09 );
+ const time_type c12 = static_cast<time_type>( 41.0 / 840.0 );
+ const time_type c13 = static_cast<time_type>( c12 );
+ const time_type c11 = static_cast<time_type>( 0. );
+
+ // the coefficients for each step
+ const time_type b02_01 = static_cast<time_type>( 2.0 / 27.0 );
+
+ const time_type b03_01 = static_cast<time_type>( 1.0 / 36.0 );
+ const time_type b03_02 = static_cast<time_type>( 1.0 / 12.0 );
+
+ const time_type b04_01 = static_cast<time_type>( 1.0 / 24.0 );
+ const time_type b04_03 = static_cast<time_type>( 1.0 / 8.0 );
+
+ const time_type b05_01 = static_cast<time_type>( 5.0 / 12.0 );
+ const time_type b05_03 = static_cast<time_type>( -25.0 / 16.0 );
+ const time_type b05_04 = static_cast<time_type>( -b05_03 );
+
+ const time_type b06_01 = static_cast<time_type>( 0.050 );
+ const time_type b06_04 = static_cast<time_type>( 0.250 );
+ const time_type b06_05 = static_cast<time_type>( 0.20 );
+
+ const time_type b07_01 = static_cast<time_type>( -25.0 / 108.0 );
+ const time_type b07_04 = static_cast<time_type>( 125.0 / 108.0 );
+ const time_type b07_05 = static_cast<time_type>( -65.0 / 27.0 );
+ const time_type b07_06 = static_cast<time_type>( 125.0 / 54.0 );
+
+ const time_type b08_01 = static_cast<time_type>( 31.0 / 300.0 );
+ const time_type b08_05 = static_cast<time_type>( 61.0 / 225.0 );
+ const time_type b08_06 = static_cast<time_type>( -2.0 / 9.0 );
+ const time_type b08_07 = static_cast<time_type>( 13.0 / 900.0 );
+
+ const time_type b09_01 = static_cast<time_type>( 2.0 );
+ const time_type b09_04 = static_cast<time_type>( -53.0 / 6.0 );
+ const time_type b09_05 = static_cast<time_type>( 704.0 / 45.0 );
+ const time_type b09_06 = static_cast<time_type>( -107.0 / 9.0 );
+ const time_type b09_07 = static_cast<time_type>( 67.0 / 90.0 );
+ const time_type b09_08 = static_cast<time_type>( 3.0 );
+
+ const time_type b10_01 = static_cast<time_type>( -91.0 / 108.0 );
+ const time_type b10_04 = static_cast<time_type>( 23.0 / 108.0 );
+ const time_type b10_05 = static_cast<time_type>( -976.0 / 135.0 );
+ const time_type b10_06 = static_cast<time_type>( 311.0 / 54.0 );
+ const time_type b10_07 = static_cast<time_type>( -19.0 / 60.0 );
+ const time_type b10_08 = static_cast<time_type>( 17.0 / 6.0 );
+ const time_type b10_09 = static_cast<time_type>( -1.0 / 12.0 );
+
+ const time_type b11_01 = static_cast<time_type>( 2383.0 / 4100.0 );
+ const time_type b11_04 = static_cast<time_type>( -341.0 / 164.0 );
+ const time_type b11_05 = static_cast<time_type>( 4496.0 / 1025.0 );
+ const time_type b11_06 = static_cast<time_type>( -301.0 / 82.0 );
+ const time_type b11_07 = static_cast<time_type>( 2133.0 / 4100.0 );
+ const time_type b11_08 = static_cast<time_type>( 45.0 / 82.0 );
+ const time_type b11_09 = static_cast<time_type>( 45.0 / 164.0 );
+ const time_type b11_10 = static_cast<time_type>( 18.0 / 41.0 );
+
+ const time_type b12_01 = static_cast<time_type>( 3.0 / 205.0 );
+ const time_type b12_06 = static_cast<time_type>( -6.0 / 41.0 );
+ const time_type b12_07 = static_cast<time_type>( -3.0 / 205.0 );
+ const time_type b12_08 = static_cast<time_type>( -3.0 / 41.0 );
+ const time_type b12_09 = static_cast<time_type>( 3.0 / 41.0 );
+ const time_type b12_10 = static_cast<time_type>( 6.0 / 41.0 );
+
+ const time_type b13_01 = static_cast<time_type>( -1777.0 / 4100.0 );
+ const time_type b13_04 = static_cast<time_type>( b11_04 );
+ const time_type b13_05 = static_cast<time_type>( b11_05 );
+ const time_type b13_06 = static_cast<time_type>( -289.0 / 82.0 );
+ const time_type b13_07 = static_cast<time_type>( 2193.0 / 4100.0 );
+ const time_type b13_08 = static_cast<time_type>( 51.0 / 82.0 );
+ const time_type b13_09 = static_cast<time_type>( 33.0 / 164.0 );
+ const time_type b13_10 = static_cast<time_type>( 12.0 / 41.0 );
+ const time_type b13_12 = static_cast<time_type>( 1.0 );
+
+ const time_type val1 = static_cast<time_type>( 1.0 );
+
+ using namespace detail::it_algebra;
+
+ // compute the times at which system is evaluated
+ m_t02 = t + a02 * dt;
+ m_t03 = t + a03 * dt;
+ m_t04 = t + a04 * dt;
+ m_t05 = t + a05 * dt;
+ m_t06 = t + a06 * dt;
+ m_t07 = t + a07 * dt;
+ m_t08 = t + a08 * dt;
+ m_t09 = t + a09 * dt;
+ m_t10 = t + a10 * dt;
+ m_t11 = t + dt;
+ m_t12 = t;
+ m_t13 = t + dt;
+
+ // resize
+ m_resizer.adjust_size( x , m_xt );
+ m_resizer.adjust_size( x , m_k02 );
+ m_resizer.adjust_size( x , m_k03 );
+ m_resizer.adjust_size( x , m_k04 );
+ m_resizer.adjust_size( x , m_k05 );
+ m_resizer.adjust_size( x , m_k06 );
+ m_resizer.adjust_size( x , m_k07 );
+ m_resizer.adjust_size( x , m_k08 );
+ m_resizer.adjust_size( x , m_k09 );
+ m_resizer.adjust_size( x , m_k10 );
+ m_resizer.adjust_size( x , m_k11 );
+ m_resizer.adjust_size( x , m_k12 );
+
+
+ // k1, the first system call has allready been evaluated
+ scale_sum( m_xt.begin() , m_xt.end() ,
+ val1 , x.begin() ,
+ dt * b02_01 , dxdt.begin() );
+
+ // k2 step
+ system( m_xt , m_k02 , m_t02 );
+ scale_sum( m_xt.begin() , m_xt.end() ,
+ val1 , x.begin() ,
+ dt * b03_01 , dxdt.begin() ,
+ dt * b03_02 , m_k02 );
+
+
+
+
+
+
+ }
+
+
+ template< class DynamicalSystem >
+ void do_step( DynamicalSystem &system ,
+ container_type &x ,
+ time_type t ,
+ time_type dt )
+ {
+ m_resizer.adjust_size( x , m_dxdt );
+ system( x , m_dxdt , t );
+ do_step( system , x , m_dxdt , t , dt );
+ }
+ };
+
+} // namespace odeint
+} // namespace numeric
+} // namespace boost
+
+
+#endif //BOOST_NUMERIC_ODEINT_STEPPER_RK78_FEHLBERG_HPP_INCLUDED
Modified: sandbox/odeint/libs/numeric/odeint/examples/lorenz_cmp_rk4_rk_generic.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/lorenz_cmp_rk4_rk_generic.cpp (original)
+++ sandbox/odeint/libs/numeric/odeint/examples/lorenz_cmp_rk4_rk_generic.cpp 2009-11-24 17:58:20 EST (Tue, 24 Nov 2009)
@@ -147,10 +147,3 @@
return 0;
}
-
-
-
-/*
- Compile with
- g++ -Wall -O3 -I$BOOST_ROOT -I../../../../ lorenz_cmp_rk4_rk_generic.cpp
-*/
Modified: sandbox/odeint/libs/numeric/odeint/examples/lorenz_controlled.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/lorenz_controlled.cpp (original)
+++ sandbox/odeint/libs/numeric/odeint/examples/lorenz_controlled.cpp 2009-11-24 17:58:20 EST (Tue, 24 Nov 2009)
@@ -36,8 +36,8 @@
const size_t olen = 10000;
-const double eps_abs = 1E-6;
-const double eps_rel = 1E-7;
+const double eps_abs = 1E-2;
+const double eps_rel = 1E-3;;
typedef array<double, 3> state_type;
@@ -67,7 +67,7 @@
cout.precision(5);
cout.setf(ios::fixed,ios::floatfield);
- size_t steps = integrate_adaptive( controlled_stepper, lorenz, x, 0.0, 10.0, 1E-4, print_state );
+ size_t steps = integrate_adaptive( controlled_stepper, lorenz, x, 0.0, 10.0, 1E-2, print_state );
clog << "Number of steps: " << steps << endl;
Modified: sandbox/odeint/libs/numeric/odeint/examples/lorenz_stepper.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/lorenz_stepper.cpp (original)
+++ sandbox/odeint/libs/numeric/odeint/examples/lorenz_stepper.cpp 2009-11-24 17:58:20 EST (Tue, 24 Nov 2009)
@@ -78,7 +78,8 @@
start= clock();
t = 0.0;
for( size_t oi=0 ; oi<olen ; ++oi,t+=dt )
- stepper1.do_step( lorenz1 , x1 , t , dt );
+ stepper1.do_step
+( lorenz1 , x1 , t , dt );
end = clock();
cout << "vector : " << double ( end - start ) / double( CLOCKS_PER_SEC ) << endl;
cout << "x: "<<x1[0]<<tab<<x1[1]<<tab<<x1[2]<<endl;