Boost logo

Boost-Commit :

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
 */


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