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

Subject: [Boost-commit] svn:boost r57683 - in sandbox/odeint: boost/numeric/odeint libs/numeric/odeint/examples
From: mario.mulansky_at_[hidden]
Date: 2009-11-15 06:53:26

Author: mariomulansky
Date: 2009-11-15 06:53:25 EST (Sun, 15 Nov 2009)
New Revision: 57683
URL: http://svn.boost.org/trac/boost/changeset/57683

Log:
added stepper_rk_generic (forgotten yesterday)
changed stepper_rk4 to use only iterator function scale_sum
Added:
   sandbox/odeint/boost/numeric/odeint/stepper_rk_generic.hpp   (contents, props changed)
Text files modified: 
   sandbox/odeint/boost/numeric/odeint/stepper_rk4.hpp                       |    40 ++++++++++++++++++++++++++++------------
   sandbox/odeint/libs/numeric/odeint/examples/lorenz_cmp_rk4_rk_generic.cpp |    23 +++++++++++++----------                 
   2 files changed, 41 insertions(+), 22 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 06:53:25 EST (Sun, 15 Nov 2009)
@@ -79,13 +79,14 @@
         template< class DynamicalSystem >
         void next_step( DynamicalSystem &system ,
                         container_type &x ,
-                        const container_type &dxdt ,
+                        container_type &dxdt ,
                         time_type t ,
                         time_type dt )
         {
             using namespace detail::it_algebra;
 
             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 );
@@ -95,22 +96,37 @@
             time_type th = t + dh;
 
             //m_xt = x + dh*dxdt
-            assign_sum( m_xt.begin() , m_xt.end() , x.begin() , dxdt.begin() , dh );
+            // 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 );
-            //m_xt = x + dh*m_dxdt
-            assign_sum( m_xt.begin() , m_xt.end() , x.begin() , m_dxt.begin() , dh );
+            //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 );
             //m_xt = x + dt*m_dxm ; m_dxm += m_dxt
-            assign_sum_increment( m_xt.begin() , m_xt.end() , x.begin() ,
-                                  m_dxm.begin() , m_dxt.begin() , dt );
-
-            system( m_xt , m_dxt , value_type( t + dt ) );
-            //x = dt/6 * ( m_dxdt + m_dxt + val2*m_dxm )
-            increment_sum_sum( x.begin() , x.end() , dxdt.begin() ,
-                               m_dxt.begin() , m_dxm.begin() ,
-                               dt /  time_type( 6.0 ) , val2 );
+            // old: assign_sum_increment( m_xt.begin() , m_xt.end() , x.begin() ,
+            //                      m_dxm.begin() , m_dxt.begin() , dt );
+            scale_sum( m_xt.begin(), m_xt.end(),
+                       val1, x.begin(),
+                       dt, m_dxm.begin() );
+
+            system( m_xt , m_xt , value_type( t + dt ) );
+            //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() ,
+            //                   dt /  time_type( 6.0 ) , val2 );
+            scale_sum( x.begin(), x.end(),
+                       val1, x.begin(),
+                       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() );
         }
 
 

