[Boost-commit] svn:boost r57376 - in sandbox/odeint: boost/numeric/odeint libs/numeric/odeint/examples

Subject: [Boost-commit] svn:boost r57376 - in sandbox/odeint: boost/numeric/odeint libs/numeric/odeint/examples
From: mario.mulansky_at_[hidden]
Date: 2009-11-04 14:17:23

Author: mariomulansky
Date: 2009-11-04 14:17:22 EST (Wed, 04 Nov 2009)
New Revision: 57376
URL: http://svn.boost.org/trac/boost/changeset/57376

Log:
added adaptive intagrate method with observer
Text files modified: 
   sandbox/odeint/boost/numeric/odeint/integrator.hpp                |   139 ++++++++++++++++++++++++++------------- 
   sandbox/odeint/boost/numeric/odeint/integrator_constant_step.hpp  |    46 ++++++-------                           
   sandbox/odeint/libs/numeric/odeint/examples/lorenz_controlled.cpp |    32 +++------                               
   sandbox/odeint/libs/numeric/odeint/examples/lorenz_integrator.cpp |    10 +-                                      
   4 files changed, 130 insertions(+), 97 deletions(-)

Modified: sandbox/odeint/boost/numeric/odeint/integrator.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/integrator.hpp	(original)
+++ sandbox/odeint/boost/numeric/odeint/integrator.hpp	2009-11-04 14:17:22 EST (Wed, 04 Nov 2009)
@@ -23,6 +23,48 @@
 namespace odeint {
 
 
+    template<
+        class Stepper,
+        class DynamicalSystem,
+        class StepController,
+        class Observer
+        >
+    size_t integrate( Stepper &stepper,
+                      DynamicalSystem &system,
+                      StepController &controller,
+                      typename Stepper::time_type start_time,
+                      typename Stepper::time_type dt,
+                      typename Stepper::container_type &state,
+                      typename Stepper::time_type end_time,
+                      Observer &observer )
+    {
+        controlled_step_result result;
+        size_t iterations = 0;
+        typename Stepper::time_type t = start_time;
+
+        observer(t, state, system);
+        
+        while( t < end_time ) {
+            result = controller.controlled_step( stepper, system, state, t, dt );
+            if( result != STEP_SIZE_DECREASED ) { // we actually did a step forward
+                observer(t, state, system);
+                iterations++;
+            }
+            while( result != SUCCESS ) {
+                result = controller.controlled_step( stepper, system, state, t, dt );
+                if( result != STEP_SIZE_DECREASED ) { // we did a step
+                    observer(t, state, system);
+                    iterations++;
+                }
+                if( !( t+dt > t) ) 
+                    throw; // we've reached machine precision with dt - no advancing in t
+            }
+        }
+        return iterations;
+    }
+
+
+
     /* Integration of an ode with adaptive stepsize methods.
        Integrates an ode give by system using the integration scheme stepper and the
        step-size controller controller.
@@ -36,10 +78,11 @@
        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.
      */
-    template< class StepType,
-	      class DynamicalSystem,
-	      class StepController,
-	      class T >
+    template< 
+        class StepType,
+        class DynamicalSystem,
+        class StepController,
+        class T >
     size_t integrate(StepType &stepper,
                      DynamicalSystem &system,
                      StepController &controller,
@@ -48,34 +91,38 @@
                      std::vector<typename StepType::container_type> &x_vec,
                      T dt)
     {
-	if( times.size() != x_vec.size() ) throw;
+        if( times.size() != x_vec.size() ) throw;
 
-	// iterators for the time and state vectors
-	typename std::vector<T>::iterator t_iter = times.begin();
-	typename std::vector<typename StepType::container_type>::iterator x_iter = x_vec.begin();
-
-	controlled_step_result result;
-	size_t iterations = 0;
-	T t = *t_iter;
-
-	while( t_iter < times.end() ) {
-
-	    if( t >= *t_iter ) { // we've reached the next time point
-		*x_iter++ = x; // save the vector
-		t_iter++; // next time point
-	    }
-
-	    result = controller.controlled_step( stepper, system, x, t, dt );
-	    while( result != SUCCESS ) {
-		result = controller.controlled_step( stepper, system, x, t, dt );
-		if( result == STEP_SIZE_INCREASED )
-		    iterations++;
-		if( !( t+dt > t) ) 
-		    throw; // we've reached machine precision with dt - no advancing in t
-	    }
-	    iterations++;
-	}
-	return iterations;
+        // iterators for the time and state vectors
+        typename std::vector<T>::iterator t_iter = times.begin();
+        typename std::vector<typename StepType::container_type>::iterator x_iter = x_vec.begin();
+
+        controlled_step_result result;
+        size_t iterations = 0;
+        T t = *t_iter;
+
+        while( true ) { // loop will break from inside
+
+            if( t >= *t_iter ) { // we've reached the next time point
+                *x_iter++ = x; // save the state
+                t_iter++; // next time point
+            }
+
+            if( t_iter >= times.end() ) // reached end of integration time
+                break; // stop loop
+
+            result = controller.controlled_step( stepper, system, x, t, dt );
+            if( result != STEP_SIZE_DECREASED )
+                iterations++;
+            while( result != SUCCESS ) {
+                result = controller.controlled_step( stepper, system, x, t, dt );
+                if( result != STEP_SIZE_DECREASED )
+                    iterations++;
+                if( !( t+dt > t) ) 
+                    throw; // we've reached machine precision with dt - no advancing in t
+            }
+        }
+        return iterations;
     }
 
 
@@ -112,23 +159,23 @@
        the increase factor to 5.0.
      */
     template< class StepType,
-	      class DynamicalSystem,
-	      class T >
+              class DynamicalSystem,
+              class T >
     size_t integrate(StepType &stepper, 
-		     DynamicalSystem &system, 
-		     typename StepType::container_type &x, 
-		     std::vector<T> &times, 
-		     std::vector<typename StepType::container_type> &x_vec,
-		     T dt = 1E-4, T eps_abs = 1E-6, 
-		     T eps_rel = 1E-7, T a_x = 1.0 , T a_dxdt = 1.0)
+                     DynamicalSystem &system, 
+                     typename StepType::container_type &x, 
+                     std::vector<T> &times, 
+                     std::vector<typename StepType::container_type> &x_vec,
+                     T dt = 1E-4, T eps_abs = 1E-6, 
+                     T eps_rel = 1E-7, T a_x = 1.0 , T a_dxdt = 1.0)
     {
-	if( times.size() != x_vec.size() ) throw;
-	// we use the standard controller for this adaptive integrator
-	step_controller_standard< typename StepType::container_type, T, typename StepType::resizer_type>
-	    controller(eps_abs, eps_rel, a_x, a_dxdt ); // initialized with values from above
-	
-	// call the normal integrator
-	return integrate(stepper, system, controller, x, times, x_vec, dt);
+        if( times.size() != x_vec.size() ) throw;
+        // we use the standard controller for this adaptive integrator
+        step_controller_standard< typename StepType::container_type, T, typename StepType::resizer_type>
+            controller(eps_abs, eps_rel, a_x, a_dxdt ); // initialized with values from above
+        
+        // call the normal integrator
+        return integrate(stepper, system, controller, x, times, x_vec, dt);
     }
     
 

