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

Subject: [Boost-commit] svn:boost r57709 - in sandbox/odeint: boost/numeric boost/numeric/odeint libs/numeric/odeint/examples
From: karsten.ahnert_at_[hidden]
Date: 2009-11-16 14:35:03

Author: karsten
Date: 2009-11-16 14:35:01 EST (Mon, 16 Nov 2009)
New Revision: 57709
URL: http://svn.boost.org/trac/boost/changeset/57709

Log:
added the classical rk4 solver
Added:
   sandbox/odeint/boost/numeric/odeint/stepper_rk4_classical.hpp   (contents, props changed)
Text files modified: 
   sandbox/odeint/boost/numeric/odeint.hpp                                   |     1 +                                       
   sandbox/odeint/boost/numeric/odeint/observer.hpp                          |    13 ++++++++-----                           
   sandbox/odeint/libs/numeric/odeint/examples/lorenz_cmp_rk4_rk_generic.cpp |    30 ++++++++++++++++++++++++------          
   3 files changed, 33 insertions(+), 11 deletions(-)

Modified: sandbox/odeint/boost/numeric/odeint.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint.hpp	(original)
+++ sandbox/odeint/boost/numeric/odeint.hpp	2009-11-16 14:35:01 EST (Mon, 16 Nov 2009)
@@ -17,6 +17,7 @@
 
 #include <boost/numeric/odeint/stepper_euler.hpp>
 #include <boost/numeric/odeint/stepper_rk4.hpp>
+#include <boost/numeric/odeint/stepper_rk4_classical.hpp>
 #include <boost/numeric/odeint/stepper_rk5_ck.hpp>
 #include <boost/numeric/odeint/stepper_rk_generic.hpp>
 #include <boost/numeric/odeint/stepper_half_step.hpp>

Modified: sandbox/odeint/boost/numeric/odeint/observer.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/observer.hpp	(original)
+++ sandbox/odeint/boost/numeric/odeint/observer.hpp	2009-11-16 14:35:01 EST (Mon, 16 Nov 2009)
@@ -26,8 +26,10 @@
     }
     
 
+
     template< class InsertIterator, class TimeContainer = std::vector<double> >
-    class state_copy_observer {
+    class state_copy_observer
+    {
         
     private:
         TimeContainer *m_times;
@@ -37,8 +39,11 @@
         typedef typename TimeContainer::value_type time_type;
 
     public:
-        state_copy_observer( TimeContainer &times, InsertIterator state_inserter ) 
-            : m_times(&times), m_state_inserter(state_inserter), m_time_iter(m_times->begin()) 
+        state_copy_observer( TimeContainer &times ,
+			     InsertIterator state_inserter ) 
+            : m_times(&times),
+	      m_state_inserter(state_inserter),
+	      m_time_iter(m_times->begin()) 
         {  }
         
         template< class Container, class System >
@@ -48,9 +53,7 @@
                 m_time_iter++; // next time point
             }
         }
-            
     };
-   
 
 } // odeint
 } // numeric

Added: sandbox/odeint/boost/numeric/odeint/stepper_rk4_classical.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/boost/numeric/odeint/stepper_rk4_classical.hpp	2009-11-16 14:35:01 EST (Mon, 16 Nov 2009)
@@ -0,0 +1,139 @@
+/* Boost odeint/stepper_rk4.hpp header file
+ 
+ Copyright 2009 Karsten Ahnert
+ Copyright 2009 Mario Mulansky
+ Copyright 2009 Andre Bergner
+ 
+ This file includes the explicit 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_RK4_CLASSICAL_HPP
+#define BOOST_NUMERIC_ODEINT_STEPPER_RK4_CLASSICAL_HPP
+
+#include <boost/concept_check.hpp>
+
+#include <boost/numeric/odeint/concepts/state_concept.hpp>
+#include <boost/numeric/odeint/resizer.hpp>
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+    template<
+        class Container ,
+        class Time = double ,
+        class Resizer = resizer< Container >
+        >
+    class stepper_rk4_classical
+    {
+
+        // 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 members
+    private:
+
+        container_type m_dxdt;
+        container_type m_dxt;
+        container_type m_dxm;
+        container_type m_dxh;
+        container_type m_xt;
+        resizer_type m_resizer;
+
+        
+
+
+
+        // public interface
+    public:
+
+        order_type order() const { return 4; }
+
+        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;
+
+            const time_type val2 = time_type( 2.0 );
+
+            m_resizer.adjust_size( x , m_dxt );
+            m_resizer.adjust_size( x , m_dxm );
+            m_resizer.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( m_xt.begin() , m_xt.end() , x.begin() , dxdt.begin() , dh );
+
+            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 );
+
+            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 , 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 );
+        }
+
+
+
+        template< class DynamicalSystem >
+        void next_step( DynamicalSystem &system ,
+                        container_type &x ,
+                        time_type t ,
+                        time_type dt )
+        {
+            m_resizer.adjust_size( x , m_dxdt );
+            system( x , m_dxdt , t );
+            next_step( system , x , m_dxdt , t , dt );
+        }
+
+
+        /* RK4 step with error estimation to 5th order according to Cash Karp */
+
+    };
+
+} // namespace odeint
+} // namespace numeric
+} // namespace boost
+
+
+#endif // BOOST_NUMERIC_ODEINT_STEPPER_RK4_CLASSICAL_HPP

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-16 14:35:01 EST (Mon, 16 Nov 2009)
@@ -21,7 +21,7 @@
 #include <tr1/array>
 
 #include <boost/numeric/odeint.hpp>
