|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r57685 - in sandbox/odeint: boost/numeric/odeint libs/numeric/odeint/examples
From: mario.mulansky_at_[hidden]
Date: 2009-11-15 12:21:59
Author: mariomulansky
Date: 2009-11-15 12:21:58 EST (Sun, 15 Nov 2009)
New Revision: 57685
URL: http://svn.boost.org/trac/boost/changeset/57685
Log:
added consistency check to generic rk stepper
Text files modified:
sandbox/odeint/boost/numeric/odeint/stepper_rk4.hpp | 4
sandbox/odeint/boost/numeric/odeint/stepper_rk_generic.hpp | 127 ++++++++++++++++++++++++++++++++++-----
sandbox/odeint/boost/numeric/odeint/stepsize_controller_standard.hpp | 24 +++---
sandbox/odeint/libs/numeric/odeint/examples/lorenz_cmp_rk4_rk_generic.cpp | 2
sandbox/odeint/libs/numeric/odeint/examples/lorenz_stepper.cpp | 9 +-
5 files changed, 131 insertions(+), 35 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 12:21:58 EST (Sun, 15 Nov 2009)
@@ -120,9 +120,9 @@
dt, m_dxm.begin() );
system( m_xt , m_dxh , value_type( t + dt ) ); // dt * m_dxh = k4
- //x += dt/6 * ( m_dxdt + m_xt + val2*m_dxm )
+ //x += dt/6 * ( m_dxdt + m_dxt + val2*m_dxm )
// old: increment_sum_sum( x.begin() , x.end() , dxdt.begin() ,
- // m_xt.begin() , m_dxm.begin() ,
+ // m_dxt.begin() , m_dxm.begin() ,
// dt / time_type( 6.0 ) , val2 );
scale_sum( x.begin(), x.end(),
val1, x.begin(),
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-15 12:21:58 EST (Sun, 15 Nov 2009)
@@ -7,7 +7,7 @@
This file includes generic Runge Kutta solver
for ordinary differential equations.
- It solves any ODE dx/dt = f(x,t).
+ It solves any ODE dx/dt = f(x,t) using a general Runge Kutta scheme.
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or
@@ -18,6 +18,9 @@
#define BOOST_NUMERIC_ODEINT_STEPPER_RK_GENERIC_HPP
#include <vector>
+#include <exception>
+#include <cmath> // for pow( ) and abs()
+#include <limits>
#include <boost/concept_check.hpp>
#include <boost/numeric/odeint/concepts/state_concept.hpp>
@@ -29,6 +32,24 @@
namespace numeric {
namespace odeint {
+
+ class butcher_tableau_consistency_exception : public std::exception {
+
+ virtual const char* what() const throw()
+ {
+ return "Consistency Check of Butcher Tableau failed!";
+ }
+ };
+
+ class butcher_tableau_order_condition_exception : public std::exception {
+
+ virtual const char* what() const throw()
+ {
+ return "Order Condition Check of Butcher Tableau failed!";
+ }
+ };
+
+
template<
class Container ,
class Time = double ,
@@ -69,8 +90,88 @@
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()
+ {
+ typename std::vector< time_type >::const_iterator a_iter = m_a.begin();
+ typename std::vector< time_type >::const_iterator c_iter = m_c.begin();
+ typename std::vector< std::vector<time_type> >::const_iterator b_iter = m_b.begin();
+ typename std::vector<time_type>::const_iterator b_iter_iter;
+
+ // check 1: a_i = sum_j b_ij
+ while( a_iter != m_a.end() ) {
+ time_type tmp = 0.0;
+ b_iter_iter = (*b_iter).begin();
+ while( b_iter_iter != (*b_iter).end() ) {
+ tmp += *b_iter_iter++;
+ }
+ b_iter++;
+ 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.begin();
+ c_iter = m_c.begin();
+ if( k == 1 ) // special treatment for 0^0 = 1
+ tmp += *c_iter++;
+ else
+ c_iter++;
+ while( a_iter != m_a.end() ) {
+ //std::clog<<pow( *a_iter , k-1 )<< ", ";
+ tmp += (*c_iter++) * pow( *a_iter++ , k-1 );
+ }
+ //std::clog << tmp << " = " << time_type(1.0)/time_type(k) << "?" << std::endl;
+ 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
+
+ a,b,c are vectors providing the butcher tableau for the Runge Kutta scheme
+
+ 0 |
+ a_1 | b_21
+ a_2 | b_31 b_32
+ ... | ... ...
+ a_q-1 | b_q1 b_q2 ... b_qq-1
+ -------------------------------
+ | c_1 c_2 ... c_q-1 c_q
+
+ The size of c is determining the order of the scheme q.
+ a is of length q-1
+ b is of length q-1, b[0] of length 1, b[1] of length 2 and so on
+ c is of length q (which defines q).
+
+ There are 2 conditions that these parameters have to fullfill:
+ Consitency:
+
+ a_i = sum_j b_ij for i = 1 ... q-1
+
+ Condition on the order of the method (all error terms dt^k , k<q+1 cancel out):
+
+ sum_i c_i * (a_(i+1))^(k-1) = 1/k k = 1 ... q
+
+ Note, that a_0 = 1 (implicitely) and 0^0 = 1
+ so this means sum_i c_i = 1 at k=1
+ */
stepper_rk_generic( std::vector<time_type> &a,
std::vector< std::vector<time_type> > &b,
std::vector< time_type > &c)
@@ -78,6 +179,8 @@
{
m_xvec.resize(m_q);
m_xiter_vec.resize(m_q);
+
+ check_consitency();
}
order_type order() const { return m_q; }
@@ -92,15 +195,17 @@
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();
- x_iter++;
+ (*xiter_iter++) = (*x_iter++).begin();
}
m_resizer.adjust_size(x, m_xtmp);
- x_iter = m_xvec.begin();
- (*x_iter++) = dxdt;
+ x_iter = m_xvec.begin()+1;
typename std::vector< time_type >::const_iterator a_iter = m_a.begin();
typename std::vector< std::vector<time_type> >::const_iterator b_iter = m_b.begin();
@@ -109,9 +214,7 @@
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++;
+ system( m_xtmp, *x_iter++ , t + dt*(*a_iter++) );
b_iter++;
}
@@ -122,14 +225,6 @@
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 ,
Modified: sandbox/odeint/boost/numeric/odeint/stepsize_controller_standard.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepsize_controller_standard.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/stepsize_controller_standard.hpp 2009-11-15 12:21:58 EST (Sun, 15 Nov 2009)
@@ -75,26 +75,26 @@
// get the maximal value of x_err/D where
// D = eps_abs + eps_rel * (a_x*|x| + a_dxdt*|dxdt|);
T err = eps_abs + eps_rel * (a_x * std::abs(*x_start++) +
- a_dxdt * dt * std::abs(*dxdt_start++));
+ a_dxdt * dt * std::abs(*dxdt_start++));
max_rel_err = max( std::abs(*x_err_start++)/err , max_rel_err );
}
//std::cout<<max_rel_err<<std::endl;
if( max_rel_err > 1.1 ) { // error too large - decrease dt
- // limit scaling factor to 0.2
- dt *= max( 0.9*pow(max_rel_err , -1.0/stepper.order()) , 0.2 );
- // reset state
- x = x_tmp;
- return STEP_SIZE_DECREASED;
+ // limit scaling factor to 0.2
+ dt *= max( 0.9*pow(max_rel_err , -1.0/stepper.order()) , 0.2 );
+ // reset state
+ x = x_tmp;
+ return STEP_SIZE_DECREASED;
} else if( max_rel_err < 0.5 ) { //error too small - increase dt
- t += dt; // we keep the evolution -> increase time
- // limit scaling factor to 5.0
- dt *= min( 0.9*pow(max_rel_err , -1.0/(stepper.order()+1.0)), 5.0 );
- return STEP_SIZE_INCREASED;
+ t += dt; // we keep the evolution -> increase time
+ // limit scaling factor to 5.0
+ dt *= min( 0.9*pow(max_rel_err , -1.0/(stepper.order()+1.0)), 5.0 );
+ return STEP_SIZE_INCREASED;
} else {
- t += dt;
- return SUCCESS;
+ t += dt;
+ return SUCCESS;
}
}
};
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 12:21:58 EST (Sun, 15 Nov 2009)
@@ -55,7 +55,7 @@
vector< double > a(3);
a[0] = 0.5;
- a[1] = 0.5;
+ a[1] = 0.5;
a[2] = 1;
vector< vector<double> > b(3);
Modified: sandbox/odeint/libs/numeric/odeint/examples/lorenz_stepper.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/lorenz_stepper.cpp (original)
+++ sandbox/odeint/libs/numeric/odeint/examples/lorenz_stepper.cpp 2009-11-15 12:21:58 EST (Sun, 15 Nov 2009)
@@ -61,7 +61,7 @@
int main( int argc , char **argv )
{
const double dt = 0.01;
- const size_t olen = 1000000000;
+ const size_t olen = 100000000;
state_type1 x1(3);
x1[0] = 1.0;
@@ -78,18 +78,19 @@
start= clock();
t = 0.0;
for( size_t oi=0 ; oi<olen ; ++oi,t+=dt )
- stepper1.next_step( lorenz1 , x1 , t , dt );
+ stepper1.next_step( lorenz1 , x1 , t , dt );
end = clock();
cout << "vector : " << double ( end - start ) / double( CLOCKS_PER_SEC ) << endl;
+ cout << "x: "<<x1[0]<<tab<<x1[1]<<tab<<x1[2]<<endl;
start= clock();
t = 0.0;
for( size_t oi=0 ; oi<olen ; ++oi,t+=dt )
- stepper2.next_step( lorenz2 , x2 , t , dt );
+ stepper2.next_step( lorenz2 , x2 , t , dt );
end = clock();
cout << "array : " << double ( end - start ) / double( CLOCKS_PER_SEC ) << endl;
-
+ cout << "x: "<<x2[0]<<tab<<x2[1]<<tab<<x2[2]<<endl;
return 0;
}
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