Modified: sandbox/odeint/boost/numeric/odeint/integrator_constant_step.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/integrator_constant_step.hpp	(original)
+++ sandbox/odeint/boost/numeric/odeint/integrator_constant_step.hpp	2009-11-04 14:17:22 EST (Wed, 04 Nov 2009)
@@ -20,37 +20,33 @@
 namespace odeint {
 
 
-
     template<
-	class Stepper ,
-	class DynamicalSystem ,
-	class Observer
-	>
+        class Stepper ,
+        class DynamicalSystem ,
+        class Observer
+        >
     void integrate(
-	Stepper stepper ,
-	DynamicalSystem system ,
-	typename Stepper::time_type start_time ,
-	typename Stepper::time_type dt ,
-	typename Stepper::container_type &state ,
-	typename Stepper::time_type end_time ,
-	Observer observer
-	)
+                   Stepper &stepper ,
+                   DynamicalSystem &system ,
+                   typename Stepper::time_type start_time ,
+                   typename Stepper::time_type dt ,
+                   typename Stepper::container_type &state ,
+                   typename Stepper::time_type end_time ,
+                   Observer &observer
+                   )
     {
-	if( start_time > end_time )
-	    throw std::invalid_argument( "integrate() : start_time > end_time" );
+        if( start_time > end_time )
+            throw std::invalid_argument( "integrate() : start_time > end_time" );
 
-	observer( start_time , state , system );
-	while( start_time < end_time )
-	{
-	    stepper.next_step( system , state , start_time , dt );
-	    start_time += dt;
-	    observer( start_time , state , system );
-	}
+        observer( start_time , state , system );
+        while( start_time < end_time ) {
+            stepper.next_step( system , state , start_time , dt );
+            start_time += dt;
+            observer( start_time , state , system );
+        }
     }
     
