Boost logo

Boost-Commit :

Subject: [Boost-commit] svn:boost r66355 - in sandbox/odeint/branches/karsten: boost/numeric/odeint/stepper boost/numeric/odeint/stepper/base libs/numeric/odeint/test
From: karsten.ahnert_at_[hidden]
Date: 2010-11-02 10:39:31


Author: karsten
Date: 2010-11-02 10:39:28 EDT (Tue, 02 Nov 2010)
New Revision: 66355
URL: http://svn.boost.org/trac/boost/changeset/66355

Log:
* interface changes for the fsal steppers
* some test routines for dense output steppers and controlled steppers
Added:
   sandbox/odeint/branches/karsten/libs/numeric/odeint/test/controlled_stepper_evolution.cpp (contents, props changed)
   sandbox/odeint/branches/karsten/libs/numeric/odeint/test/dense_output_stepper_evolution.cpp (contents, props changed)
Text files modified:
   sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/explicit_stepper_and_error_stepper_fsal_base.hpp | 50 +++++++++++++----------
   sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/controlled_error_stepper.hpp | 46 +++++++++++----------
   sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/dense_output_dopri5.hpp | 84 ++++++++++++++++++++++++---------------
   sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_error_dopri5.hpp | 66 +++++++-----------------------
   sandbox/odeint/branches/karsten/libs/numeric/odeint/test/Jamfile | 23 ++++------
   5 files changed, 131 insertions(+), 138 deletions(-)

Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/explicit_stepper_and_error_stepper_fsal_base.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/explicit_stepper_and_error_stepper_fsal_base.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/explicit_stepper_and_error_stepper_fsal_base.hpp 2010-11-02 10:39:28 EDT (Tue, 02 Nov 2010)
@@ -95,32 +95,41 @@
         template< class System >
         void do_step( System &system , state_type &x , const time_type t , const time_type dt )
         {
- m_size_adjuster.adjust_size_by_policy( x , adjust_size_policy() );
- system( x , m_dxdt ,t );
- this->stepper().do_step_impl( system , x , m_dxdt , t , x , dt );
+ if( m_size_adjuster.adjust_size_by_policy( x , adjust_size_policy() ) || m_first_call )
+ {
+ system( x , m_dxdt ,t );
+ m_first_call = false;
+ }
+ this->stepper().do_step_impl( system , x , m_dxdt , t , x , m_dxdt , dt );
         }
 
         // do_step( sys , x , dxdt , t , dt )
         template< class System >
- void do_step( System &system , state_type &x , const state_type &dxdt , const time_type t , const time_type dt )
+ void do_step( System &system , state_type &x , state_type &dxdt , const time_type t , const time_type dt )
         {
- this->stepper().do_step_impl( system , x , dxdt , t , x , dt );
+ m_first_call = true;
+ this->stepper().do_step_impl( system , x , dxdt , t , x , dxdt , dt );
         }
 
         // do_step( sys , in , t , out , dt )
         template< class System >
         void do_step( System &system , const state_type &in , const time_type t , state_type &out , const time_type dt )
         {
- m_size_adjuster.adjust_size_by_policy( in , adjust_size_policy() );
- system( in , m_dxdt ,t );
- this->stepper().do_step_impl( system , in , m_dxdt , t , out , dt );
+ if( m_size_adjuster.adjust_size_by_policy( in , adjust_size_policy() ) || m_first_call )
+ {
+ system( in , m_dxdt ,t );
+ m_first_call = false;
+ }
+ this->stepper().do_step_impl( system , in , m_dxdt , t , out , m_dxdt , dt );
         }
 
- // do_step( sys , in , dxdt , t , out , dt )
+ // do_step( sys , in , dxdt_in , t , out , dxdt_out , dt )
         template< class System >
- void do_step( System &system , const state_type &in , const state_type &dxdt , const time_type t , state_type &out , const time_type dt )
+ void do_step( System &system , const state_type &in , const state_type &dxdt_in , const time_type t ,
+ state_type &out , state_type &dxdt_out , const time_type dt )
         {
- this->stepper().do_step_impl( system , in , dxdt , t , out , dt );
+ m_first_call = true;
+ this->stepper().do_step_impl( system , in , dxdt_in , t , out , dxdt_out , dt );
         }
 
 
@@ -137,14 +146,15 @@
                 system( x , m_dxdt ,t );
                 m_first_call = false;
             }
