|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r65710 - in sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper: . base
From: mario.mulansky_at_[hidden]
Date: 2010-10-01 09:43:48
Author: mariomulansky
Date: 2010-10-01 09:43:46 EDT (Fri, 01 Oct 2010)
New Revision: 65710
URL: http://svn.boost.org/trac/boost/changeset/65710
Log:
controlled stepper now works properly with fsal steppers
Added:
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/stepper_categories.hpp (contents, props changed)
Text files modified:
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/explicit_stepper_and_error_stepper_fsal_base.hpp | 7
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/controlled_error_stepper.hpp | 171 +++++++++++++++++++++++++++++++++++++--
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/error_checker.hpp | 4
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_error_dopri5.hpp | 5
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_error_rk54_ck.hpp | 3
5 files changed, 175 insertions(+), 15 deletions(-)
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/explicit_stepper_and_error_stepper_fsal_base.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/explicit_stepper_and_error_stepper_fsal_base.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/explicit_stepper_and_error_stepper_fsal_base.hpp 2010-10-01 09:43:46 EDT (Fri, 01 Oct 2010)
@@ -141,14 +141,15 @@
this->stepper().do_step_impl( system , in , dxdt , t , out , dt , xerr );
}
- void reset_dxdt( state_type &dxdt )
+ void reset_step( state_type &dxdt )
{
- this->stepper().reset_dxdt_impl( dxdt );
+ this->stepper().reset_step_impl( dxdt );
}
void adjust_size( const state_type &x )
{
- m_size_adjuster.adjust_size( x );
+ if( m_size_adjuster.adjust_size( x ) )
+ m_first_call = true;
}
Added: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/stepper_categories.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/stepper_categories.hpp 2010-10-01 09:43:46 EDT (Fri, 01 Oct 2010)
@@ -0,0 +1,33 @@
+/*
+ boost header: numeric/odeint/stepper_categories.hpp
+
+ Copyright 2009 Karsten Ahnert
+ Copyright 2009 Mario Mulansky
+
+ 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_CATEGORIES_HPP_
+#define BOOST_NUMERIC_ODEINT_STEPPER_CATEGORIES_HPP_
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+/*
+ * Tags to specify stepper types
+ */
+
+struct explicit_error_stepper_tag {};
+struct explicit_stepper_and_error_stepper_tag {};
+struct explicit_error_stepper_fsal_tag {};
+struct explicit_stepper_and_error_stepper_fsal_tag {};
+
+} // odeint
+} // numeric
+} // boost
+
+
+#endif /* BOOST_NUMERIC_ODEINT_STEPPER_CATEGORIES_HPP_ */
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/controlled_error_stepper.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/controlled_error_stepper.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/controlled_error_stepper.hpp 2010-10-01 09:43:46 EDT (Fri, 01 Oct 2010)
@@ -19,6 +19,7 @@
#include <boost/numeric/odeint/stepper/adjust_size.hpp>
#include <boost/numeric/odeint/stepper/error_checker.hpp>
+#include <boost/numeric/odeint/stepper/base/stepper_categories.hpp>
namespace boost {
namespace numeric {
@@ -32,14 +33,27 @@
} controlled_step_result;
+/* error stepper categroy dispatcher */
+template<
+ class ErrorStepper ,
+ class ErrorChecker = error_checker_standard< typename ErrorStepper::state_type ,
+ typename ErrorStepper::time_type ,
+ typename ErrorStepper::algebra_type ,
+ typename ErrorStepper::operations_type > ,
+ class ErrorStepperCategory = typename ErrorStepper::error_stepper_category
+>
+class controlled_error_stepper { };
+
+
+
+/****************************/
+/* explicit stepper version */
+
template<
class ErrorStepper ,
- class ErrorChecker = error_checker_standard< typename ErrorStepper::state_type ,
- typename ErrorStepper::time_type ,
- typename ErrorStepper::algebra_type ,
- typename ErrorStepper::operations_type >
+ class ErrorChecker
>
-class controlled_error_stepper : boost::noncopyable
+class controlled_error_stepper< ErrorStepper , ErrorChecker , explicit_error_stepper_tag > : boost::noncopyable
{
public:
@@ -79,7 +93,8 @@
template< class System >
- controlled_step_result try_step( System &sys , state_type &x , state_type &dxdt , time_type &t , time_type &dt )
+ controlled_step_result try_step( System &sys , state_type &x , const state_type &dxdt ,
+ time_type &t , time_type &dt )
{
using std::max;
using std::min;
@@ -87,6 +102,7 @@
m_xerr_size_adjuster.adjust_size_by_policy( x , adjust_size_policy() );
boost::numeric::odeint::copy( x , m_x_old );
+
m_stepper.do_step( sys , x , dxdt , t , dt , m_x_err );
// this potentially overwrites m_x_err! (standard_error_checker does, at least)
@@ -121,15 +137,15 @@
controlled_step_result try_step( System &sys , state_type &x , time_type &t , time_type &dt )
{
m_dxdt_size_adjuster.adjust_size_by_policy( x , adjust_size_policy() );
- sys( x , m_dxdt , t );
+ sys( x , m_dxdt ,t );
return try_step( sys , x , m_dxdt , t , dt );
}
void adjust_size( const state_type &x )
{
- m_dxdt_size_adjuster.adjust_size( x );
- m_xerr_size_adjuster.adjust_size( x );
- m_stepper.adjust_size( x );
+ m_dxdt_size_adjuster.adjust_size( x );
+ m_xerr_size_adjuster.adjust_size( x );
+ m_stepper.adjust_size( x );
}
@@ -150,6 +166,141 @@
+
+
+
+
+
+
+/*********************************/
+/* explicit stepper fsal version */
+
+template<
+ class ErrorStepper ,
+ class ErrorChecker
+ >
+class controlled_error_stepper< ErrorStepper , ErrorChecker , explicit_error_stepper_fsal_tag > : boost::noncopyable
+{
+public:
+
+ typedef ErrorStepper error_stepper_type;
+ typedef typename error_stepper_type::state_type state_type;
+ typedef typename error_stepper_type::time_type time_type;
+ typedef typename error_stepper_type::order_type order_type;
+
+ // ToDo : check if the next line can be avoided
+ typedef typename error_stepper_type::adjust_size_policy adjust_size_policy;
+
+ typedef ErrorChecker error_checker_type;
+
+ // ToDo : check if stepper could be constructed by the controlled stepper
+ controlled_error_stepper(
+ error_stepper_type &stepper ,
+ const error_checker_type &error_checker = error_checker_type()
+ )
+ : m_stepper( stepper ) , m_error_checker( error_checker ) ,
+ m_dxdt_size_adjuster() , m_xerr_size_adjuster() ,
+ m_dxdt() , m_x_old() , m_x_err() ,
+ m_first_call( true )
+ {
+ boost::numeric::odeint::construct( m_dxdt );
+ boost::numeric::odeint::construct( m_x_err );
+ boost::numeric::odeint::construct( m_x_old );
+ m_dxdt_size_adjuster.register_state( 0 , m_dxdt );
+ m_xerr_size_adjuster.register_state( 0 , m_x_err );
+ }
+
+ ~controlled_error_stepper( void )
+ {
+ boost::numeric::odeint::destruct( m_dxdt );
+ boost::numeric::odeint::destruct( m_x_err );
+ boost::numeric::odeint::destruct( m_x_old );
+ }
+
+
+
+ template< class System >
+ controlled_step_result try_step( System &sys , state_type &x , state_type &dxdt ,
+ time_type &t , time_type &dt )
+ {
+ using std::max;
+ using std::min;
+ using std::pow;
+
+ m_xerr_size_adjuster.adjust_size_by_policy( x , adjust_size_policy() );
+ boost::numeric::odeint::copy( x , m_x_old );
+ //fsal: m_stepper.get_dxdt( dxdt );
+ //fsal: m_stepper.do_step( sys , x , dxdt , t , dt , m_x_err );
+ m_stepper.do_step( sys , x , dxdt , t , dt , m_x_err );
+
+ // this potentially overwrites m_x_err! (standard_error_checker does, at least)
+ time_type max_rel_err = m_error_checker.error( m_x_old , dxdt , m_x_err , dt );
+
+ if( max_rel_err > 1.1 )
+ {
+ // error too large - decrease dt ,limit scaling factor to 0.2 and reset state
+ /* ToDo: for fsal steppers we have to do some resetting of dxdt */
+ dt *= max( 0.9 * pow( max_rel_err , -1.0 / ( m_stepper.error_order() - 1.0 ) ) , 0.2 );
+ boost::numeric::odeint::copy( m_x_old , x );
+ m_stepper.reset_step( dxdt );
+ return step_size_decreased;
+ }
+ else
+ {
+ if( max_rel_err < 0.5 )
+ {
+ //error too small - increase dt and keep the evolution and limit scaling factor to 5.0
+ t += dt;
+ dt *= min( 0.9 * pow( max_rel_err , -1.0 / m_stepper.stepper_order() ) , 5.0 );
+ return success_step_size_increased;
+ }
+ else
+ {
+ t += dt;
+ return success_step_size_unchanged;
+ }
+ }
+ }
+
+ template< class System >
+ controlled_step_result try_step( System &sys , state_type &x , time_type &t , time_type &dt )
+ {
+ if( m_dxdt_size_adjuster.adjust_size_by_policy( x , adjust_size_policy() ) || m_first_call )
+ {
+ sys( x , m_dxdt ,t );
+ m_first_call = false;
+ }
+ return try_step( sys , x , m_dxdt , t , dt );
+ }
+
+ void adjust_size( const state_type &x )
+ {
+ bool changed = false;
+ changed |= m_dxdt_size_adjuster.adjust_size( x );
+ changed |= m_xerr_size_adjuster.adjust_size( x );
+ changed |= m_stepper.adjust_size( x );
+ if( changed )
+ m_first_call = true;
+ }
+
+
+private:
+
+ error_stepper_type &m_stepper;
+ error_checker_type m_error_checker;
+
+ size_adjuster< state_type , 1 > m_dxdt_size_adjuster;
+ size_adjuster< state_type , 1 > m_xerr_size_adjuster;
+
+ state_type m_dxdt;
+ state_type m_x_old;
+ state_type m_x_err;
+ bool m_first_call;
+};
+
+
+
+
} // odeint
} // numeric
} // boost
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/error_checker.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/error_checker.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/error_checker.hpp 2010-10-01 09:43:46 EDT (Fri, 01 Oct 2010)
@@ -38,7 +38,9 @@
error_checker_standard( void ) : m_eps_abs( 1E-6 ) , m_eps_rel( 1E-6 ) , m_a_x( 1.0 ) , m_a_dxdt( 1.0 )
{}
- time_type error( const state_type &x_old , const state_type &dxdt_old , state_type &x_err , time_type dt )
+ /* ToDo: implement constructor with epsilons */
+
+ time_type error( const state_type &x_old , const state_type &dxdt_old , state_type &x_err , const time_type &dt )
{ // this overwrites x_err !
algebra_type::for_each3( x_old , dxdt_old , x_err ,
typename operations_type::rel_error( m_eps_abs , m_eps_rel , m_a_x , m_a_dxdt*dt ) );
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_error_dopri5.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_error_dopri5.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_error_dopri5.hpp 2010-10-01 09:43:46 EDT (Fri, 01 Oct 2010)
@@ -17,6 +17,7 @@
#include <boost/numeric/odeint/algebra/standard_resize.hpp>
#include <boost/numeric/odeint/stepper/base/explicit_stepper_and_error_stepper_fsal_base.hpp>
+#include <boost/numeric/odeint/stepper/base/stepper_categories.hpp>
#include <boost/numeric/odeint/stepper/detail/macros.hpp>
namespace boost {
@@ -43,6 +44,8 @@
BOOST_ODEINT_EXPLICIT_STEPPERS_AND_ERROR_STEPPERS_TYPEDEFS( explicit_error_dopri5 , 5 , 5 , 4 );
+ typedef explicit_error_stepper_fsal_tag error_stepper_category;
+
explicit_error_dopri5( void )
: m_size_adjuster() , m_x1() , m_x2() , m_x3() , m_x4() , m_x5() , m_x6()
{
@@ -168,7 +171,7 @@
typename operations_type::scale_sum6( 1.0 , dt*c1 , dt*c3 , dt*c4 , dt*c5 , dt*c6 ));
}
- void reset_dxdt_impl( state_type &dxdt )
+ void reset_step_impl( state_type &dxdt )
{
boost::numeric::odeint::copy( m_x1 , dxdt );
}
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_error_rk54_ck.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_error_rk54_ck.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_error_rk54_ck.hpp 2010-10-01 09:43:46 EDT (Fri, 01 Oct 2010)
@@ -19,6 +19,7 @@
#include <boost/numeric/odeint/algebra/standard_resize.hpp>
#include <boost/numeric/odeint/stepper/base/explicit_stepper_and_error_stepper_base.hpp>
+#include <boost/numeric/odeint/stepper/base/stepper_categories.hpp>
#include <boost/numeric/odeint/stepper/detail/macros.hpp>
namespace boost {
@@ -49,6 +50,8 @@
BOOST_ODEINT_EXPLICIT_STEPPERS_AND_ERROR_STEPPERS_TYPEDEFS( explicit_error_rk54_ck , 5 , 5 , 4);
+ typedef explicit_error_stepper_tag error_stepper_category;
+
explicit_error_rk54_ck( void ) : m_size_adjuster() , m_x1() , m_x2() , m_x3() , m_x4() , m_x5() , m_x6()
{
boost::numeric::odeint::construct( m_x1 );
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