boost/odeint
2
0
Fork 0
mirror of https://github.com/boostorg/odeint.git synced 2026-01-26 06:42:23 +00:00
Code Issues Projects Releases Wiki Activity

* interface changes for the fsal steppers

Browse Source 
* some test routines for dense output steppers and controlled steppers
This commit is contained in:
karsten
2010-11-02 14:39:28 +00:00
parent 0cf7ab86e5
commit 34575ca069
7 changed files with 300 additions and 138 deletions
Download Patch File Download Diff File Expand all files Collapse all files

50
boost/numeric/odeint/stepper/base/explicit_stepper_and_error_stepper_fsal_base.hpp
View File

@@ -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 );
}

46
boost/numeric/odeint/stepper/controlled_error_stepper.hpp
View File

@@ -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;
};

84
boost/numeric/odeint/stepper/dense_output_dopri5.hpp
View File

@@ -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 @@ 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<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 @@ 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;
};

66
boost/numeric/odeint/stepper/explicit_error_dopri5.hpp
View File

@@ -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<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 @@ public :
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 );

23
libs/numeric/odeint/test/Jamfile
View File

@@ -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 @@ unit-test basic_test
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
;

77
libs/numeric/odeint/test/controlled_stepper_evolution.cpp Normal file
View File

@@ -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;
}

92
libs/numeric/odeint/test/dense_output_stepper_evolution.cpp Normal file
View File

@@ -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" );
}