Boost logo

Boost-Commit :

Subject: [Boost-commit] svn:boost r60376 - in sandbox/odeint: . boost/numeric/odeint libs/numeric/odeint/test
From: karsten.ahnert_at_[hidden]
Date: 2010-03-09 03:09:38


Author: karsten
Date: 2010-03-09 03:09:36 EST (Tue, 09 Mar 2010)
New Revision: 60376
URL: http://svn.boost.org/trac/boost/changeset/60376

Log:
check the stepper concepts, order->order_step, adjust_size introduced...
Text files modified:
   sandbox/odeint/ToDo | 8 ++--
   sandbox/odeint/boost/numeric/odeint/controlled_stepper_bs.hpp | 8 ++--
   sandbox/odeint/boost/numeric/odeint/controlled_stepper_standard.hpp | 4 +-
   sandbox/odeint/boost/numeric/odeint/stepper_euler.hpp | 4 +-
   sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp | 4 +-
   sandbox/odeint/boost/numeric/odeint/stepper_midpoint.hpp | 74 ++++++++++++++++++++++++---------------
   sandbox/odeint/boost/numeric/odeint/stepper_rk4.hpp | 39 ++++++++++++++-------
   sandbox/odeint/boost/numeric/odeint/stepper_rk4_classical.hpp | 44 +++++++++++++++--------
   sandbox/odeint/libs/numeric/odeint/test/check_stepper_concepts.cpp | 27 ++++++++++++++
   9 files changed, 141 insertions(+), 71 deletions(-)

Modified: sandbox/odeint/ToDo
==============================================================================
--- sandbox/odeint/ToDo (original)
+++ sandbox/odeint/ToDo 2010-03-09 03:09:36 EST (Tue, 09 Mar 2010)
@@ -19,13 +19,13 @@
 * check: orders and adjust_size() in :
   * stepper_euler.hpp - DONE
   * stepper_half_step.hpp - DONE
- * stepper_midpoint.hpp
- * stepper_rk4_classical.hpp
- * stepper_rk4.hpp
+ * stepper_midpoint.hpp - DONE, changed stepcount -> step_number
+ * stepper_rk4_classical.hpp - DONE
+ * stepper_rk4.hpp - DONE
   * stepper_rk5_ck.hpp
   * stepper_rk78_fehlberg.hpp
   * stepper_rk_generic.hpp
-
+* die auskommentierten Iteratoren typedefs entfernen
 
 
 Mario:

