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;
}