Subject: [Boost-commit] svn:boost r60376 - in sandbox/odeint: . boost/numeric/odeint libs/numeric/odeint/test
From: karsten.ahnert_at_[hidden]
Date: 2010-03-09 03:09:38
Author: karsten
Date: 2010-03-09 03:09:36 EST (Tue, 09 Mar 2010)
New Revision: 60376
URL: http://svn.boost.org/trac/boost/changeset/60376
Log:
check the stepper concepts, order->order_step, adjust_size introduced...
Text files modified:
sandbox/odeint/ToDo | 8 ++--
sandbox/odeint/boost/numeric/odeint/controlled_stepper_bs.hpp | 8 ++--
sandbox/odeint/boost/numeric/odeint/controlled_stepper_standard.hpp | 4 +-
sandbox/odeint/boost/numeric/odeint/stepper_euler.hpp | 4 +-
sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp | 4 +-
sandbox/odeint/boost/numeric/odeint/stepper_midpoint.hpp | 74 ++++++++++++++++++++++++---------------
sandbox/odeint/boost/numeric/odeint/stepper_rk4.hpp | 39 ++++++++++++++-------
sandbox/odeint/boost/numeric/odeint/stepper_rk4_classical.hpp | 44 +++++++++++++++--------
sandbox/odeint/libs/numeric/odeint/test/check_stepper_concepts.cpp | 27 ++++++++++++++
9 files changed, 141 insertions(+), 71 deletions(-)
Modified: sandbox/odeint/ToDo
==============================================================================
--- sandbox/odeint/ToDo (original)
+++ sandbox/odeint/ToDo 2010-03-09 03:09:36 EST (Tue, 09 Mar 2010)
@@ -19,13 +19,13 @@
* check: orders and adjust_size() in :
* stepper_euler.hpp - DONE
* stepper_half_step.hpp - DONE
- * stepper_midpoint.hpp
- * stepper_rk4_classical.hpp
- * stepper_rk4.hpp
+ * stepper_midpoint.hpp - DONE, changed stepcount -> step_number
+ * stepper_rk4_classical.hpp - DONE
+ * stepper_rk4.hpp - DONE
* stepper_rk5_ck.hpp
* stepper_rk78_fehlberg.hpp
* stepper_rk_generic.hpp
-
+* die auskommentierten Iteratoren typedefs entfernen
Mario:
Modified: sandbox/odeint/boost/numeric/odeint/controlled_stepper_bs.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/controlled_stepper_bs.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/controlled_stepper_bs.hpp 2010-03-09 03:09:36 EST (Tue, 09 Mar 2010)
@@ -162,12 +162,12 @@
for( unsigned short k=0; k<=m_current_k_max; k++ )
{ // loop through interval numbers
- unsigned short stepcount = m_interval_sequence[k];
+ unsigned short step_number = m_interval_sequence[k];
//out-of-place midpoint step
- m_stepper_mp.set_stepcount(stepcount);
- m_stepper_mp.do_step(system, m_x0, dxdt, t, dt, m_x_mp);
+ m_stepper_mp.set_step_number(step_number);
+ m_stepper_mp.midpoint_step(system, m_x0, dxdt, t, dt, m_x_mp);
//std::clog << "x_mp: " << k << '\t' << m_x_mp[0] << '\t' << m_x_mp[1] << std::endl;
- time_type t_est = (dt/stepcount)*(dt/stepcount);
+ time_type t_est = (dt/step_number)*(dt/step_number);
extrapolate(k, t_est, m_x_mp, x, m_xerr);
//std::clog << "Error: " << k << '\t' << m_xerr[0] << '\t' << m_xerr[1] << std::endl;
if( k != 0 )
Modified: sandbox/odeint/boost/numeric/odeint/controlled_stepper_standard.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/controlled_stepper_standard.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/controlled_stepper_standard.hpp 2010-03-09 03:09:36 EST (Tue, 09 Mar 2010)
@@ -71,8 +71,8 @@
typedef typename stepper_type::time_type time_type;
typedef typename stepper_type::traits_type traits_type;
typedef typename stepper_type::value_type value_type;
- typedef typename stepper_type::iterator iterator;
- typedef typename stepper_type::const_iterator const_iterator;
+// typedef typename stepper_type::iterator iterator;
+// typedef typename stepper_type::const_iterator const_iterator;
// private members
Modified: sandbox/odeint/boost/numeric/odeint/stepper_euler.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepper_euler.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/stepper_euler.hpp 2010-03-09 03:09:36 EST (Tue, 09 Mar 2010)
@@ -45,8 +45,8 @@
typedef Traits traits_type;
typedef typename traits_type::container_type container_type;
typedef typename traits_type::value_type value_type;
- typedef typename traits_type::iterator iterator;
- typedef typename traits_type::const_iterator const_iterator;
+// typedef typename traits_type::iterator iterator;
+// typedef typename traits_type::const_iterator const_iterator;
//
Modified: sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp 2010-03-09 03:09:36 EST (Tue, 09 Mar 2010)
@@ -39,8 +39,8 @@
typedef typename stepper_type::time_type time_type;
typedef typename stepper_type::order_type order_type;
typedef typename stepper_type::value_type value_type;
- typedef typename stepper_type::iterator iterator;
- typedef typename stepper_type::const_iterator const_iterator;
+// typedef typename stepper_type::iterator iterator;
+// typedef typename stepper_type::const_iterator const_iterator;
Modified: sandbox/odeint/boost/numeric/odeint/stepper_midpoint.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepper_midpoint.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/stepper_midpoint.hpp 2010-03-09 03:09:36 EST (Tue, 09 Mar 2010)
@@ -44,8 +44,8 @@
typedef Traits traits_type;
typedef typename traits_type::container_type container_type;
typedef typename traits_type::value_type value_type;
- typedef typename traits_type::iterator iterator;
- typedef typename traits_type::const_iterator const_iterator;
+// typedef typename traits_type::iterator iterator;
+// typedef typename traits_type::const_iterator const_iterator;
@@ -54,7 +54,7 @@
// private memebers
private:
- unsigned short m_stepcount;
+ unsigned short m_step_number;
container_type m_x0;
container_type m_x1;
@@ -65,23 +65,45 @@
public:
- stepper_midpoint( unsigned short stepcount = 2 ) { }
-
- order_type order() const { return 2; }
+ order_type order_step() const { return 2; }
+
+ // standard constructor, the size of the internal container is not set
+ stepper_midpoint( unsigned short step_number = 2 )
+ {
+ set_step_number( step_number );
+ }
+
+ // constructor, which adjusts the size of the internal containers
+ stepper_midpoint( const container_type &x , unsigned short step_number = 2 )
+ {
+ adjust_size( x );
+ set_step_number( step_number );
+ }
+
+ // adjusts the size of the internal containers
+ void adjust_size( const container_type &x )
+ {
+ traits_type::adjust_size( x , m_x0 );
+ traits_type::adjust_size( x , m_x1 );
+ traits_type::adjust_size( x , m_dxdt );
+ }
- void set_stepcount( unsigned short stepcount )
+ void set_step_number( unsigned short step_number )
{
- if( stepcount > 1 )
- m_stepcount = stepcount;
+ if( step_number > 1 )
+ m_step_number = step_number;
}
- unsigned short get_stepcount() const { return m_stepcount; }
-
+ unsigned short get_step_number() const
+ {
+ return m_step_number;
+ }
+ // performs a midpoint step
template< class DynamicalSystem >
- void do_step(
+ void midpoint_step(
DynamicalSystem &system ,
container_type &x ,
container_type &dxdt ,
@@ -89,34 +111,33 @@
time_type dt ,
container_type &x_out )
{
+ using namespace detail::it_algebra;
+
const time_type t_1 = static_cast<time_type>( 1.0 );
const time_type t_05 = static_cast<time_type>( 0.5 );
- const time_type h = dt / static_cast<time_type>( m_stepcount );
+ const time_type h = dt / static_cast<time_type>( m_step_number );
const time_type h2 = static_cast<time_type>( 2.0 ) * h;
- time_type th = t + h;
- traits_type::adjust_size(x, m_x0);
- traits_type::adjust_size(x, m_x1);
- traits_type::adjust_size(x, m_dxdt);
- using namespace detail::it_algebra;
+ time_type th = t + h;
// m_x1 = x + h*dxdt
- scale_sum( traits_type::begin(m_x1), traits_type::end(m_x1),
+ scale_sum( traits_type::begin(m_x1),
+ traits_type::end(m_x1),
t_1, traits_type::begin(x),
h, traits_type::begin(dxdt) );
system( m_x1, m_dxdt, th );
m_x0 = x;
- //std::clog<< "mp: " << 0 << '\t' << m_x1[0] << '\t' << m_x1[1] << '\t' << m_dxdt[0] << '\t' << m_dxdt[1] << std::endl;
-
unsigned short i = 1;
- while( i != m_stepcount )
- { // general step
+ while( i != m_step_number )
+ {
+ // general step
//tmp = m_x1; m_x1 = m_x0 + h2*m_dxdt; m_x0 = tmp
- scale_sum_swap( traits_type::begin(m_x1), traits_type::end(m_x1),
+ scale_sum_swap( traits_type::begin(m_x1),
+ traits_type::end(m_x1),
traits_type::begin(m_x0),
h2, traits_type::begin(m_dxdt) );
th += h;
@@ -124,8 +145,6 @@
i++;
}
- //std::clog<< "mp: " << i << '\t' << m_x1[0] << '\t' << m_x1[1] << '\t' << m_dxdt[0] << '\t' << m_dxdt[1] << std::endl;
-
// last step
// x = 0.5*( m_x0 + m_x1 + h*m_dxdt )
scale_sum( traits_type::begin(x_out), traits_type::end(x_out),
@@ -146,7 +165,7 @@
time_type t ,
time_type dt )
{
- do_step( system, x, dxdt, t, dt, x );
+ midpoint_step( system, x, dxdt, t, dt, x );
}
@@ -160,7 +179,6 @@
time_type t ,
time_type dt )
{
- traits_type::adjust_size(x, m_dxdt);
system( x, m_dxdt, t );
do_step( system , x, m_dxdt, t, dt );
}
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 2010-03-09 03:09:36 EST (Tue, 09 Mar 2010)
@@ -39,8 +39,8 @@
typedef Traits traits_type;
typedef typename traits_type::container_type container_type;
typedef typename traits_type::value_type value_type;
- typedef typename traits_type::iterator iterator;
- typedef typename traits_type::const_iterator const_iterator;
+// typedef typename traits_type::iterator iterator;
+// typedef typename traits_type::const_iterator const_iterator;
@@ -65,7 +65,26 @@
public:
- order_type order() const { return 4; }
+ order_type order_step() const { return 4; }
+
+ // standard constructor, internal containers are not initialized
+ stepper_rk4( void )
+ {
+ }
+
+ // constructor, which adjusts the internal containers
+ stepper_rk4( const container_type &x )
+ {
+ }
+
+ void adjust_size( const container_type &x )
+ {
+ traits_type::adjust_size( x , m_dxdt );
+ traits_type::adjust_size( x , m_dxt );
+ traits_type::adjust_size( x , m_dxm );
+ traits_type::adjust_size( x , m_xt );
+ traits_type::adjust_size( x , m_dxh );
+ }
@@ -80,11 +99,6 @@
const time_type val1 = static_cast<time_type>( 1.0 );
- traits_type::adjust_size( x , m_dxt );
- traits_type::adjust_size( x , m_dxm );
- traits_type::adjust_size( x , m_xt );
- traits_type::adjust_size( x , m_dxh );
-
time_type dh = static_cast<time_type>( 0.5 ) * dt;
time_type th = t + dh;
@@ -95,14 +109,16 @@
val1, traits_type::begin(x),
dh, traits_type::begin(dxdt) );
+
// dt * m_dxt = k2
system( m_xt , m_dxt , th );
- //m_xt = x + dh*m_dxt
+ // m_xt = x + dh*m_dxt
scale_sum( traits_type::begin(m_xt) ,
traits_type::end(m_xt) ,
val1, traits_type::begin(x) ,
dh, traits_type::begin(m_dxt) );
+
// dt * m_dxm = k3
system( m_xt , m_dxm , th );
//m_xt = x + dt*m_dxm
@@ -110,6 +126,7 @@
val1, traits_type::begin(x) ,
dt, traits_type::begin(m_dxm) );
+
// dt * m_dxh = k4
system( m_xt , m_dxh , t + dt );
//x += dt/6 * ( m_dxdt + m_dxt + val2*m_dxm )
@@ -132,14 +149,10 @@
time_type t ,
time_type dt )
{
- traits_type::adjust_size( x , m_dxdt );
system( x , m_dxdt , t );
do_step( system , x , m_dxdt , t , dt );
}
-
- /* RK4 step with error estimation to 5th order according to Cash Karp */
-
};
} // namespace odeint
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 2010-03-09 03:09:36 EST (Tue, 09 Mar 2010)
@@ -39,8 +39,8 @@
typedef Traits traits_type;
typedef typename traits_type::container_type container_type;
typedef typename traits_type::value_type value_type;
- typedef typename traits_type::iterator iterator;
- typedef typename traits_type::const_iterator const_iterator;
+// typedef typename traits_type::iterator iterator;
+// typedef typename traits_type::const_iterator const_iterator;
@@ -53,7 +53,6 @@
container_type m_dxdt;
container_type m_dxt;
container_type m_dxm;
- container_type m_dxh;
container_type m_xt;
@@ -63,7 +62,26 @@
// public interface
public:
- order_type order() const { return 4; }
+ order_type order_step() const { return 4; }
+
+ // standard constructor, internal containers are not initialized
+ stepper_rk4_classical( void )
+ {
+ }
+
+ // constructor, which adjusts the internal containers
+ stepper_rk4_classical( const container_type &x )
+ {
+ adjust_size( x );
+ }
+
+ void adjust_size( const container_type &x )
+ {
+ traits_type::adjust_size( x , m_dxdt );
+ traits_type::adjust_size( x , m_dxt );
+ traits_type::adjust_size( x , m_dxm );
+ traits_type::adjust_size( x , m_xt );
+ }
template< class DynamicalSystem >
void do_step( DynamicalSystem &system ,
@@ -76,41 +94,36 @@
const time_type val2 = time_type( 2.0 );
- traits_type::adjust_size( x , m_dxt );
- traits_type::adjust_size( x , m_dxm );
- traits_type::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( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
- traits_type::begin(x) , traits_type::begin(dxdt) , dh );*/
scale_sum( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
static_cast< time_type >(1.0) ,
traits_type::begin(x) ,
dh,
traits_type::begin(dxdt) );
- system( m_xt , m_dxt , th );
+
//m_xt = x + dh*m_dxdt
- /*assign_sum( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
- traits_type::begin(x) , traits_type::begin(m_dxt) , dh );
- */
+ system( m_xt , m_dxt , th );
scale_sum( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
static_cast< time_type >(1.0),
traits_type::begin(x) ,
dh ,
traits_type::begin(m_dxt) );
- system( m_xt , m_dxm , th );
+
//m_xt = x + dt*m_dxm ; m_dxm += m_dxt
+ system( m_xt , m_dxm , th );
assign_sum_increment( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
traits_type::begin(x) , traits_type::begin(m_dxm) ,
traits_type::begin(m_dxt) , dt );
- system( m_xt , m_dxt , t + dt );
+
//x = dt/6 * ( m_dxdt + m_dxt + val2*m_dxm )
+ system( m_xt , m_dxt , t + dt );
increment_sum_sum( traits_type::begin(x) , traits_type::end(x) ,
traits_type::begin(dxdt) ,
traits_type::begin(m_dxt) ,
@@ -126,7 +139,6 @@
time_type t ,
time_type dt )
{
- traits_type::adjust_size( x , m_dxdt );
system( x , m_dxdt , t );
do_step( system , x , m_dxdt , t , dt );
}
Modified: sandbox/odeint/libs/numeric/odeint/test/check_stepper_concepts.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/test/check_stepper_concepts.cpp (original)
+++ sandbox/odeint/libs/numeric/odeint/test/check_stepper_concepts.cpp 2010-03-09 03:09:36 EST (Tue, 09 Mar 2010)
@@ -17,6 +17,9 @@
#include <boost/numeric/odeint/stepper_euler.hpp>
#include <boost/numeric/odeint/stepper_half_step.hpp>
+#include <boost/numeric/odeint/stepper_midpoint.hpp>
+#include <boost/numeric/odeint/stepper_rk4_classical.hpp>
+#include <boost/numeric/odeint/stepper_rk4.hpp>
using namespace boost::unit_test;
using namespace boost::numeric::odeint;
@@ -119,6 +122,27 @@
check_error_stepper_concept( stepper , 1 , 2 );
}
+void test_midpoint_concept()
+{
+ stepper_midpoint< std::vector< double > > stepper;
+ stepper.set_step_number( 4 );
+ unsigned short step_number = stepper.get_step_number();
+ step_number = 5; // no warnings
+ check_stepper_concept( stepper , 2 );
+}
+
+void test_rk4_classical_concept()
+{
+ stepper_rk4_classical< std::vector<double> > stepper;
+ check_stepper_concept( stepper , 4 );
+}
+
+void test_rk4_concept()
+{
+ stepper_rk4< std::vector<double> > stepper;
+ check_stepper_concept( stepper , 4 );
+}
+
@@ -128,6 +152,9 @@
test->add( BOOST_TEST_CASE( &test_euler_concept ) );
test->add( BOOST_TEST_CASE( &test_half_step_euler_concept ) );
+ test->add( BOOST_TEST_CASE( &test_midpoint_concept ) );
+ test->add( BOOST_TEST_CASE( &test_rk4_classical_concept ) );
+ test->add( BOOST_TEST_CASE( &test_rk4_concept ) );
return test;
}