- this->stepper().do_step_impl( system , x , m_dxdt , t , x , dt , xerr );
+ this->stepper().do_step_impl( system , x , m_dxdt , t , x , m_dxdt , dt , xerr );
         }
 
         // do_step( sys , x , dxdt , t , dt , xerr )
         template< class System >
         void do_step( System &system , state_type &x , state_type &dxdt , const time_type t , const time_type dt , state_type &xerr )
         {
- this->stepper().do_step_impl( system , x , dxdt , t , x , dt , xerr );
+ m_first_call = true;
+ this->stepper().do_step_impl( system , x , dxdt , t , x , dxdt , dt , xerr );
         }
 
         // do_step( sys , in , t , out , dt , xerr )
@@ -156,22 +166,20 @@
                 system( in , m_dxdt ,t );
                 m_first_call = false;
             }
- this->stepper().do_step_impl( system , in , m_dxdt , t , out , dt , xerr );
+ this->stepper().do_step_impl( system , in , m_dxdt , t , out , m_dxdt , dt , xerr );
         }
 
- // do_step( sys , in , dxdt , t , out , dt , xerr )
+ // do_step( sys , in , dxdt_in , t , out , dxdt_out , dt )
         template< class System >
- void do_step( System &system , const state_type &in , state_type &dxdt , const time_type t , state_type &out , const time_type dt , state_type &xerr )
+ void do_step( System &system , const state_type &in , const state_type &dxdt_in , const time_type t ,
+ state_type &out , state_type &dxdt_out , const time_type dt , state_type &xerr )
         {
- this->stepper().do_step_impl( system , in , dxdt , t , out , dt , xerr );
+ m_first_call = true;
+ this->stepper().do_step_impl( system , in , dxdt_in , t , out , dxdt_out , dt , xerr );
         }
 
 
 
- void reset_step( state_type &dxdt )
- {
- this->stepper().reset_step_impl( dxdt );
- }
 
 
 

Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/controlled_error_stepper.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/controlled_error_stepper.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/controlled_error_stepper.hpp 2010-11-02 10:39:28 EDT (Tue, 02 Nov 2010)
@@ -238,16 +238,18 @@
             const error_checker_type &error_checker = error_checker_type()
             )
     : m_stepper( stepper ) , m_error_checker( error_checker ) ,
- m_dxdt_size_adjuster() , m_xerr_size_adjuster() , m_xnew_size_adjuster() ,
- m_dxdt() , m_xerr() , m_xnew() ,
+ m_dxdt_size_adjuster() , m_xerr_size_adjuster() , m_new_size_adjuster() ,
+ m_dxdt() , m_xerr() , m_xnew() , m_dxdtnew() ,
       m_first_call( true )
     {
         boost::numeric::odeint::construct( m_dxdt );
         boost::numeric::odeint::construct( m_xerr );
         boost::numeric::odeint::construct( m_xnew );
+ boost::numeric::odeint::construct( m_dxdtnew );
         m_dxdt_size_adjuster.register_state( 0 , m_dxdt );
         m_xerr_size_adjuster.register_state( 0 , m_xerr );
- m_xnew_size_adjuster.register_state( 0 , m_xnew );
+ m_new_size_adjuster.register_state( 0 , m_xnew );
+ m_new_size_adjuster.register_state( 1 , m_dxdtnew );
     }
 
     ~controlled_error_stepper( void )
@@ -255,6 +257,7 @@
         boost::numeric::odeint::destruct( m_dxdt );
         boost::numeric::odeint::destruct( m_xerr );
         boost::numeric::odeint::destruct( m_xnew );
+ boost::numeric::odeint::destruct( m_dxdtnew );
     }
 
 
@@ -273,19 +276,6 @@
         return try_step( sys , x , m_dxdt , t , dt );
     }
 
- // try_step( sys , x , dxdt , t , dt )
- template< class System >
- controlled_step_result try_step( System &sys , state_type &x , state_type &dxdt , time_type &t , time_type &dt )
- {
- m_xnew_size_adjuster.adjust_size_by_policy( x , adjust_size_policy() );
- controlled_step_result res = try_step( sys , x , dxdt , t , m_xnew , dt );
- if( ( res == success_step_size_increased ) || ( res == success_step_size_unchanged) )
- {
- boost::numeric::odeint::copy( m_xnew , x );
- }
- return res;
- }
-
     // try_step( sys , in , t , out , dt );
     template< class System >
     controlled_step_result try_step( System &sys , const state_type &in , time_type &t , state_type &out , time_type &dt )
