diff --git a/boost/numeric/odeint/stepper/base/explicit_stepper_and_error_stepper_fsal_base.hpp b/boost/numeric/odeint/stepper/base/explicit_stepper_and_error_stepper_fsal_base.hpp index abbada37..be2ff1db 100644 --- a/boost/numeric/odeint/stepper/base/explicit_stepper_and_error_stepper_fsal_base.hpp +++ b/boost/numeric/odeint/stepper/base/explicit_stepper_and_error_stepper_fsal_base.hpp @@ -95,32 +95,41 @@ public: template< class System > void do_step( System &system , state_type &x , const time_type t , const time_type dt ) { - m_size_adjuster.adjust_size_by_policy( x , adjust_size_policy() ); - system( x , m_dxdt ,t ); - this->stepper().do_step_impl( system , x , m_dxdt , t , x , dt ); + if( m_size_adjuster.adjust_size_by_policy( x , adjust_size_policy() ) || m_first_call ) + { + system( x , m_dxdt ,t ); + m_first_call = false; + } + this->stepper().do_step_impl( system , x , m_dxdt , t , x , m_dxdt , 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 ) + void do_step( System &system , state_type &x , state_type &dxdt , const time_type t , const time_type dt ) { - this->stepper().do_step_impl( system , x , dxdt , t , x , dt ); + m_first_call = true; + this->stepper().do_step_impl( system , x , dxdt , t , x , dxdt , 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 ) { - m_size_adjuster.adjust_size_by_policy( in , adjust_size_policy() ); - system( in , m_dxdt ,t ); - this->stepper().do_step_impl( system , in , m_dxdt , t , out , dt ); + if( m_size_adjuster.adjust_size_by_policy( in , adjust_size_policy() ) || m_first_call ) + { + system( in , m_dxdt ,t ); + m_first_call = false; + } + this->stepper().do_step_impl( system , in , m_dxdt , t , out , m_dxdt , dt ); } - // do_step( sys , in , dxdt , t , out , dt ) + // do_step( sys , in , dxdt_in , t , out , dxdt_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 ) + void do_step( 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 ) { - this->stepper().do_step_impl( system , in , dxdt , t , out , dt ); + m_first_call = true; + this->stepper().do_step_impl( system , in , dxdt_in , t , out , dxdt_out , dt ); } @@ -137,14 +146,15 @@ public: system( x , m_dxdt ,t ); m_first_call = false; } - this->stepper().do_step_impl( system , x , m_dxdt , t , x , dt , xerr ); + this->stepper().do_step_impl( system , x , m_dxdt , t , x , m_dxdt , 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 ); + m_first_call = true; + this->stepper().do_step_impl( system , x , dxdt , t , x , dxdt , dt , xerr ); } // do_step( sys , in , t , out , dt , xerr ) @@ -156,22 +166,20 @@ public: system( in , m_dxdt ,t ); m_first_call = false; } - this->stepper().do_step_impl( system , in , m_dxdt , t , out , dt , xerr ); + this->stepper().do_step_impl( system , in , m_dxdt , t , out , m_dxdt , dt , xerr ); } - // do_step( sys , in , dxdt , t , out , dt , xerr ) + // do_step( sys , in , dxdt_in , t , out , dxdt_out , dt ) 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 ) + void do_step( 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 ) { - this->stepper().do_step_impl( system , in , dxdt , t , out , dt , xerr ); + m_first_call = true; + this->stepper().do_step_impl( system , in , dxdt_in , t , out , dxdt_out , dt , xerr ); } - void reset_step( state_type &dxdt ) - { - this->stepper().reset_step_impl( dxdt ); - } diff --git a/boost/numeric/odeint/stepper/controlled_error_stepper.hpp b/boost/numeric/odeint/stepper/controlled_error_stepper.hpp index 6c1426f8..2830143f 100644 --- a/boost/numeric/odeint/stepper/controlled_error_stepper.hpp +++ b/boost/numeric/odeint/stepper/controlled_error_stepper.hpp @@ -238,16 +238,18 @@ public: 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_xnew_size_adjuster() , - m_dxdt() , m_xerr() , m_xnew() , + m_dxdt_size_adjuster() , m_xerr_size_adjuster() , m_new_size_adjuster() , + m_dxdt() , m_xerr() , m_xnew() , m_dxdtnew() , m_first_call( true ) { boost::numeric::odeint::construct( m_dxdt ); boost::numeric::odeint::construct( m_xerr ); boost::numeric::odeint::construct( m_xnew ); + boost::numeric::odeint::construct( m_dxdtnew ); m_dxdt_size_adjuster.register_state( 0 , m_dxdt ); m_xerr_size_adjuster.register_state( 0 , m_xerr ); - m_xnew_size_adjuster.register_state( 0 , m_xnew ); + m_new_size_adjuster.register_state( 0 , m_xnew ); + m_new_size_adjuster.register_state( 1 , m_dxdtnew ); } ~controlled_error_stepper( void ) @@ -255,6 +257,7 @@ public: boost::numeric::odeint::destruct( m_dxdt ); boost::numeric::odeint::destruct( m_xerr ); boost::numeric::odeint::destruct( m_xnew ); + boost::numeric::odeint::destruct( m_dxdtnew ); } @@ -273,19 +276,6 @@ public: 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 ) @@ -298,10 +288,24 @@ public: return try_step( sys , in , m_dxdt , t , out , 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_new_size_adjuster.adjust_size_by_policy( x , adjust_size_policy() ); + controlled_step_result res = try_step( sys , x , dxdt , t , m_xnew , m_dxdtnew , dt ); + if( ( res == success_step_size_increased ) || ( res == success_step_size_unchanged) ) + { + boost::numeric::odeint::copy( m_xnew , x ); + boost::numeric::odeint::copy( m_dxdtnew , dxdt ); + } + return res; + } // 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 ) + controlled_step_result try_step( System &sys , const state_type &in , const state_type &dxdt_in , time_type &t , + state_type &out , state_type &dxdt_out , time_type &dt ) { using std::max; using std::min; @@ -311,16 +315,15 @@ public: //fsal: m_stepper.get_dxdt( dxdt ); //fsal: m_stepper.do_step( sys , x , dxdt , t , dt , m_x_err ); - m_stepper.do_step( sys , in , dxdt , t , out , dt , m_xerr ); + m_stepper.do_step( sys , in , dxdt_in , t , out , dxdt_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( in , dxdt , m_xerr , dt ); + time_type max_rel_err = m_error_checker.error( in , dxdt_in , m_xerr , dt ); if( max_rel_err > 1.1 ) { // error too large - decrease dt ,limit scaling factor to 0.2 and reset state dt *= max( 0.9 * pow( max_rel_err , -1.0 / ( m_stepper.error_order() - 1.0 ) ) , 0.2 ); - m_stepper.reset_step( dxdt ); return step_size_decreased; } else @@ -375,11 +378,12 @@ private: 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; + size_adjuster< state_type , 2 > m_new_size_adjuster; state_type m_dxdt; state_type m_xerr; state_type m_xnew; + state_type m_dxdtnew; bool m_first_call; }; diff --git a/boost/numeric/odeint/stepper/dense_output_dopri5.hpp b/boost/numeric/odeint/stepper/dense_output_dopri5.hpp index d920e109..9e87dfdf 100644 --- a/boost/numeric/odeint/stepper/dense_output_dopri5.hpp +++ b/boost/numeric/odeint/stepper/dense_output_dopri5.hpp @@ -12,12 +12,12 @@ #ifndef BOOST_BOOST_NUMERIC_ODEINT_DENSE_OUTPUT_DOPRI5_HPP_INCLUDED #define BOOST_BOOST_NUMERIC_ODEINT_DENSE_OUTPUT_DOPRI5_HPP_INCLUDED -#include -#define tab "\t" -using std::cout; -using std::cerr; -using std::clog; -using std::endl; +//#include +//#define tab "\t" +//using std::cout; +//using std::cerr; +//using std::clog; +//using std::endl; #include @@ -52,25 +52,35 @@ public: 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 ) + m_x1() , m_x2() , m_dxdt1() , m_dxdt2() , + m_current_state( &m_x1 ) , m_old_state( &m_x2 ) , + m_current_deriv( &m_dxdt1 ) , m_old_deriv( &m_dxdt2 ) , + m_t( 0.0 ) , m_t_old( 0.0 ) , m_dt( 1.0 ) , + m_is_deriv_initialized( false ) { boost::numeric::odeint::construct( m_x1 ); boost::numeric::odeint::construct( m_x2 ); + boost::numeric::odeint::construct( m_dxdt1 ); + boost::numeric::odeint::construct( m_dxdt2 ); m_size_adjuster.register_state( 0 , m_x1 ); m_size_adjuster.register_state( 1 , m_x2 ); + m_size_adjuster.register_state( 1 , m_dxdt1); + m_size_adjuster.register_state( 1 , m_dxdt2 ); } ~dense_output_dopri5( void ) { boost::numeric::odeint::destruct( m_x1 ); boost::numeric::odeint::destruct( m_x2 ); + boost::numeric::odeint::destruct( m_dxdt1 ); + boost::numeric::odeint::destruct( m_dxdt2 ); } void adjust_size( const state_type &x ) { m_size_adjuster.adjust_size( x ); m_stepper.adjust_size( x ); + m_is_deriv_initialized = false; } void initialize( const state_type &x0 , const time_type t0 , const time_type dt0 ) @@ -78,23 +88,29 @@ public: boost::numeric::odeint::copy( x0 , *m_current_state ); m_t = t0; m_dt = dt0; + m_is_deriv_initialized = false; } template< class System > std::pair< time_type , time_type > do_step( System &system ) { const size_t max_count = 1000; - controlled_step_result res; + + if( !m_is_deriv_initialized ) + system( *m_current_state , *m_current_deriv , m_t ); + + 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 ); + res = m_stepper.try_step( system , *m_current_state , *m_current_deriv , m_t , *m_old_state , *m_old_deriv , 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 ); + std::swap( m_current_deriv , m_old_deriv ); return std::make_pair( m_t_old , m_t ); } @@ -122,27 +138,27 @@ public: void calc_state( time_type t , state_type &x ) { - const time_type b1 = static_cast ( 35.0 ) / static_cast( 384.0 ); - const time_type b3 = static_cast ( 500.0 ) / static_cast( 1113.0 ); - const time_type b4 = static_cast ( 125.0 ) / static_cast( 192.0 ); - const time_type b5 = static_cast ( -2187.0 ) / static_cast( 6784.0 ); - const time_type b6 = static_cast ( 11.0 ) / static_cast( 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; +// const time_type b1 = static_cast ( 35.0 ) / static_cast( 384.0 ); +// const time_type b3 = static_cast ( 500.0 ) / static_cast( 1113.0 ); +// const time_type b4 = static_cast ( 125.0 ) / static_cast( 192.0 ); +// const time_type b5 = static_cast ( -2187.0 ) / static_cast( 6784.0 ); +// const time_type b6 = static_cast ( 11.0 ) / static_cast( 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 * @@ -164,10 +180,12 @@ private: stepper_type &m_stepper; - size_adjuster< state_type , 2 > m_size_adjuster; - state_type m_x1 , m_x2; + size_adjuster< state_type , 4 > m_size_adjuster; + state_type m_x1 , m_x2 , m_dxdt1 , m_dxdt2; state_type *m_current_state , *m_old_state; + state_type *m_current_deriv , *m_old_deriv; time_type m_t , m_t_old , m_dt; + bool m_is_deriv_initialized; }; diff --git a/boost/numeric/odeint/stepper/explicit_error_dopri5.hpp b/boost/numeric/odeint/stepper/explicit_error_dopri5.hpp index ec990c54..005a043d 100644 --- a/boost/numeric/odeint/stepper/explicit_error_dopri5.hpp +++ b/boost/numeric/odeint/stepper/explicit_error_dopri5.hpp @@ -78,7 +78,8 @@ public : template< class System > - void do_step_impl( System system , const state_type &in , const state_type &dxdt , const time_type t , state_type &out , const time_type dt ) + 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 ) { m_size_adjuster.adjust_size_by_policy( in , adjust_size_policy() ); @@ -114,34 +115,39 @@ public : const time_type c6 = static_cast ( 11.0 ) / static_cast( 84.0 ); //m_x1 = x + dt*b21*dxdt - algebra_type::for_each3( m_x1 , in , dxdt , + algebra_type::for_each3( m_x1 , in , dxdt_in , typename operations_type::scale_sum2( 1.0 , dt*b21 ) ); system( m_x1 , m_x2 , t + dt*a2 ); // m_x1 = x + dt*b31*dxdt + dt*b32*m_x2 - algebra_type::for_each4( m_x1 , in , dxdt , m_x2 , + algebra_type::for_each4( m_x1 , in , dxdt_in , m_x2 , typename operations_type::scale_sum3( 1.0 , dt*b31 , dt*b32 )); system( m_x1 , m_x3 , t + dt*a3 ); // m_x1 = x + dt * (b41*dxdt + b42*m_x2 + b43*m_x3) - algebra_type::for_each5( m_x1 , in , dxdt , m_x2 , m_x3 , + algebra_type::for_each5( m_x1 , in , dxdt_in , m_x2 , m_x3 , typename operations_type::scale_sum4( 1.0 , dt*b41 , dt*b42 , dt*b43 )); system( m_x1, m_x4 , t + dt*a4 ); - algebra_type::for_each6( m_x1 , in , dxdt , m_x2 , m_x3 , m_x4 , + algebra_type::for_each6( m_x1 , in , dxdt_in , m_x2 , m_x3 , m_x4 , typename operations_type::scale_sum5( 1.0 , dt*b51 , dt*b52 , dt*b53 , dt*b54 )); system( m_x1 , m_x5 , t + dt*a5 ); - algebra_type::for_each7( m_x1 , in , dxdt , m_x2 , m_x3 , m_x4 , m_x5 , + algebra_type::for_each7( m_x1 , in , dxdt_in , m_x2 , m_x3 , m_x4 , m_x5 , typename operations_type::scale_sum6( 1.0 , dt*b61 , dt*b62 , dt*b63 , dt*b64 , dt*b65 )); system( m_x1 , m_x6 , t + dt ); - algebra_type::for_each7( out , in , dxdt , m_x3 , m_x4 , m_x5 , m_x6 , + algebra_type::for_each7( out , in , dxdt_in , m_x3 , m_x4 , m_x5 , m_x6 , typename operations_type::scale_sum6( 1.0 , dt*c1 , dt*c3 , dt*c4 , dt*c5 , dt*c6 )); + + // the new derivative + system( out , dxdt_out , t + dt ); } + 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 ) + 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 ( 35.0 ) / static_cast( 384.0 ); @@ -157,56 +163,16 @@ public : const time_type dc6 = c6 - static_cast ( 187.0 ) / static_cast( 2100.0 ); const time_type dc7 = static_cast( -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 ); + do_step_impl( system , in , dxdt_in , t , out , dxdt_out , dt ); //error estimate - algebra_type::for_each7( xerr , m_x1 , m_x3 , m_x4 , m_x5 , m_x6 , dxdt , + algebra_type::for_each7( xerr , dxdt_in , 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 ) ); } -// 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 ( 35.0 ) / static_cast( 384.0 ); -// const time_type c3 = static_cast ( 500.0 ) / static_cast( 1113.0 ); -// const time_type c4 = static_cast ( 125.0 ) / static_cast( 192.0 ); -// const time_type c5 = static_cast ( -2187.0 ) / static_cast( 6784.0 ); -// const time_type c6 = static_cast ( 11.0 ) / static_cast( 84.0 ); -// -// const time_type dc1 = c1 - static_cast ( 5179.0 ) / static_cast( 57600.0 ); -// const time_type dc3 = c3 - static_cast ( 7571.0 ) / static_cast( 16695.0 ); -// const time_type dc4 = c4 - static_cast ( 393.0 ) / static_cast( 640.0 ); -// const time_type dc5 = c5 - static_cast ( -92097.0 ) / static_cast( 339200.0 ); -// const time_type dc6 = c6 - static_cast ( 187.0 ) / static_cast( 2100.0 ); -// const time_type dc7 = static_cast( -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 ); - } - - void adjust_size( const state_type &x ) { m_size_adjuster.adjust_size( x ); diff --git a/libs/numeric/odeint/test/Jamfile b/libs/numeric/odeint/test/Jamfile index 92e9fcd4..52e958ff 100644 --- a/libs/numeric/odeint/test/Jamfile +++ b/libs/numeric/odeint/test/Jamfile @@ -6,18 +6,10 @@ import testing ; -#project -# : requirements -# ../../../.. -# $(BOOST_ROOT) -# BOOST_ALL_NO_LIB=1 -# : build-dir . -# ; project : requirements BOOST_ALL_NO_LIB=1 - /boost/test//boost_unit_test_framework ../../../.. $(BOOST_ROOT) ; @@ -31,9 +23,14 @@ unit-test basic_test check_dense_output_explicit_euler.cpp check_dense_output_dopri5.cpp : - : static + : /boost/test//boost_unit_test_framework + static ; - - - - + +exe controlled_stepper_evolution + : controlled_stepper_evolution.cpp + ; + +exe dense_output_stepper_evolution + : dense_output_stepper_evolution.cpp + ; \ No newline at end of file diff --git a/libs/numeric/odeint/test/controlled_stepper_evolution.cpp b/libs/numeric/odeint/test/controlled_stepper_evolution.cpp new file mode 100644 index 00000000..ae09c693 --- /dev/null +++ b/libs/numeric/odeint/test/controlled_stepper_evolution.cpp @@ -0,0 +1,77 @@ +/* + * controlled_stepper_evolution.cpp + * + * Created on: Nov 2, 2010 + * Author: karsten + */ + +#include +#include +#include +#include + +#include + +#define tab "\t" + +using namespace std; +using namespace boost::numeric::odeint; + + +typedef std::tr1::array< double , 2 > state_type; + +ostream& operator<<( ostream &out , const state_type &x ) +{ + out << x[0] << tab << x[1]; + return out; +} + +void sys( const state_type &x , state_type &dxdt , double t ) +{ + dxdt[0] = x[1]; + dxdt[1] = -x[0]; +} + +state_type sys_solution( double t , const state_type &x0 ) +{ + state_type sol = {{ 0.0 , 0.0 }}; + sol[0] = x0[0] * cos( t ) + x0[1] * sin( t ); + sol[1] = -x0[0] * sin( t ) + x0[1] * cos( t ); + return sol; +} + + + +template< class Stepper > +void evolution( Stepper &stepper , double t_end , const state_type &x0 , const string &filename ) +{ + ofstream fout( filename.c_str() ); + state_type x = x0; + double t = 0.0 , dt = 0.01; + while( t < t_end ) + { + state_type orig = sys_solution( t , x0 ); + state_type diff = {{ orig[0] - x[0] , orig[1] - x[1] }}; + double diff_abs = sqrt( diff[0] * diff[0] + diff[1] * diff[1] ); + fout << t << tab << orig << tab << x << tab << diff << tab << diff_abs << endl; + stepper.try_step( sys , x , t , dt ); + } +} + + +int main( int argc , char **argv ) +{ + typedef explicit_error_rk54_ck< state_type > rk54_type; + typedef explicit_error_dopri5< state_type > dopri5_type; + + rk54_type rk54; + dopri5_type dopri5; + controlled_error_stepper< rk54_type > rk54_controlled( rk54 ); + controlled_error_stepper< dopri5_type > dopri5_controlled( dopri5 ); + + state_type x0 = {{ 1.25 , 0.43 }}; + evolution( rk54_controlled , 100.0 , x0 , "rk54.dat" ); + evolution( dopri5_controlled , 100.0 , x0 , "dopri5.dat" ); + + return 0; +} diff --git a/libs/numeric/odeint/test/dense_output_stepper_evolution.cpp b/libs/numeric/odeint/test/dense_output_stepper_evolution.cpp new file mode 100644 index 00000000..7ec74c30 --- /dev/null +++ b/libs/numeric/odeint/test/dense_output_stepper_evolution.cpp @@ -0,0 +1,92 @@ +/* + * dense_output_stepper_evolution.cpp + * + * Created on: Nov 2, 2010 + * Author: karsten + */ + +#include +#include +#include +#include + +#include + +#define tab "\t" + +using namespace std; +using namespace boost::numeric::odeint; + +typedef std::tr1::array< double , 2 > state_type; + +ostream& operator<<( ostream &out , const state_type &x ) +{ + out << x[0] << tab << x[1]; + return out; +} + +void sys( const state_type &x , state_type &dxdt , double t ) +{ + dxdt[0] = x[1]; + dxdt[1] = -x[0]; +} + +state_type sys_solution( double t , const state_type &x0 ) +{ + state_type sol = {{ 0.0 , 0.0 }}; + sol[0] = x0[0] * cos( t ) + x0[1] * sin( t ); + sol[1] = -x0[0] * sin( t ) + x0[1] * cos( t ); + return sol; +} + + + +template< class Stepper > +void evolution( Stepper &stepper , double t_end , const state_type &x0 , const string &stepper_filename , const string &state_filename ) +{ + ofstream state_out( state_filename.c_str() ); + ofstream stepper_out( stepper_filename.c_str() ); + + stepper.initialize( x0 , 0.0 , 0.01 ); + + state_type x = x0; + double t = 0.0; + double dt = 0.001; + while( t < t_end ) + { + if( t < stepper.current_time() ) + { + stepper.calc_state( t , x ); + state_type orig = sys_solution( t , x0 ); + state_type diff = {{ orig[0] - x[0] , orig[1] - x[1] }}; + double diff_abs = sqrt( diff[0] * diff[0] + diff[1] * diff[1] ); + state_out << t << tab << x << tab << orig << tab << diff << diff_abs << endl; + } + else + { + stepper.do_step( sys ); + stepper_out << stepper.current_time() << "\t" << stepper.current_state() << std::endl; + continue; + } + t += dt; + + } +} + + +int main( int argc , char **argv ) +{ + typedef explicit_error_dopri5< state_type > dopri5_type; + typedef controlled_error_stepper< dopri5_type > controlled_dopri5_type; + + dopri5_type dopri5; + controlled_dopri5_type controlled_dopri5( dopri5 ); + + dense_output_explicit_euler< state_type > dense_euler; + dense_output_dopri5< controlled_dopri5_type > dense_dopri5( controlled_dopri5 ); + + state_type x0 = {{ 1.25 , 0.43 }}; + evolution( dense_euler , 100.0 , x0 , "dense_euler_stepper.dat" , "dense_euler_state.dat" ); + evolution( dense_dopri5 , 100.0 , x0 , "dense_dopri5_stepper.dat" , "dense_dopri5_state.dat" ); + +}