|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r66355 - 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-02 10:39:31
Author: karsten
Date: 2010-11-02 10:39:28 EDT (Tue, 02 Nov 2010)
New Revision: 66355
URL: http://svn.boost.org/trac/boost/changeset/66355
Log:
* interface changes for the fsal steppers
* some test routines for dense output steppers and controlled steppers
Added:
sandbox/odeint/branches/karsten/libs/numeric/odeint/test/controlled_stepper_evolution.cpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/test/dense_output_stepper_evolution.cpp (contents, props changed)
Text files modified:
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/explicit_stepper_and_error_stepper_fsal_base.hpp | 50 +++++++++++++----------
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/controlled_error_stepper.hpp | 46 +++++++++++----------
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/dense_output_dopri5.hpp | 84 ++++++++++++++++++++++++---------------
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_error_dopri5.hpp | 66 +++++++-----------------------
sandbox/odeint/branches/karsten/libs/numeric/odeint/test/Jamfile | 23 ++++------
5 files changed, 131 insertions(+), 138 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-11-02 10:39:28 EDT (Tue, 02 Nov 2010)
@@ -95,32 +95,41 @@
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 @@
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 @@
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 );
- }
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-02 10:39:28 EDT (Tue, 02 Nov 2010)
@@ -238,16 +238,18 @@
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 @@
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 @@
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 @@
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 @@
//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 @@
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;
};
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-02 10:39:28 EDT (Tue, 02 Nov 2010)
@@ -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 <iostream>
-#define tab "\t"
-using std::cout;
-using std::cerr;
-using std::clog;
-using std::endl;
+//#include <iostream>
+//#define tab "\t"
+//using std::cout;
+//using std::cerr;
+//using std::clog;
+//using std::endl;
#include <stdexcept>
@@ -52,25 +52,35 @@
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 @@
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 @@
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;
+// 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 *
@@ -164,10 +180,12 @@
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;
};
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-02 10:39:28 EDT (Tue, 02 Nov 2010)
@@ -78,7 +78,8 @@
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 @@
const time_type c6 = static_cast<time_type> ( 11.0 ) / static_cast<time_type>( 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<time_type> ( 35.0 ) / static_cast<time_type>( 384.0 );
@@ -157,56 +163,16 @@
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 );
+ 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<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 );
- }
-
-
void adjust_size( const state_type &x )
{
m_size_adjuster.adjust_size( x );
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-02 10:39:28 EDT (Tue, 02 Nov 2010)
@@ -6,18 +6,10 @@
import testing ;
-#project
-# : requirements
-# <include>../../../..
-# <include>$(BOOST_ROOT)
-# <define>BOOST_ALL_NO_LIB=1
-# : build-dir .
-# ;
project
: requirements
<define>BOOST_ALL_NO_LIB=1
- <library>/boost/test//boost_unit_test_framework
<include>../../../..
<include>$(BOOST_ROOT)
;
@@ -31,9 +23,14 @@
check_dense_output_explicit_euler.cpp
check_dense_output_dopri5.cpp
:
- : <link>static
+ : <library>/boost/test//boost_unit_test_framework
+ <link>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
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/test/controlled_stepper_evolution.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/test/controlled_stepper_evolution.cpp 2010-11-02 10:39:28 EDT (Tue, 02 Nov 2010)
@@ -0,0 +1,77 @@
+/*
+ * controlled_stepper_evolution.cpp
+ *
+ * Created on: Nov 2, 2010
+ * Author: karsten
+ */
+
+#include <string>
+#include <fstream>
+#include <iostream>
+#include <tr1/array>
+
+#include <boost/numeric/odeint.hpp>
+
+#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;
+}
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/test/dense_output_stepper_evolution.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/test/dense_output_stepper_evolution.cpp 2010-11-02 10:39:28 EDT (Tue, 02 Nov 2010)
@@ -0,0 +1,92 @@
+/*
+ * dense_output_stepper_evolution.cpp
+ *
+ * Created on: Nov 2, 2010
+ * Author: karsten
+ */
+
+#include <string>
+#include <fstream>
+#include <iostream>
+#include <tr1/array>
+
+#include <boost/numeric/odeint.hpp>
+
+#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" );
+
+}
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