|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r66318 - in sandbox/odeint/branches/karsten: boost/numeric/odeint/stepper boost/numeric/odeint/stepper/base libs/numeric/odeint/test
From: karsten.ahnert_at_[hidden]
Date: 2010-11-01 13:22:16
Author: karsten
Date: 2010-11-01 13:22:14 EDT (Mon, 01 Nov 2010)
New Revision: 66318
URL: http://svn.boost.org/trac/boost/changeset/66318
Log:
* dense output dopri 5
* some formating and comments
Added:
sandbox/odeint/branches/karsten/libs/numeric/odeint/test/check_dense_output_dopri5.cpp (contents, props changed)
Text files modified:
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp | 34 +++++-
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/explicit_stepper_and_error_stepper_base.hpp | 52 +++++++++-
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/explicit_stepper_and_error_stepper_fsal_base.hpp | 52 +++++++++-
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp | 25 ++++
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/controlled_error_stepper.hpp | 191 +++++++++++++++++++++++++++------------
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/dense_output_dopri5.hpp | 148 ++++++++++++++++++++++++++++++
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/dense_output_explicit_euler.hpp | 18 ++-
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_error_dopri5.hpp | 92 ++++++++++++------
sandbox/odeint/branches/karsten/libs/numeric/odeint/test/Jamfile | 1
sandbox/odeint/branches/karsten/libs/numeric/odeint/test/check_dense_output_explicit_euler.cpp | 4
10 files changed, 492 insertions(+), 125 deletions(-)
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp 2010-11-01 13:22:14 EDT (Mon, 01 Nov 2010)
@@ -54,9 +54,17 @@
static const order_type stepper_order_value = StepperOrder;
static const order_type error_order_value = ErrorOrder;
- order_type stepper_order( void ) const { return stepper_order_value; }
+ order_type stepper_order( void ) const
+ {
+ return stepper_order_value;
+ }
+
+ order_type error_order( void ) const
+ {
+ return error_order_value;
+ }
+
- order_type error_order( void ) const { return error_order_value; }
explicit_error_stepper_base( void ) : m_size_adjuster() , m_dxdt()
@@ -72,11 +80,8 @@
- stepper_type& stepper( void ) { return *static_cast< stepper_type* >( this ); }
- const stepper_type& stepper( void ) const {return *static_cast< const stepper_type* >( this );}
-
-
+ // do_step( system , x , t , dt , xerr )
template< class System >
void do_step( System &system , state_type &x , const time_type t , const time_type dt , state_type &xerr )
{
@@ -85,12 +90,14 @@
this->stepper().do_step_impl( system , x , m_dxdt , t , x , dt , xerr );
}
+ // do_step( system , x , dxdt , t , dt , xerr )
template< class System >
void do_step( System &system , state_type &x , const state_type &dxdt , const time_type t , const time_type dt , state_type &xerr )
{
this->stepper().do_step_impl( system , x , dxdt , t , x , dt , xerr );
}
+ // do_step( system , in , t , out , dt , xerr )
template< class System >
void do_step( System &system , const state_type &in , const time_type t , state_type &out , const time_type dt , state_type &xerr )
{
@@ -99,6 +106,7 @@
this->stepper().do_step_impl( system , in , m_dxdt , t , out , dt , xerr );
}
+ // do_step( system , in , dxdt , t , out , dt , xerr )
template< class System >
void do_step( System &system , const state_type &in , const state_type &dxdt , const time_type t , state_type &out , const time_type dt , state_type &xerr )
{
@@ -106,6 +114,8 @@
}
+
+
void adjust_size( const state_type &x )
{
m_size_adjuster.adjust_size( x );
@@ -114,6 +124,18 @@
private:
+ // ToDo : make the next two methods private?
+ stepper_type& stepper( void )
+ {
+ return *static_cast< stepper_type* >( this );
+ }
+
+ const stepper_type& stepper( void ) const
+ {
+ return *static_cast< const stepper_type* >( this );
+ }
+
+
size_adjuster< state_type , 1 > m_size_adjuster;
state_type m_dxdt;
};
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/explicit_stepper_and_error_stepper_base.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/explicit_stepper_and_error_stepper_base.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/explicit_stepper_and_error_stepper_base.hpp 2010-11-01 13:22:14 EDT (Mon, 01 Nov 2010)
@@ -23,6 +23,7 @@
namespace numeric {
namespace odeint {
+
/*
* base class for explicit stepper and error steppers
* models the stepper AND the error stepper concept
@@ -55,9 +56,25 @@
static const order_type stepper_order_value = StepperOrder;
static const order_type error_order_value = ErrorOrder;
- order_type order( void ) const { return order_value; }
- order_type stepper_order( void ) const { return stepper_order_value; }
- order_type error_order( void ) const { return error_order_value; }
+
+
+
+ order_type order( void ) const
+ {
+ return order_value;
+ }
+
+ order_type stepper_order( void ) const
+ {
+ return stepper_order_value;
+ }
+
+ order_type error_order( void ) const
+ {
+ return error_order_value;
+ }
+
+
explicit_stepper_and_error_stepper_base( void ) : m_size_adjuster() , m_dxdt()
@@ -73,10 +90,8 @@
- stepper_type& stepper( void ) { return *static_cast< stepper_type* >( this ); }
- const stepper_type& stepper( void ) const {return *static_cast< const stepper_type* >( this );}
-
+ // do_step( sys , x , t , dt )
template< class System >
void do_step( System &system , state_type &x , const time_type t , const time_type dt )
{
@@ -85,12 +100,14 @@
this->stepper().do_step_impl( system , x , m_dxdt , t , x , dt );
}
+ // do_step( sys , x , dxdt , t , dt )
template< class System >
void do_step( System &system , state_type &x , const state_type &dxdt , const time_type t , const time_type dt )
{
this->stepper().do_step_impl( system , x , dxdt , t , x , dt );
}
+ // do_step( sys , in , t , out , dt )
template< class System >
void do_step( System &system , const state_type &in , const time_type t , state_type &out , const time_type dt )
{
@@ -99,6 +116,7 @@
this->stepper().do_step_impl( system , in , m_dxdt , t , out , dt );
}
+ // do_step( sys , in , dxdt , t , out , dt )
template< class System >
void do_step( System &system , const state_type &in , const state_type &dxdt , const time_type t , state_type &out , const time_type dt )
{
@@ -107,6 +125,10 @@
+
+
+
+ // do_step( sys , x , t , dt , xerr )
template< class System >
void do_step( System &system , state_type &x , const time_type t , const time_type dt , state_type &xerr )
{
@@ -115,13 +137,14 @@
this->stepper().do_step_impl( system , x , m_dxdt , t , x , dt , xerr );
}
-
+ // do_step( sys , x , dxdt , t , dt , xerr )
template< class System >
void do_step( System &system , state_type &x , const state_type &dxdt , const time_type t , const time_type dt , state_type &xerr )
{
this->stepper().do_step_impl( system , x , dxdt , t , x , dt , xerr );
}
+ // do_step( sys , in , t , out , dt , xerr )
template< class System >
void do_step( System &system , const state_type &in , const time_type t , state_type &out , const time_type dt , state_type &xerr )
{
@@ -130,12 +153,16 @@
this->stepper().do_step_impl( system , in , m_dxdt , t , out , dt , xerr );
}
+ // do_step( sys , in , dxdt , t , out , dt , xerr )
template< class System >
void do_step( System &system , const state_type &in , const state_type &dxdt , const time_type t , state_type &out , const time_type dt , state_type &xerr )
{
this->stepper().do_step_impl( system , in , dxdt , t , out , dt , xerr );
}
+
+
+
void adjust_size( const state_type &x )
{
m_size_adjuster.adjust_size( x );
@@ -144,6 +171,17 @@
private:
+ stepper_type& stepper( void )
+ {
+ return *static_cast< stepper_type* >( this );
+ }
+
+ const stepper_type& stepper( void ) const
+ {
+ return *static_cast< const stepper_type* >( this );
+ }
+
+
size_adjuster< state_type , 1 > m_size_adjuster;
state_type m_dxdt;
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-11-01 13:22:14 EDT (Mon, 01 Nov 2010)
@@ -55,12 +55,28 @@
static const order_type stepper_order_value = StepperOrder;
static const order_type error_order_value = ErrorOrder;
- order_type order( void ) const { return order_value; }
- order_type stepper_order( void ) const { return stepper_order_value; }
- order_type error_order( void ) const { return error_order_value; }
- explicit_stepper_and_error_stepper_fsal_base( void ) : m_size_adjuster() , m_dxdt() , m_first_call( true )
+
+ order_type order( void ) const
+ {
+ return order_value;
+ }
+
+ order_type stepper_order( void ) const
+ {
+ return stepper_order_value;
+ }
+
+ order_type error_order( void ) const
+ {
+ return error_order_value;
+ }
+
+
+
+
+ explicit_stepper_and_error_stepper_fsal_base( void ) : m_size_adjuster() , m_dxdt() , m_first_call( true )
{
boost::numeric::odeint::construct( m_dxdt );
m_size_adjuster.register_state( 0 , m_dxdt );
@@ -73,10 +89,9 @@
- stepper_type& stepper( void ) { return *static_cast< stepper_type* >( this ); }
- const stepper_type& stepper( void ) const {return *static_cast< const stepper_type* >( this );}
+ // do_step( sys , x , t , dt )
template< class System >
void do_step( System &system , state_type &x , const time_type t , const time_type dt )
{
@@ -85,12 +100,14 @@
this->stepper().do_step_impl( system , x , m_dxdt , t , x , dt );
}
+ // do_step( sys , x , dxdt , t , dt )
template< class System >
void do_step( System &system , state_type &x , const state_type &dxdt , const time_type t , const time_type dt )
{
this->stepper().do_step_impl( system , x , dxdt , t , x , dt );
}
+ // do_step( sys , in , t , out , dt )
template< class System >
void do_step( System &system , const state_type &in , const time_type t , state_type &out , const time_type dt )
{
@@ -99,6 +116,7 @@
this->stepper().do_step_impl( system , in , m_dxdt , t , out , dt );
}
+ // do_step( sys , in , dxdt , t , out , dt )
template< class System >
void do_step( System &system , const state_type &in , const state_type &dxdt , const time_type t , state_type &out , const time_type dt )
{
@@ -107,6 +125,10 @@
+
+
+
+ // do_step( sys , x , t , dt , xerr )
template< class System >
void do_step( System &system , state_type &x , const time_type t , const time_type dt , state_type &xerr )
{
@@ -118,12 +140,14 @@
this->stepper().do_step_impl( system , x , m_dxdt , t , x , dt , xerr );
}
+ // do_step( sys , x , dxdt , t , dt , xerr )
template< class System >
void do_step( System &system , state_type &x , state_type &dxdt , const time_type t , const time_type dt , state_type &xerr )
{
this->stepper().do_step_impl( system , x , dxdt , t , x , dt , xerr );
}
+ // do_step( sys , in , t , out , dt , xerr )
template< class System >
void do_step( System &system , const state_type &in , const time_type t , state_type &out , const time_type dt , state_type &xerr )
{
@@ -135,17 +159,22 @@
this->stepper().do_step_impl( system , in , m_dxdt , t , out , dt , xerr );
}
+ // do_step( sys , in , dxdt , t , out , dt , xerr )
template< class System >
void do_step( System &system , const state_type &in , state_type &dxdt , const time_type t , state_type &out , const time_type dt , state_type &xerr )
{
this->stepper().do_step_impl( system , in , dxdt , t , out , dt , xerr );
}
+
+
void reset_step( state_type &dxdt )
{
this->stepper().reset_step_impl( dxdt );
}
+
+
void adjust_size( const state_type &x )
{
if( m_size_adjuster.adjust_size( x ) )
@@ -155,6 +184,17 @@
private:
+ stepper_type& stepper( void )
+ {
+ return *static_cast< stepper_type* >( this );
+ }
+
+ const stepper_type& stepper( void ) const
+ {
+ return *static_cast< const stepper_type* >( this );
+ }
+
+
size_adjuster< state_type , 1 > m_size_adjuster;
state_type m_dxdt;
bool m_first_call;
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp 2010-11-01 13:22:14 EDT (Mon, 01 Nov 2010)
@@ -53,7 +53,11 @@
typedef unsigned short order_type;
static const order_type order_value = Order;
- order_type order( void ) const { return order_value; }
+
+ order_type order( void ) const
+ {
+ return order_value;
+ }
explicit_stepper_base( void ) : m_size_adjuster() , m_dxdt()
@@ -68,11 +72,9 @@
}
- stepper_type& stepper( void ) { return *static_cast< stepper_type* >( this ); }
- const stepper_type& stepper( void ) const {return *static_cast< const stepper_type* >( this );}
-
+ // do_step( sys , x , t , dt )
template< class System >
void do_step( System &system , state_type &x , const time_type t , const time_type dt )
{
@@ -81,12 +83,14 @@
this->stepper().do_step_impl( system , x , m_dxdt , t , x , dt );
}
+ // do_step( sys , x , dxdt , t , dt )
template< class System >
void do_step( System &system , state_type &x , const state_type dxdt , const time_type t , const time_type dt )
{
this->stepper().do_step_impl( system , x , dxdt , t , x , dt );
}
+ // do_step( sys , in , t , out , dt )
template< class System >
void do_step( System &system , const state_type &in , const time_type t , state_type &out , const time_type dt )
{
@@ -95,6 +99,7 @@
this->stepper().do_step_impl( system , in , m_dxdt , t , out , dt );
}
+ // do_step( sys , in , dxdt , t , out , dt )
template< class System >
void do_step( System &system , const state_type &in , const state_type &dxdt , const time_type t , state_type &out , const time_type dt )
{
@@ -102,6 +107,7 @@
}
+
void adjust_size( const state_type &x )
{
m_size_adjuster.adjust_size( x );
@@ -110,6 +116,17 @@
protected:
+ stepper_type& stepper( void )
+ {
+ return *static_cast< stepper_type* >( this );
+ }
+
+ const stepper_type& stepper( void ) const
+ {
+ return *static_cast< const stepper_type* >( this );
+ }
+
+
size_adjuster< state_type , 1 > m_size_adjuster;
state_type m_dxdt;
};
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-11-01 13:22:14 EDT (Mon, 01 Nov 2010)
@@ -33,7 +33,9 @@
} controlled_step_result;
-/* error stepper category dispatcher */
+/*
+ * error stepper category dispatcher
+ */
template<
class ErrorStepper ,
class ErrorChecker = error_checker_standard< typename ErrorStepper::state_type ,
@@ -46,9 +48,10 @@
-/****************************/
-/* explicit stepper version */
+/*
+ * explicit stepper version
+ */
template<
class ErrorStepper ,
class ErrorChecker
@@ -65,53 +68,84 @@
typedef ErrorChecker error_checker_type;
- // ToDo : check if stepper could be constructed by the controlled stepper
+
controlled_error_stepper(
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_dxdt_size_adjuster() , m_xerr_size_adjuster() , m_xnew_size_adjuster() ,
+ m_dxdt() , m_xerr() , m_xnew()
{
boost::numeric::odeint::construct( m_dxdt );
- boost::numeric::odeint::construct( m_x_err );
- boost::numeric::odeint::construct( m_x_old );
+ boost::numeric::odeint::construct( m_xerr );
+ boost::numeric::odeint::construct( m_xnew );
m_dxdt_size_adjuster.register_state( 0 , m_dxdt );
- m_xerr_size_adjuster.register_state( 0 , m_x_err );
+ m_xerr_size_adjuster.register_state( 0 , m_xerr );
+ m_xnew_size_adjuster.register_state( 0 , m_xnew );
}
~controlled_error_stepper( void )
{
boost::numeric::odeint::destruct( m_dxdt );
- boost::numeric::odeint::destruct( m_x_err );
- boost::numeric::odeint::destruct( m_x_old );
+ boost::numeric::odeint::destruct( m_xerr );
+ boost::numeric::odeint::destruct( m_xnew );
+ }
+
+
+
+
+ // try_step( sys , x , t , dt )
+ template< class System >
+ 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 );
+ return try_step( sys , x , m_dxdt , t , dt );
}
+ // try_step( sys , x , dxdt , t , dt )
+ template< class System >
+ controlled_step_result try_step( System &sys , state_type &x , const state_type &dxdt , time_type &t , time_type &dt )
+ {
+ m_xnew_size_adjuster.adjust_size_by_policy( x , adjust_size_policy() );
+ controlled_step_result res = try_step( sys , x , dxdt , t , m_xnew , dt );
+ if( ( res == success_step_size_increased ) || ( res == success_step_size_unchanged ) )
+ {
+ boost::numeric::odeint::copy( m_xnew , x );
+ }
+ return res;
+ }
+
+ // try_step( sys , in , t , out , dt )
+ template< class System >
+ controlled_step_result try_step( System &sys , const state_type &in , time_type &t , state_type &out , time_type &dt )
+ {
+ m_dxdt_size_adjuster.adjust_size_by_policy( in , adjust_size_policy() );
+ sys( in , m_dxdt , t );
+ return try_step( sys , in , m_dxdt , t , out , dt );
+ }
+ // try_step( sys , in , dxdt , t , out , dt )
template< class System >
- controlled_step_result try_step( System &sys , state_type &x , const state_type &dxdt ,
- time_type &t , time_type &dt )
+ controlled_step_result try_step( System &sys , const state_type &in , const state_type &dxdt , time_type &t , state_type &out , 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 );
+ m_xerr_size_adjuster.adjust_size_by_policy( in , adjust_size_policy() );
- m_stepper.do_step( sys , x , dxdt , t , dt , m_x_err );
+ // do one step with error calculation
+ m_stepper.do_step( sys , in , dxdt , t , out , dt , m_xerr );
- // 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 );
+ time_type max_rel_err = m_error_checker.error( in , dxdt , m_xerr , 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 );
return step_size_decreased;
}
else
@@ -131,21 +165,20 @@
}
}
- template< class System >
- 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 );
- 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_xnew_size_adjuster.adjust_size( x );
m_stepper.adjust_size( x );
}
+
+
+
stepper_type& stepper( void )
{
return m_stepper;
@@ -164,10 +197,11 @@
size_adjuster< state_type , 1 > m_dxdt_size_adjuster;
size_adjuster< state_type , 1 > m_xerr_size_adjuster;
+ size_adjuster< state_type , 1 > m_xnew_size_adjuster;
state_type m_dxdt;
- state_type m_x_old;
- state_type m_x_err;
+ state_type m_xerr;
+ state_type m_xnew;
};
@@ -179,10 +213,11 @@
-
-/*********************************/
-/* explicit stepper fsal version */
-
+/*
+ * explicit stepper fsal version
+ *
+ * ToDo : introduce the same functions as for the above stepper
+ */
template<
class ErrorStepper ,
class ErrorChecker
@@ -198,55 +233,93 @@
typedef typename 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(
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_dxdt_size_adjuster() , m_xerr_size_adjuster() , m_xnew_size_adjuster() ,
+ m_dxdt() , m_xerr() , m_xnew() ,
m_first_call( true )
{
boost::numeric::odeint::construct( m_dxdt );
- boost::numeric::odeint::construct( m_x_err );
- boost::numeric::odeint::construct( m_x_old );
+ boost::numeric::odeint::construct( m_xerr );
+ boost::numeric::odeint::construct( m_xnew );
m_dxdt_size_adjuster.register_state( 0 , m_dxdt );
- m_xerr_size_adjuster.register_state( 0 , m_x_err );
+ m_xerr_size_adjuster.register_state( 0 , m_xerr );
+ m_xnew_size_adjuster.register_state( 0 , m_xnew );
}
~controlled_error_stepper( void )
{
boost::numeric::odeint::destruct( m_dxdt );
- boost::numeric::odeint::destruct( m_x_err );
- boost::numeric::odeint::destruct( m_x_old );
+ boost::numeric::odeint::destruct( m_xerr );
+ boost::numeric::odeint::destruct( m_xnew );
}
+
+
+ // try_step( sys , x , t , dt )
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 , 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 );
+ }
+
+ // try_step( sys , x , dxdt , t , dt )
+ template< class System >
+ controlled_step_result try_step( System &sys , state_type &x , state_type &dxdt , time_type &t , time_type &dt )
+ {
+ m_xnew_size_adjuster.adjust_size_by_policy( x , adjust_size_policy() );
+ controlled_step_result res = try_step( sys , x , dxdt , t , m_xnew , dt );
+ if( ( res == success_step_size_increased ) || ( res == success_step_size_unchanged) )
+ {
+ boost::numeric::odeint::copy( m_xnew , x );
+ }
+ return res;
+ }
+
+ // try_step( sys , in , t , out , dt );
+ template< class System >
+ controlled_step_result try_step( System &sys , const state_type &in , time_type &t , state_type &out , time_type &dt )
+ {
+ if( m_dxdt_size_adjuster.adjust_size_by_policy( in , adjust_size_policy() ) || m_first_call )
+ {
+ sys( in , m_dxdt ,t );
+ m_first_call = false;
+ }
+ return try_step( sys , in , m_dxdt , t , out , dt );
+ }
+
+
+ // try_step( sys , in , dxdt , t , out , dt )
+ template< class System >
+ controlled_step_result try_step( System &sys , const state_type &in , state_type &dxdt , time_type &t , state_type &out , 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 );
+ m_xerr_size_adjuster.adjust_size_by_policy( in , adjust_size_policy() );
+
//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 );
+ m_stepper.do_step( sys , in , dxdt , t , out , dt , m_xerr );
// 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 );
+ time_type max_rel_err = m_error_checker.error( in , dxdt , m_xerr , 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;
}
@@ -267,16 +340,9 @@
}
}
- 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 )
{
@@ -288,6 +354,8 @@
m_first_call = true;
}
+
+
stepper_type& stepper( void )
{
return m_stepper;
@@ -307,10 +375,11 @@
size_adjuster< state_type , 1 > m_dxdt_size_adjuster;
size_adjuster< state_type , 1 > m_xerr_size_adjuster;
+ size_adjuster< state_type , 1 > m_xnew_size_adjuster;
state_type m_dxdt;
- state_type m_x_old;
- state_type m_x_err;
+ state_type m_xerr;
+ state_type m_xnew;
bool m_first_call;
};
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/dense_output_dopri5.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/dense_output_dopri5.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/dense_output_dopri5.hpp 2010-11-01 13:22:14 EDT (Mon, 01 Nov 2010)
@@ -12,7 +12,19 @@
#ifndef BOOST_BOOST_NUMERIC_ODEINT_DENSE_OUTPUT_DOPRI5_HPP_INCLUDED
#define BOOST_BOOST_NUMERIC_ODEINT_DENSE_OUTPUT_DOPRI5_HPP_INCLUDED
-#include <boost/numeric/odeint/stepper/dense_output_dopri5.hpp>
+#include <iostream>
+#define tab "\t"
+using std::cout;
+using std::cerr;
+using std::clog;
+using std::endl;
+
+#include <stdexcept>
+
+#include <boost/static_assert.hpp>
+#include <boost/type_traits.hpp>
+
+#include <boost/numeric/odeint/stepper/explicit_error_dopri5.hpp>
namespace boost {
namespace numeric {
@@ -22,10 +34,142 @@
template< class ControlledStepper >
class dense_output_dopri5
{
-private:
+public:
typedef ControlledStepper stepper_type;
+ typedef typename stepper_type::stepper_type dopri5_type;
+ typedef typename dopri5_type::state_type state_type;
+ typedef typename dopri5_type::time_type time_type;
+ typedef typename dopri5_type::algebra_type algebra_type;
+ typedef typename dopri5_type::operations_type operations_type;
+ typedef typename dopri5_type::adjust_size_policy adjust_size_policy;
+
+ BOOST_STATIC_ASSERT(( boost::is_same<
+ dopri5_type ,
+ explicit_error_dopri5< state_type , time_type , algebra_type , operations_type , adjust_size_policy >
+ >::value ));
+
+ dense_output_dopri5( stepper_type &stepper )
+ : m_stepper( stepper ) , m_size_adjuster() ,
+ m_x1() , m_x2() , m_current_state( &m_x1 ) , m_old_state( &m_x2 ) ,
+ m_t( 0.0 ) , m_t_old( 0.0 ) , m_dt( 1.0 )
+ {
+ boost::numeric::odeint::construct( m_x1 );
+ boost::numeric::odeint::construct( m_x2 );
+ m_size_adjuster.register_state( 0 , m_x1 );
+ m_size_adjuster.register_state( 1 , m_x2 );
+ }
+
+ ~dense_output_dopri5( void )
+ {
+ boost::numeric::odeint::destruct( m_x1 );
+ boost::numeric::odeint::destruct( m_x2 );
+ }
+
+ void adjust_size( const state_type &x )
+ {
+ m_size_adjuster.adjust_size( x );
+ m_stepper.adjust_size( x );
+ }
+
+ void initialize( const state_type &x0 , const time_type t0 , const time_type dt0 )
+ {
+ boost::numeric::odeint::copy( x0 , *m_current_state );
+ m_t = t0;
+ m_dt = dt0;
+ }
+
+ template< class System >
+ std::pair< time_type , time_type > do_step( System &system )
+ {
+ const size_t max_count = 1000;
+ controlled_step_result res;
+ m_t_old = m_t;
+ size_t count = 0;
+ do
+ {
+ res = m_stepper.try_step( system , *m_current_state , m_t , *m_old_state , m_dt );
+ if( count++ == max_count )
+ throw std::overflow_error( "dense_output_dopri5 : too much iterations!");
+ }
+ while( res == step_size_decreased );
+ std::swap( m_current_state , m_old_state );
+ return std::make_pair( m_t_old , m_t );
+
+ }
+
+ /*
+ * Calculates Dense-Output for Dopri5
+ *
+ * See Hairer, Norsett, Wanner: Solving Ordinary Differential Equations, Nonstiff Problems. I, p.191/192
+ *
+ * y(t+theta) = y(t) + h * sum_i^7 b_i(theta) * k_i
+ *
+ * A = theta^2 * ( 3 - 2 theta )
+ * B = theta^2 * ( theta - 1 )
+ * C = theta^2 * ( theta - 1 )^2
+ * D = theta * ( theta - 1 )^2
+ *
+ * b_1( theta ) = A * b_1 + C * X1( theta ) + D
+ * b_2( theta ) = 0
+ * b_3( theta ) = A * b_3 + C * X3( theta )
+ * b_4( theta ) = A * b_4 + C * X4( theta )
+ * b_5( theta ) = A * b_5 + C * X5( theta )
+ * b_6( theta ) = A * b_6 + C * X6( theta )
+ * b_7( theta ) = B + C * X7( theta )
+ */
+ void calc_state( time_type t , state_type &x )
+ {
+
+ const time_type b1 = static_cast<time_type> ( 35.0 ) / static_cast<time_type>( 384.0 );
+ const time_type b3 = static_cast<time_type> ( 500.0 ) / static_cast<time_type>( 1113.0 );
+ const time_type b4 = static_cast<time_type> ( 125.0 ) / static_cast<time_type>( 192.0 );
+ const time_type b5 = static_cast<time_type> ( -2187.0 ) / static_cast<time_type>( 6784.0 );
+ const time_type b6 = static_cast<time_type> ( 11.0 ) / static_cast<time_type>( 84.0 );
+
+
+ time_type theta = t - m_t_old;
+ time_type X1 = static_cast< time_type >( 5.0 ) * ( static_cast< time_type >( 2558722523.0 ) - static_cast< time_type >( 31403016.0 ) * theta ) / static_cast< time_type >( 11282082432.0 );
+ time_type X3 = static_cast< time_type >( 100.0 ) * ( static_cast< time_type >( 882725551.0 ) - static_cast< time_type >( 15701508.0 ) * theta ) / static_cast< time_type >( 32700410799.0 );
+ time_type X4 = static_cast< time_type >( 25.0 ) * ( static_cast< time_type >( 443332067.0 ) - static_cast< time_type >( 31403016.0 ) * theta ) / static_cast< time_type >( 1880347072.0 ) ;
+ time_type X5 = static_cast< time_type >( 32805.0 ) * ( static_cast< time_type >( 23143187.0 ) - static_cast< time_type >( 3489224.0 ) * theta ) / static_cast< time_type >( 199316789632.0 );
+ time_type X6 = static_cast< time_type >( 55.0 ) * ( static_cast< time_type >( 29972135.0 ) - static_cast< time_type >( 7076736.0 ) * theta ) / static_cast< time_type >( 822651844.0 );
+ time_type X7 = static_cast< time_type >( 10.0 ) * ( static_cast< time_type >( 7414447.0 ) - static_cast< time_type >( 829305.0 ) * theta ) / static_cast< time_type >( 29380423.0 );
+
+ time_type theta_m_1 = theta - static_cast< time_type >( 1.0 );
+ time_type theta_sq = theta * theta;
+ time_type A = theta_sq * ( static_cast< time_type >( 3.0 ) - static_cast< time_type >( 2.0 ) * theta );
+ time_type B = theta_sq * theta_m_1;
+ time_type C = theta_sq * theta_m_1 * theta_m_1;
+ time_type D = theta * theta_m_1 * theta_m_1;
+
+// time_type b1_theta = A *
+
+
+
+
+
+// algebra_type::for_each3( x , *m_old_state , m_euler.m_dxdt , typename operations_type::scale_sum2( 1.0 , delta ) );
+ }
+
+ const state_type& current_state( void ) const { return *m_current_state; }
+ const time_type& current_time( void ) const { return m_t; }
+ const time_type& previous_state( void ) const { return *m_old_state; }
+ const time_type& previous_time( void ) const { return m_t_old; }
+
+
+private:
+
+
+
+ stepper_type &m_stepper;
+ size_adjuster< state_type , 2 > m_size_adjuster;
+ state_type m_x1 , m_x2;
+ state_type *m_current_state , *m_old_state;
+ time_type m_t , m_t_old , m_dt;
+
+
};
} // namespace odeint
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/dense_output_explicit_euler.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/dense_output_explicit_euler.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/dense_output_explicit_euler.hpp 2010-11-01 13:22:14 EDT (Mon, 01 Nov 2010)
@@ -55,19 +55,18 @@
boost::numeric::odeint::destruct( m_x2 );
}
- void adjust_size( const state_type &x )
- {
- m_size_adjuster.adjust_size( x );
- m_euler.adjust_size( x );
- }
- void initialize( const state_type x0 , const time_type t0 , const time_type dt0 )
+
+
+ void initialize( const state_type &x0 , const time_type t0 , const time_type dt0 )
{
boost::numeric::odeint::copy( x0 , *m_current_state );
m_t = t0;
m_dt = dt0;
}
+
+
template< class System >
std::pair< time_type , time_type > do_step( System &system )
{
@@ -84,6 +83,13 @@
algebra_type::for_each3( x , *m_old_state , m_euler.m_dxdt , typename operations_type::scale_sum2( 1.0 , delta ) );
}
+ void adjust_size( const state_type &x )
+ {
+ m_size_adjuster.adjust_size( x );
+ m_euler.adjust_size( x );
+ }
+
+
const state_type& current_state( void ) const { return *m_current_state; }
const time_type& current_time( void ) const { return m_t; }
const time_type& previous_state( void ) const { return *m_old_state; }
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-11-01 13:22:14 EDT (Mon, 01 Nov 2010)
@@ -74,37 +74,6 @@
}
- template< class System >
- void do_step_impl( System system , const state_type &in , state_type &dxdt , const time_type t , state_type &out , const time_type dt , state_type &xerr )
- {
-
- const time_type c1 = static_cast<time_type> ( 35.0 ) / static_cast<time_type>( 384.0 );
- const time_type c3 = static_cast<time_type> ( 500.0 ) / static_cast<time_type>( 1113.0 );
- const time_type c4 = static_cast<time_type> ( 125.0 ) / static_cast<time_type>( 192.0 );
- const time_type c5 = static_cast<time_type> ( -2187.0 ) / static_cast<time_type>( 6784.0 );
- const time_type c6 = static_cast<time_type> ( 11.0 ) / static_cast<time_type>( 84.0 );
-
- const time_type dc1 = c1 - static_cast<time_type> ( 5179.0 ) / static_cast<time_type>( 57600.0 );
- const time_type dc3 = c3 - static_cast<time_type> ( 7571.0 ) / static_cast<time_type>( 16695.0 );
- const time_type dc4 = c4 - static_cast<time_type> ( 393.0 ) / static_cast<time_type>( 640.0 );
- const time_type dc5 = c5 - static_cast<time_type> ( -92097.0 ) / static_cast<time_type>( 339200.0 );
- const time_type dc6 = c6 - static_cast<time_type> ( 187.0 ) / static_cast<time_type>( 2100.0 );
- const time_type dc7 = static_cast<time_type>( -0.025 );
-
- do_step_impl( system , in , dxdt , t , out , dt );
-
- // we need to copy the old dxdt
- boost::numeric::odeint::copy( dxdt , m_x1 );
-
- // store the new result in dxdt
- system( out , dxdt , t + dt );
-
- //error estimate
- algebra_type::for_each7( xerr , m_x1 , m_x3 , m_x4 , m_x5 , m_x6 , dxdt ,
- typename operations_type::scale_sum6( dt*dc1 , dt*dc3 , dt*dc4 , dt*dc5 , dt*dc6 , dt*dc7 ) );
-
-
- }
@@ -171,6 +140,67 @@
typename operations_type::scale_sum6( 1.0 , dt*c1 , dt*c3 , dt*c4 , dt*c5 , dt*c6 ));
}
+ template< class System >
+ void do_step_impl( System system , const state_type &in , state_type &dxdt , const time_type t , state_type &out , const time_type dt , state_type &xerr )
+ {
+
+ const time_type c1 = static_cast<time_type> ( 35.0 ) / static_cast<time_type>( 384.0 );
+ const time_type c3 = static_cast<time_type> ( 500.0 ) / static_cast<time_type>( 1113.0 );
+ const time_type c4 = static_cast<time_type> ( 125.0 ) / static_cast<time_type>( 192.0 );
+ const time_type c5 = static_cast<time_type> ( -2187.0 ) / static_cast<time_type>( 6784.0 );
+ const time_type c6 = static_cast<time_type> ( 11.0 ) / static_cast<time_type>( 84.0 );
+
+ const time_type dc1 = c1 - static_cast<time_type> ( 5179.0 ) / static_cast<time_type>( 57600.0 );
+ const time_type dc3 = c3 - static_cast<time_type> ( 7571.0 ) / static_cast<time_type>( 16695.0 );
+ const time_type dc4 = c4 - static_cast<time_type> ( 393.0 ) / static_cast<time_type>( 640.0 );
+ const time_type dc5 = c5 - static_cast<time_type> ( -92097.0 ) / static_cast<time_type>( 339200.0 );
+ const time_type dc6 = c6 - static_cast<time_type> ( 187.0 ) / static_cast<time_type>( 2100.0 );
+ const time_type dc7 = static_cast<time_type>( -0.025 );
+
+ do_step_impl( system , in , dxdt , t , out , dt );
+
+ // we need to copy the old dxdt
+ boost::numeric::odeint::copy( dxdt , m_x1 );
+
+ // store the new result in dxdt
+ system( out , dxdt , t + dt );
+
+ //error estimate
+ algebra_type::for_each7( xerr , m_x1 , m_x3 , m_x4 , m_x5 , m_x6 , dxdt ,
+ typename operations_type::scale_sum6( dt*dc1 , dt*dc3 , dt*dc4 , dt*dc5 , dt*dc6 , dt*dc7 ) );
+ }
+
+// template< class System >
+// void do_step_impl( System system , const state_type &in , const state_type &dxdt_in , const time_type t , state_type &out , state_type &dxdt_out , const time_type dt , state_type &xerr )
+// {
+//
+// const time_type c1 = static_cast<time_type> ( 35.0 ) / static_cast<time_type>( 384.0 );
+// const time_type c3 = static_cast<time_type> ( 500.0 ) / static_cast<time_type>( 1113.0 );
+// const time_type c4 = static_cast<time_type> ( 125.0 ) / static_cast<time_type>( 192.0 );
+// const time_type c5 = static_cast<time_type> ( -2187.0 ) / static_cast<time_type>( 6784.0 );
+// const time_type c6 = static_cast<time_type> ( 11.0 ) / static_cast<time_type>( 84.0 );
+//
+// const time_type dc1 = c1 - static_cast<time_type> ( 5179.0 ) / static_cast<time_type>( 57600.0 );
+// const time_type dc3 = c3 - static_cast<time_type> ( 7571.0 ) / static_cast<time_type>( 16695.0 );
+// const time_type dc4 = c4 - static_cast<time_type> ( 393.0 ) / static_cast<time_type>( 640.0 );
+// const time_type dc5 = c5 - static_cast<time_type> ( -92097.0 ) / static_cast<time_type>( 339200.0 );
+// const time_type dc6 = c6 - static_cast<time_type> ( 187.0 ) / static_cast<time_type>( 2100.0 );
+// const time_type dc7 = static_cast<time_type>( -0.025 );
+//
+// do_step_impl( system , in , dxdt_in , t , out , dt );
+//
+// // store the new result in dxdt
+// system( out , dxdt_out , t + dt );
+//
+// //error estimate
+// algebra_type::for_each7( xerr , m_x1 , m_x3 , m_x4 , m_x5 , m_x6 , dxdt_out ,
+// typename operations_type::scale_sum6( dt*dc1 , dt*dc3 , dt*dc4 , dt*dc5 , dt*dc6 , dt*dc7 ) );
+// }
+
+
+
+
+
void reset_step_impl( state_type &dxdt )
{
boost::numeric::odeint::copy( m_x1 , dxdt );
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/test/Jamfile
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/test/Jamfile (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/test/Jamfile 2010-11-01 13:22:14 EDT (Mon, 01 Nov 2010)
@@ -29,6 +29,7 @@
check_resize.cpp
check_implicit_euler.cpp
check_dense_output_explicit_euler.cpp
+ check_dense_output_dopri5.cpp
:
: <link>static
;
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/test/check_dense_output_dopri5.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/test/check_dense_output_dopri5.cpp 2010-11-01 13:22:14 EDT (Mon, 01 Nov 2010)
@@ -0,0 +1,85 @@
+/*
+ * check_dense_output_dopri5.cpp
+ *
+ * Created on: Nov 1, 2010
+ * Author: karsten
+ */
+
+#include <tr1/array>
+#include <fstream>
+#include <iostream>
+
+#include <boost/test/unit_test.hpp>
+
+#include <boost/numeric/odeint.hpp>
+#include <boost/numeric/odeint/stepper/implicit_euler.hpp>
+
+using namespace boost::unit_test;
+using namespace boost::numeric::odeint;
+
+typedef double value_type;
+typedef std::tr1::array< double , 2 > state_type;
+
+inline std::ostream& operator<<( std::ostream &out , const state_type &x )
+{
+ out << x[0] << "\t" << x[1];
+ return out;
+}
+
+inline void sys( const state_type &x , state_type &dxdt , const value_type t )
+{
+ dxdt[0] = x[1];
+ dxdt[1] = -x[0];
+}
+
+
+BOOST_AUTO_TEST_SUITE( dense_output_dopri5_test )
+
+BOOST_AUTO_TEST_CASE( test_dopri5 )
+{
+ using std::abs;
+
+ typedef explicit_error_dopri5< state_type > dopri5_type;
+ typedef controlled_error_stepper< dopri5_type > controlled_stepper_type;
+ dopri5_type dopri5;
+ controlled_stepper_type controlled_stepper( dopri5 );
+ dense_output_dopri5< controlled_stepper_type > stepper( controlled_stepper );
+
+ state_type x0;
+ x0[0] = 0.0;
+ x0[1] = 1.0;
+
+ stepper.initialize( x0 , 0.0 , 0.1 );
+// stepper.do_step( sys );
+
+ std::ofstream stepper_out( "dopri5_stepper_states.dat" );
+ std::ofstream states_out( "dopri5_states.dat" );
+
+
+ double t = stepper.current_time();
+ double t_end = 10.0;
+ double dt = 0.02;
+ state_type x;
+ while( t < t_end )
+ {
+ if( t < stepper.current_time() )
+ {
+ stepper.calc_state( t , x );
+ states_out << t << "\t" << x << std::endl;
+ }
+ else
+ {
+ stepper.do_step( sys );
+ stepper_out << stepper.current_time() << "\t" << stepper.current_state() << std::endl;
+ continue;
+ }
+ t += dt;
+ }
+
+// // compare with analytic solution of above system
+// BOOST_CHECK_MESSAGE( abs( x(0) - 20.0/81.0 ) < eps , x(0) - 20.0/81.0 );
+// BOOST_CHECK_MESSAGE( abs( x(1) - 10.0/9.0 ) < eps , x(0) - 10.0/9.0 );
+
+}
+
+BOOST_AUTO_TEST_SUITE_END()
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/test/check_dense_output_explicit_euler.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/test/check_dense_output_explicit_euler.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/test/check_dense_output_explicit_euler.cpp 2010-11-01 13:22:14 EDT (Mon, 01 Nov 2010)
@@ -25,13 +25,13 @@
typedef double value_type;
typedef std::tr1::array< double , 2 > state_type;
-std::ostream& operator<<( std::ostream &out , const state_type &x )
+inline std::ostream& operator<<( std::ostream &out , const state_type &x )
{
out << x[0] << "\t" << x[1];
return out;
}
-void sys( const state_type &x , state_type &dxdt , const value_type t )
+inline void sys( const state_type &x , state_type &dxdt , const value_type t )
{
dxdt[0] = x[1];
dxdt[1] = -x[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