-
-
-    
+   
 
 } // odeint
 } // numeric

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-04 14:17:22 EST (Wed, 04 Nov 2009)
@@ -35,10 +35,8 @@
 
 const size_t olen = 10000;
 
-const double eps_abs = 1E-7;
-const double eps_rel = 1E-8;
-
-const double min_dt = 1E-10;
+const double eps_abs = 1E-6;
+const double eps_rel = 1E-7;
 
 typedef array<double, 3> state_type;
 
@@ -49,6 +47,11 @@
     dxdt[2] = x[0]*x[1] - b * x[2];
 }
 
+void print_state( double t, state_type &x, ... )
+{
+    cout << t << tab << x[0] << tab << x[1] << tab << x[2] << endl;
+}
+
 int main( int argc , char **argv )
 {
     state_type x;
@@ -57,28 +60,15 @@
     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);
+    step_controller_standard< state_type, double > 
+        controller(eps_abs, eps_rel, 1.0, 1.0);
     
-    double t = 0.0;
-    double dt = 1E-4;
-    controlled_step_result result;
-
     cout.precision(5);
     cout.setf(ios::fixed,ios::floatfield);
     
+    size_t steps = integrate( euler, lorenz, controller, 0.0, 1E-4, x, 10.0, print_state );
 
-    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;
-    }
+    clog << "Number of steps: " << steps << endl;
 
     return 0;
 }

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-04 14:17:22 EST (Wed, 04 Nov 2009)
@@ -37,7 +37,7 @@
 const double eps_abs = 1E-3;
 const double eps_rel = 1E-3;
 
-const size_t time_points = 100;
+const size_t time_points = 101;
 
 typedef array<double, 3> state_type;
 
@@ -58,7 +58,7 @@
     vector<state_type> x_t_vec(time_points);
     vector<double> times(time_points);
     for( size_t i=0; i<time_points; i++ ) {
-	times[i] = 0.1*i;
+        times[i] = 0.1*i;
     }
 
     ode_step_euler< state_type > euler;
@@ -71,9 +71,9 @@
     
 
     for( size_t i=0; i<time_points; i++ ) {
-	//cout << "current state: " ;
-	cout << times[i] << tab;
-	cout << x_t_vec[i][0] << tab << x_t_vec[i][1] << tab << x_t_vec[i][2] << endl;
+        //cout << "current state: " ;
+        cout << times[i] << tab;
+        cout << x_t_vec[i][0] << tab << x_t_vec[i][1] << tab << x_t_vec[i][2] << endl;
     }
 
     return 0;