@@ -298,10 +288,24 @@
         return try_step( sys , in , m_dxdt , t , out , dt );
     }
 
+ // try_step( sys , x , dxdt , t , dt )
+ template< class System >
+ controlled_step_result try_step( System &sys , state_type &x , state_type &dxdt , time_type &t , time_type &dt )
+ {
+ m_new_size_adjuster.adjust_size_by_policy( x , adjust_size_policy() );
+ controlled_step_result res = try_step( sys , x , dxdt , t , m_xnew , m_dxdtnew , dt );
+ if( ( res == success_step_size_increased ) || ( res == success_step_size_unchanged) )
+ {
+ boost::numeric::odeint::copy( m_xnew , x );
+ boost::numeric::odeint::copy( m_dxdtnew , dxdt );
+ }
+ return res;
+ }
 
     // try_step( sys , in , dxdt , t , out , dt )
     template< class System >
- controlled_step_result try_step( System &sys , const state_type &in , state_type &dxdt , time_type &t , state_type &out , time_type &dt )
+ controlled_step_result try_step( System &sys , const state_type &in , const state_type &dxdt_in , time_type &t ,
+ state_type &out , state_type &dxdt_out , time_type &dt )
     {
         using std::max;
         using std::min;
@@ -311,16 +315,15 @@
 
         //fsal: m_stepper.get_dxdt( dxdt );
         //fsal: m_stepper.do_step( sys , x , dxdt , t , dt , m_x_err );
- m_stepper.do_step( sys , in , dxdt , t , out , dt , m_xerr );
+ m_stepper.do_step( sys , in , dxdt_in , t , out , dxdt_out , dt , m_xerr );
 
         // this potentially overwrites m_x_err! (standard_error_checker does, at least)
- time_type max_rel_err = m_error_checker.error( in , dxdt , m_xerr , dt );
+ time_type max_rel_err = m_error_checker.error( in , dxdt_in , m_xerr , dt );
 
         if( max_rel_err > 1.1 )
         {
             // error too large - decrease dt ,limit scaling factor to 0.2 and reset state
             dt *= max( 0.9 * pow( max_rel_err , -1.0 / ( m_stepper.error_order() - 1.0 ) ) , 0.2 );
- m_stepper.reset_step( dxdt );
             return step_size_decreased;
         }
         else
@@ -375,11 +378,12 @@
 
     size_adjuster< state_type , 1 > m_dxdt_size_adjuster;
     size_adjuster< state_type , 1 > m_xerr_size_adjuster;
- size_adjuster< state_type , 1 > m_xnew_size_adjuster;
+ size_adjuster< state_type , 2 > m_new_size_adjuster;
 
     state_type m_dxdt;
     state_type m_xerr;
     state_type m_xnew;
+ state_type m_dxdtnew;
     bool m_first_call;
 };
 

Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/dense_output_dopri5.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/dense_output_dopri5.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/dense_output_dopri5.hpp 2010-11-02 10:39:28 EDT (Tue, 02 Nov 2010)
@@ -12,12 +12,12 @@
 #ifndef BOOST_BOOST_NUMERIC_ODEINT_DENSE_OUTPUT_DOPRI5_HPP_INCLUDED
 #define BOOST_BOOST_NUMERIC_ODEINT_DENSE_OUTPUT_DOPRI5_HPP_INCLUDED
 
-#include <iostream>
-#define tab "\t"
-using std::cout;
-using std::cerr;
-using std::clog;
-using std::endl;
+//#include <iostream>
+//#define tab "\t"
+//using std::cout;
+//using std::cerr;
+//using std::clog;
+//using std::endl;
 
 #include <stdexcept>
 
@@ -52,25 +52,35 @@
 
         dense_output_dopri5( stepper_type &stepper )
         : m_stepper( stepper ) , m_size_adjuster() ,
