[Boost-commit] svn:boost r57742 - in sandbox/odeint: boost/numeric/odeint libs/numeric/odeint/examples libs/numeric/odeint/stuff/gsl_compare

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;