[Boost-commit] svn:boost r57760 - in sandbox/odeint: boost/numeric/odeint boost/numeric/odeint/detail libs/numeric/odeint/stuff/test

Subject: [Boost-commit] svn:boost r57760 - in sandbox/odeint: boost/numeric/odeint boost/numeric/odeint/detail libs/numeric/odeint/stuff/test
From: mario.mulansky_at_[hidden]
Date: 2009-11-18 16:53:02

Author: mariomulansky
Date: 2009-11-18 16:53:01 EST (Wed, 18 Nov 2009)
New Revision: 57760
URL: http://svn.boost.org/trac/boost/changeset/57760

Log:
midpoint method complete
Text files modified: 
   sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp    |    23 +++++++++++++++++++++++                 
   sandbox/odeint/boost/numeric/odeint/stepper_midpoint.hpp           |    24 +++++++++++++-----------                
   sandbox/odeint/libs/numeric/odeint/stuff/test/stepper_midpoint.cpp |    16 ++++++++++------                        
   3 files changed, 46 insertions(+), 17 deletions(-)

Modified: sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp	(original)
+++ sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp	2009-11-18 16:53:01 EST (Wed, 18 Nov 2009)
@@ -260,6 +260,29 @@
     }
 
 
+    /* calculates tmp = y, y = x1 + alpha*x2, x1 = tmp */
+    template <
+        class OutputIterator,
+        class InputIterator,
+        class T
+        >
+    inline void scale_sum_swap(
+            OutputIterator y_begin,
+            OutputIterator y_end,
+            OutputIterator x1_begin,
+            T alpha,
+            InputIterator x2_begin )
+    {
+        T swap;
+        while( y_begin != y_end )
+        {
+            swap = (*x1_begin) + alpha*(*x2_begin++);
+            *x1_begin++ = *y_begin;
+            *y_begin++ = swap;
+        }
+    }
+
+
     // computes y = x1 + alpha2 * x2 ; x2 += x3
     template<
         class OutputIterator ,

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	2009-11-18 16:53:01 EST (Wed, 18 Nov 2009)
@@ -62,7 +62,7 @@
 
         container_type m_x0;
         container_type m_x1;
-        container_type m_x2;
+        //container_type m_x2;
         container_type m_dxdt;
 
     public:
@@ -88,34 +88,36 @@
 
             m_resizer.adjust_size(x, m_x0);
             m_resizer.adjust_size(x, m_x1);
-            m_resizer.adjust_size(x, m_x2);
+            m_resizer.adjust_size(x, m_dxdt);
 
             using namespace detail::it_algebra;
 
-            scale_sum( m_x1.begin(), m_x1.end(),
+            // m_x0 = x + h*dxdt
+            scale_sum( m_x0.begin(), m_x0.end(),
                        t_1, x.begin(),
                        h, dxdt.begin() );
-            system( m_x1, m_x2, th );
+            system( m_x0, m_dxdt, th );
 
             m_x1 = x;
 
-            unsigned short i = 2;
+            unsigned short i = 1;
             while( i != n )
             {   // general step
-                m_x0 = m_x1;
-                scale_sum( m_x1.begin(), m_x1.end(),
-                           t_1, m_x0.begin(),
-                           h2, m_x2.begin() );
+                //tmp = m_x1; m_x1 = m_x0 + h2*m_dxdt; m_x0 = tmp
+                scale_sum_swap( m_x1.begin(), m_x1.end(), 
+                                m_x0.begin(),
+                                h2, m_dxdt.begin() );
                 th += h;
-                system( m_x1, m_x2, th);
+                system( m_x1, m_dxdt, th);
                 i++;
             }
             
             // last step
+            // x = 0.5*( m_x0 + m_x1 + h*m_dxdt
             scale_sum( x.begin(), x.end(),
                        t_05, m_x0.begin(),
                        t_05, m_x1.begin(),
-                       t_05*h, m_x2.begin() );
+                       t_05*h, m_dxdt.begin() );
         }
 
 

Modified: sandbox/odeint/libs/numeric/odeint/stuff/test/stepper_midpoint.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/stuff/test/stepper_midpoint.cpp	(original)
+++ sandbox/odeint/libs/numeric/odeint/stuff/test/stepper_midpoint.cpp	2009-11-18 16:53:01 EST (Wed, 18 Nov 2009)
@@ -32,25 +32,29 @@
     state_type x2 = {{ 1.0, 0.0 }};
     state_type x3 = {{ 1.0, 0.0 }};
     state_type x4 = {{ 1.0, 0.0 }};
+    state_type x5 = {{ 1.0, 0.0 }};
 
+    stepper_euler< state_type > stepper_euler;
     stepper_midpoint< state_type > stepper_mp;
     stepper_rk4< state_type > stepper_rk4;
 
     double t = 0.0;
     for( size_t oi=0 ; oi<olen ; ++oi,t+=dt )
     {
-        stepper_mp.next_step( harmonic_oscillator , x1 , t , dt , 2 );
-        stepper_mp.next_step( harmonic_oscillator , x2 , t , dt , 4 );
-        stepper_mp.next_step( harmonic_oscillator , x3 , t , dt , 8 );
-        stepper_rk4.next_step( harmonic_oscillator , x4 , t , dt );
+        stepper_euler.next_step( harmonic_oscillator , x1 , t , dt );
+        stepper_mp.next_step( harmonic_oscillator , x2 , t , dt , 2 );
+        stepper_mp.next_step( harmonic_oscillator , x3 , t , dt , 4 );
+        stepper_mp.next_step( harmonic_oscillator , x4 , t , dt , 8 );
+        stepper_rk4.next_step( harmonic_oscillator , x5 , t , dt );
         cout<< t << tab << x1[0]*x1[0] + x1[1]*x1[1];
         cout<< tab << x2[0]*x2[0] + x2[1]*x2[1];
         cout<< tab << x3[0]*x3[0] + x3[1]*x3[1];
-        cout<< tab << x4[0]*x4[0] + x4[1]*x4[1] << endl;
+        cout<< tab << x4[0]*x4[0] + x4[1]*x4[1];
+        cout<< tab << x5[0]*x5[0] + x5[1]*x5[1] << endl;
     }
 }
 
 /*
   Compile with
-  g++ -Wall -O3 -I$BOOST_ROOT -I../../../../../ lorenz_stepper_cmp.cpp
+  g++ -Wall -O3 -I$BOOST_ROOT -I../../../../../ stepper_midpoint.cpp
 */