- m_x1() , m_x2() , m_current_state( &m_x1 ) , m_old_state( &m_x2 ) ,
- m_t( 0.0 ) , m_t_old( 0.0 ) , m_dt( 1.0 )
+ m_x1() , m_x2() , m_dxdt1() , m_dxdt2() ,
+ m_current_state( &m_x1 ) , m_old_state( &m_x2 ) ,
+ m_current_deriv( &m_dxdt1 ) , m_old_deriv( &m_dxdt2 ) ,
+ m_t( 0.0 ) , m_t_old( 0.0 ) , m_dt( 1.0 ) ,
+ m_is_deriv_initialized( false )
         {
                 boost::numeric::odeint::construct( m_x1 );
                 boost::numeric::odeint::construct( m_x2 );
+ boost::numeric::odeint::construct( m_dxdt1 );
+ boost::numeric::odeint::construct( m_dxdt2 );
                 m_size_adjuster.register_state( 0 , m_x1 );
                 m_size_adjuster.register_state( 1 , m_x2 );
+ m_size_adjuster.register_state( 1 , m_dxdt1);
+ m_size_adjuster.register_state( 1 , m_dxdt2 );
         }
 
         ~dense_output_dopri5( void )
         {
                 boost::numeric::odeint::destruct( m_x1 );
                 boost::numeric::odeint::destruct( m_x2 );
+ boost::numeric::odeint::destruct( m_dxdt1 );
+ boost::numeric::odeint::destruct( m_dxdt2 );
         }
 
         void adjust_size( const state_type &x )
         {
                 m_size_adjuster.adjust_size( x );
                 m_stepper.adjust_size( x );
+ m_is_deriv_initialized = false;
         }
 
         void initialize( const state_type &x0 , const time_type t0 , const time_type dt0 )
@@ -78,23 +88,29 @@
                 boost::numeric::odeint::copy( x0 , *m_current_state );
                 m_t = t0;
                 m_dt = dt0;
+ m_is_deriv_initialized = false;
         }
 
         template< class System >
         std::pair< time_type , time_type > do_step( System &system )
         {
                 const size_t max_count = 1000;
- controlled_step_result res;
+
+ if( !m_is_deriv_initialized )
+ system( *m_current_state , *m_current_deriv , m_t );
+
+ controlled_step_result res = step_size_decreased;
                 m_t_old = m_t;
                 size_t count = 0;
                 do
                 {
- res = m_stepper.try_step( system , *m_current_state , m_t , *m_old_state , m_dt );
+ res = m_stepper.try_step( system , *m_current_state , *m_current_deriv , m_t , *m_old_state , *m_old_deriv , m_dt );
                         if( count++ == max_count )
                                 throw std::overflow_error( "dense_output_dopri5 : too much iterations!");
                 }
                 while( res == step_size_decreased );
                 std::swap( m_current_state , m_old_state );
+ std::swap( m_current_deriv , m_old_deriv );
                 return std::make_pair( m_t_old , m_t );
 
         }
