[Boost-commit] svn:boost r57684 - sandbox/odeint/boost/numeric/odeint

Subject: [Boost-commit] svn:boost r57684 - sandbox/odeint/boost/numeric/odeint
From: mario.mulansky_at_[hidden]
Date: 2009-11-15 07:17:38

Author: mariomulansky
Date: 2009-11-15 07:17:38 EST (Sun, 15 Nov 2009)
New Revision: 57684
URL: http://svn.boost.org/trac/boost/changeset/57684

Log:
fixed small bug in rk4
Text files modified: 
   sandbox/odeint/boost/numeric/odeint/stepper_rk4.hpp        |    15 +++++++++------                         
   sandbox/odeint/boost/numeric/odeint/stepper_rk_generic.hpp |    12 ++++++------                            
   2 files changed, 15 insertions(+), 12 deletions(-)

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	2009-11-15 07:17:38 EST (Sun, 15 Nov 2009)
@@ -63,8 +63,8 @@
         container_type m_dxdt;
         container_type m_dxt;
         container_type m_dxm;
+        container_type m_dxh;
         container_type m_xt;
-        container_type m_x4, m_x5, m_x6;
         resizer_type m_resizer;
 
         
@@ -85,30 +85,33 @@
         {
             using namespace detail::it_algebra;
 
-            const time_type val2 = time_type( 2.0 );
+            //const time_type val2 = time_type( 2.0 );
             const time_type val1 = time_type( 1.0 );
 
             m_resizer.adjust_size( x , m_dxt );
             m_resizer.adjust_size( x , m_dxm );
             m_resizer.adjust_size( x , m_xt );
+            m_resizer.adjust_size( x , m_dxh );
 
             time_type  dh = time_type( 0.5 ) * dt;
             time_type th = t + dh;
 
+            // dt * dxdt = k1
+
             //m_xt = x + dh*dxdt
             // old: assign_sum( m_xt.begin() , m_xt.end() , x.begin() , dxdt.begin() , dh );
             scale_sum( m_xt.begin(), m_xt.end(),
                        val1, x.begin(),
                        dh, dxdt.begin() );
 
-            system( m_xt , m_dxt , th );
+            system( m_xt , m_dxt , th ); // dt * m_dxt = k2
             //m_xt = x + dh*m_dxt
             // old: assign_sum( m_xt.begin() , m_xt.end() , x.begin() , m_dxt.begin() , dh );
             scale_sum( m_xt.begin(), m_xt.end(),
                        val1, x.begin(),
                        dh, m_dxt.begin() );
 
-            system( m_xt , m_dxm , th );
+            system( m_xt , m_dxm , th ); // dt * m_dxm = k3
             //m_xt = x + dt*m_dxm ; m_dxm += m_dxt
             // old: assign_sum_increment( m_xt.begin() , m_xt.end() , x.begin() ,
             //                      m_dxm.begin() , m_dxt.begin() , dt );
@@ -116,7 +119,7 @@
                        val1, x.begin(),
                        dt, m_dxm.begin() );
 
-            system( m_xt , m_xt , value_type( t + dt ) );
+            system( m_xt , m_dxh , value_type( t + dt ) );  // dt * m_dxh = k4
             //x += dt/6 * ( m_dxdt + m_xt + val2*m_dxm )
             // old: increment_sum_sum( x.begin() , x.end() , dxdt.begin() ,
             //                   m_xt.begin() , m_dxm.begin() ,
@@ -126,7 +129,7 @@
                        dt / time_type( 6.0 ), dxdt.begin(),
                        dt / time_type( 3.0 ), m_dxt.begin(),
                        dt / time_type( 3.0 ), m_dxm.begin(),
-                       dt / time_type( 6.0 ), m_xt.begin() );
+                       dt / time_type( 6.0 ), m_dxh.begin() );
         }
 
 

Modified: sandbox/odeint/boost/numeric/odeint/stepper_rk_generic.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepper_rk_generic.hpp	(original)
+++ sandbox/odeint/boost/numeric/odeint/stepper_rk_generic.hpp	2009-11-15 07:17:38 EST (Sun, 15 Nov 2009)
@@ -61,9 +61,9 @@
         container_vector m_xvec;
         container_iterator_vector m_xiter_vec;
         container_type m_xtmp;
-        std::vector< time_type > m_a;
-        std::vector< std::vector<time_type> > m_b;
-        std::vector< time_type > m_c;
+        const std::vector< time_type > m_a;
+        const std::vector< std::vector<time_type> > m_b;
+        const std::vector< time_type > m_c;
 
         order_type m_q;
 
@@ -102,8 +102,8 @@
             x_iter = m_xvec.begin(); 
             (*x_iter++) = dxdt;
             
-            typename std::vector< time_type >::iterator a_iter = m_a.begin();
-            typename std::vector< std::vector<time_type> >::iterator b_iter = m_b.begin();
+            typename std::vector< time_type >::const_iterator a_iter = m_a.begin();
+            typename std::vector< std::vector<time_type> >::const_iterator b_iter = m_b.begin();
             while( x_iter != m_xvec.end() ) {
                 reset_iter(m_xiter_vec.begin());
                 scale_sum_generic( m_xtmp.begin(), m_xtmp.end(),
@@ -116,7 +116,7 @@
             }
 
             reset_iter(m_xiter_vec.begin());
-            typename std::vector< time_type >::iterator c_iter = m_c.begin();
+            typename std::vector< time_type >::const_iterator c_iter = m_c.begin();
             scale_sum_generic( x.begin(), x.end(),
                                m_c.begin(), m_c.end(), dt,
                                x.begin(), m_xiter_vec.begin() );