Modified: sandbox/odeint/boost/numeric/odeint/controlled_stepper_bs.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/controlled_stepper_bs.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/controlled_stepper_bs.hpp 2010-03-09 03:09:36 EST (Tue, 09 Mar 2010)
@@ -162,12 +162,12 @@
             
             for( unsigned short k=0; k<=m_current_k_max; k++ )
             { // loop through interval numbers
- unsigned short stepcount = m_interval_sequence[k];
+ unsigned short step_number = m_interval_sequence[k];
                 //out-of-place midpoint step
- m_stepper_mp.set_stepcount(stepcount);
- m_stepper_mp.do_step(system, m_x0, dxdt, t, dt, m_x_mp);
+ m_stepper_mp.set_step_number(step_number);
+ m_stepper_mp.midpoint_step(system, m_x0, dxdt, t, dt, m_x_mp);
                 //std::clog << "x_mp: " << k << '\t' << m_x_mp[0] << '\t' << m_x_mp[1] << std::endl;
- time_type t_est = (dt/stepcount)*(dt/stepcount);
+ time_type t_est = (dt/step_number)*(dt/step_number);
                 extrapolate(k, t_est, m_x_mp, x, m_xerr);
                 //std::clog << "Error: " << k << '\t' << m_xerr[0] << '\t' << m_xerr[1] << std::endl;
                 if( k != 0 )

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 2010-03-09 03:09:36 EST (Tue, 09 Mar 2010)
@@ -71,8 +71,8 @@
         typedef typename stepper_type::time_type time_type;
         typedef typename stepper_type::traits_type traits_type;
         typedef typename stepper_type::value_type value_type;
- typedef typename stepper_type::iterator iterator;
- typedef typename stepper_type::const_iterator const_iterator;
+// typedef typename stepper_type::iterator iterator;
+// typedef typename stepper_type::const_iterator const_iterator;
 
 
         // private members

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 2010-03-09 03:09:36 EST (Tue, 09 Mar 2010)
@@ -45,8 +45,8 @@
         typedef Traits traits_type;
         typedef typename traits_type::container_type container_type;
         typedef typename traits_type::value_type value_type;
- typedef typename traits_type::iterator iterator;
- typedef typename traits_type::const_iterator const_iterator;
+// typedef typename traits_type::iterator iterator;
+// typedef typename traits_type::const_iterator const_iterator;
 
 
         //

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 2010-03-09 03:09:36 EST (Tue, 09 Mar 2010)
@@ -39,8 +39,8 @@
         typedef typename stepper_type::time_type time_type;
         typedef typename stepper_type::order_type order_type;
         typedef typename stepper_type::value_type value_type;
- typedef typename stepper_type::iterator iterator;
- typedef typename stepper_type::const_iterator const_iterator;
+// typedef typename stepper_type::iterator iterator;
+// typedef typename stepper_type::const_iterator const_iterator;
 
 
 

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 2010-03-09 03:09:36 EST (Tue, 09 Mar 2010)
@@ -44,8 +44,8 @@
         typedef Traits traits_type;
         typedef typename traits_type::container_type container_type;
         typedef typename traits_type::value_type value_type;
- typedef typename traits_type::iterator iterator;
- typedef typename traits_type::const_iterator const_iterator;
+// typedef typename traits_type::iterator iterator;
+// typedef typename traits_type::const_iterator const_iterator;
 
 
 
@@ -54,7 +54,7 @@
         // private memebers
     private:
 
- unsigned short m_stepcount;
+ unsigned short m_step_number;
 
         container_type m_x0;
         container_type m_x1;
@@ -65,23 +65,45 @@
 
     public:
 
- stepper_midpoint( unsigned short stepcount = 2 ) { }
-
- order_type order() const { return 2; }
+ order_type order_step() const { return 2; }
+
+ // standard constructor, the size of the internal container is not set
+ stepper_midpoint( unsigned short step_number = 2 )
+ {
+ set_step_number( step_number );
+ }
+
+ // constructor, which adjusts the size of the internal containers
+ stepper_midpoint( const container_type &x , unsigned short step_number = 2 )
+ {
+ adjust_size( x );
+ set_step_number( step_number );
+ }
+
+ // adjusts the size of the internal containers
+ void adjust_size( const container_type &x )
+ {
+ traits_type::adjust_size( x , m_x0 );
+ traits_type::adjust_size( x , m_x1 );
+ traits_type::adjust_size( x , m_dxdt );
+ }
 
- void set_stepcount( unsigned short stepcount )
+ void set_step_number( unsigned short step_number )
         {
- if( stepcount > 1 )
- m_stepcount = stepcount;
+ if( step_number > 1 )
+ m_step_number = step_number;
         }
 
- unsigned short get_stepcount() const { return m_stepcount; }
-
+ unsigned short get_step_number() const
+ {
+ return m_step_number;
+ }
 
 
 
+ // performs a midpoint step
         template< class DynamicalSystem >
- void do_step(
+ void midpoint_step(
                 DynamicalSystem &system ,
                 container_type &x ,
                 container_type &dxdt ,
@@ -89,34 +111,33 @@
                 time_type dt ,
                 container_type &x_out )
         {
+ using namespace detail::it_algebra;
+
             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 h = dt / static_cast<time_type>( m_step_number );
             const time_type h2 = static_cast<time_type>( 2.0 ) * h;
- time_type th = t + h;
 
- traits_type::adjust_size(x, m_x0);
- traits_type::adjust_size(x, m_x1);
- traits_type::adjust_size(x, m_dxdt);
 
- using namespace detail::it_algebra;
+ time_type th = t + h;
 
             // m_x1 = x + h*dxdt
- scale_sum( traits_type::begin(m_x1), traits_type::end(m_x1),
+ scale_sum( traits_type::begin(m_x1),
+ traits_type::end(m_x1),
                        t_1, traits_type::begin(x),
                        h, traits_type::begin(dxdt) );
             system( m_x1, m_dxdt, th );
 
             m_x0 = x;
             
- //std::clog<< "mp: " << 0 << '\t' << m_x1[0] << '\t' << m_x1[1] << '\t' << m_dxdt[0] << '\t' << m_dxdt[1] << std::endl;
-
             unsigned short i = 1;
- while( i != m_stepcount )
- { // general step
+ while( i != m_step_number )
+ {
+ // general step
                 //tmp = m_x1; m_x1 = m_x0 + h2*m_dxdt; m_x0 = tmp
- scale_sum_swap( traits_type::begin(m_x1), traits_type::end(m_x1),
+ scale_sum_swap( traits_type::begin(m_x1),
+ traits_type::end(m_x1),
                                 traits_type::begin(m_x0),
                                 h2, traits_type::begin(m_dxdt) );
                 th += h;
@@ -124,8 +145,6 @@
                 i++;
             }
 
- //std::clog<< "mp: " << i << '\t' << m_x1[0] << '\t' << m_x1[1] << '\t' << m_dxdt[0] << '\t' << m_dxdt[1] << std::endl;
-
             // last step
             // x = 0.5*( m_x0 + m_x1 + h*m_dxdt )
             scale_sum( traits_type::begin(x_out), traits_type::end(x_out),
@@ -146,7 +165,7 @@
                 time_type t ,
                 time_type dt )
         {
- do_step( system, x, dxdt, t, dt, x );
+ midpoint_step( system, x, dxdt, t, dt, x );
         }
 
 
@@ -160,7 +179,6 @@
                 time_type t ,
                 time_type dt )
         {
- traits_type::adjust_size(x, m_dxdt);
             system( x, m_dxdt, t );
             do_step( system , x, m_dxdt, t, dt );
         }

Modified: sandbox/odeint/boost/numeric/odeint/stepper_rk4.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepper_rk4.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/stepper_rk4.hpp 2010-03-09 03:09:36 EST (Tue, 09 Mar 2010)
@@ -39,8 +39,8 @@
         typedef Traits traits_type;
         typedef typename traits_type::container_type container_type;
         typedef typename traits_type::value_type value_type;
- typedef typename traits_type::iterator iterator;
- typedef typename traits_type::const_iterator const_iterator;
+// typedef typename traits_type::iterator iterator;
+// typedef typename traits_type::const_iterator const_iterator;
 
 
 
@@ -65,7 +65,26 @@
     public:
 
 
- order_type order() const { return 4; }
+ order_type order_step() const { return 4; }
+
+ // standard constructor, internal containers are not initialized
+ stepper_rk4( void )
+ {
+ }
+
+ // constructor, which adjusts the internal containers
+ stepper_rk4( const container_type &x )
+ {
+ }
+
+ void adjust_size( const container_type &x )
+ {
+ traits_type::adjust_size( x , m_dxdt );
+ traits_type::adjust_size( x , m_dxt );
+ traits_type::adjust_size( x , m_dxm );
+ traits_type::adjust_size( x , m_xt );
+ traits_type::adjust_size( x , m_dxh );
+ }
 
 
 
@@ -80,11 +99,6 @@
 
             const time_type val1 = static_cast<time_type>( 1.0 );
 
- traits_type::adjust_size( x , m_dxt );
- traits_type::adjust_size( x , m_dxm );
- traits_type::adjust_size( x , m_xt );
- traits_type::adjust_size( x , m_dxh );
-
             time_type dh = static_cast<time_type>( 0.5 ) * dt;
             time_type th = t + dh;
 
@@ -95,14 +109,16 @@
                        val1, traits_type::begin(x),
                        dh, traits_type::begin(dxdt) );
 
+
             // dt * m_dxt = k2
             system( m_xt , m_dxt , th );
- //m_xt = x + dh*m_dxt
+ // m_xt = x + dh*m_dxt
             scale_sum( traits_type::begin(m_xt) ,
                        traits_type::end(m_xt) ,
                        val1, traits_type::begin(x) ,
                        dh, traits_type::begin(m_dxt) );
 
+
             // dt * m_dxm = k3
             system( m_xt , m_dxm , th );
             //m_xt = x + dt*m_dxm
@@ -110,6 +126,7 @@
                        val1, traits_type::begin(x) ,
                        dt, traits_type::begin(m_dxm) );
 
+
             // dt * m_dxh = k4
             system( m_xt , m_dxh , t + dt );
             //x += dt/6 * ( m_dxdt + m_dxt + val2*m_dxm )
@@ -132,14 +149,10 @@
                         time_type t ,
                         time_type dt )
         {
- traits_type::adjust_size( x , m_dxdt );
             system( x , m_dxdt , t );
             do_step( system , x , m_dxdt , t , dt );
         }
 
-
- /* RK4 step with error estimation to 5th order according to Cash Karp */
-
     };
 
 } // namespace odeint