Added: sandbox/odeint/boost/numeric/odeint/stepper_rk_generic.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/boost/numeric/odeint/stepper_rk_generic.hpp	2009-11-15 06:53:25 EST (Sun, 15 Nov 2009)
@@ -0,0 +1,150 @@
+/* Boost odeint/stepper_rk_generic.hpp header file
+ 
+ Copyright 2009 Karsten Ahnert
+ Copyright 2009 Mario Mulansky
+ Copyright 2009 Andre Bergner
+ 
+ This file includes generic Runge Kutta solver
+ for ordinary differential equations.
+
+ It solves any ODE dx/dt = f(x,t).
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+*/
+
+#ifndef BOOST_NUMERIC_ODEINT_STEPPER_RK_GENERIC_HPP
+#define BOOST_NUMERIC_ODEINT_STEPPER_RK_GENERIC_HPP
+
+#include <vector>
+#include <boost/concept_check.hpp>
+
+#include <boost/numeric/odeint/concepts/state_concept.hpp>
+#include <boost/numeric/odeint/resizer.hpp>
+
+#include <iostream>
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+    template<
+        class Container ,
+        class Time = double ,
+        class Resizer = resizer< Container >
+        >
+    class stepper_rk_generic {
+
+        // provide basic typedefs
+    public:
+
+        typedef Container container_type;
+        typedef Resizer resizer_type;
+        typedef Time time_type;
+        typedef const unsigned short order_type;
+        typedef typename container_type::value_type value_type;
+        typedef typename container_type::iterator iterator;
+
+
+        // check the concept of the ContainerType
+    private:
+
+        BOOST_CLASS_REQUIRE( container_type ,
+                             boost::numeric::odeint, Container );
+
+        // private variables
+    private:
+        typedef std::vector< container_type > container_vector;
+        typedef std::vector< iterator > container_iterator_vector;
+
+        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;
+
+        order_type m_q;
+
+        resizer_type m_resizer;
+
+    public:
+
+        stepper_rk_generic( std::vector<time_type> &a,
+                            std::vector< std::vector<time_type> > &b,
+                            std::vector< time_type > &c)
+            : m_a(a), m_b(b), m_c(c), m_q(c.size())
+        {
+            m_xvec.resize(m_q);
+            m_xiter_vec.resize(m_q);
+        }
+
+        order_type order() const { return m_q; }
+
+        template< class DynamicalSystem >
+        void next_step( DynamicalSystem &system ,
+                        container_type &x ,
+                        container_type &dxdt ,
+                        time_type t ,
+                        time_type dt )
+        {
+            using namespace detail::it_algebra;
+            typename container_vector::iterator x_iter = m_xvec.begin();
+            typename container_iterator_vector::iterator xiter_iter = m_xiter_vec.begin();
+            while( x_iter != m_xvec.end() ) {
+                m_resizer.adjust_size(x, (*x_iter));
+                (*xiter_iter++) = (*x_iter).begin();
+                x_iter++;
+            }
+            m_resizer.adjust_size(x, m_xtmp);
+            
+            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();
+            while( x_iter != m_xvec.end() ) {
+                reset_iter(m_xiter_vec.begin());
+                scale_sum_generic( m_xtmp.begin(), m_xtmp.end(),
+                                   (*b_iter).begin(), (*b_iter).end(), dt,
+                                   x.begin(), m_xiter_vec.begin() );
+                system( m_xtmp, *x_iter , t + dt*(*a_iter) );
+                x_iter++;
+                a_iter++;
+                b_iter++;
+            }
+
+            reset_iter(m_xiter_vec.begin());
+            typename std::vector< time_type >::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() );
+        }
+
+        void reset_iter(typename container_iterator_vector::iterator xiter_iter)
+        {
+            typename container_vector::iterator x_iter = m_xvec.begin();
+            while( x_iter != m_xvec.end() ) {
+                (*xiter_iter++) = (*x_iter++).begin();
+            }
+        }
+
+        template< class DynamicalSystem >
+        void next_step( DynamicalSystem &system ,
+                        container_type &x ,
+                        time_type t ,
+                        time_type dt )
+        {
+            m_resizer.adjust_size(x, m_xtmp);
+            system(x, m_xtmp, t);
+            next_step( system, x, m_xtmp, t, dt);
+        }
+
+    };
+
+}
+}
+}
+
+#endif

Modified: sandbox/odeint/libs/numeric/odeint/examples/lorenz_cmp_rk4_rk_generic.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/lorenz_cmp_rk4_rk_generic.cpp	(original)
+++ sandbox/odeint/libs/numeric/odeint/examples/lorenz_cmp_rk4_rk_generic.cpp	2009-11-15 06:53:25 EST (Sun, 15 Nov 2009)
@@ -80,29 +80,32 @@
 
     cout.precision(16);
 
+    cout << "Hand written Runge Kutta 4"<<endl;
+
     t = 0.0;
-    stepper_rk4.next_step( lorenz , x1 , t , dt );
-    cout << "x after one step: "<<x1[0]<<tab<<x1[1]<<tab<<x1[2]<<endl;
-    t += dt;
     start= clock();
-    for( size_t oi=0 ; oi<olen ; ++oi,t+=dt ) {
+    for( size_t oi=1 ; oi<olen ; ++oi,t+=dt ) {
         stepper_rk4.next_step( lorenz , x1 , t , dt );
+        if( oi < 5 )
+            cout << "x after step "<<oi<<": "<<x1[0]<<tab<<x1[1]<<tab<<x1[2]<<endl;
+
     }
     end = clock();
-    cout << "RK4 : " << double ( end - start ) / double( CLOCKS_PER_SEC ) << endl;
     cout << "x after "<<olen<<" steps: "<<x1[0]<<tab<<x1[1]<<tab<<x1[2]<<endl;
+    cout << "Time for "<<olen<<" steps: " << double ( end - start ) / double( CLOCKS_PER_SEC ) <<"s"<< endl;
+    cout << endl << "###########################" << endl << endl;
+    cout << "Runge Kutta 4 via generic Runge Kutta implementation"<<endl;
 
     t = 0.0;
-    stepper_generic4.next_step( lorenz , x2 , t , dt );
-    cout << "x after one step: "<<x2[0]<<tab<<x2[1]<<tab<<x2[2]<<endl;
-    t += dt;
     start= clock();
-    for( size_t oi=0 ; oi<olen ; ++oi,t+=dt ) {
+    for( size_t oi=1 ; oi<olen ; ++oi,t+=dt ) {
         stepper_generic4.next_step( lorenz , x2 , t , dt );
+        if( oi < 5 )
+            cout << "x after step "<<oi<<":  "<<x2[0]<<tab<<x2[1]<<tab<<x2[2]<<endl;        
     }
     end = clock();
-    cout << "Generic RK4 : " << double ( end - start ) / double( CLOCKS_PER_SEC ) << endl;
     cout << "x after "<<olen<<" steps: "<<x2[0]<<tab<<x2[1]<<tab<<x2[2]<<endl;
+    cout << "Time for "<<olen<<" steps: " << double ( end - start ) / double( CLOCKS_PER_SEC ) << endl;
 
     return 0;
 }