Subject: [Boost-commit] svn:boost r57742 - in sandbox/odeint: boost/numeric/odeint libs/numeric/odeint/examples libs/numeric/odeint/stuff/gsl_compare
From: mario.mulansky_at_[hidden]
Date: 2009-11-18 07:50:35
Author: mariomulansky
Date: 2009-11-18 07:50:34 EST (Wed, 18 Nov 2009)
New Revision: 57742
URL: http://svn.boost.org/trac/boost/changeset/57742
Log:
added correct compuation of the order of error term
Binary files modified:
sandbox/odeint/boost/numeric/odeint/stepper_rk5_ck.hpp
Text files modified:
sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp | 5 ++
sandbox/odeint/boost/numeric/odeint/stepsize_controller_standard.hpp | 62 +++++++++++++++++++++++++++++----------
sandbox/odeint/libs/numeric/odeint/examples/lorenz_integrator.cpp | 19 +++++++++---
sandbox/odeint/libs/numeric/odeint/stuff/gsl_compare/lorenz_stepper_cmp.cpp | 3 +
4 files changed, 66 insertions(+), 23 deletions(-)
Modified: sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp 2009-11-18 07:50:34 EST (Wed, 18 Nov 2009)
@@ -65,7 +65,10 @@
order_type order() const { return m_stepper.order(); }
-
+ order_type order_error() const
+ { /* Order of the error term is the order of the underlying stepper + 1 */
+ return m_stepper.order() + 1;
+ }
template< class DynamicalSystem >
void next_step( DynamicalSystem &system ,
Modified: sandbox/odeint/boost/numeric/odeint/stepper_rk5_ck.hpp
==============================================================================
Binary files. No diff available.
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-11-18 07:50:34 EST (Wed, 18 Nov 2009)
@@ -28,6 +28,36 @@
typedef enum{SUCCESS, STEP_SIZE_DECREASED, STEP_SIZE_INCREASED} controlled_step_result;
+ /*
+ The initial state is given in x.
+ t is an vector including the times at which the state will be written into
+ the vector x_vec.
+ x_vec must provide enough space to hold times.size() states.
+ dt is the initial step size (will be adjusted according to the errors).
+ This function returns the total number of steps required to integrate the
+ whole intervale times.begin() - times.end().
+ Note that the values in times don't influence the stepsize, but only the
+ time points at which the state is stored into x_vec.
+
+ The stepsize is adjust such that the following maximal relative error is
+ small enough for each step:
+ R = max( x_err_n / [eps_abs + eps_rel*( a_x * |x_n| + a_dxdt * |dxdt_n| )] )
+ where the max refers to the componentwise maximum the expression.
+
+ if R > 1.1 the stepsize is decreased:
+ dt = dt*S*R^(-1/(q-1))
+
+ if R < 0.5 the stepsize is increased:
+ dt = dt*S*R^(-1/q))
+
+ q is the order of the error term provided by the stepper.order_error() function
+ (e.g. 2 for simple euler and 5 for rk5_ck) and S is a safety factor set to
+ S = 0.9. See e.g. Numerical Recipes for more details on this strategy.
+
+ To avoid extensive chages in dt, the decrease factor is limited to 0.2 and
+ the increase factor to 5.0.
+ */
+
template<
class ContainerType,
class T,
@@ -72,29 +102,29 @@
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 );
+ // 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;
+ // limit scaling factor to 0.2
+ dt *= max( 0.9*pow(max_rel_err , -1.0/(stepper.order_error()-1.0)) , 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;
+ 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()), 5.0 );
+ return STEP_SIZE_INCREASED;
} else {
- t += dt;
- return SUCCESS;
+ t += dt;
+ return SUCCESS;
}
}
};
Modified: sandbox/odeint/libs/numeric/odeint/examples/lorenz_integrator.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/lorenz_integrator.cpp (original)
+++ sandbox/odeint/libs/numeric/odeint/examples/lorenz_integrator.cpp 2009-11-18 07:50:34 EST (Wed, 18 Nov 2009)
@@ -56,10 +56,12 @@
x1[1] = 0.0;
x1[2] = 20.0;
state_type x2(x1);
+ state_type x3(x1);
vector<state_type> x1_t_vec;
vector<state_type> x2_t_vec;
+ vector<state_type> x3_t_vec;
vector<double> times(time_points);
for( size_t i=0; i<time_points; i++ ) {
times[i] = 0.1*i;
@@ -68,14 +70,20 @@
stepper_half_step< stepper_euler< state_type > > euler;
size_t steps = integrate( euler, lorenz, x1, times, back_inserter(x1_t_vec));
- clog << "Euler Steps: " << steps << endl;
+ clog << "Euler Half Step: #steps " << steps << endl;
+
+ stepper_half_step< stepper_rk4< state_type > > rk4;
+ steps = integrate( rk4, lorenz, x2, times, back_inserter(x2_t_vec));
+
+ clog << "RK4 Half Step: #steps " << steps << endl;
+
stepper_rk5_ck< state_type > rk5;
- steps = integrate( rk5, lorenz, x2, times, back_inserter(x2_t_vec));
+ steps = integrate( rk5, lorenz, x3, times, back_inserter(x3_t_vec));
- clog << "RK5 Steps: " << steps << endl;
+ clog << "RK5 Cash-Karp: #steps: " << steps << endl;
- cout.precision(5);
+ cout.precision(16);
cout.setf(ios::fixed,ios::floatfield);
@@ -83,7 +91,8 @@
//cout << "current state: " ;
cout << times[i] << tab;
cout << x1_t_vec[i][0] << tab << x1_t_vec[i][1] << tab << x1_t_vec[i][2] << tab;
- cout << x2_t_vec[i][0] << tab << x2_t_vec[i][1] << tab << x2_t_vec[i][2] << endl;
+ cout << x2_t_vec[i][0] << tab << x2_t_vec[i][1] << tab << x2_t_vec[i][2] << tab;
+ cout << x3_t_vec[i][0] << tab << x3_t_vec[i][1] << tab << x3_t_vec[i][2] << endl;
}
return 0;
Modified: sandbox/odeint/libs/numeric/odeint/stuff/gsl_compare/lorenz_stepper_cmp.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/stuff/gsl_compare/lorenz_stepper_cmp.cpp (original)
+++ sandbox/odeint/libs/numeric/odeint/stuff/gsl_compare/lorenz_stepper_cmp.cpp 2009-11-18 07:50:34 EST (Wed, 18 Nov 2009)
@@ -142,6 +142,7 @@
// gsl method rk4 cash-karp
gsl_odeiv_step *s_rkck = gsl_odeiv_step_alloc( gsl_odeiv_step_rkck , 3);
+ clog << "Order of gsl " << gsl_odeiv_step_name(s_rkck) << ": " << gsl_odeiv_step_order(s_rkck) << endl;
// sys = { lorenz_gsl , 0 , 3 , 0 };
copy( x_start.begin() , x_start.end() , x2 );
@@ -150,7 +151,7 @@
for( size_t oi=0 ; oi<olen ; ++oi,t+=dt )
gsl_odeiv_step_apply ( s_rkck , t , dt , x2 , x2_err , 0 , 0 , &sys );
end = clock();
- clog << "gsl rk4 Cash-Karp : " << double ( end - start ) / double( CLOCKS_PER_SEC ) << endl;
+ clog << "gsl rk Cash-Karp : " << double ( end - start ) / double( CLOCKS_PER_SEC ) << endl;