Subject: [Boost-commit] svn:boost r57716 - in sandbox/odeint: boost/numeric/odeint libs/numeric/odeint/examples
From: mario.mulansky_at_[hidden]
Date: 2009-11-16 17:18:45
Author: mariomulansky
Date: 2009-11-16 17:18:44 EST (Mon, 16 Nov 2009)
New Revision: 57716
URL: http://svn.boost.org/trac/boost/changeset/57716
Log:
+ generic RK with c-styled parameter handling
unfortunately not faster...
Text files modified:
sandbox/odeint/boost/numeric/odeint/stepper_rk4_classical.hpp | 2
sandbox/odeint/boost/numeric/odeint/stepper_rk_generic.hpp | 186 ++++++++++++++++++++++++++++++++++++++++
sandbox/odeint/libs/numeric/odeint/examples/lorenz_cmp_rk4_rk_generic.cpp | 21 +++
3 files changed, 203 insertions(+), 6 deletions(-)
Modified: sandbox/odeint/boost/numeric/odeint/stepper_rk4_classical.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepper_rk4_classical.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/stepper_rk4_classical.hpp 2009-11-16 17:18:44 EST (Mon, 16 Nov 2009)
@@ -127,8 +127,6 @@
}
- /* RK4 step with error estimation to 5th order according to Cash Karp */
-
};
} // namespace odeint
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-16 17:18:44 EST (Mon, 16 Nov 2009)
@@ -238,6 +238,192 @@
};
+
+
+ /* ############################################################
+ ############################################################
+
+ C-Array Version of a,b,c handling
+
+ ############################################################
+ ############################################################
+ */
+
+
+
+
+ template<
+ class Container ,
+ class Time = double ,
+ class Resizer = resizer< Container >
+ >
+ class stepper_rk_generic_ptr {
+
+ // 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;
+ const time_type* m_a;
+ const time_type* m_b;
+ const time_type* m_c;
+
+ order_type m_q;
+
+ resizer_type m_resizer;
+
+ private:
+ 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();
+ }
+ }
+
+ void check_consitency()
+ {
+ const time_type* a_iter = &m_a[0];
+ const time_type* b_iter = &m_b[0];
+ const time_type* c_iter = &m_c[0];
+
+ unsigned short b_len = 1;
+ // check 1: a_i = sum_j b_ij
+ while( a_iter != &m_a[+m_q-1] ) {
+ time_type tmp = 0.0;
+ const time_type* b_end = b_iter + b_len;
+ while( b_iter != b_end ) {
+ tmp += *b_iter++;
+ }
+ b_len++;
+ if( *a_iter++ != tmp )
+ throw butcher_tableau_consistency_exception();
+ }
+
+ // check 2: sum_i c_i * (a_i)^(k-1) = 1/k k = 1..q
+ for( unsigned short k = 1; k <= m_q; k++ ) {
+ time_type tmp = 0.0;
+ a_iter = &m_a[0];
+ c_iter = &m_c[0];
+ if( k == 1 ) // special treatment for 0^0 = 1
+ tmp += *c_iter++;
+ else
+ c_iter++;
+ while( a_iter != &m_a[m_q-1] ) {
+ tmp += (*c_iter++) * pow( *a_iter++ , k-1 );
+ }
+ if( std::abs(time_type(tmp - time_type(1.0)/time_type(k))) >
+ std::numeric_limits<time_type>::epsilon() ) {
+ //std::clog << std::abs(time_type(tmp - time_type(1.0)/time_type(k))) << std::endl;
+ throw butcher_tableau_order_condition_exception();
+ }
+ }
+ }
+
+
+ public:
+
+ /* Constructor
+
+ Same as for the stepper_rk_generic class, but with a, b and c provided
+ as time_type*. The order q of the integration scheme is given separately
+ as parameter. a is a pointer to an array of time_type[] of length q-1,
+ b is of length 1 + 2 + ... q-1 = (q-1)(q-2)/2 and ordered as follows:
+ b[0] = b_21, b[1] = b_31, b[2] = b_32, b[3] = b_41, b[4] = b42 ...
+ c has length q.
+
+ */
+ stepper_rk_generic_ptr( const time_type* a,
+ const time_type* b,
+ const time_type* c,
+ const unsigned short q)
+ : m_a(a), m_b(b), m_c(c), m_q(q)
+ {
+ m_xvec.resize(m_q);
+ m_xiter_vec.resize(m_q);
+
+ check_consitency();
+ }
+
+ 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();
+
+ (*x_iter) = dxdt;
+ (*xiter_iter++) = (*x_iter++).begin();
+
+ while( x_iter != m_xvec.end() ) {
+ m_resizer.adjust_size(x, (*x_iter));
+ (*xiter_iter++) = (*x_iter++).begin();
+ }
+ m_resizer.adjust_size(x, m_xtmp);
+
+ x_iter = m_xvec.begin()+1;
+
+ const time_type* a_iter = &m_a[0];
+ const time_type* b_iter = &m_b[0];
+ unsigned short b_len= 1;
+ while( x_iter != m_xvec.end() ) {
+ reset_iter(m_xiter_vec.begin());
+ const time_type* b_end = b_iter + b_len;
+ scale_sum_generic( m_xtmp.begin(), m_xtmp.end(),
+ b_iter, b_end, dt,
+ x.begin(), m_xiter_vec.begin() );
+ system( m_xtmp, *x_iter++ , t + dt*(*a_iter++) );
+ b_iter = b_end;
+ b_len++;
+ }
+
+ reset_iter(m_xiter_vec.begin());
+ scale_sum_generic( x.begin(), x.end(),
+ &m_c[0], &m_c[m_q], dt,
+ x.begin(), m_xiter_vec.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);
+ }
+
+ };
+
}
}
}
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 17:18:44 EST (Mon, 16 Nov 2009)
@@ -21,7 +21,6 @@
#include <tr1/array>
#include <boost/numeric/odeint.hpp>
- // #include <boost/numeric/odeint/stepper_rk_generic.hpp>
#define tab "\t"
@@ -76,13 +75,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_ptr< state_type > stepper_generic4_ptr(a2, b2, c2, 4);
clock_t start , end;
double t;
@@ -117,6 +116,20 @@
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 via C-Array styled generic Runge Kutta implementation"<<endl;
+
+ t = 0.0;
+ start= clock();
+ for( size_t oi=1 ; oi<olen ; ++oi,t+=dt ) {
+ stepper_generic4_ptr.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;
+
cout << endl << "###########################" << endl << endl;