Modified: sandbox/odeint/boost/numeric/odeint/stepper_rk4_classical.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepper_rk4_classical.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/stepper_rk4_classical.hpp 2010-03-09 03:09:36 EST (Tue, 09 Mar 2010)
@@ -39,8 +39,8 @@
         typedef Traits traits_type;
         typedef typename traits_type::container_type container_type;
         typedef typename traits_type::value_type value_type;
- typedef typename traits_type::iterator iterator;
- typedef typename traits_type::const_iterator const_iterator;
+// typedef typename traits_type::iterator iterator;
+// typedef typename traits_type::const_iterator const_iterator;
 
 
 
@@ -53,7 +53,6 @@
         container_type m_dxdt;
         container_type m_dxt;
         container_type m_dxm;
- container_type m_dxh;
         container_type m_xt;
 
         
@@ -63,7 +62,26 @@
         // public interface
     public:
 
- order_type order() const { return 4; }
+ order_type order_step() const { return 4; }
+
+ // standard constructor, internal containers are not initialized
+ stepper_rk4_classical( void )
+ {
+ }
+
+ // constructor, which adjusts the internal containers
+ stepper_rk4_classical( const container_type &x )
+ {
+ adjust_size( x );
+ }
+
+ void adjust_size( const container_type &x )
+ {
+ traits_type::adjust_size( x , m_dxdt );
+ traits_type::adjust_size( x , m_dxt );
+ traits_type::adjust_size( x , m_dxm );
+ traits_type::adjust_size( x , m_xt );
+ }
 
         template< class DynamicalSystem >
         void do_step( DynamicalSystem &system ,
@@ -76,41 +94,36 @@
 
             const time_type val2 = time_type( 2.0 );
 
- traits_type::adjust_size( x , m_dxt );
- traits_type::adjust_size( x , m_dxm );
- traits_type::adjust_size( x , m_xt );
 
             time_type dh = time_type( 0.5 ) * dt;
             time_type th = t + dh;
 
             //m_xt = x + dh*dxdt
- /*assign_sum( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
- traits_type::begin(x) , traits_type::begin(dxdt) , dh );*/
             scale_sum( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
                        static_cast< time_type >(1.0) ,
                        traits_type::begin(x) ,
                        dh,
                        traits_type::begin(dxdt) );
 
- system( m_xt , m_dxt , th );
+
             //m_xt = x + dh*m_dxdt
- /*assign_sum( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
- traits_type::begin(x) , traits_type::begin(m_dxt) , dh );
- */
+ system( m_xt , m_dxt , th );
             scale_sum( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
                         static_cast< time_type >(1.0),
                         traits_type::begin(x) ,
                         dh ,
                         traits_type::begin(m_dxt) );
 
- system( m_xt , m_dxm , th );
+
             //m_xt = x + dt*m_dxm ; m_dxm += m_dxt
+ system( m_xt , m_dxm , th );
             assign_sum_increment( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
                                   traits_type::begin(x) , traits_type::begin(m_dxm) ,
                                   traits_type::begin(m_dxt) , dt );
 
- system( m_xt , m_dxt , t + dt );
+
             //x = dt/6 * ( m_dxdt + m_dxt + val2*m_dxm )
+ system( m_xt , m_dxt , t + dt );
             increment_sum_sum( traits_type::begin(x) , traits_type::end(x) ,
                                traits_type::begin(dxdt) ,
                                traits_type::begin(m_dxt) ,
@@ -126,7 +139,6 @@
                         time_type t ,
                         time_type dt )
         {
- traits_type::adjust_size( x , m_dxdt );
             system( x , m_dxdt , t );
             do_step( system , x , m_dxdt , t , dt );
         }

Modified: sandbox/odeint/libs/numeric/odeint/test/check_stepper_concepts.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/test/check_stepper_concepts.cpp (original)
+++ sandbox/odeint/libs/numeric/odeint/test/check_stepper_concepts.cpp 2010-03-09 03:09:36 EST (Tue, 09 Mar 2010)
@@ -17,6 +17,9 @@
 
 #include <boost/numeric/odeint/stepper_euler.hpp>
 #include <boost/numeric/odeint/stepper_half_step.hpp>
+#include <boost/numeric/odeint/stepper_midpoint.hpp>
+#include <boost/numeric/odeint/stepper_rk4_classical.hpp>
+#include <boost/numeric/odeint/stepper_rk4.hpp>
 
 using namespace boost::unit_test;
 using namespace boost::numeric::odeint;
@@ -119,6 +122,27 @@
     check_error_stepper_concept( stepper , 1 , 2 );
 }
 
+void test_midpoint_concept()
+{
+ stepper_midpoint< std::vector< double > > stepper;
+ stepper.set_step_number( 4 );
+ unsigned short step_number = stepper.get_step_number();
+ step_number = 5; // no warnings
+ check_stepper_concept( stepper , 2 );
+}
+
+void test_rk4_classical_concept()
+{
+ stepper_rk4_classical< std::vector<double> > stepper;
+ check_stepper_concept( stepper , 4 );
+}
+
+void test_rk4_concept()
+{
+ stepper_rk4< std::vector<double> > stepper;
+ check_stepper_concept( stepper , 4 );
+}
+
 
 
 
@@ -128,6 +152,9 @@
 
     test->add( BOOST_TEST_CASE( &test_euler_concept ) );
     test->add( BOOST_TEST_CASE( &test_half_step_euler_concept ) );
+ test->add( BOOST_TEST_CASE( &test_midpoint_concept ) );
+ test->add( BOOST_TEST_CASE( &test_rk4_classical_concept ) );
+ test->add( BOOST_TEST_CASE( &test_rk4_concept ) );
 
     return test;
 }


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