-#include <boost/numeric/odeint/stepper_rk_generic_test.hpp>
+ // #include <boost/numeric/odeint/stepper_rk_generic.hpp>
 
 #define tab "\t"
 
@@ -50,8 +50,10 @@
     
     state_type x1 = {{1.0, 0.0, 0.0}};
     state_type x2 = {{1.0, 0.0, 0.0}};
+    state_type x3 = {{1.0, 0.0, 0.0}};
 
     stepper_rk4< state_type > stepper_rk4;
+    stepper_rk4_classical< state_type > stepper_rk4_classical;
 
     vector< double > a(3);
 
@@ -74,13 +76,13 @@
     c[2] = 1.0/3.0;
     c[3] = 1.0/6.0;
 
-    const double a2[3] = {0.5, 0.5, 1.0 };
+/*    const double a2[3] = {0.5, 0.5, 1.0 };
     const double b2[6] = {0.5, 0.0, 0.5, 0.0, 0.0, 1.0};
-    const double c2[4] = {1.0/6.0, 1.0/3.0, 1.0/3.0, 1.0/6.0};
+    const double c2[4] = {1.0/6.0, 1.0/3.0, 1.0/3.0, 1.0/6.0};*/
 
     stepper_rk_generic< state_type > stepper_generic4(a, b, c);
 
-    stepper_rk_generic_test< state_type > stepper_generic4_test(a2, b2, c2, 4);
+//    stepper_rk_generic_test< state_type > stepper_generic4_test(a2, b2, c2, 4);
 
     clock_t start , end;
     double t;
@@ -100,9 +102,10 @@
     end = clock();
     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;
 
+    cout << "Runge Kutta 4 via generic Runge Kutta implementation"<<endl;
     t = 0.0;
     start= clock();
     for( size_t oi=1 ; oi<olen ; ++oi,t+=dt ) {
@@ -114,6 +117,21 @@
     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;
 
+
+    cout << endl << "###########################" << endl << endl;
+
+    cout << "Runge Kutta 4 with classical NR-style Runge Kutta implementation"<<endl;
+    t = 0.0;
+    start= clock();
+    for( size_t oi=1 ; oi<olen ; ++oi,t+=dt ) {
+        stepper_rk4_classical.next_step( lorenz , x3 , t , dt );
+        if( oi < 5 )
+            cout << "x after step "<<oi<<":  "<<x3[0]<<tab<<x3[1]<<tab<<x3[2]<<endl;        
+    }
+    end = clock();
+    cout << "x after "<<olen<<" steps: "<<x3[0]<<tab<<x3[1]<<tab<<x3[2]<<endl;
+    cout << "Time for "<<olen<<" steps: " << double ( end - start ) / double( CLOCKS_PER_SEC ) << endl;
+
     return 0;
 }
 
@@ -121,5 +139,5 @@
 
 /*
   Compile with
-  g++ -Wall -O3 -I$BOOST_ROOT -I../../../../ lorenz_stepper.cpp
+  g++ -Wall -O3 -I$BOOST_ROOT -I../../../../ lorenz_cmp_rk4_rk_generic.cpp
 */