@@ -122,27 +138,27 @@
         void calc_state( time_type t , state_type &x )
         {
 
- const time_type b1 = static_cast<time_type> ( 35.0 ) / static_cast<time_type>( 384.0 );
- const time_type b3 = static_cast<time_type> ( 500.0 ) / static_cast<time_type>( 1113.0 );
- const time_type b4 = static_cast<time_type> ( 125.0 ) / static_cast<time_type>( 192.0 );
- const time_type b5 = static_cast<time_type> ( -2187.0 ) / static_cast<time_type>( 6784.0 );
- const time_type b6 = static_cast<time_type> ( 11.0 ) / static_cast<time_type>( 84.0 );
-
-
- time_type theta = t - m_t_old;
- time_type X1 = static_cast< time_type >( 5.0 ) * ( static_cast< time_type >( 2558722523.0 ) - static_cast< time_type >( 31403016.0 ) * theta ) / static_cast< time_type >( 11282082432.0 );
- time_type X3 = static_cast< time_type >( 100.0 ) * ( static_cast< time_type >( 882725551.0 ) - static_cast< time_type >( 15701508.0 ) * theta ) / static_cast< time_type >( 32700410799.0 );
- time_type X4 = static_cast< time_type >( 25.0 ) * ( static_cast< time_type >( 443332067.0 ) - static_cast< time_type >( 31403016.0 ) * theta ) / static_cast< time_type >( 1880347072.0 ) ;
- time_type X5 = static_cast< time_type >( 32805.0 ) * ( static_cast< time_type >( 23143187.0 ) - static_cast< time_type >( 3489224.0 ) * theta ) / static_cast< time_type >( 199316789632.0 );
- time_type X6 = static_cast< time_type >( 55.0 ) * ( static_cast< time_type >( 29972135.0 ) - static_cast< time_type >( 7076736.0 ) * theta ) / static_cast< time_type >( 822651844.0 );
- time_type X7 = static_cast< time_type >( 10.0 ) * ( static_cast< time_type >( 7414447.0 ) - static_cast< time_type >( 829305.0 ) * theta ) / static_cast< time_type >( 29380423.0 );
-
- time_type theta_m_1 = theta - static_cast< time_type >( 1.0 );
- time_type theta_sq = theta * theta;
- time_type A = theta_sq * ( static_cast< time_type >( 3.0 ) - static_cast< time_type >( 2.0 ) * theta );
- time_type B = theta_sq * theta_m_1;
- time_type C = theta_sq * theta_m_1 * theta_m_1;
- time_type D = theta * theta_m_1 * theta_m_1;
+// const time_type b1 = static_cast<time_type> ( 35.0 ) / static_cast<time_type>( 384.0 );
+// const time_type b3 = static_cast<time_type> ( 500.0 ) / static_cast<time_type>( 1113.0 );
+// const time_type b4 = static_cast<time_type> ( 125.0 ) / static_cast<time_type>( 192.0 );
+// const time_type b5 = static_cast<time_type> ( -2187.0 ) / static_cast<time_type>( 6784.0 );
+// const time_type b6 = static_cast<time_type> ( 11.0 ) / static_cast<time_type>( 84.0 );
+//
+//
+// time_type theta = t - m_t_old;
+// time_type X1 = static_cast< time_type >( 5.0 ) * ( static_cast< time_type >( 2558722523.0 ) - static_cast< time_type >( 31403016.0 ) * theta ) / static_cast< time_type >( 11282082432.0 );
+// time_type X3 = static_cast< time_type >( 100.0 ) * ( static_cast< time_type >( 882725551.0 ) - static_cast< time_type >( 15701508.0 ) * theta ) / static_cast< time_type >( 32700410799.0 );
+// time_type X4 = static_cast< time_type >( 25.0 ) * ( static_cast< time_type >( 443332067.0 ) - static_cast< time_type >( 31403016.0 ) * theta ) / static_cast< time_type >( 1880347072.0 ) ;
+// time_type X5 = static_cast< time_type >( 32805.0 ) * ( static_cast< time_type >( 23143187.0 ) - static_cast< time_type >( 3489224.0 ) * theta ) / static_cast< time_type >( 199316789632.0 );
+// time_type X6 = static_cast< time_type >( 55.0 ) * ( static_cast< time_type >( 29972135.0 ) - static_cast< time_type >( 7076736.0 ) * theta ) / static_cast< time_type >( 822651844.0 );
+// time_type X7 = static_cast< time_type >( 10.0 ) * ( static_cast< time_type >( 7414447.0 ) - static_cast< time_type >( 829305.0 ) * theta ) / static_cast< time_type >( 29380423.0 );
+//
+// time_type theta_m_1 = theta - static_cast< time_type >( 1.0 );
+// time_type theta_sq = theta * theta;
+// time_type A = theta_sq * ( static_cast< time_type >( 3.0 ) - static_cast< time_type >( 2.0 ) * theta );
+// time_type B = theta_sq * theta_m_1;
+// time_type C = theta_sq * theta_m_1 * theta_m_1;
+// time_type D = theta * theta_m_1 * theta_m_1;
 
 // time_type b1_theta = A *
 
@@ -164,10 +180,12 @@
 
 
         stepper_type &m_stepper;
- size_adjuster< state_type , 2 > m_size_adjuster;
- state_type m_x1 , m_x2;
+ size_adjuster< state_type , 4 > m_size_adjuster;
+ state_type m_x1 , m_x2 , m_dxdt1 , m_dxdt2;
         state_type *m_current_state , *m_old_state;
+ state_type *m_current_deriv , *m_old_deriv;
         time_type m_t , m_t_old , m_dt;
+ bool m_is_deriv_initialized;
 
 
 };

Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_error_dopri5.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_error_dopri5.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_error_dopri5.hpp 2010-11-02 10:39:28 EDT (Tue, 02 Nov 2010)
@@ -78,7 +78,8 @@
 
 
         template< class System >
