|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r70758 - sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper
From: karsten.ahnert_at_[hidden]
Date: 2011-03-30 16:24:47
Author: karsten
Date: 2011-03-30 16:24:46 EDT (Wed, 30 Mar 2011)
New Revision: 70758
URL: http://svn.boost.org/trac/boost/changeset/70758
Log:
rosenbrock4 dense output
Added:
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/rosenbrock4_dense_output.hpp (contents, props changed)
Text files modified:
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/rosenbrock4.hpp | 45 +++++++++++-
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp | 138 +++++++++++----------------------------
2 files changed, 79 insertions(+), 104 deletions(-)
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/rosenbrock4.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/rosenbrock4.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/rosenbrock4.hpp 2011-03-30 16:24:46 EDT (Wed, 30 Mar 2011)
@@ -58,8 +58,10 @@
a41 ( 0.3314825187068521e+01 ) , a42 ( 0.2896124015972201e+01 ) , a43 ( 0.9986419139977817e+00 ) ,
c51 ( 0.7496443313967647e+01 ) , c52 ( -0.1024680431464352e+02 ) , c53 ( -0.3399990352819905e+02 ) , c54 ( 0.1170890893206160e+02 ) ,
a51 ( 0.1221224509226641e+01 ) , a52 ( 0.6019134481288629e+01 ) , a53 ( 0.1253708332932087e+02 ) , a54 ( -0.6878860361058950e+00 ) ,
- c61 ( 0.8083246795921522e+01 ) , c62 ( -0.7981132988064893e+01 ) , c63 ( -0.3152159432874371e+02 ) , c64 ( 0.1631930543123136e+02 ) , c65 ( -0.6058818238834054e+01 )
- {}
+ c61 ( 0.8083246795921522e+01 ) , c62 ( -0.7981132988064893e+01 ) , c63 ( -0.3152159432874371e+02 ) , c64 ( 0.1631930543123136e+02 ) , c65 ( -0.6058818238834054e+01 ) ,
+ d21 ( 0.1012623508344586e+02 ) , d22 ( 0.7487995877610167e+01 ) , d23 ( -0.3480091861555747e+02 ) , d24 ( -0.7992771707568823e+01 ) , d25 ( 0.1025137723295662e+01 ) ,
+ d31 ( -0.6762803392801253e+00 ) , d32 ( 0.6087714651680015e+01 ) , d33 ( 0.1643084320892478e+02 ) , d34 ( 0.2476722511418386e+02 ) , d35 ( -0.6594389125716872e+01 )
+ {}
const value_type gamma;
const value_type d1 , d2 , d3 , d4;
@@ -73,6 +75,8 @@
const value_type c51 , c52 , c53 , c54;
const value_type a51 , a52 , a53 , a54;
const value_type c61 , c62 , c63 , c64 , c65;
+ const value_type d21 , d22 , d23 , d24 , d25;
+ const value_type d31 , d32 , d33 , d34 , d35;
};
@@ -95,8 +99,10 @@
m_state_adjuster.register_state( 4 , m_g3 );
m_state_adjuster.register_state( 5 , m_g4 );
m_state_adjuster.register_state( 6 , m_g5 );
- m_state_adjuster.register_state( 7 , m_xtmp );
- m_state_adjuster.register_state( 8 , m_dxdtnew );
+ m_state_adjuster.register_state( 7 , m_cont3 );
+ m_state_adjuster.register_state( 8 , m_cont4 );
+ m_state_adjuster.register_state( 9 , m_xtmp );
+ m_state_adjuster.register_state( 10 , m_dxdtnew );
}
void copy( const rosenbrock4 &rb )
@@ -110,6 +116,8 @@
m_g3 = rb.m_g3;
m_g4 = rb.m_g4;
m_g5 = rb.m_g5;
+ m_cont3 = rb.m_cont3;
+ m_cont4 = rb.m_cont4;
m_xtmp = rb.m_xtmp;
m_dxdtnew = rb.m_dxdtnew;
}
@@ -131,6 +139,7 @@
m_jac() , m_pm( 1 ) ,
m_dfdt() , m_dxdt() ,
m_g1() , m_g2() , m_g3() , m_g4() , m_g5() ,
+ m_cont3() , m_cont4() ,
m_xtmp() , m_dxdtnew() ,
m_coef()
{
@@ -142,6 +151,7 @@
m_jac() , m_pm( 1 ) ,
m_dfdt() , m_dxdt() ,
m_g1() , m_g2() , m_g3() , m_g4() , m_g5() ,
+ m_cont3() , m_cont4() ,
m_xtmp() , m_dxdtnew() ,
m_coef()
{
@@ -235,6 +245,30 @@
do_step( system , x , t , x , dt , xerr );
}
+ void prepare_dense_output()
+ {
+ const size_t n = m_g1.size();
+ for( size_t i=0 ; i<n ; ++i )
+ {
+ m_cont3[i] = m_coef.d21 * m_g1[i] + m_coef.d22 * m_g2[i] + m_coef.d23 * m_g3[i] + m_coef.d24 * m_g4[i] + m_coef.d25 * m_g5[i];
+ m_cont4[i] = m_coef.d31 * m_g1[i] + m_coef.d32 * m_g2[i] + m_coef.d33 * m_g3[i] + m_coef.d34 * m_g4[i] + m_coef.d35 * m_g5[i];
+ }
+ }
+
+
+ void calc_state( const time_type &t , state_type &x ,
+ const state_type &x_old , const time_type &t_old ,
+ const state_type &x_new , const time_type &t_new )
+ {
+ const size_t n = m_g1.size();
+ time_type dt = t_new - t_old;
+ time_type s = ( t - t_old ) / dt;
+ time_type s1 = 1.0 - s;
+ for( size_t i=0 ; i<n ; ++i )
+ x[i] = x_old[i] * s1 + s * ( x_new[i] + s1 * ( m_cont3[i] + s * m_cont4[i] ) );
+ }
+
+
template< class StateType >
@@ -248,7 +282,7 @@
private:
- size_adjuster< state_type , 9 > m_state_adjuster;
+ size_adjuster< state_type , 11 > m_state_adjuster;
size_adjuster< matrix_type , 1 , matrix_vector_adjust_size > m_matrix_adjuster;
size_adjuster< pmatrix_type , 1 > m_pmatrix_adjuster;
@@ -256,6 +290,7 @@
pmatrix_type m_pm;
state_type m_dfdt , m_dxdt;
state_type m_g1 , m_g2 , m_g3 , m_g4 , m_g5;
+ state_type m_cont3 , m_cont4;
state_type m_xtmp , m_dxdtnew;
rosenbrock_coefficients m_coef;
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp 2011-03-30 16:24:46 EDT (Wed, 30 Mar 2011)
@@ -32,7 +32,7 @@
rosenbrock4_controller( value_type atol = 1.0e-6 , value_type rtol = 1.0e-6 , const stepper_type &stepper = stepper_type() )
- : m_rb4() , m_atol( atol ) , m_rtol( rtol ) ,
+ : m_stepper() , m_atol( atol ) , m_rtol( rtol ) ,
m_first_step( true ) , m_err_old( 0.0 ) , m_dt_old( 0.0 ) ,
m_last_rejected( false )
{
@@ -56,16 +56,36 @@
}
+
+
+ /*
+ * ToDo : xout in size_adjuster
+ * Check
+ */
template< class System >
boost::numeric::odeint::controlled_step_result
try_step( System sys , state_type &x , value_type &t , value_type &dt )
{
+ state_type xout( x.size() );
+ boost::numeric::odeint::controlled_step_result res = try_step( sys , x , t , xout , dt );
+ x = xout;
+ return res;
+ }
+
+
+ /*
+ * ToDo : xerr in size_adjuster
+ */
+ template< class System >
+ boost::numeric::odeint::controlled_step_result
+ try_step( System sys , const state_type &x , value_type &t , state_type &xout , value_type &dt )
+ {
static const value_type safe = 0.9 , fac1 = 5.0 , fac2 = 1.0 / 6.0;
const size_t n = x.size();
- state_type xnew( n ) , xerr( n );
- m_rb4.do_step( sys , x , t , xnew , dt , xerr );
- value_type err = error( xnew , x , xerr );
+ state_type xerr( n );
+ m_stepper.do_step( sys , x , t , xout , dt , xerr );
+ value_type err = error( xout , x , xerr );
value_type fac = std::max( fac2 ,std::min( fac1 , std::pow( err , 0.25 ) / safe ) );
value_type dt_new = dt / fac;
@@ -90,7 +110,6 @@
t += dt;
dt = dt_new;
m_last_rejected = false;
- x = xnew;
return success_step_size_increased;
}
else
@@ -103,9 +122,23 @@
+
+ stepper_type& stepper( void )
+ {
+ return m_stepper;
+ }
+
+ const stepper_type& stepper( void ) const
+ {
+ return m_stepper;
+ }
+
+
+
+
private:
- stepper_type m_rb4;
+ stepper_type m_stepper;
value_type m_atol , m_rtol;
bool m_first_step;
value_type m_err_old , m_dt_old;
@@ -115,99 +148,6 @@
-//template< class Value >
-//class rosenbrock4_controller
-//{
-//public:
-//
-// typedef Value value_type;
-// typedef boost::numeric::ublas::vector< value_type > state_type;
-// typedef state_type deriv_type;
-// typedef value_type time_type;
-// typedef boost::numeric::ublas::matrix< value_type > matrix_type;
-// typedef boost::numeric::ublas::permutation_matrix< size_t > pmatrix_type;
-//
-// rosenbrock4_controller( value_type atol = 1.0e-6 , value_type rtol = 1.0e-6 )
-// : m_rb4() , m_atol( atol ) , m_rtol( rtol ) ,
-// m_first_step( true ) , m_err_old( 0.0 ) , m_dt_old( 0.0 ) ,
-// m_last_rejected( false )
-// {
-// }
-//
-// value_type error( const state_type &x , const state_type &xold , const state_type &xerr )
-// {
-// const size_t n = x.size();
-// value_type err = 0.0 , sk = 0.0;
-// for( size_t i=0 ; i<n ; ++i )
-// {
-// sk = m_atol + m_rtol * std::max( std::abs( xold[i] ) , std::abs( x[i] ) );
-// err += xerr[i] * xerr[i] / sk / sk;
-// }
-// return sqrt( err / value_type( n ) );
-// }
-//
-// value_type last_error( void ) const
-// {
-// return m_err_old;
-// }
-//
-//
-// template< class System >
-// boost::numeric::odeint::controlled_step_result
-// try_step( System sys , state_type &x , value_type &t , value_type &dt )
-// {
-// static const value_type safe = 0.9 , fac1 = 5.0 , fac2 = 1.0 / 6.0;
-//
-// const size_t n = x.size();
-// state_type xnew( n ) , xerr( n );
-// m_rb4.do_step( sys , x , t , xnew , dt , xerr );
-// value_type err = error( xnew , x , xerr );
-//
-// value_type fac = std::max( fac2 ,std::min( fac1 , std::pow( err , 0.25 ) / safe ) );
-// value_type dt_new = dt / fac;
-// if ( err <= 1.0 )
-// {
-// if( m_first_step )
-// {
-// m_first_step = false;
-// }
-// else
-// {
-// value_type fac_pred = ( m_dt_old / dt ) * pow( err * err / m_err_old , 0.25 ) / safe;
-// fac_pred = std::max( fac2 , std::min( fac1 , fac_pred ) );
-// fac = std::max( fac , fac_pred );
-// dt_new = dt / fac;
-// }
-//
-// m_dt_old = dt;
-// m_err_old = std::max( 0.01 , err );
-// if( m_last_rejected )
-// dt_new = ( dt >= 0.0 ? std::min( dt_new , dt ) : std::max( dt_new , dt ) );
-// t += dt;
-// dt = dt_new;
-// m_last_rejected = false;
-// x = xnew;
-// return success_step_size_increased;
-// }
-// else
-// {
-// dt = dt_new;
-// m_last_rejected = true;
-// return step_size_decreased;
-// }
-// }
-//
-//
-//
-//private:
-//
-// rosenbrock4< value_type > m_rb4;
-// value_type m_atol , m_rtol;
-// bool m_first_step;
-// value_type m_err_old , m_dt_old;
-// bool m_last_rejected;
-//};
-
} // namespace odeint
Added: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/rosenbrock4_dense_output.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/rosenbrock4_dense_output.hpp 2011-03-30 16:24:46 EDT (Wed, 30 Mar 2011)
@@ -0,0 +1,150 @@
+/*
+ * rosenbrock4_dense_output.hpp
+ *
+ * Created on: Mar 29, 2011
+ * Author: karsten
+ */
+
+#ifndef BOOST_NUMERIC_ODEINT_STEPPER_ROSENBROCK4_DENSE_OUTPUT_HPP_
+#define STEPPER_ROSENBROCK4_DENSE_OUTPUT_HPP_
+
+#include <utility>
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+template< class ControlledStepper >
+class rosenbrock4_dense_output
+{
+ void initialize_variables( void )
+ {
+ 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 );
+ }
+
+ void copy_variables( const rosenbrock4_dense_output &rb )
+ {
+ m_stepper = rb.m_stepper;
+ boost::numeric::copy( rb.m_x1 , m_x1 );
+ boost::numeric::copy( rb.m_x2 , m_x2 );
+ if( rb.m_current_state == ( & ( rb.m_x1 ) ) )
+ {
+ m_current_state = &m_x1;
+ m_old_state = &m_x2;
+ }
+ else
+ {
+ m_current_state = &m_x2;
+ m_old_state = &m_x1;
+ }
+ m_t = rb.m_t;
+ m_t_old = rb.m_t_old;
+ m_dt = rb.m_dt;
+ }
+
+public:
+
+ typedef ControlledStepper controlled_stepper_type;
+ typedef typename stepper_type::value_type value_type;
+ typedef typename stepper_type::state_type state_type;
+ typedef typename stepper_type::time_type time_type;
+ typedef typename stepper_type::deriv_type deriv_type;
+ typedef dense_output_stepper_tag stepper_category;
+
+ rosenbrock4_dense_output( const controlled_stepper_type &stepper = controlled_stepper_type() )
+ : m_stepper( stepper ) , m_state_adjuster() ,
+ m_x1() , m_x2() , m_current_state( &m_x1 ) , m_old_state( &m_x2 ) ,
+ m_t() , m_t_old() , m_dt()
+ {
+ initialize_variables();
+ }
+
+ rosenbrock4_dense_output( const rosenbrock4_dense_output &rb )
+ : m_stepper( stepper ) , m_state_adjuster() ,
+ m_x1() , m_x2() , m_current_state( &m_x1 ) , m_old_state( &m_x2 ) ,
+ m_t() , m_t_old() , m_dt()
+ {
+ initialize_variables();
+ copy_variables( rb );
+ }
+
+ rosenbrock4_dense_output& operator=( const rosenbrock4_dense_output &rb )
+ {
+ copy_variables( rb );
+ return *this;
+ }
+
+ ~rosenbrock4_dense_output( void )
+ {
+ boost::numeric::odeint::destruct( m_x1 );
+ boost::numeric::odeint::destruct( m_x2 );
+ }
+
+
+
+ template< class StateType >
+ void initialize( const StateType &x0 , const time_type &t0 , const time_type &dt0 )
+ {
+ m_state_adjuster.adjust_size_by_policy( x0 , adjust_size_policy() );
+ 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 = step_size_decreased;
+ 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_controlled_explicit_fsal : too much iterations!");
+ }
+ while( res == step_size_decreased );
+ m_stepper.stepper().prepare_dense_output();
+ std::swap( m_current_state , m_old_state );
+ return std::make_pair( m_t_old , m_t );
+ }
+
+
+ /*
+ * The two overloads are needed in order to solve the forwarding problem.
+ */
+ template< class StateOut >
+ void calc_state( const time_type &t , StateOut &x )
+ {
+ m_stepper.stepper().calc_state( t , x , *m_old_state , m_t_old , *m_current_state , m_t );
+ }
+
+ template< class StateOut >
+ void calc_state( const time_type &t , const StateOut &x )
+ {
+ m_stepper.stepper().calc_state( t , x , *m_old_state , m_t_old , *m_current_state , m_t );
+ }
+
+
+private:
+
+ controlled_stepper_type m_stepper;
+ size_adjuster< state_type , 2 > m_state_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
+} // namespace numeric
+} // namespace boost
+
+
+#endif /* STEPPER_ROSENBROCK4_DENSE_OUTPUT_HPP_ */
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