- void do_step_impl( System system , const state_type &in , const state_type &dxdt , const time_type t , state_type &out , const time_type dt )
+ void do_step_impl( System system , const state_type &in , const state_type &dxdt_in , const time_type t ,
+ state_type &out , state_type &dxdt_out , const time_type dt )
         {
             m_size_adjuster.adjust_size_by_policy( in , adjust_size_policy() );
 
@@ -114,34 +115,39 @@
         const time_type c6 = static_cast<time_type> ( 11.0 ) / static_cast<time_type>( 84.0 );
 
         //m_x1 = x + dt*b21*dxdt
- algebra_type::for_each3( m_x1 , in , dxdt ,
+ algebra_type::for_each3( m_x1 , in , dxdt_in ,
                     typename operations_type::scale_sum2( 1.0 , dt*b21 ) );
 
         system( m_x1 , m_x2 , t + dt*a2 );
         // m_x1 = x + dt*b31*dxdt + dt*b32*m_x2
- algebra_type::for_each4( m_x1 , in , dxdt , m_x2 ,
+ algebra_type::for_each4( m_x1 , in , dxdt_in , m_x2 ,
                     typename operations_type::scale_sum3( 1.0 , dt*b31 , dt*b32 ));
 
         system( m_x1 , m_x3 , t + dt*a3 );
         // m_x1 = x + dt * (b41*dxdt + b42*m_x2 + b43*m_x3)
- algebra_type::for_each5( m_x1 , in , dxdt , m_x2 , m_x3 ,
+ algebra_type::for_each5( m_x1 , in , dxdt_in , m_x2 , m_x3 ,
                     typename operations_type::scale_sum4( 1.0 , dt*b41 , dt*b42 , dt*b43 ));
 
         system( m_x1, m_x4 , t + dt*a4 );
- algebra_type::for_each6( m_x1 , in , dxdt , m_x2 , m_x3 , m_x4 ,
+ algebra_type::for_each6( m_x1 , in , dxdt_in , m_x2 , m_x3 , m_x4 ,
                     typename operations_type::scale_sum5( 1.0 , dt*b51 , dt*b52 , dt*b53 , dt*b54 ));
 
         system( m_x1 , m_x5 , t + dt*a5 );
- algebra_type::for_each7( m_x1 , in , dxdt , m_x2 , m_x3 , m_x4 , m_x5 ,
+ algebra_type::for_each7( m_x1 , in , dxdt_in , m_x2 , m_x3 , m_x4 , m_x5 ,
                             typename operations_type::scale_sum6( 1.0 , dt*b61 , dt*b62 , dt*b63 , dt*b64 , dt*b65 ));
 
         system( m_x1 , m_x6 , t + dt );
- algebra_type::for_each7( out , in , dxdt , m_x3 , m_x4 , m_x5 , m_x6 ,
+ algebra_type::for_each7( out , in , dxdt_in , m_x3 , m_x4 , m_x5 , m_x6 ,
                     typename operations_type::scale_sum6( 1.0 , dt*c1 , dt*c3 , dt*c4 , dt*c5 , dt*c6 ));
+
+ // the new derivative
+ system( out , dxdt_out , t + dt );
         }
 
+
         template< class System >
- void do_step_impl( System system , const state_type &in , state_type &dxdt , const time_type t , state_type &out , const time_type dt , state_type &xerr )
+ void do_step_impl( System system , const state_type &in , const state_type &dxdt_in , const time_type t ,
+ state_type &out , state_type &dxdt_out , const time_type dt , state_type &xerr )
         {
 
         const time_type c1 = static_cast<time_type> ( 35.0 ) / static_cast<time_type>( 384.0 );
@@ -157,56 +163,16 @@
         const time_type dc6 = c6 - static_cast<time_type> ( 187.0 ) / static_cast<time_type>( 2100.0 );
         const time_type dc7 = static_cast<time_type>( -0.025 );
 
- do_step_impl( system , in , dxdt , t , out , dt );
-
- // we need to copy the old dxdt
- boost::numeric::odeint::copy( dxdt , m_x1 );
-
- // store the new result in dxdt
- system( out , dxdt , t + dt );
+ do_step_impl( system , in , dxdt_in , t , out , dxdt_out , dt );
 
         //error estimate
- algebra_type::for_each7( xerr , m_x1 , m_x3 , m_x4 , m_x5 , m_x6 , dxdt ,
+ algebra_type::for_each7( xerr , dxdt_in , m_x3 , m_x4 , m_x5 , m_x6 , dxdt_out ,
                     typename operations_type::scale_sum6( dt*dc1 , dt*dc3 , dt*dc4 , dt*dc5 , dt*dc6 , dt*dc7 ) );
         }
 
-// template< class System >
-// void do_step_impl( System system , const state_type &in , const state_type &dxdt_in , const time_type t , state_type &out , state_type &dxdt_out , const time_type dt , state_type &xerr )
-// {
-//
-// const time_type c1 = static_cast<time_type> ( 35.0 ) / static_cast<time_type>( 384.0 );
-// const time_type c3 = static_cast<time_type> ( 500.0 ) / static_cast<time_type>( 1113.0 );
-// const time_type c4 = static_cast<time_type> ( 125.0 ) / static_cast<time_type>( 192.0 );
-// const time_type c5 = static_cast<time_type> ( -2187.0 ) / static_cast<time_type>( 6784.0 );
-// const time_type c6 = static_cast<time_type> ( 11.0 ) / static_cast<time_type>( 84.0 );
-//
-// const time_type dc1 = c1 - static_cast<time_type> ( 5179.0 ) / static_cast<time_type>( 57600.0 );
-// const time_type dc3 = c3 - static_cast<time_type> ( 7571.0 ) / static_cast<time_type>( 16695.0 );
-// const time_type dc4 = c4 - static_cast<time_type> ( 393.0 ) / static_cast<time_type>( 640.0 );
-// const time_type dc5 = c5 - static_cast<time_type> ( -92097.0 ) / static_cast<time_type>( 339200.0 );
-// const time_type dc6 = c6 - static_cast<time_type> ( 187.0 ) / static_cast<time_type>( 2100.0 );
-// const time_type dc7 = static_cast<time_type>( -0.025 );
-//
-// do_step_impl( system , in , dxdt_in , t , out , dt );
-//
-// // store the new result in dxdt
-// system( out , dxdt_out , t + dt );
-//
-// //error estimate
-// algebra_type::for_each7( xerr , m_x1 , m_x3 , m_x4 , m_x5 , m_x6 , dxdt_out ,
-// typename operations_type::scale_sum6( dt*dc1 , dt*dc3 , dt*dc4 , dt*dc5 , dt*dc6 , dt*dc7 ) );
-// }
-
 
 
 
-
- void reset_step_impl( state_type &dxdt )
- {
- boost::numeric::odeint::copy( m_x1 , dxdt );
- }
-
-
         void adjust_size( const state_type &x )
         {
                 m_size_adjuster.adjust_size( x );

Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/test/Jamfile
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/test/Jamfile (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/test/Jamfile 2010-11-02 10:39:28 EDT (Tue, 02 Nov 2010)
@@ -6,18 +6,10 @@
 
 import testing ;
 
-#project
-# : requirements
-# <include>../../../..
-# <include>$(BOOST_ROOT)
-# <define>BOOST_ALL_NO_LIB=1
-# : build-dir .
-# ;
 
 project
     : requirements
       <define>BOOST_ALL_NO_LIB=1
- <library>/boost/test//boost_unit_test_framework
       <include>../../../..
       <include>$(BOOST_ROOT)
     ;
@@ -31,9 +23,14 @@
           check_dense_output_explicit_euler.cpp
           check_dense_output_dopri5.cpp
         :
- : <link>static
+ : <library>/boost/test//boost_unit_test_framework
+ <link>static
         ;
-
-
-
-
+
+exe controlled_stepper_evolution
+ : controlled_stepper_evolution.cpp
+ ;
+
+exe dense_output_stepper_evolution
+ : dense_output_stepper_evolution.cpp
+ ;
\ No newline at end of file

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/test/controlled_stepper_evolution.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/test/controlled_stepper_evolution.cpp 2010-11-02 10:39:28 EDT (Tue, 02 Nov 2010)
@@ -0,0 +1,77 @@
+/*
+ * controlled_stepper_evolution.cpp
+ *
+ * Created on: Nov 2, 2010
+ * Author: karsten
+ */
+
+#include <string>
+#include <fstream>
+#include <iostream>
+#include <tr1/array>
+
+#include <boost/numeric/odeint.hpp>
+
+#define tab "\t"
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+
+typedef std::tr1::array< double , 2 > state_type;
+
+ostream& operator<<( ostream &out , const state_type &x )
+{
+ out << x[0] << tab << x[1];
+ return out;
+}
+
+void sys( const state_type &x , state_type &dxdt , double t )
+{
+ dxdt[0] = x[1];
+ dxdt[1] = -x[0];
+}
+
+state_type sys_solution( double t , const state_type &x0 )
+{
+ state_type sol = {{ 0.0 , 0.0 }};
+ sol[0] = x0[0] * cos( t ) + x0[1] * sin( t );
+ sol[1] = -x0[0] * sin( t ) + x0[1] * cos( t );
+ return sol;
+}
+
+
+
+template< class Stepper >
+void evolution( Stepper &stepper , double t_end , const state_type &x0 , const string &filename )
+{
+ ofstream fout( filename.c_str() );
+ state_type x = x0;
+ double t = 0.0 , dt = 0.01;
+ while( t < t_end )
+ {
+ state_type orig = sys_solution( t , x0 );
+ state_type diff = {{ orig[0] - x[0] , orig[1] - x[1] }};
+ double diff_abs = sqrt( diff[0] * diff[0] + diff[1] * diff[1] );
+ fout << t << tab << orig << tab << x << tab << diff << tab << diff_abs << endl;
+ stepper.try_step( sys , x , t , dt );
+ }
+}
+
+
+int main( int argc , char **argv )
+{
+ typedef explicit_error_rk54_ck< state_type > rk54_type;
+ typedef explicit_error_dopri5< state_type > dopri5_type;
+
+ rk54_type rk54;
+ dopri5_type dopri5;
+ controlled_error_stepper< rk54_type > rk54_controlled( rk54 );
+ controlled_error_stepper< dopri5_type > dopri5_controlled( dopri5 );
+
+ state_type x0 = {{ 1.25 , 0.43 }};
+ evolution( rk54_controlled , 100.0 , x0 , "rk54.dat" );
+ evolution( dopri5_controlled , 100.0 , x0 , "dopri5.dat" );
+
+ return 0;
+}

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/test/dense_output_stepper_evolution.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/test/dense_output_stepper_evolution.cpp 2010-11-02 10:39:28 EDT (Tue, 02 Nov 2010)
@@ -0,0 +1,92 @@
+/*
+ * dense_output_stepper_evolution.cpp
+ *
+ * Created on: Nov 2, 2010
+ * Author: karsten
+ */
+
+#include <string>
+#include <fstream>
+#include <iostream>
+#include <tr1/array>
+
+#include <boost/numeric/odeint.hpp>
+
+#define tab "\t"
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+typedef std::tr1::array< double , 2 > state_type;
+
+ostream& operator<<( ostream &out , const state_type &x )
+{
+ out << x[0] << tab << x[1];
+ return out;
+}
+
+void sys( const state_type &x , state_type &dxdt , double t )
+{
+ dxdt[0] = x[1];
+ dxdt[1] = -x[0];
+}
+
+state_type sys_solution( double t , const state_type &x0 )
+{
+ state_type sol = {{ 0.0 , 0.0 }};
+ sol[0] = x0[0] * cos( t ) + x0[1] * sin( t );
+ sol[1] = -x0[0] * sin( t ) + x0[1] * cos( t );
+ return sol;
+}
+
+
+
+template< class Stepper >
+void evolution( Stepper &stepper , double t_end , const state_type &x0 , const string &stepper_filename , const string &state_filename )
+{
+ ofstream state_out( state_filename.c_str() );
+ ofstream stepper_out( stepper_filename.c_str() );
+
+ stepper.initialize( x0 , 0.0 , 0.01 );
+
+ state_type x = x0;
+ double t = 0.0;
+ double dt = 0.001;
+ while( t < t_end )
+ {
+ if( t < stepper.current_time() )
+ {
+ stepper.calc_state( t , x );
+ state_type orig = sys_solution( t , x0 );
+ state_type diff = {{ orig[0] - x[0] , orig[1] - x[1] }};
+ double diff_abs = sqrt( diff[0] * diff[0] + diff[1] * diff[1] );
+ state_out << t << tab << x << tab << orig << tab << diff << diff_abs << endl;
+ }
+ else
+ {
+ stepper.do_step( sys );
+ stepper_out << stepper.current_time() << "\t" << stepper.current_state() << std::endl;
+ continue;
+ }
+ t += dt;
+
+ }
+}
+
+
+int main( int argc , char **argv )
+{
+ typedef explicit_error_dopri5< state_type > dopri5_type;
+ typedef controlled_error_stepper< dopri5_type > controlled_dopri5_type;
+
+ dopri5_type dopri5;
+ controlled_dopri5_type controlled_dopri5( dopri5 );
+
+ dense_output_explicit_euler< state_type > dense_euler;
+ dense_output_dopri5< controlled_dopri5_type > dense_dopri5( controlled_dopri5 );
+
+ state_type x0 = {{ 1.25 , 0.43 }};
+ evolution( dense_euler , 100.0 , x0 , "dense_euler_stepper.dat" , "dense_euler_state.dat" );
+ evolution( dense_dopri5 , 100.0 , x0 , "dense_dopri5_stepper.dat" , "dense_dopri5_state.dat" );
+
+}


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