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

manual merge

Browse Source
This commit is contained in:
Karsten Ahnert
2012-07-09 10:10:28 +02:00
parent 3b0d56a05c 735071f11e
commit bcabbe61a5
181 changed files with 1095 additions and 18021 deletions
Download Patch File Download Diff File Expand all files Collapse all files

9
Jamroot
View File

@@ -19,7 +19,6 @@ project
# tests, regression tests and examples
build-project libs/numeric/odeint/test ;
build-project libs/numeric/odeint/examples ;
build-project libs/numeric/odeint/regression_test ;
# additional tests with external libraries :
@@ -28,14 +27,6 @@ build-project libs/numeric/odeint/regression_test ;
# build-project libs/numeric/odeint/test_external/gsl ;
# ideas
# build-project libs/numeric/odeint/ideas/butcher ;
# build-project libs/numeric/odeint/ideas/generic_stepper ;
#build-project libs/numeric/odeint/ideas/rosenbrock4 ;
#build-project libs/numeric/odeint/ideas/units ;
#build-project libs/numeric/odeint/ideas/algebra ;
# docs:
# build-project libs/numeric/odeint/doc ;

162
boost/numeric/odeint/external/mtl4/implicit_euler_mtl4.hpp vendored Normal file
View File

@@ -0,0 +1,162 @@
/*
[begin_description]
Modification of the implicit Euler method, works with the MTL4 matrix library only.
[end_description]
Copyright 2009-2011 Karsten Ahnert
Copyright 2009-2011 Mario Mulansky
Copyright 2012 Andreas Angelopoulos
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or
copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#ifndef BOOST_NUMERIC_ODEINT_EXTERNAL_IMPLICIT_EULER_MTL4_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_EXTERNAL_IMPLICIT_EULER_MTL4_HPP_INCLUDED
#include <utility>
#include <boost/numeric/odeint/util/bind.hpp>
#include <boost/numeric/odeint/util/unwrap_reference.hpp>
#include <boost/numeric/odeint/stepper/stepper_categories.hpp>
#include <boost/numeric/odeint/external/mtl4/mtl4_resize.hpp>
#include <boost/numeric/mtl/mtl.hpp>
#include <boost/numeric/itl/itl.hpp>
namespace boost {
namespace numeric {
namespace odeint {
template< class ValueType , class Resizer = initially_resizer >
class implicit_euler_mtl4
{
public:
typedef ValueType value_type;
typedef value_type time_type;
typedef mtl::dense_vector<value_type> state_type;
typedef state_wrapper< state_type > wrapped_state_type;
typedef state_type deriv_type;
typedef state_wrapper< deriv_type > wrapped_deriv_type;
typedef mtl::compressed2D< value_type > matrix_type;
typedef state_wrapper< matrix_type > wrapped_matrix_type;
typedef Resizer resizer_type;
typedef stepper_tag stepper_category;
typedef implicit_euler_mtl4< ValueType , Resizer > stepper_type;
implicit_euler_mtl4( const value_type epsilon = 1E-6 )
: m_epsilon( epsilon ) , m_resizer() ,
m_dxdt() , m_x() ,
m_identity() , m_jacobi()
{ }
template< class System >
void do_step( System system , state_type &x , time_type t , time_type dt )
{
typedef typename odeint::unwrap_reference< System >::type system_type;
typedef typename odeint::unwrap_reference< typename system_type::first_type >::type deriv_func_type;
typedef typename odeint::unwrap_reference< typename system_type::second_type >::type jacobi_func_type;
system_type &sys = system;
deriv_func_type &deriv_func = sys.first;
jacobi_func_type &jacobi_func = sys.second;
m_resizer.adjust_size( x , detail::bind(
&stepper_type::template resize_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
m_identity.m_v = 1;
t += dt;
m_x.m_v = x;
deriv_func( x , m_dxdt.m_v , t );
jacobi_func( x , m_jacobi.m_v , t );
m_dxdt.m_v *= -dt;
m_jacobi.m_v *= dt;
m_jacobi.m_v -= m_identity.m_v ;
// using ilu_0 preconditioning -incomplete LU factorisation
// itl::pc::diagonal<matrix_type,double> L(m_jacobi.m_v);
itl::pc::ilu_0<matrix_type> L( m_jacobi.m_v );
solve( m_jacobi.m_v , m_x.m_v , m_dxdt.m_v , L );
x+= m_x.m_v;
}
template< class StateType >
void adjust_size( const StateType &x )
{
resize_impl( x );
}
private:
/*
Applying approximate iteractive linear solvers
default solver is Biconjugate gradient stabilized method
itl::bicgstab(A, x, b, L, iter);
*/
template < class LinearOperator, class HilbertSpaceX, class HilbertSpaceB, class Preconditioner>
void solve(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
const Preconditioner& L, int max_iteractions =500)
{
// Termination criterion: r < 1e-6 * b or N iterations
itl::basic_iteration< double > iter( b , max_iteractions , 1e-6 );
itl::bicgstab( A , x , b , L , iter );
}
template< class StateIn >
bool resize_impl( const StateIn &x )
{
bool resized = false;
resized |= adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() );
resized |= adjust_size_by_resizeability( m_x , x , typename is_resizeable<state_type>::type() );
resized |= adjust_size_by_resizeability( m_identity , x , typename is_resizeable<matrix_type>::type() );
resized |= adjust_size_by_resizeability( m_jacobi , x , typename is_resizeable<matrix_type>::type() );
return resized;
}
private:
value_type m_epsilon;
resizer_type m_resizer;
wrapped_deriv_type m_dxdt;
wrapped_state_type m_x;
wrapped_matrix_type m_identity;
wrapped_matrix_type m_jacobi;
};
} // odeint
} // numeric
} // boost
#endif // BOOST_NUMERIC_ODEINT_EXTERNAL_IMPLICIT_EULER_MTL4_HPP_INCLUDED

133
boost/numeric/odeint/external/mtl4/mtl4_resize.hpp vendored Normal file
View File

@@ -0,0 +1,133 @@
/*
[begin_description]
Modification of the implicit Euler method, works with the MTL4 matrix library only.
[end_description]
Copyright 2009-2011 Karsten Ahnert
Copyright 2009-2011 Mario Mulansky
Copyright 2012 Andreas Angelopoulos
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or
copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#ifndef BOOST_NUMERIC_ODEINT_EXTERNAL_MTL4_RESIZE_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_EXTERNAL_MTL4_RESIZE_HPP_INCLUDED
#include <boost/numeric/odeint/util/is_resizeable.hpp>
#include <boost/numeric/odeint/util/resize.hpp>
#include <boost/numeric/odeint/util/same_size.hpp>
#include <boost/numeric/mtl/vector/dense_vector.hpp>
#include <boost/numeric/mtl/matrix/dense2D.hpp>
#include <boost/numeric/mtl/matrix/compressed2D.hpp>
namespace boost {
namespace numeric {
namespace odeint {
template< class Value , class Parameters >
struct is_resizeable< mtl::dense_vector< Value , Parameters > >
{
typedef boost::true_type type;
const static bool value = type::value;
};
template< class Value , class Parameters >
struct is_resizeable< mtl::dense2D< Value , Parameters > >
{
typedef boost::true_type type;
const static bool value = type::value;
};
template< class Value , class Parameters >
struct is_resizeable< mtl::compressed2D< Value , Parameters > >
{
typedef boost::true_type type;
const static bool value = type::value;
};
template< class Value , class Parameters >
struct same_size_impl< mtl::dense_vector< Value , Parameters > , mtl::dense_vector< Value , Parameters > >
{
static bool same_size( const mtl::dense_vector< Value , Parameters > &v1 ,
const mtl::dense_vector< Value , Parameters > &v2 )
{
return mtl::size( v1 ) == mtl::size( v2 );
}
};
template< class Value , class Parameters >
struct resize_impl< mtl::dense_vector< Value , Parameters > , mtl::dense_vector< Value , Parameters > >
{
static void resize( mtl::dense_vector< Value , Parameters > &v1 ,
const mtl::dense_vector< Value , Parameters > &v2 )
{
v1.change_dim( mtl::size( v2 ) );
}
};
template< class Value , class MatrixParameters , class VectorParameters >
struct same_size_impl< mtl::dense2D< Value , MatrixParameters > , mtl::dense_vector< Value , VectorParameters > >
{
static bool same_size( const mtl::dense2D< Value , MatrixParameters > &m ,
const mtl::dense_vector< Value , VectorParameters > &v )
{
return ( ( mtl::size( v ) == m.num_cols() ) && ( mtl::size( v ) == m.num_rows() ) );
}
};
template< class Value , class MatrixParameters , class VectorParameters >
struct resize_impl< mtl::dense2D< Value , MatrixParameters > , mtl::dense_vector< Value , VectorParameters > >
{
static void resize( mtl::dense2D< Value , MatrixParameters > &m ,
const mtl::dense_vector< Value , VectorParameters > &v )
{
m.change_dim( mtl::size( v ) , mtl::size( v ) , false );
}
};
template< class Value , class MatrixParameters , class VectorParameters >
struct same_size_impl< mtl::compressed2D< Value , MatrixParameters > , mtl::dense_vector< Value , VectorParameters > >
{
static bool same_size( const mtl::compressed2D< Value , MatrixParameters > &m ,
const mtl::dense_vector< Value , VectorParameters > &v )
{
return ( ( mtl::size( v ) == m.num_cols() ) && ( mtl::size( v ) == m.num_rows() ) );
}
};
template< class Value , class MatrixParameters , class VectorParameters >
struct resize_impl< mtl::compressed2D< Value , MatrixParameters > , mtl::dense_vector< Value , VectorParameters > >
{
static void resize( mtl::compressed2D< Value , MatrixParameters > &m ,
const mtl::dense_vector< Value , VectorParameters > &v )
{
m.change_dim( mtl::size( v ) , mtl::size( v ) );
}
};
} // namespace odeint
} // namesapce numeric
} // namespace boost
#endif // BOOST_NUMERIC_ODEINT_EXTERNAL_MTL4_RESIZE_HPP_INCLUDED

19
boost/numeric/odeint/external/thrust/thrust_operations.hpp vendored
View File

@@ -46,6 +46,25 @@ struct thrust_operations
}
};
template< class Fac1 = double , class Fac2 = Fac1 >
struct scale_sum_swap2
{
const Fac1 m_alpha1;
const Fac2 m_alpha2;
scale_sum_swap2( const Fac1 alpha1 , const Fac2 alpha2 )
: m_alpha1( alpha1 ) , m_alpha2( alpha2 ) { }
template< class Tuple >
__host__ __device__
void operator()( Tuple t ) const
{
typename thrust::tuple_element<0,Tuple>::type tmp = thrust::get<0>(t);
thrust::get<0>(t) = m_alpha1 * thrust::get<1>(t) + m_alpha2 * thrust::get<2>(t);
thrust::get<1>(t) = tmp;
}
};
template< class Fac1 = double , class Fac2 = Fac1 , class Fac3 = Fac2 >
struct scale_sum3
{

159
boost/numeric/odeint/integrate/detail/integrate_const.hpp Normal file
View File

@@ -0,0 +1,159 @@
/*
[auto_generated]
boost/numeric/odeint/integrate/detail/integrate_const.hpp
[begin_description]
integrate const implementation
[end_description]
Copyright 2009-2012 Karsten Ahnert
Copyright 2009-2012 Mario Mulansky
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or
copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#ifndef BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_CONST_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_CONST_HPP_INCLUDED
#include <boost/numeric/odeint/util/unwrap_reference.hpp>
#include <boost/numeric/odeint/stepper/stepper_categories.hpp>
#include <boost/numeric/odeint/util/unit_helper.hpp>
#include <boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp>
#include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
namespace boost {
namespace numeric {
namespace odeint {
namespace detail {
// forward declaration
template< class Stepper , class System , class State , class Time , class Observer >
size_t integrate_adaptive(
Stepper stepper , System system , State &start_state ,
Time &start_time , Time end_time , Time &dt ,
Observer observer , controlled_stepper_tag
);
template< class Stepper , class System , class State , class Time , class Observer >
size_t integrate_const(
Stepper stepper , System system , State &start_state ,
Time start_time , Time end_time , Time dt ,
Observer observer , stepper_tag
)
{
typename odeint::unwrap_reference< Observer >::type &obs = observer;
Time time = start_time;
int step = 0;
while( less_eq_with_sign( time+dt , end_time , dt ) )
{
obs( start_state , time );
stepper.do_step( system , start_state , time , dt );
// direct computation of the time avoids error propagation happening when using time += dt
// we need clumsy type analysis to get boost units working here
++step;
time = start_time + static_cast< typename unit_value_type<Time>::type >(step) * dt;
}
obs( start_state , time );
return step;
}
template< class Stepper , class System , class State , class Time , class Observer >
size_t integrate_const(
Stepper stepper , System system , State &start_state ,
Time start_time , Time end_time , Time dt ,
Observer observer , controlled_stepper_tag
)
{
typename odeint::unwrap_reference< Observer >::type &obs = observer;
Time time = start_time;
const Time time_step = dt;
int step = 0;
while( less_eq_with_sign( time+time_step , end_time , dt ) )
{
obs( start_state , time );
detail::integrate_adaptive( stepper , system , start_state , time , time+time_step , dt ,
null_observer() , controlled_stepper_tag() );
// direct computation of the time avoids error propagation happening when using time += dt
// we need clumsy type analysis to get boost units working here
++step;
time = start_time + static_cast< typename unit_value_type<Time>::type >(step) * time_step;
}
obs( start_state , time );
return step;
}
template< class Stepper , class System , class State , class Time , class Observer >
size_t integrate_const(
Stepper stepper , System system , State &start_state ,
Time start_time , Time end_time , Time dt ,
Observer observer , dense_output_stepper_tag
)
{
typename odeint::unwrap_reference< Observer >::type &obs = observer;
Time time = start_time;
stepper.initialize( start_state , time , dt );
obs( start_state , time );
time += dt;
int obs_step( 1 );
int real_step( 0 );
while( less_with_sign( time+dt , end_time , dt ) )
{
while( less_eq_with_sign( time , stepper.current_time() , dt ) )
{
stepper.calc_state( time , start_state );
obs( start_state , time );
++obs_step;
// direct computation of the time avoids error propagation happening when using time += dt
// we need clumsy type analysis to get boost units working here
time = start_time + static_cast< typename unit_value_type<Time>::type >(obs_step) * dt;
}
// we have not reached the end, do another real step
if( less_with_sign( stepper.current_time()+stepper.current_time_step() ,
end_time ,
stepper.current_time_step() ) )
{
while( less_eq_with_sign( stepper.current_time() , time , dt ) )
{
stepper.do_step( system );
++real_step;
}
}
else if( less_with_sign( stepper.current_time() , end_time , stepper.current_time_step() ) )
{ // do the last step ending exactly on the end point
stepper.initialize( stepper.current_state() , stepper.current_time() , end_time - stepper.current_time() );
stepper.do_step( system );
++real_step;
}
}
// last observation, if we are still in observation interval
if( less_eq_with_sign( time , end_time , dt ) )
{
stepper.calc_state( time , start_state );
obs( start_state , time );
}
return real_step;
}
} } } }
#endif

1
boost/numeric/odeint/integrate/detail/integrate_n_steps.hpp
View File

@@ -105,6 +105,7 @@ Time integrate_n_steps(
stepper.initialize( start_state , time , dt );
size_t step = 0;
while( step < num_of_steps )
{
while( less_with_sign( time , stepper.current_time() , stepper.current_time_step() ) )

17
boost/numeric/odeint/integrate/integrate_const.hpp
View File

@@ -23,7 +23,8 @@
#include <boost/numeric/odeint/stepper/stepper_categories.hpp>
#include <boost/numeric/odeint/integrate/null_observer.hpp>
#include <boost/numeric/odeint/integrate/integrate_n_steps.hpp>
#include <boost/numeric/odeint/integrate/detail/integrate_const.hpp>
#include <boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp>
namespace boost {
namespace numeric {
@@ -53,10 +54,10 @@ size_t integrate_const(
}
else
{
const size_t steps = static_cast<size_t>( (end_time-start_time)/dt );
integrate_n_steps( stepper , system , start_state , start_time , dt , steps , observer );
return steps;
}
return detail::integrate_const( stepper , system , start_state ,
start_time , end_time , dt ,
observer , typename Stepper::stepper_category() );
}
}
@@ -77,9 +78,9 @@ size_t integrate_const(
}
else
{
const size_t steps = static_cast<size_t>( (end_time-start_time)/dt );
integrate_n_steps( stepper , system , start_state , start_time , dt , steps , observer );
return steps;
return detail::integrate_const( stepper , system , start_state ,
start_time , end_time , dt ,
observer , typename Stepper::stepper_category() );
}
}

67
boost/numeric/odeint/stepper/bulirsch_stoer.hpp
View File

@@ -38,13 +38,19 @@
#include <boost/numeric/odeint/util/state_wrapper.hpp>
#include <boost/numeric/odeint/util/is_resizeable.hpp>
#include <boost/numeric/odeint/util/resizer.hpp>
#include <boost/numeric/odeint/util/unit_helper.hpp>
namespace boost {
namespace numeric {
namespace odeint {
/** ToDo try_step stepsize changed return values doesn't make too much sense here as we have order control as well */
template<
class State ,
class Value = double ,
@@ -71,7 +77,11 @@ public:
typedef bulirsch_stoer< State , Value , Deriv , Time , Algebra , Operations , Resizer > controlled_error_bs_type;
typedef std::vector< time_type > value_vector;
typedef typename inverse_time< time_type >::type inv_time_type;
typedef std::vector< value_type > value_vector;
typedef std::vector< time_type > time_vector;
typedef std::vector< inv_time_type > inv_time_vector; //should be 1/time_type for boost.units
typedef std::vector< value_vector > value_matrix;
typedef std::vector< size_t > int_vector;
typedef std::vector< wrapped_state_type > state_table_type;
@@ -80,21 +90,22 @@ public:
bulirsch_stoer(
time_type eps_abs = 1E-6 , time_type eps_rel = 1E-6 ,
time_type factor_x = 1.0 , time_type factor_dxdt = 1.0 )
value_type eps_abs = 1E-6 , value_type eps_rel = 1E-6 ,
value_type factor_x = 1.0 , value_type factor_dxdt = 1.0 )
: m_error_checker( eps_abs , eps_rel , factor_x, factor_dxdt ) , m_midpoint() ,
m_last_step_rejected( false ) , m_first( true ) ,
m_dt_last( 1.0E30 ) , m_t_last() ,
/* , m_t_last() ,
m_current_k_opt() ,
m_algebra() ,
m_dxdt_resizer() , m_xnew_resizer() , m_resizer() ,
m_xnew() , m_err() , m_dxdt() ,
m_xnew() , m_err() , m_dxdt() ,*/
m_interval_sequence( m_k_max+1 ) ,
m_coeff( m_k_max+1 ) ,
m_cost( m_k_max+1 ) ,
m_table( m_k_max ) ,
STEPFAC1( 0.65 ) , STEPFAC2( 0.94 ) , STEPFAC3( 0.02 ) , STEPFAC4( 4.0 ) , KFAC1( 0.8 ) , KFAC2( 0.9 )
{
//m_dt_last = 1.0E30;
for( unsigned short i = 0; i < m_k_max+1; i++ )
{
m_interval_sequence[i] = 2 * (i+1);
@@ -105,13 +116,13 @@ public:
m_coeff[i].resize(i);
for( size_t k = 0 ; k < i ; ++k )
{
const time_type r = static_cast< time_type >( m_interval_sequence[i] ) / static_cast< time_type >( m_interval_sequence[k] );
m_coeff[i][k] = 1.0 / ( r*r - static_cast< time_type >( 1.0 ) ); // coefficients for extrapolation
const value_type r = static_cast< value_type >( m_interval_sequence[i] ) / static_cast< value_type >( m_interval_sequence[k] );
m_coeff[i][k] = 1.0 / ( r*r - static_cast< value_type >( 1.0 ) ); // coefficients for extrapolation
//std::cout << i << "," << k << " " << m_coeff[i][k] << '\t' ;
}
//std ::cout << std::endl;
// crude estimate of optimal order
const time_type logfact( -log10( std::max( eps_rel , 1.0E-12 ) ) * 0.6 + 0.5 );
const value_type logfact( -log10( std::max( eps_rel , 1.0E-12 ) ) * 0.6 + 0.5 );
m_current_k_opt = std::max( 1 , std::min( static_cast<int>( m_k_max-1 ) , static_cast<int>( logfact ) ));
//m_current_k_opt = m_k_max - 1;
//std::cout << m_cost[i] << std::endl;
@@ -172,7 +183,7 @@ public:
template< class System , class StateIn , class DerivIn , class StateOut >
controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , time_type &dt )
{
static const time_type val1( static_cast< time_type >( 1.0 ) );
static const value_type val1( 1.0 );
typename odeint::unwrap_reference< System >::type &sys = system;
if( m_resizer.adjust_size( in , detail::bind( &controlled_error_bs_type::template resize_impl< StateIn > , detail::ref( *this ) , detail::_1 ) ) )
@@ -187,8 +198,8 @@ public:
bool reject( true );
value_vector h_opt( m_k_max+1 );
value_vector work( m_k_max+1 );
time_vector h_opt( m_k_max+1 );
inv_time_vector work( m_k_max+1 );
//std::cout << "t=" << t <<", dt=" << dt << "(" << m_dt_last << ")" << ", k_opt=" << m_current_k_opt << std::endl;
@@ -208,10 +219,10 @@ public:
extrapolate( k , m_table , m_coeff , out );
// get error estimate
m_algebra.for_each3( m_err.m_v , out , m_table[0].m_v ,
typename operations_type::template scale_sum2< time_type , time_type >( val1 , -val1 ) );
const time_type error = m_error_checker.error( m_algebra , in , dxdt , m_err.m_v , dt );
typename operations_type::template scale_sum2< value_type , value_type >( val1 , -val1 ) );
const value_type error = m_error_checker.error( m_algebra , in , dxdt , m_err.m_v , dt );
h_opt[k] = calc_h_opt( dt , error , k );
work[k] = m_cost[k]/h_opt[k];
work[k] = static_cast<value_type>( m_cost[k] ) / h_opt[k];
//std::cout << '\t' << "h_opt=" << h_opt[k] << ", work=" << work[k] << std::endl;
//std::cout << '\t' << "error: " << error << std::endl;
@@ -225,7 +236,8 @@ public:
{
// leave order as is (except we were in first round)
m_current_k_opt = std::min( static_cast<int>(m_k_max)-1 , std::max( 2 , static_cast<int>(k)+1 ) );
new_h = h_opt[k] * m_cost[k+1]/m_cost[k];
new_h = h_opt[k];
new_h *= static_cast<value_type>( m_cost[k+1] ) / static_cast<value_type>( m_cost[k] );
} else {
m_current_k_opt = std::min( static_cast<int>(m_k_max)-1 , std::max( 2 , static_cast<int>(k) ) );
new_h = h_opt[k];
@@ -253,7 +265,8 @@ public:
else if( (work[k] < KFAC2*work[k-1]) && !m_last_step_rejected )
{
m_current_k_opt = std::min( static_cast<int>(m_k_max-1) , static_cast<int>(m_current_k_opt)+1 );
new_h = h_opt[k]*m_cost[m_current_k_opt]/m_cost[k];
new_h = h_opt[k];
new_h *= m_cost[m_current_k_opt]/m_cost[k];
//std::cout << new_h << std::endl;
} else
new_h = h_opt[m_current_k_opt];
@@ -369,23 +382,23 @@ private:
//polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf
{
//std::cout << "extrapolate k=" << k << ":" << std::endl;
static const time_type val1 = static_cast< time_type >( 1.0 );
static const value_type val1 = static_cast< value_type >( 1.0 );
for( int j=k-1 ; j>0 ; --j )
{
//std::cout << '\t' << m_coeff[k][j];
m_algebra.for_each3( table[j-1].m_v , table[j].m_v , table[j-1].m_v ,
typename operations_type::template scale_sum2< time_type , time_type >( val1 + coeff[k][j] , -coeff[k][j] ) );
typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k][j] , -coeff[k][j] ) );
}
//std::cout << std::endl << m_coeff[k][0] << std::endl;
m_algebra.for_each3( xest , table[0].m_v , xest ,
typename operations_type::template scale_sum2< time_type , time_type >( val1 + coeff[k][0] , -coeff[k][0]) );
typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k][0] , -coeff[k][0]) );
}
time_type calc_h_opt( time_type h , value_type error , size_t k ) const
{
time_type expo=1.0/(2*k+1);
time_type facmin = std::pow( STEPFAC3 , expo );
time_type fac;
value_type expo=1.0/(2*k+1);
value_type facmin = std::pow( STEPFAC3 , expo );
value_type fac;
if (error == 0.0)
fac=1.0/facmin;
else
@@ -397,7 +410,7 @@ private:
return h*fac;
}
controlled_step_result set_k_opt( size_t k , const value_vector &work , const value_vector &h_opt , time_type &dt )
controlled_step_result set_k_opt( size_t k , const inv_time_vector &work , const time_vector &h_opt , time_type &dt )
{
//std::cout << "finding k_opt..." << std::endl;
if( k == 1 )
@@ -433,18 +446,18 @@ private:
return ( (k == m_current_k_opt) || (k == m_current_k_opt+1) );
}
bool should_reject( time_type error , size_t k ) const
bool should_reject( value_type error , size_t k ) const
{
if( (k == m_current_k_opt-1) )
{
const time_type d = m_interval_sequence[m_current_k_opt] * m_interval_sequence[m_current_k_opt+1] /
const value_type d = m_interval_sequence[m_current_k_opt] * m_interval_sequence[m_current_k_opt+1] /
(m_interval_sequence[0]*m_interval_sequence[0]);
//step will fail, criterion 17.3.17 in NR
return ( error > d*d );
}
else if( k == m_current_k_opt )
{
const time_type d = m_interval_sequence[m_current_k_opt] / m_interval_sequence[0];
const value_type d = m_interval_sequence[m_current_k_opt] / m_interval_sequence[0];
return ( error > d*d );
} else
return error > 1.0;
@@ -477,7 +490,7 @@ private:
state_table_type m_table; // sequence of states for extrapolation
const time_type STEPFAC1 , STEPFAC2 , STEPFAC3 , STEPFAC4 , KFAC1 , KFAC2;
const value_type STEPFAC1 , STEPFAC2 , STEPFAC3 , STEPFAC4 , KFAC1 , KFAC2;
};

92
boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp
View File

@@ -36,6 +36,7 @@
#include <boost/numeric/odeint/util/state_wrapper.hpp>
#include <boost/numeric/odeint/util/is_resizeable.hpp>
#include <boost/numeric/odeint/util/resizer.hpp>
#include <boost/numeric/odeint/util/unit_helper.hpp>
#include <boost/type_traits.hpp>
@@ -87,7 +88,11 @@ public:
typedef bulirsch_stoer_dense_out< State , Value , Deriv , Time , Algebra , Operations , Resizer > controlled_error_bs_type;
typedef std::vector< time_type > value_vector;
typedef typename inverse_time< time_type >::type inv_time_type;
typedef std::vector< value_type > value_vector;
typedef std::vector< time_type > time_vector;
typedef std::vector< inv_time_type > inv_time_vector; //should be 1/time_type for boost.units
typedef std::vector< value_vector > value_matrix;
typedef std::vector< size_t > int_vector;
typedef std::vector< wrapped_state_type > state_vector_type;
@@ -100,18 +105,12 @@ public:
bulirsch_stoer_dense_out(
time_type eps_abs = 1E-6 , time_type eps_rel = 1E-6 ,
time_type factor_x = 1.0 , time_type factor_dxdt = 1.0 ,
value_type eps_abs = 1E-6 , value_type eps_rel = 1E-6 ,
value_type factor_x = 1.0 , value_type factor_dxdt = 1.0 ,
bool control_interpolation = false )
: m_error_checker( eps_abs , eps_rel , factor_x, factor_dxdt ) , m_midpoint() ,
: m_error_checker( eps_abs , eps_rel , factor_x, factor_dxdt ) ,
m_control_interpolation( control_interpolation) ,
m_last_step_rejected( false ) , m_first( true ) ,
m_t() , m_dt() ,
m_dt_last( 1.0E30 ) , m_t_last() ,
m_current_k_opt() ,
m_k_final(0) ,
m_algebra() , m_resizer() ,
m_x1() , m_x2() , m_dxdt1() , m_dxdt2() , m_err() ,
m_current_state_x1( true ) ,
m_error( m_k_max ) ,
m_interval_sequence( m_k_max+1 ) ,
@@ -134,13 +133,13 @@ public:
m_coeff[i].resize(i);
for( size_t k = 0 ; k < i ; ++k )
{
const time_type r = static_cast< time_type >( m_interval_sequence[i] ) / static_cast< time_type >( m_interval_sequence[k] );
m_coeff[i][k] = 1.0 / ( r*r - static_cast< time_type >( 1.0 ) ); // coefficients for extrapolation
const value_type r = static_cast< value_type >( m_interval_sequence[i] ) / static_cast< value_type >( m_interval_sequence[k] );
m_coeff[i][k] = 1.0 / ( r*r - static_cast< value_type >( 1.0 ) ); // coefficients for extrapolation
//std::cout << i << "," << k << " " << m_coeff[i][k] << '\t' ;
}
//std ::cout << std::endl;
// crude estimate of optimal order
const time_type logfact( -log10( std::max( eps_rel , 1.0E-12 ) ) * 0.6 + 0.5 );
const value_type logfact( -log10( std::max( eps_rel , 1.0E-12 ) ) * 0.6 + 0.5 );
m_current_k_opt = std::max( 1 , std::min( static_cast<int>( m_k_max-1 ) , static_cast<int>( logfact ) ));
//m_current_k_opt = m_k_max - 1;
//std::cout << m_cost[i] << std::endl;
@@ -173,7 +172,7 @@ public:
template< class System , class StateIn , class DerivIn , class StateOut , class DerivOut >
controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , DerivOut &dxdt_new , time_type &dt )
{
static const time_type val1( static_cast< time_type >( 1.0 ) );
static const value_type val1( 1.0 );
typename odeint::unwrap_reference< System >::type &sys = system;
// if( m_resizer.adjust_size( in , detail::bind( &controlled_error_bs_type::template resize_impl< StateIn > , detail::ref( *this ) , detail::_1 ) ) )
@@ -183,8 +182,8 @@ public:
bool reject( true );
value_vector h_opt( m_k_max+1 );
value_vector work( m_k_max+1 );
time_vector h_opt( m_k_max+1 );
inv_time_vector work( m_k_max+1 );
m_k_final = 0;
time_type new_h = dt;
@@ -205,10 +204,10 @@ public:
extrapolate( k , m_table , m_coeff , out );
// get error estimate
m_algebra.for_each3( m_err.m_v , out , m_table[0].m_v ,
typename operations_type::template scale_sum2< time_type , time_type >( val1 , -val1 ) );
const time_type error = m_error_checker.error( m_algebra , in , dxdt , m_err.m_v , dt );
typename operations_type::template scale_sum2< value_type , value_type >( val1 , -val1 ) );
const value_type error = m_error_checker.error( m_algebra , in , dxdt , m_err.m_v , dt );
h_opt[k] = calc_h_opt( dt , error , k );
work[k] = m_cost[k]/h_opt[k];
work[k] = static_cast<value_type>( m_cost[k] ) / h_opt[k];
//std::cout << '\t' << "h_opt=" << h_opt[k] << ", work=" << work[k] << std::endl;
//std::cout << '\t' << "error: " << error << std::endl;
@@ -224,7 +223,7 @@ public:
{
// leave order as is (except we were in first round)
m_current_k_opt = std::min( static_cast<int>(m_k_max)-1 , std::max( 2 , static_cast<int>(k)+1 ) );
new_h = h_opt[k] * m_cost[k+1]/m_cost[k];
new_h = h_opt[k] * static_cast<value_type>( m_cost[k+1] ) / static_cast<value_type>( m_cost[k] );
} else {
m_current_k_opt = std::min( static_cast<int>(m_k_max)-1 , std::max( 2 , static_cast<int>(k) ) );
new_h = h_opt[k];
@@ -252,7 +251,7 @@ public:
else if( (work[k] < KFAC2*work[k-1]) && !m_last_step_rejected )
{
m_current_k_opt = std::min( static_cast<int>(m_k_max)-1 , static_cast<int>(m_current_k_opt)+1 );
new_h = h_opt[k]*m_cost[m_current_k_opt]/m_cost[k];
new_h = h_opt[k]*static_cast<value_type>( m_cost[m_current_k_opt] ) / static_cast<value_type>( m_cost[k] );
} else
new_h = h_opt[m_current_k_opt];
break;
@@ -412,15 +411,15 @@ private:
void extrapolate( size_t k , StateVector &table , const value_matrix &coeff , StateInOut &xest , size_t order_start_index = 0 )
//polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf
{
static const time_type val1 = static_cast< time_type >( 1.0 );
static const value_type val1( 1.0 );
for( int j=k-1 ; j>0 ; --j )
{
m_algebra.for_each3( table[j-1].m_v , table[j].m_v , table[j-1].m_v ,
typename operations_type::template scale_sum2< time_type , time_type >( val1 + coeff[k + order_start_index][j + order_start_index] ,
typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k + order_start_index][j + order_start_index] ,
-coeff[k + order_start_index][j + order_start_index] ) );
}
m_algebra.for_each3( xest , table[0].m_v , xest ,
typename operations_type::template scale_sum2< time_type , time_type >( val1 + coeff[k + order_start_index][0 + order_start_index] ,
typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k + order_start_index][0 + order_start_index] ,
-coeff[k + order_start_index][0 + order_start_index]) );
}
@@ -431,25 +430,25 @@ private:
{
// result is written into table[0]
//std::cout << "extrapolate k=" << k << ":" << std::endl;
static const time_type val1 = static_cast< time_type >( 1.0 );
static const value_type val1( 1.0 );
for( int j=k ; j>1 ; --j )
{
//std::cout << '\t' << coeff[k + order_start_index][j + order_start_index - 1];
m_algebra.for_each3( table[j-1].m_v , table[j].m_v , table[j-1].m_v ,
typename operations_type::template scale_sum2< time_type , time_type >( val1 + coeff[k + order_start_index][j + order_start_index - 1] ,
typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k + order_start_index][j + order_start_index - 1] ,
-coeff[k + order_start_index][j + order_start_index - 1] ) );
}
//std::cout << std::endl << coeff[k + order_start_index][order_start_index] << std::endl;
m_algebra.for_each3( table[0].m_v , table[1].m_v , table[0].m_v ,
typename operations_type::template scale_sum2< time_type , time_type >( val1 + coeff[k + order_start_index][order_start_index] ,
typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k + order_start_index][order_start_index] ,
-coeff[k + order_start_index][order_start_index]) );
}
time_type calc_h_opt( time_type h , value_type error , size_t k ) const
{
time_type expo=1.0/(m_interval_sequence[k-1]);
time_type facmin = std::pow( STEPFAC3 , expo );
time_type fac;
value_type expo=1.0/(m_interval_sequence[k-1]);
value_type facmin = std::pow( STEPFAC3 , expo );
value_type fac;
if (error == 0.0)
fac=1.0/facmin;
else
@@ -468,18 +467,18 @@ private:
return ( (k == m_current_k_opt) || (k == m_current_k_opt+1) );
}
bool should_reject( time_type error , size_t k ) const
bool should_reject( value_type error , size_t k ) const
{
if( (k == m_current_k_opt-1) )
{
const time_type d = m_interval_sequence[m_current_k_opt] * m_interval_sequence[m_current_k_opt+1] /
const value_type d = m_interval_sequence[m_current_k_opt] * m_interval_sequence[m_current_k_opt+1] /
(m_interval_sequence[0]*m_interval_sequence[0]);
//step will fail, criterion 17.3.17 in NR
return ( error > d*d );
}
else if( k == m_current_k_opt )
{
const time_type d = m_interval_sequence[m_current_k_opt+1] / m_interval_sequence[0];
const value_type d = m_interval_sequence[m_current_k_opt+1] / m_interval_sequence[0];
return ( error > d*d );
} else
return error > 1.0;
@@ -500,8 +499,9 @@ private:
// calculate finite difference approximations to derivatives at the midpoint
for( int j = 0 ; j<=k ; j++ )
{
const time_type d = m_interval_sequence[j] / static_cast<time_type>(2*dt);
time_type f = 1.0; //factor 1/2 here because our interpolation interval has length 2 !!!
/* not working with boost units... */
const value_type d = m_interval_sequence[j] / ( static_cast<value_type>(2) * dt );
value_type f = 1.0; //factor 1/2 here because our interpolation interval has length 2 !!!
for( int kappa = 0 ; kappa <= 2*j+1 ; ++kappa )
{
calculate_finite_difference( j , kappa , f , dxdt_start );
@@ -527,7 +527,7 @@ private:
//std::cout << "extrapolation result: " << m_diffs[kappa][0].m_v << std::endl;
// divide kappa-th derivative by kappa because we need these terms for dense output interpolation
m_algebra.for_each1( m_diffs[kappa][0].m_v , typename operations_type::template scale< value_type >( static_cast<value_type>(d) ) );
m_algebra.for_each1( m_diffs[kappa][0].m_v , typename operations_type::template scale< time_type >( static_cast<time_type>(d) ) );
d *= dt/(2*(kappa+2));
@@ -608,13 +608,13 @@ private:
}
template< class DerivIn >
void calculate_finite_difference( size_t j , size_t kappa , time_type fac , const DerivIn &dxdt )
void calculate_finite_difference( size_t j , size_t kappa , value_type fac , const DerivIn &dxdt )
{
const int m = m_interval_sequence[j]/2-1;
if( kappa == 0) // no calculation required for 0th derivative of f
{
m_algebra.for_each2( m_diffs[0][j].m_v , m_derivs[j][m].m_v ,
typename operations_type::template scale_sum1< time_type >( fac ) );
typename operations_type::template scale_sum1< value_type >( fac ) );
//std::cout << "j=" << j << ", kappa=" << kappa << ", m=" << m;
//std::cout << ": m_diffs[" << kappa << "][" << j << "] = " << fac << " * f[" << m << "] ";
//std::cout << "(size(f)=" << m_derivs[j].size() << ") = " << m_diffs[0][j].m_v << std::endl;
@@ -630,8 +630,8 @@ private:
//std::cout << "j=" << j << ", kappa=" << kappa << ", m=" << m << ": m_diffs[" << kappa << "][" << j_diffs << "] = " << fac << " ( 1*f[" << m+kappa << "]";
m_algebra.for_each2( m_diffs[kappa][j_diffs].m_v , m_derivs[j][m+kappa].m_v ,
typename operations_type::template scale_sum1< time_type >( fac ) );
time_type sign = -1.0;
typename operations_type::template scale_sum1< value_type >( fac ) );
value_type sign = -1.0;
int c = 1;
//computes the j-th order finite difference for the kappa-th derivative of f at t+dt/2 using function evaluations stored in m_derivs
for( int i = m+static_cast<int>(kappa)-2 ; i >= m-static_cast<int>(kappa) ; i -= 2 )
@@ -639,7 +639,7 @@ private:
if( i >= 0 )
{
m_algebra.for_each3( m_diffs[kappa][j_diffs].m_v , m_diffs[kappa][j_diffs].m_v , m_derivs[j][i].m_v ,
typename operations_type::template scale_sum2< time_type , time_type >( 1.0 ,
typename operations_type::template scale_sum2< value_type , value_type >( 1.0 ,
sign * fac * boost::math::binomial_coefficient< double >( kappa , c ) ) );
//std::cout << ( (sign > 0.0) ? " + " : " - " ) <<
// boost::math::binomial_coefficient< double >( kappa , c ) << "*f[" << i << "]";
@@ -647,7 +647,7 @@ private:
else
{
m_algebra.for_each3( m_diffs[kappa][j_diffs].m_v , m_diffs[kappa][j_diffs].m_v , dxdt ,
typename operations_type::template scale_sum2< time_type , time_type >( 1.0 , sign * fac ) );
typename operations_type::template scale_sum2< value_type , value_type >( 1.0 , sign * fac ) );
//std::cout << ( (sign > 0.0) ? " + " : " - " ) << "dxdt";
}
sign *= -1;
@@ -662,7 +662,7 @@ private:
{
// interpolation polynomial is defined for theta = -1 ... 1
// m_k_final is the number of order-iterations done for the last step - it governs the order of the interpolation polynomial
const time_type theta = 2*(t - m_t_last)/(m_t - m_t_last) - 1;
const value_type theta = 2 * get_unit_value( (t - m_t_last) / (m_t - m_t_last) ) - 1;
//std::cout << "theta=" << theta << std::endl;
//start with x = a0 + a_{2k+1} theta^{2k+1} + a_{2k+2} theta^{2k+2} + a_{2k+3} theta^{2k+3} + a_{2k+4} theta^{2k+4}
//std::cout << "x = a_0 + ";
@@ -681,12 +681,12 @@ private:
boost::numeric::odeint::copy( m_mp_states[0].m_v , out );
// add remaining terms: x += a_1 theta + a2 theta^2 + ... + a_{2k} theta^{2k}
time_type theta_pow( theta );
value_type theta_pow( theta );
for( size_t i=0 ; i<=2*m_k_final+1 ; ++i )
{
//std::cout << "a_" << i+1 << " theta^" << i+1 << " = " << m_diffs[i][0].m_v[0] * std::pow( theta , i+1 ) << std::endl;
m_algebra.for_each3( out , out , m_diffs[i][0].m_v ,
typename operations_type::template scale_sum2< time_type >( static_cast<time_type>(1) , theta_pow ) );
typename operations_type::template scale_sum2< value_type >( static_cast<value_type>(1) , theta_pow ) );
theta_pow *= theta;
}
}
@@ -807,7 +807,7 @@ private:
//wrapped_state_type m_a1 , m_a2 , m_a3 , m_a4;
const time_type STEPFAC1 , STEPFAC2 , STEPFAC3 , STEPFAC4 , KFAC1 , KFAC2;
const value_type STEPFAC1 , STEPFAC2 , STEPFAC3 , STEPFAC4 , KFAC1 , KFAC2;
};
}

4
boost/numeric/odeint/stepper/controlled_runge_kutta.hpp
View File

@@ -475,7 +475,7 @@ public:
if( max_rel_err > 1.0 )
{
// error too large - decrease dt ,limit scaling factor to 0.2 and reset state
dt *= max( static_cast<value_type>(9)/static_cast<value_type>(10) * pow( max_rel_err , static_cast<value_type>(-1) / ( m_stepper.error_order() - 1 ) ) , static_cast<value_type>(1)/static_cast<value_type> (5) );
dt *= max( static_cast<value_type>( static_cast<value_type>(9)/static_cast<value_type>(10) * pow( max_rel_err , static_cast<value_type>(-1) / ( m_stepper.error_order() - 1 ) ) ) , static_cast<value_type>( static_cast<value_type>(1)/static_cast<value_type> (5)) );
return fail;
}
else
@@ -483,7 +483,7 @@ public:
if( max_rel_err < 0.5 )
{ //error too small - increase dt and keep the evolution and limit scaling factor to 5.0
t += dt;
dt *= min( static_cast<value_type>(9)/static_cast<value_type>(10) * pow( max_rel_err , static_cast<value_type>(-1) / m_stepper.stepper_order() ) , static_cast<value_type>(5) );
dt *= min( static_cast<value_type>( static_cast<value_type>(9)/static_cast<value_type>(10) * pow( max_rel_err , static_cast<value_type>(-1) / m_stepper.stepper_order() ) ) , static_cast<value_type>(5) );
return success;
}
else

8
boost/numeric/odeint/stepper/modified_midpoint.hpp
View File

@@ -76,8 +76,8 @@ public :
m_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_impl< StateIn > , detail::ref( *this ) , detail::_1 ) );
const time_type h = dt / static_cast<time_type>( m_steps );
const time_type h2 = static_cast<time_type>( 2 ) * h;
const time_type h = dt / static_cast<value_type>( m_steps );
const time_type h2 = static_cast<value_type>(2) * h;
typename odeint::unwrap_reference< System >::type &sys = system;
@@ -200,8 +200,8 @@ public :
m_resizer.adjust_size( in , detail::bind( &stepper_type::template resize< StateIn > , detail::ref( *this ) , detail::_1 ) );
const time_type h = dt / static_cast<time_type>( m_steps );
const time_type h2 = static_cast<time_type>( 2 ) * h;
const time_type h = dt / static_cast<value_type>( m_steps );
const time_type h2 = static_cast<value_type>( 2 ) * h;
typename odeint::unwrap_reference< System >::type &sys = system;

28
boost/numeric/odeint/util/unit_helper.hpp
View File

@@ -21,6 +21,8 @@
#ifndef __CUDACC__
#include <boost/units/quantity.hpp>
#include <boost/units/get_dimension.hpp>
#include <boost/units/get_system.hpp>
#endif
@@ -115,6 +117,32 @@ namespace detail {
template< typename Time >
struct inverse_time
{
typedef Time type;
};
template< typename Unit , typename Value >
struct inverse_time< boost::units::quantity< Unit , Value > >
{
typedef boost::units::quantity< Unit , Value > time_type;
typedef typename boost::units::get_dimension< time_type >::type dimension;
typedef typename boost::units::get_system< time_type >::type system;
typedef typename boost::mpl::divides< boost::units::dimensionless_type , dimension >::type inv_dimension;
typedef boost::units::unit< inv_dimension , system > inv_unit;
typedef boost::units::quantity< inv_unit , Value > type;
};
} // namespace odeint
} // namespace numeric
} // namespace boost

9
doc/boostbook.css
View File

@@ -73,7 +73,14 @@
}
h1 { font-size: 140%; }
h2 { font-weight: bold; font-size: 140%; }
h2 {
/* changed by Mario Mulansky to show logo instead of title */
height: 127px;
background: white url(logo.jpg) left no-repeat;
text-indent: -999px;
/* end change */
font-weight: bold; font-size: 140%;
}
h3 { font-weight: bold; font-size: 130%; }
h4 { font-weight: bold; font-size: 120%; }
h5 { font-weight: normal; font-style: italic; font-size: 110%; }

2
doc/index.html
View File

@@ -24,7 +24,7 @@
<div><p class="copyright">Copyright &#169; 2009-2011 Karsten Ahnert
and Mario Mulansky</p></div>
<div><div class="legalnotice">
<a name="odeint.legal"></a><p>
<a name="id533358"></a><p>
Distributed under the Boost Software License, Version 1.0. (See accompanying
file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>)
</p>

48
libs/numeric/odeint/code_analyzer/check_includes.py
View File

@@ -1,48 +0,0 @@
#! /usr/bin/python
import os
import glob
import subprocess
topleveldir = "../../../"
filename = []
filename += glob.glob( topleveldir + "boost/numeric/odeint/*.hpp" )
filename += glob.glob( topleveldir + "boost/numeric/odeint/algebra/*.hpp" )
filename += glob.glob( topleveldir + "boost/numeric/odeint/algebra/detail/*.hpp" )
filename += glob.glob( topleveldir + "boost/numeric/odeint/algebra/external/*.hpp" )
filename += glob.glob( topleveldir + "boost/numeric/odeint/stepper/*.hpp" )
filename += glob.glob( topleveldir + "boost/numeric/odeint/stepper/base/*.hpp" )
filename += glob.glob( topleveldir + "boost/numeric/odeint/stepper/detail/*.hpp" )
test_prog = '#include <INCLUDE>\n\
\n\
int main( void )\n\
{\n\
return 0;\n\
}\n\
'
count = 0
for fn in filename :
fn = fn[len(topleveldir) : ]
prog = test_prog.replace( "INCLUDE" , fn )
file = open( "test.cc" , "w" )
file.write( prog )
file.close
cmd = ["g++" , "-I"+topleveldir , "-I$BOOST_ROOT" , "-c" , "test.cc" ]
print "Test " + fn + " : "
# print prog
print "Commandline : " + subprocess.list2cmdline( cmd )
p = subprocess.Popen( subprocess.list2cmdline( cmd ) , shell=True)
sts = os.waitpid(p.pid, 0)[1]
print sts
print ""

15
libs/numeric/odeint/code_analyzer/glob_dir.py
View File

@@ -1,15 +0,0 @@
import fnmatch
import os
def glob_dir( dir , pattern , exclude = [ ] ):
matches = []
for root , dirnames, filenames in os.walk( dir ):
exclude_dir = False
for e in exclude:
if root.find( e ) != -1:
exclude_dir = True
if exclude_dir == False:
for filename in fnmatch.filter( filenames , pattern ):
matches.append( [ root , filename , os.path.join(root, filename) ] )
return matches

48
libs/numeric/odeint/code_analyzer/header_preamble.py
View File

@@ -1,48 +0,0 @@
#! /usr/bin/python
from glob_dir import *
preamble_template = "/*\n\
[auto_generated]\n\
FILENAME\n\
\n\
[begin_description]\n\
\n\
[end_description]\n\
\n\
Copyright 2009-2011 Karsten Ahnert\n\
Copyright 2009-2011 Mario Mulansky\n\
\n\
Distributed under the Boost Software License, Version 1.0.\n\
(See accompanying file LICENSE_1_0.txt or\n\
copy at http://www.boost.org/LICENSE_1_0.txt)\n\
*/\n\
\n\
"
def header_guard( file_desc ):
return file_desc[2].upper().replace( "/" , "_" ).replace( "." , "_" ) + "_INCLUDED"
def insert_preamble( content , file_desc ):
preamble = preamble_template.replace( "FILENAME" , file_desc[2] )
ret = preamble + "\n"
ret += "#ifndef " + header_guard( file_desc ) + "\n"
ret += "#define " + header_guard( file_desc ) + "\n\n"
ret += content
ret += "#endif // " + header_guard( file_desc ) + "\n"
return ret
files = glob_dir( "../../../../boost" , "*.hpp" , [ "algebra" , "external" , "integrate" , "util" ] )
for f in files:
print str( f[2] ) #+ " " + header_guard( f )
file = open( f[2] )
content = file.read()
content = insert_preamble( content , f )
file.close()
file = open( f[2] , "w" )
file.write( content )
file.close()

9
libs/numeric/odeint/doc/html/boostbook.css
View File

@@ -73,7 +73,14 @@
}
h1 { font-size: 140%; }
h2 { font-weight: bold; font-size: 140%; }
h2 {
/* changed by Mario Mulansky to show logo instead of title */
height: 127px;
background: white url(logo.jpg) left no-repeat;
text-indent: -999px;
/* end change */
font-weight: bold; font-size: 140%;
}
h3 { font-weight: bold; font-size: 130%; }
h4 { font-weight: bold; font-size: 120%; }
h5 { font-weight: normal; font-style: italic; font-size: 110%; }

10
libs/numeric/odeint/examples/generation_functions.cpp
View File

@@ -11,11 +11,11 @@
copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#include <array>
#include <boost/array.hpp>
#include <boost/numeric/odeint.hpp>
typedef std::array< double , 1 > state_type;
typedef boost::array< double , 1 > state_type;
using namespace boost::numeric::odeint;
@@ -73,10 +73,12 @@ int main( int argc , char **argv )
{
typedef runge_kutta_dopri5< state_type > stepper_type;
/*
//[ generation_functions_syntax_auto
auto stepper1 = make_controlled( 1.0e-6 , 1.0e-6 , stepper_type() );
auto stepper2 = make_dense_output( 1.0e-6 , 1.0e-6 , stepper_type() );
//]
*/
//[ generation_functions_syntax_result_of
result_of::make_controlled< stepper_type >::type stepper3 = make_controlled( 1.0e-6 , 1.0e-6 , stepper_type() );
@@ -85,9 +87,13 @@ int main( int argc , char **argv )
}
{
/*
//[ generation_functions_example_custom_controller
auto stepper5 = make_controlled( 1.0e-6 , 1.0e-6 , custom_stepper() );
//]
*/
result_of::make_controlled< custom_stepper >::type stepper5 = make_controlled( 1.0e-6 , 1.0e-6 , custom_stepper() );
}
return 0;
}

21
libs/numeric/odeint/examples/heun.cpp
View File

@@ -15,7 +15,6 @@
*/
#include <iostream>
#include <array>
#include <boost/fusion/container/vector.hpp>
@@ -136,6 +135,19 @@ struct lorenz
}
};
struct streaming_observer
{
std::ostream &m_out;
streaming_observer( std::ostream &out ) : m_out( out ) { }
template< typename State , typename Value >
void operator()( const State &x , Value t ) const
{
m_out << t;
for( size_t i=0 ; i<x.size() ; ++i ) m_out << "\t" << x[i];
m_out << "\n";
}
};
int main( int argc , char **argv )
@@ -145,14 +157,13 @@ int main( int argc , char **argv )
//[ heun_example
typedef array< double , 3 > state_type;
typedef boost::array< double , 3 > state_type;
heun< state_type > h;
state_type x = {{ 10.0 , 10.0 , 10.0 }};
integrate_const( h , lorenz() , x , 0.0 , 100.0 , 0.01 ,
[]( const state_type &_x , double _t ) {
std::cout << _t << "\t" << _x[0] << "\t" << _x[1] << "\t" << _x[2] << "\n";
});
streaming_observer( std::cout ) );
//]
return 0;

5
libs/numeric/odeint/examples/mtl/Jamfile
View File

@@ -3,13 +3,16 @@
# accompanying file LICENSE_1_0.txt or copy at
# http://www.boost.org/LICENSE_1_0.txt)
MTL4_INCLUDE = /home/karsten/boost/MTL-4.0.8862-Linux/usr/include ;
project
: requirements
<include>../../../../..
<include>$(MTL4_INCLUDE)
<define>BOOST_ALL_NO_LIB=1
: build-dir .
: default-build release
;
exe gauss_packet : gauss_packet.cpp ;
exe gauss_packet : gauss_packet.cpp ;
exe implicit_euler_mtl : implicit_euler_mtl.cpp ;

4
libs/numeric/odeint/examples/mtl/gauss_packet.cpp
View File

@@ -16,9 +16,9 @@
#include <complex>
#include <boost/numeric/odeint.hpp>
#include "mtl_bindings.hpp"
#include <boost/numeric/odeint/external/mtl4/mtl4_resize.hpp>
#include <boost/numeric/mtl/mtl.hpp>
using namespace std;

314
libs/numeric/odeint/examples/mtl/implicit_euler_mtl.cpp Normal file
View File

@@ -0,0 +1,314 @@
#include <iostream>
#include <fstream>
#include <utility>
#include "time.h"
#include <boost/numeric/odeint.hpp>
#include <boost/phoenix/phoenix.hpp>
#include <boost/numeric/mtl/mtl.hpp>
#include <boost/numeric/odeint/external/mtl4/implicit_euler_mtl4.hpp>
using namespace std;
using namespace boost::numeric::odeint;
namespace phoenix = boost::phoenix;
typedef mtl::dense_vector< double > vec_mtl4;
typedef mtl::compressed2D< double > mat_mtl4;
typedef boost::numeric::ublas::vector< double > vec_ublas;
typedef boost::numeric::ublas::matrix< double > mat_ublas;
// two systems defined 1 & 2 both are mostly sparse with the number of element variable
struct system1_mtl4
{
void operator()( const vec_mtl4 &x , vec_mtl4 &dxdt , double t )
{
int size = mtl::size(x);
dxdt[ 0 ] = -0.06*x[0];
for (int i =1; i< size ; ++i){
dxdt[ i ] = 4.2*x[i-1]-2.2*x[i]*x[i];
}
}
};
struct jacobi1_mtl4
{
void operator()( const vec_mtl4 &x , mat_mtl4 &J , const double &t )
{
int size = mtl::size(x);
mtl::matrix::inserter<mat_mtl4> ins(J);
ins[0][0]=-0.06;
for (int i =1; i< size ; ++i)
{
ins[i][i-1] = + 4.2;
ins[i][i] = -4.2*x[i];
}
}
};
struct system1_ublas
{
void operator()( const vec_ublas &x , vec_ublas &dxdt , double t )
{
int size = x.size();
dxdt[ 0 ] = -0.06*x[0];
for (int i =1; i< size ; ++i){
dxdt[ i ] = 4.2*x[i-1]-2.2*x[i]*x[i];
}
}
};
struct jacobi1_ublas
{
void operator()( const vec_ublas &x , mat_ublas &J , const double &t )
{
int size = x.size();
// mtl::matrix::inserter<mat_mtl4> ins(J);
J(0,0)=-0.06;
for (int i =1; i< size ; ++i){
//ins[i][0]=120.0*x[i];
J(i,i-1) = + 4.2;
J(i,i) = -4.2*x[i];
}
}
};
struct system2_mtl4
{
void operator()( const vec_mtl4 &x , vec_mtl4 &dxdt , double t )
{
int size = mtl::size(x);
for (int i =0; i< size/5 ; i+=5){
dxdt[ i ] = -0.5*x[i];
dxdt[i+1]= +25*x[i+1]*x[i+2]-740*x[i+3]*x[i+3]+4.2e-2*x[i];
dxdt[i+2]= +25*x[i]*x[i]-740*x[i+3]*x[i+3];
dxdt[i+3]= -25*x[i+1]*x[i+2]+740*x[i+3]*x[i+3];
dxdt[i+4] = 0.250*x[i]*x[i+1]-44.5*x[i+3];
}
}
};
struct jacobi2_mtl4
{
void operator()( const vec_mtl4 &x , mat_mtl4 &J , const double &t )
{
int size = mtl::size(x);
mtl::matrix::inserter<mat_mtl4> ins(J);
for (int i =0; i< size/5 ; i+=5){
ins[ i ][i] = -0.5;
ins[i+1][i+1]=25*x[i+2];
ins[i+1][i+2] = 25*x[i+1];
ins[i+1][i+3] = -740*2*x[i+3];
ins[i+1][i] =+4.2e-2;
ins[i+2][i]= 50*x[i];
ins[i+2][i+3]= -740*2*x[i+3];
ins[i+3][i+1] = -25*x[i+2];
ins[i+3][i+2] = -25*x[i+1];
ins[i+3][i+3] = +740*2*x[i+3];
ins[i+4][i] = 0.25*x[i+1];
ins[i+4][i+1] =0.25*x[i];
ins[i+4][i+3]=-44.5;
}
}
};
struct system2_ublas
{
void operator()( const vec_ublas &x , vec_ublas &dxdt , double t )
{
int size = x.size();
for (int i =0; i< size/5 ; i+=5){
dxdt[ i ] = -4.2e-2*x[i];
dxdt[i+1]= +25*x[i+1]*x[i+2]-740*x[i+3]*x[i+3]+4.2e-2*x[i];
dxdt[i+2]= +25*x[i]*x[i]-740*x[i+3]*x[i+3];
dxdt[i+3]= -25*x[i+1]*x[i+2]+740*x[i+3]*x[i+3];
dxdt[i+4] = 0.250*x[i]*x[i+1]-44.5*x[i+3];
}
}
};
struct jacobi2_ublas
{
void operator()( const vec_ublas &x , mat_ublas &J , const double &t )
{
int size = x.size();
for (int i =0; i< size/5 ; i+=5){
J(i ,i) = -4.2e-2;
J(i+1,i+1)=25*x[i+2];
J(i+1,i+2) = 25*x[i+1];
J(i+1,i+3) = -740*2*x[i+3];
J(i+1,i) =+4.2e-2;
J(i+2,i)= 50*x[i];
J(i+2,i+3)= -740*2*x[i+3];
J(i+3,i+1) = -25*x[i+2];
J(i+3,i+2) = -25*x[i+1];
J(i+3,i+3) = +740*2*x[i+3];
J(i+4,i) = 0.25*x[i+1];
J(i+4,i+1) =0.25*x[i];
J(i+4,i+3)=-44.5;
}
}
};
void testRidiculouslyMassiveArray( int size )
{
typedef boost::numeric::odeint::implicit_euler_mtl4 < double > mtl4stepper;
typedef boost::numeric::odeint::implicit_euler< double > booststepper;
vec_mtl4 x(size , 0.0);
x[0]=1;
double dt = 0.02;
double endtime = 10.0;
clock_t tstart_mtl4 = clock();
size_t num_of_steps_mtl4 = integrate_const(
mtl4stepper() ,
make_pair( system1_mtl4() , jacobi1_mtl4() ) ,
x , 0.0 , endtime , dt );
clock_t tend_mtl4 = clock() ;
clog << x[0] << endl;
clog << num_of_steps_mtl4 << " time elapsed: " << (double)(tend_mtl4-tstart_mtl4 )/CLOCKS_PER_SEC << endl;
vec_ublas x_ublas(size , 0.0);
x_ublas[0]=1;
clock_t tstart_boost = clock();
size_t num_of_steps_ublas = integrate_const(
booststepper() ,
make_pair( system1_ublas() , jacobi1_ublas() ) ,
x_ublas , 0.0 , endtime , dt );
clock_t tend_boost = clock() ;
clog << x_ublas[0] << endl;
clog << num_of_steps_ublas << " time elapsed: " << (double)(tend_boost-tstart_boost)/CLOCKS_PER_SEC<< endl;
clog << "dt_ublas/dt_mtl4 = " << (double)( tend_boost-tstart_boost )/( tend_mtl4-tstart_mtl4 ) << endl << endl;
return ;
}
void testRidiculouslyMassiveArray2( int size )
{
typedef boost::numeric::odeint::implicit_euler_mtl4 < double > mtl4stepper;
typedef boost::numeric::odeint::implicit_euler< double > booststepper;
vec_mtl4 x(size , 0.0);
x[0]=100;
double dt = 0.01;
double endtime = 10.0;
clock_t tstart_mtl4 = clock();
size_t num_of_steps_mtl4 = integrate_const(
mtl4stepper() ,
make_pair( system1_mtl4() , jacobi1_mtl4() ) ,
x , 0.0 , endtime , dt );
clock_t tend_mtl4 = clock() ;
clog << x[0] << endl;
clog << num_of_steps_mtl4 << " time elapsed: " << (double)(tend_mtl4-tstart_mtl4 )/CLOCKS_PER_SEC << endl;
vec_ublas x_ublas(size , 0.0);
x_ublas[0]=100;
clock_t tstart_boost = clock();
size_t num_of_steps_ublas = integrate_const(
booststepper() ,
make_pair( system1_ublas() , jacobi1_ublas() ) ,
x_ublas , 0.0 , endtime , dt );
clock_t tend_boost = clock() ;
clog << x_ublas[0] << endl;
clog << num_of_steps_ublas << " time elapsed: " << (double)(tend_boost-tstart_boost)/CLOCKS_PER_SEC<< endl;
clog << "dt_ublas/dt_mtl4 = " << (double)( tend_boost-tstart_boost )/( tend_mtl4-tstart_mtl4 ) << endl << endl;
return ;
}
int main( int argc , char **argv )
{
std::vector< size_t > length;
length.push_back( 8 );
length.push_back( 16 );
length.push_back( 32 );
length.push_back( 64 );
length.push_back( 128 );
length.push_back( 256 );
for( size_t i=0 ; i<length.size() ; ++i )
{
clog << "Testing with size " << length[i] << endl;
testRidiculouslyMassiveArray( length[i] );
}
clog << endl << endl;
for( size_t i=0 ; i<length.size() ; ++i )
{
clog << "Testing with size " << length[i] << endl;
testRidiculouslyMassiveArray2( length[i] );
}
}

49
libs/numeric/odeint/examples/mtl/mtl_bindings.hpp
View File

@@ -1,49 +0,0 @@
/*
* mtl_bindings.hpp
*
* Created on: Oct 26, 2011
* Author: mario
*/
#ifndef MTL_BINDINGS_HPP_
#define MTL_BINDINGS_HPP_
//[ mtl_bindings
#include <boost/numeric/mtl/mtl.hpp>
#include <boost/numeric/odeint/util/state_wrapper.hpp>
namespace boost { namespace numeric { namespace odeint {
template< class Value , class Parameters >
struct is_resizeable< mtl::dense_vector< Value , Parameters > >
{ // declare resizeablility
typedef boost::true_type type;
const static bool value = type::value;
};
template< class Value , class Parameters > //, class Value2 , class Parameters2 >
struct same_size_impl< mtl::dense_vector< Value , Parameters > , mtl::dense_vector< Value , Parameters > >
{ // define how to check size
static bool same_size( const mtl::dense_vector< Value , Parameters > &v1 ,
const mtl::dense_vector< Value , Parameters > &v2 )
{
return mtl::size( v1 ) == mtl::size( v2 );
}
};
template< class Value , class Parameters >
struct resize_impl< mtl::dense_vector< Value , Parameters > , mtl::dense_vector< Value , Parameters > >
{ // define how to resize
static void resize( mtl::dense_vector< Value , Parameters > &v1 ,
const mtl::dense_vector< Value , Parameters > &v2 )
{
v1.change_dim( mtl::size( v2 ) );
}
};
} } }
//]
#endif /* MTL_BINDINGS_HPP_ */

18
libs/numeric/odeint/examples/stochastic_euler.cpp
View File

@@ -14,8 +14,8 @@
#include <vector>
#include <iostream>
#include <random>
#include <array>
#include <boost/random.hpp>
#include <boost/array.hpp>
#include <boost/numeric/odeint.hpp>
@@ -26,8 +26,8 @@ template< size_t N > class stochastic_euler
{
public:
typedef std::array< double , N > state_type;
typedef std::array< double , N > deriv_type;
typedef boost::array< double , N > state_type;
typedef boost::array< double , N > deriv_type;
typedef double value_type;
typedef double time_type;
typedef unsigned short order_type;
@@ -71,8 +71,8 @@ class stochastic_euler
{
public:
typedef std::array< double , N > state_type;
typedef std::array< double , N > deriv_type;
typedef boost::array< double , N > state_type;
typedef boost::array< double , N > deriv_type;
typedef double value_type;
typedef double time_type;
typedef unsigned short order_type;
@@ -97,7 +97,7 @@ public:
//[ stochastic_euler_ornstein_uhlenbeck_def
const static size_t N = 1;
typedef std::array< double , N > state_type;
typedef boost::array< double , N > state_type;
struct ornstein_det
{
@@ -109,8 +109,8 @@ struct ornstein_det
struct ornstein_stoch
{
std::mt19937 m_rng;
std::normal_distribution<> m_dist;
boost::mt19937 m_rng;
boost::normal_distribution<> m_dist;
ornstein_stoch( double sigma ) : m_rng() , m_dist( 0.0 , sigma ) { }

1237
libs/numeric/odeint/ideas/adam_bashforrth_moulton/coefficients.nb
View File

File diff suppressed because it is too large Load Diff

13
libs/numeric/odeint/ideas/algebra/Jamfile
View File

@@ -1,13 +0,0 @@
# (C) Copyright 2010 : Karsten Ahnert, Mario Mulansky
# Distributed under the Boost Software License, Version 1.0.
# (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#project
# : requirements
# <define>BOOST_ALL_NO_LIB=1
# <include>../../../../..
# ;
exe test_algebra
: test_algebra.cpp
;

65
libs/numeric/odeint/ideas/algebra/algebra.hpp
View File

@@ -1,65 +0,0 @@
/*
* range_algebra.hpp
*
* Created on: Feb 12, 2011
* Author: karsten
*/
#ifndef RANGE_ALGEBRA_HPP_
#define RANGE_ALGEBRA_HPP_
#include <iostream>
#include <algorithm>
#include <boost/range.hpp>
#include <boost/functional/forward_adapter.hpp>
struct algebra1
{
struct for_each1_impl
{
template< class S1 , class Op >
void operator()( S1 &s1 , Op op ) const
{
std::for_each( boost::begin( s1 ) , boost::end( s1 ) , op );
tmp2 = 0.11;
}
typedef void result_type;
};
struct for_each2_impl
{
template< class S1 , class Op >
void operator()( S1 &s1 , Op op ) const
{
std::for_each( boost::begin( s1 ) , boost::end( s1 ) , op );
std::cout << tmp2 << std::endl;
}
typedef void result_type;
};
typedef boost::forward_adapter< for_each1_impl , 2 > for_each1;
double tmp2;
};
// alternativ
struct algebra2
{
template< class S1 , class Op >
void for_each1( S1 &s1 , Op op )
{
tmp = 123.0;
std::for_each( boost::begin( s1 ) , boost::end( s1 ) , op );
};
double tmp;
};
#endif /* RANGE_ALGEBRA_HPP_ */

41
libs/numeric/odeint/ideas/algebra/test_algebra.cpp
View File

@@ -1,41 +0,0 @@
/*
* test_algebra.cpp
*
* Created on: Feb 12, 2011
* Author: karsten
*/
#include <vector>
#include <iostream>
#include "algebra.hpp"
using namespace std;
struct print
{
template< class T >
void operator()( T t )
{
cout << t << endl;
}
};
int main( int argc , char **argv )
{
algebra2 al;
vector< double > vec( 3 );
vec[0] = 0.1;
vec[1] = 0.95;
vec[2] = 1.23;
al.for_each1( vec , print() );
cout << al.tmp << endl;
algebra1 al2;
// al2.for_each1()( vec , print() );
// al.for_each2()( vec , print() );
return 0;
}

39
libs/numeric/odeint/ideas/algebra_dispatcher.hpp
View File

@@ -1,39 +0,0 @@
/*
boost header: BOOST_NUMERIC_ODEINT/algebra_dispatcher.hpp
Copyright 2009 Karsten Ahnert
Copyright 2009 Mario Mulansky
Copyright 2009 Andre Bergner
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or
copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#ifndef BOOST_BOOST_NUMERIC_ODEINT_ALGEBRA_DISPATCHER_HPP_INCLUDED
#define BOOST_BOOST_NUMERIC_ODEINT_ALGEBRA_DISPATCHER_HPP_INCLUDED
namespace boost {
namespace numeric {
namespace odeint {
/*
* ToDo : Implement
*/
template< class Container , class Time >
struct algebra_dispatcher
{
typedef Container container_type;
typedef Time time_type;
typedef range_algebra< container_type > type;
};
} // odeint
} // numeric
} // boost
#endif //BOOST_BOOST_NUMERIC_ODEINT_ALGEBRA_DISPATCHER_HPP_INCLUDED

29
libs/numeric/odeint/ideas/butcher/Jamfile
View File

@@ -1,29 +0,0 @@
# (C) Copyright 2010 : Karsten Ahnert, Mario Mulansky
# Distributed under the Boost Software License, Version 1.0.
# (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
import os ;
import modules ;
import path ;
path-constant HOME : [ os.environ HOME ] ;
path-constant CHRONO_ROOT : [ os.environ CHRONO_ROOT ] ;
project
: requirements
<define>BOOST_ALL_NO_LIB=1
<include>../../../../..
<include>$(CHRONO_ROOT)
;
exe butcher_performance
: butcher_performance.cpp
;
exe butcher_lorenz
: butcher_lorenz.cpp
;
exe butcher_vdp
: butcher_vdp.cpp
;

199
libs/numeric/odeint/ideas/butcher/algebra.hpp
View File

@@ -1,199 +0,0 @@
/*
* algebra.hpp
*
* Created on: Nov 7, 2010
* Author: karsten
*/
#ifndef ALGEBRA_HPP_
#define ALGEBRA_HPP_
#include <boost/mpl/int.hpp>
#include <boost/mpl/at.hpp>
#include <boost/numeric/odeint/algebra/range_algebra.hpp>
#include <boost/numeric/odeint/algebra/default_operations.hpp>
#include "convert_value.hpp"
namespace mpl = boost::mpl;
template< class state_type , class coef_tye , class index >
struct algebra
{
static void do_step( state_type &x_tmp , const state_type &x , const state_type *k_vector , const double dt )
{}
};
template< class state_type , class coef_type >
struct algebra< state_type , coef_type , mpl::int_< 0 > >
{
typedef boost::numeric::odeint::range_algebra std_algebra;
typedef boost::numeric::odeint::default_operations std_op;
typedef typename state_type::iterator iterator;
typedef typename state_type::const_iterator const_iterator;
static void do_step( state_type &x_tmp , const state_type &x , const state_type *k_vector , const double dt )
{
const static double a1 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value();
// iterator first1 = x_tmp.begin() , last1 = x_tmp.end();
// const_iterator first2 = x.begin() , first3 = k_vector[0].begin();
// while( first1 != last1 )
// *first1++ = *first2++ + a1 * *first3++ ;
std_algebra::for_each3( x_tmp , x , k_vector[0] ,
std_op::scale_sum2<>( 1.0 , a1 ) );
// for( size_t i=0 ; i<x.size() ; ++i )
// {
// x_tmp[i] = x[i] + a1 * k_vector[0][i];
// }
}
};
template< class state_type , class coef_type >
struct algebra< state_type , coef_type , mpl::int_< 1 > >
{
typedef boost::numeric::odeint::range_algebra std_algebra;
typedef boost::numeric::odeint::default_operations std_op;
typedef typename state_type::iterator iterator;
typedef typename state_type::const_iterator const_iterator;
static void do_step( state_type &x_tmp , const state_type &x , const state_type *k_vector , const double dt )
{
const static double a1 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value();
const static double a2 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 1 > >::type >::get_value();
// iterator first1 = x_tmp.begin() , last1 = x_tmp.end();
// const_iterator first2 = x.begin() , first3 = k_vector[0].begin() , first4 = k_vector[1].begin();
// while( first1 != last1 )
// *first1++ = *first2++ + a1 * *first3++ + a2 * *first4++;
std_algebra::for_each4( x_tmp , x , k_vector[0] , k_vector[1] ,
std_op::scale_sum3<>( 1.0 , a1 , a2 ) );
// for( size_t i=0 ; i<x.size() ; ++i )
// {
// x_tmp[i] = x[i] +
// dt * convert_value< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value() * k_vector[0][i] +
// dt * convert_value< typename mpl::at< coef_type , mpl::int_< 1 > >::type >::get_value() * k_vector[1][i];
// }
}
};
template< class state_type , class coef_type >
struct algebra< state_type , coef_type , mpl::int_< 2 > >
{
typedef boost::numeric::odeint::range_algebra std_algebra;
typedef boost::numeric::odeint::default_operations std_op;
typedef typename state_type::iterator iterator;
typedef typename state_type::const_iterator const_iterator;
static void do_step( state_type &x_tmp , const state_type &x , const state_type *k_vector , const double dt )
{
const static double a1 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value();
const static double a2 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 1 > >::type >::get_value();
const static double a3 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 2 > >::type >::get_value();
//
// iterator first1 = x_tmp.begin() , last1 = x_tmp.end();
// const_iterator first2 = x.begin() , first3 = k_vector[0].begin() , first4 = k_vector[1].begin() , first5 = k_vector[2].begin();
// while( first1 != last1 )
// *first1++ = *first2++ + a1 * *first3++ + a2 * *first4++ + a3 * *first5++;
std_algebra::for_each5( x_tmp , x , k_vector[0] , k_vector[1] , k_vector[2] ,
std_op::scale_sum4<>( 1.0 , a1 , a2 , a3 ) );
// for( size_t i=0 ; i<x.size() ; ++i )
// {
// x_tmp[i] = x[i] +
// dt * convert_value< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value() * k_vector[0][i] +
// dt * convert_value< typename mpl::at< coef_type , mpl::int_< 1 > >::type >::get_value() * k_vector[1][i] +
// dt * convert_value< typename mpl::at< coef_type , mpl::int_< 2 > >::type >::get_value() * k_vector[2][i];
// }
}
};
template< class state_type , class coef_type >
struct algebra< state_type , coef_type , mpl::int_< 3 > >
{
typedef boost::numeric::odeint::range_algebra std_algebra;
typedef boost::numeric::odeint::default_operations std_op;
typedef typename state_type::iterator iterator;
typedef typename state_type::const_iterator const_iterator;
static void do_step( state_type &x_tmp , const state_type &x , const state_type *k_vector , const double dt )
{
const static double a1 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value();
const static double a2 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 1 > >::type >::get_value();
const static double a3 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 2 > >::type >::get_value();
const static double a4 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 3 > >::type >::get_value();
//
// iterator first1 = x_tmp.begin() , last1 = x_tmp.end();
// const_iterator first2 = x.begin() , first3 = k_vector[0].begin() , first4 = k_vector[1].begin() , first5 = k_vector[2].begin();
// const_iterator first6 = k_vector[3].begin();
// while( first1 != last1 )
// *first1++ = *first2++ + a1 * *first3++ + a2 * *first4++ + a3 * *first5++ + a4 * *first6++;
std_algebra::for_each6( x_tmp , x , k_vector[0] , k_vector[1] , k_vector[2] , k_vector[3] ,
std_op::scale_sum5<>( 1.0 , a1 , a2 , a3 , a4 ) );
// for( size_t i=0 ; i<x.size() ; ++i )
// {
// x_tmp[i] = x[i] +
// dt * convert_value< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value() * k_vector[0][i] +
// dt * convert_value< typename mpl::at< coef_type , mpl::int_< 1 > >::type >::get_value() * k_vector[1][i] +
// dt * convert_value< typename mpl::at< coef_type , mpl::int_< 2 > >::type >::get_value() * k_vector[2][i] +
// dt * convert_value< typename mpl::at< coef_type , mpl::int_< 3 > >::type >::get_value() * k_vector[3][i];
// }
}
};
template< class state_type , class coef_type >
struct algebra< state_type , coef_type , mpl::int_< 4 > >
{
typedef boost::numeric::odeint::range_algebra std_algebra;
typedef boost::numeric::odeint::default_operations std_op;
typedef typename state_type::iterator iterator;
typedef typename state_type::const_iterator const_iterator;
static void do_step( state_type &x_tmp , const state_type &x , const state_type *k_vector , const double dt )
{
const static double a1 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value();
const static double a2 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 1 > >::type >::get_value();
const static double a3 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 2 > >::type >::get_value();
const static double a4 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 3 > >::type >::get_value();
const static double a5 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 4 > >::type >::get_value();
//
// iterator first1 = x_tmp.begin() , last1 = x_tmp.end();
// const_iterator first2 = x.begin() , first3 = k_vector[0].begin() , first4 = k_vector[1].begin() , first5 = k_vector[2].begin();
// const_iterator first6 = k_vector[3].begin() , first7 = k_vector[4].begin();
// while( first1 != last1 )
// *first1++ = *first2++ + a1 * *first3++ + a2 * *first4++ + a3 * *first5++ + a4 * *first6++ + a5 * *first7++;
std_algebra::for_each7( x_tmp , x , k_vector[0] , k_vector[1] , k_vector[2] , k_vector[3] , k_vector[4] ,
std_op::scale_sum6<>( 1.0 , a1 , a2 , a3 , a4 , a5 ) );
// for( size_t i=0 ; i<x.size() ; ++i )
// {
// x_tmp[i] = x[i] +
// dt * convert_value< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value() * k_vector[0][i] +
// dt * convert_value< typename mpl::at< coef_type , mpl::int_< 1 > >::type >::get_value() * k_vector[1][i] +
// dt * convert_value< typename mpl::at< coef_type , mpl::int_< 2 > >::type >::get_value() * k_vector[2][i] +
// dt * convert_value< typename mpl::at< coef_type , mpl::int_< 3 > >::type >::get_value() * k_vector[3][i] +
// dt * convert_value< typename mpl::at< coef_type , mpl::int_< 4 > >::type >::get_value() * k_vector[4][i];
// }
}
};
#endif /* ALGEBRA_HPP_ */

66
libs/numeric/odeint/ideas/butcher/butcher_lorenz.cpp
View File

@@ -1,66 +0,0 @@
/*
* butcher_test.cpp
*
* Created on: Nov 5, 2010
* Author: karsten
*/
#include <iostream>
#include <fstream>
#include <tr1/array>
#include <boost/numeric/odeint/stepper/runge_kutta4_classic.hpp>
#include "predefined_steppers.hpp"
#define tab "\t"
using namespace std;
typedef std::tr1::array< double , 3 > state_type;
typedef mpl_rk4_stepper< state_type > stepper_type;
typedef boost::numeric::odeint::runge_kutta4< state_type > stepper_std_type;
const double sigma = 10.0;
const double R = 28.0;
const double b = 8.0 / 3.0;
void lorenz( const state_type &x , state_type &dxdt , double t )
{
dxdt[0] = sigma * ( x[1] - x[0] );
dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
dxdt[2] = x[0]*x[1] - b * x[2];
}
int main( int argc , char **argv )
{
stepper_type stepper;
stepper_std_type stepper_std;
state_type x = {{ 1.0 , 1.0 , 2.0 }};
state_type x_std = x;
cout << "Tableau : " << endl;
stepper.print_tableau();
cout << endl;
double t = 0.0 , dt = 0.01;
ofstream fout( "lorenz.dat" );
for( size_t i=0 ; i<10000 ; ++i )
{
fout << t << tab;
fout << x[0] << tab << x[1] << tab << x[2] << tab;
fout << x_std[0] << tab << x_std[1] << tab << x_std[2] << "\n";
stepper.do_step( lorenz , x , t , dt );
stepper_std.do_step( lorenz , x_std , t , dt );
t += dt;
}
return 0;
}

100
libs/numeric/odeint/ideas/butcher/butcher_performance.cpp
View File

@@ -1,100 +0,0 @@
/*
* butcher_test.cpp
*
* Created on: Nov 5, 2010
* Author: karsten
*/
#include <iostream>
#include <fstream>
#include <tr1/array>
#include <boost/numeric/odeint/stepper/euler.hpp>
#include <boost/accumulators/accumulators.hpp>
#include <boost/accumulators/statistics.hpp>
#include <boost/timer.hpp>
#include "predefined_steppers.hpp"
#define tab "\t"
using namespace std;
using namespace boost::accumulators;
typedef accumulator_set<
double , stats< tag::mean , tag::variance >
> accumulator_type;
ostream& operator<<( ostream& out , accumulator_type &acc )
{
out << boost::accumulators::mean( acc ) << tab;
// out << boost::accumulators::variance( acc ) << tab;
return out;
}
typedef boost::timer timer_type;
typedef std::tr1::array< double , 3 > state_type;
typedef mpl_euler_stepper< state_type > stepper_type;
typedef boost::numeric::odeint::euler< state_type > stepper_std_type;
void lorenz( const state_type &x , state_type &dxdt , double t )
{
const double sigma = 10.0;
const double R = 28.0;
const double b = 8.0 / 3.0;
dxdt[0] = sigma * ( x[1] - x[0] );
dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
dxdt[2] = x[0]*x[1] - b * x[2];
}
int main( int argc , char **argv )
{
stepper_type stepper;
stepper_std_type stepper_std;
const size_t num_of_steps = 20000000;
const size_t dt = 0.01;
accumulator_type acc1 , acc2;
timer_type timer;
srand48( 12312354 );
while( true )
{
state_type x = {{ 10.0 * drand48() , 10.0 * drand48() , 10.0 * drand48() }};
state_type x_std = x;
double t = 0.0 ;
double t_std = 0.0;
timer.restart();
for( size_t i=0 ; i<num_of_steps ; ++i,t+=dt )
stepper.do_step( lorenz , x , t , dt );
acc1( timer.elapsed() );
timer.restart();
for( size_t i=0 ; i<num_of_steps ; ++i,t_std+=dt )
stepper_std.do_step( lorenz , x_std , t_std , dt );
acc2( timer.elapsed() );
clog.precision( 3 );
clog.width( 5 );
clog << acc1 << " " << acc2 << " " << x[0] << " " << x_std[0] << endl;
}
return 0;
}
/*
* Compile with :
* g++ -O3 -I$BOOST_ROOT -I$HOME/boost/chrono -I$ODEINT_ROOT butcher_performance.cpp
*/

67
libs/numeric/odeint/ideas/butcher/butcher_vdp.cpp
View File

@@ -1,67 +0,0 @@
/*
* butcher_test.cpp
*
* Created on: Nov 5, 2010
* Author: karsten
*/
#include <iostream>
#include <fstream>
#include <tr1/array>
#include <boost/numeric/odeint/stepper/runge_kutta4_classic.hpp>
#include "predefined_steppers.hpp"
#define tab "\t"
using namespace std;
typedef std::tr1::array< double , 2 > state_type;
typedef mpl_rk4_stepper< state_type > stepper_type;
typedef boost::numeric::odeint::runge_kutta4< state_type > stepper_std_type;
void vdp( const state_type &x , state_type &dxdt , double t )
{
const double mu = 0.2;
dxdt[0] = x[1];
dxdt[1] = -x[0] + mu * ( 1.0 - x[0]*x[0] ) * x[1];
}
int main( int argc , char **argv )
{
stepper_type midpoint;
stepper_std_type midpoint_std;
state_type x = {{ 1.0 , 1.0 }};
state_type x_std = x;
cout << "Tableau : " << endl;
midpoint.print_tableau();
cout << endl;
double t = 0.0 , dt = 0.01;
ofstream fout( "vdp.dat" );
for( size_t i=0 ; i<10000 ; ++i )
{
fout << t << tab;
fout << x[0] << tab << x[1] << tab ;
fout << x_std[0] << tab << x_std[1] << "\n";
midpoint.do_step( vdp , x , t , dt );
midpoint_std.do_step( vdp , x_std , t , dt );
t += dt;
}
// cout << double( one_half::num ) / double( one_half::den ) << endl;
// cout << conv< one_half >::get_value() << endl;
// cout << conv< one >::get_value() << endl;
return 0;
}

36
libs/numeric/odeint/ideas/butcher/convert_value.hpp
View File

@@ -1,36 +0,0 @@
/*
* convert_values.hpp
*
* Created on: Nov 7, 2010
* Author: karsten
*/
#ifndef CONVERT_VALUE_HPP_
#define CONVERT_VALUE_HPP_
#include <boost/ratio.hpp>
template< class T >
struct convert_value
{
static double get_value( void )
{
return double( T::value );
}
};
template< long N , long D >
struct convert_value< boost::ratio< N , D > >
{
typedef typename boost::ratio< N , D > number;
static double get_value( void )
{
return double( number::num ) / double( number::den );
}
};
#endif /* CONVERT_VALUE_HPP_ */

59
libs/numeric/odeint/ideas/butcher/diagnostic.hpp
View File

@@ -1,59 +0,0 @@
/*
* diagnostic.hpp
*
* Created on: Nov 7, 2010
* Author: karsten
*/
#ifndef DIAGNOSTIC_HPP_
#define DIAGNOSTIC_HPP_
#include <iostream>
#include <boost/mpl/for_each.hpp>
#include <boost/mpl/at.hpp>
#include <boost/mpl/int.hpp>
#include <boost/mpl/size.hpp>
#include "convert_value.hpp"
namespace mpl = boost::mpl;
struct print_value
{
template< class T >
void operator()( T )
{
std::cout << convert_value< T >::get_value() << " ";
}
};
struct print_vector
{
template< class Vec >
void operator()( Vec )
{
mpl::for_each< Vec >( print_value() );
std::cout << "\n";
}
};
struct print_tableau_vector
{
template< class Vec >
void operator()( Vec )
{
typedef typename mpl::at< Vec , mpl::int_< 0 > >::type index;
typedef typename mpl::at< Vec , mpl::int_< 1 > >::type c;
typedef typename mpl::at< Vec , mpl::int_< 2 > >::type coef;
std::cout << "Size : " << mpl::size< Vec >::value << "\t";
std::cout << convert_value< index >::get_value() << " ";
std::cout << convert_value< c >::get_value() << " ";
mpl::for_each< coef >( print_value() );
std::cout << "\n";
}
};
#endif /* DIAGNOSTIC_HPP_ */

159
libs/numeric/odeint/ideas/butcher/explicit_runge_kutta.hpp
View File

@@ -1,159 +0,0 @@
/*
* explicit_runge_kutta.hpp
*
* Created on: Nov 5, 2010
* Author: karsten
*/
#ifndef EXPLICIT_RUNGE_KUTTA_HPP_
#define EXPLICIT_RUNGE_KUTTA_HPP_
#include <boost/static_assert.hpp>
#include <boost/ratio.hpp>
#include <boost/type_traits/is_same.hpp>
#include <boost/mpl/vector.hpp>
#include <boost/mpl/size.hpp>
#include <boost/mpl/push_back.hpp>
#include <boost/mpl/for_each.hpp>
#include <boost/mpl/zip_view.hpp>
#include <boost/mpl/range_c.hpp>
#include <boost/mpl/int.hpp>
#include <boost/mpl/at.hpp>
#include "convert_value.hpp"
#include "diagnostic.hpp"
#include "algebra.hpp"
namespace mpl = boost::mpl;
/*
* c_1 |
* c_2 | a_{2,1}
* c_3 | a_{3,1} a_{3,2}
* ...
* c_s | a_{s,1} a_{s,2} ... a_{s,s-1}
* -----------------------------------
* | b_1 b_2 b_{s-1} b_s
*
* ButcherA : the coefficients a_ij
*
* ButcherC : the time coefficients c_i
*
* ButcherC : the sum coefficients b_j
*/
template< class State , class ButcherA , class ButcherC , class ButcherB >
class explicit_runge_kutta
{
public:
typedef ButcherA butcher_a;
typedef ButcherB butcher_b;
typedef ButcherC butcher_c;
static const size_t dim = mpl::size< ButcherC >::value;
typedef mpl::range_c< size_t , 0 , dim > indices;
typedef typename mpl::push_back< butcher_a , butcher_b >::type butcher_right;
typedef mpl::zip_view< mpl::vector< indices , butcher_c , butcher_right > > butcher_tableau;
typedef double time_type;
typedef State state_type;
/*
* t , dt , system, x, k[]
*/
template< class System >
struct calculate_stage
{
System &system;
state_type &x , &x_tmp;
state_type *k_vector;
const double t;
const double dt;
calculate_stage( System &_system , state_type &_x , state_type &_x_tmp , state_type *_k_vector , const double _t , const double _dt )
: system( _system ) , x( _x ) , x_tmp( _x_tmp ) , k_vector( _k_vector ) , t( _t ) , dt( _dt )
{}
/*
* for( size_t j=0 ; j<dim ; ++j )
* {
* system( x[j] , k[j] , c[j] )
* x[j+1] = sum( i , a[j][i] , k[j] )
* }
*/
template< class Stage > void operator()( Stage v )
{
typedef typename mpl::at< Stage , mpl::int_< 0 > >::type index_type;
typedef typename mpl::at< Stage , mpl::int_< 1 > >::type time_value_type;
typedef typename mpl::at< Stage , mpl::int_< 2 > >::type coef_type;
const static size_t index = index_type::value;
const static double time_value = convert_value< time_value_type >::get_value();
BOOST_STATIC_ASSERT(( ( mpl::size< coef_type >::value - 1 ) == index ));
// cout << index << " " << time_value << " " << endl;
if( index == 0 ) system( x , k_vector[index] , t + time_value * dt );
else system( x_tmp , k_vector[index] , t + time_value * dt );
if( index == ( dim - 1 ) )
algebra< state_type , coef_type , mpl::int_< index > >::do_step( x , x , k_vector , dt );
else
algebra< state_type , coef_type , mpl::int_< index > >::do_step( x_tmp , x , k_vector , dt );
};
};
explicit_runge_kutta( void )
{
}
template< class System >
void do_step( System &system , state_type &x , double t , const double dt )
{
mpl::for_each< butcher_tableau >( calculate_stage< System >( system , x , m_x_tmp , m_k_vector , t , dt ) );
}
void print_butcher_a( void )
{
mpl::for_each< butcher_a >( print_vector() );
}
void print_butcher_b( void )
{
mpl::for_each< butcher_b >( print_value() );
}
void print_butcher_c( void )
{
mpl::for_each< butcher_c >( print_value() );
}
void print_tableau( void )
{
mpl::for_each< butcher_tableau >( print_tableau_vector() );
}
private:
state_type m_k_vector[dim];
state_type m_x_tmp;
BOOST_STATIC_ASSERT(( dim > 0 ));
BOOST_STATIC_ASSERT(( dim == size_t( mpl::size< butcher_b >::value ) ));
BOOST_STATIC_ASSERT(( ( dim - 1 ) == mpl::size< butcher_a >::value ));
};
#endif /* EXPLICIT_RUNGE_KUTTA_HPP_ */

77
libs/numeric/odeint/ideas/butcher/predefined_steppers.hpp
View File

@@ -1,77 +0,0 @@
/*
* predefined_steppers.hpp
*
* Created on: Nov 7, 2010
* Author: karsten
*/
#ifndef PREDEFINED_STEPPERS_HPP_
#define PREDEFINED_STEPPERS_HPP_
#include <boost/mpl/int.hpp>
#include <boost/ratio.hpp>
#include "explicit_runge_kutta.hpp"
typedef mpl::int_< 0 > null;
typedef mpl::int_< 1 > one;
typedef boost::ratio< 1 , 2 > one_half;
typedef boost::ratio< 1 , 3 > one_third;
typedef boost::ratio< 1 , 6 > one_sixth;
/*
* euler :
* 0 |
* | 1
*/
typedef mpl::vector< null > euler_c;
typedef mpl::vector<> euler_a;
typedef mpl::vector< one > euler_b;
template< class state_type >
struct mpl_euler_stepper : public explicit_runge_kutta< state_type , euler_a , euler_c , euler_b > { };
/*
* midpoint :
* 0 |
* 1/2 | 1/2
* -----------
* | 0 1
*/
typedef mpl::vector< null , one_half > midpoint_c;
typedef mpl::vector< mpl::vector< one_half > > midpoint_a;
typedef mpl::vector< null , one > midpoint_b;
template< class state_type >
struct mpl_midpoint_stepper : public explicit_runge_kutta< state_type , midpoint_a , midpoint_c , midpoint_b > { };
/*
* classical Runge Kutta
* 0 |
* 1/2 | 1/2
* 1/2 | 0 1/2
* 1 | 0 0 1
* ------------------
* | 1/6 1/3 1/6 1/6
*
*/
typedef mpl::vector< null , one_half , one_half , one > rk4_c;
typedef mpl::vector<
mpl::vector< one_half > ,
mpl::vector< null , one_half > ,
mpl::vector< null , null , one >
> rk4_a;
typedef mpl::vector< one_sixth , one_third , one_third , one_sixth > rk4_b;
template< class state_type >
struct mpl_rk4_stepper : public explicit_runge_kutta< state_type , rk4_a , rk4_c , rk4_b > { };
#endif /* PREDEFINED_STEPPERS_HPP_ */

22
libs/numeric/odeint/ideas/controller/Jamfile
View File

@@ -1,22 +0,0 @@
# (C) Copyright 2010 : Karsten Ahnert, Mario Mulansky
# Distributed under the Boost Software License, Version 1.0.
# (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
import os ;
import modules ;
import path ;
project
: requirements
<define>BOOST_ALL_NO_LIB=1
<include>../../../../..
<include>.
;
exe controller : controller.cpp ;
exe rb_controller : rb_controller.cpp ;
exe dopri5_controller : dopri5_controller.cpp ;
exe dense_out_dopri5_controller : dense_out_dopri5_controller.cpp ;
exe rb_dense_output : rb_dense_output.cpp ;
build-project libs/numeric/odeint/test ;

70
libs/numeric/odeint/ideas/controller/boost/numeric/odeint/stepper/controller/default_controller.hpp
View File

@@ -1,70 +0,0 @@
/*
[auto_generated]
boost/numeric/odeint/stepper/controller/default_controller.hpp
[begin_description]
Default controller for the use in the generic steppers.
[end_description]
Copyright 2009-2011 Karsten Ahnert
Copyright 2009-2011 Mario Mulansky
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or
copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#ifndef BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLER_DEFAULT_CONTROLLER_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLER_DEFAULT_CONTROLLER_HPP_INCLUDED
#include <cmath>
#include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
namespace boost {
namespace numeric {
namespace odeint {
class default_controller
{
public:
template< class Value , class Time , class Order >
controlled_step_result operator()( Value error , Time &t , Time &dt , Order stepper_order , Order error_order )
{
using std::max;
using std::min;
using std::pow;
if( error > 1.0 )
{
// error too large - decrease dt ,limit scaling factor to 0.2 and reset state
dt *= max( 0.9 * pow( error , -1.0 / ( Value( error_order ) - 1.0 ) ) , 0.2 );
return fail;
}
else
{
if( error < 0.5 )
{
//error too small - increase dt and keep the evolution and limit scaling factor to 5.0
t += dt;
dt *= min( 0.9 * pow( error , -1.0 / Value( stepper_order ) ) , 5.0 );
return success;
}
else
{
t += dt;
return success;
}
}
}
};
}
}
}
#endif // BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLER_DEFAULT_CONTROLLER_HPP_INCLUDED

114
libs/numeric/odeint/ideas/controller/boost/numeric/odeint/stepper/controller/pi_controller.hpp
View File

@@ -1,114 +0,0 @@
/*
[auto_generated]
boost/numeric/odeint/stepper/controller/pi_controller.hpp
[begin_description]
PI-controller for the use in the generic controller steppers.
[end_description]
Copyright 2009-2011 Karsten Ahnert
Copyright 2009-2011 Mario Mulansky
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or
copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#ifndef BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLER_PI_CONTROLLER_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLER_PI_CONTROLLER_HPP_INCLUDED
#include <cmath>
#include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
namespace boost {
namespace numeric {
namespace odeint {
template< class Value >
class pi_controller
{
public:
typedef Value value_type;
pi_controller( void )
: m_err_old( static_cast< value_type >( 1.0e-4 ) ) , m_last_rejected( false ) { }
/*
* Returns true if err Ä 1, false otherwise. If step was successful, sets hnext to the estimated
* optimal stepsize for the next step. If the step failed, reduces h appropriately for another try.
*
* Set beta to a nonzero value for PI control. beta D 0:04–0.08 is a good default.
*/
template< class Time , class Order >
controlled_step_result operator()( value_type error , Time &t , Time &dt , Order stepper_order , Order error_order )
{
using std::max;
using std::min;
using std::pow;
// ToDo : adapt for values of stepper_order or error_order
const value_type beta = 0.02;
const value_type alpha = 0.2 - beta * 0.75;
const value_type safe = 0.9;
const value_type min_scale = 0.2;
const value_type max_scale = 10.0;
if( error <= 1.0 )
{
value_type scale = 1.0;
if( error == 0.0 )
{
scale = max_scale;
}
else
{
// PI control if beta != 0.
scale = safe * pow( error ,-alpha ) * pow( m_err_old , beta );
if( scale < min_scale )
scale = min_scale;
if( scale > max_scale )
scale = max_scale;
}
if( m_last_rejected )
{
dt *= min( scale , 1.0 );
}
else
{
dt *= scale;
}
m_err_old = std::max( error , 1.0e-4 );
m_last_rejected = false;
t += dt;
return success;
}
else
{
value_type scale = max( safe * pow( error , -alpha ) , min_scale );
dt *= scale;
m_last_rejected = true;
return fail;
}
}
private:
value_type m_err_old;
bool m_last_rejected;
};
}
}
}
#endif // BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLER_PI_CONTROLLER_HPP_INCLUDED

99
libs/numeric/odeint/ideas/controller/boost/numeric/odeint/stepper/controller/stiff_controller.hpp
View File

@@ -1,99 +0,0 @@
/*
[auto_generated]
boost/numeric/odeint/stepper/controller/stiff_controller.hpp
[begin_description]
Controller for the use in the generic steppers in combination with solvers for stiff systems.
[end_description]
Copyright 2009-2011 Karsten Ahnert
Copyright 2009-2011 Mario Mulansky
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or
copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#ifndef BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLER_STIFF_CONTROLLER_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLER_STIFF_CONTROLLER_HPP_INCLUDED
#include <cmath>
#include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
namespace boost {
namespace numeric {
namespace odeint {
template< class Value >
class stiff_controller
{
public:
typedef Value value_type;
stiff_controller( void )
: m_first_step( true ) , m_last_rejected( false ) ,
m_err_old( 0.0 ) , m_dt_old( 0.0 )
{ }
template< class Order >
controlled_step_result operator()( value_type error , value_type &t , value_type &dt , Order stepper_order , Order error_order )
{
typedef Value value_type;
static const value_type safe = 0.9 , fac1 = 5.0 , fac2 = 1.0 / 6.0;
value_type fac = std::max( fac2 ,std::min( fac1 , std::pow( error , 0.25 ) / safe ) );
value_type dt_new = dt / fac;
if ( error <= 1.0 )
{
if( m_first_step )
{
m_first_step = false;
}
else
{
value_type fac_pred = ( m_dt_old / dt ) * pow( error * error / 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 , error );
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;
return success;
}
else
{
dt = dt_new;
m_last_rejected = true;
return fail;
}
}
private:
bool m_first_step;
bool m_last_rejected;
value_type m_err_old;
value_type m_dt_old;
};
}
}
}
#endif // BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLER_STIFF_CONTROLLER_HPP_INCLUDED

137
libs/numeric/odeint/ideas/controller/boost/numeric/odeint/stepper/error_checker_explicit.hpp
View File

@@ -1,137 +0,0 @@
/*
[auto_generated]
boost/numeric/odeint/stepper/error_checker_explicit.hpp
[begin_description]
Error checker for the use in the generic_controlled_steppers in combination with explicit_error_steppers.
[end_description]
Copyright 2009-2011 Karsten Ahnert
Copyright 2009-2011 Mario Mulansky
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or
copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#ifndef BOOST_NUMERIC_ODEINT_STEPPER_ERROR_CHECKER_EXPLICIT_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_STEPPER_ERROR_CHECKER_EXPLICIT_HPP_INCLUDED
#include <boost/numeric/odeint/algebra/default_operations.hpp>
namespace boost {
namespace numeric {
namespace odeint {
/*
* Error checker for controlled_error_stepper
*
* ToDo: remove, is redundant now
*/
template< class Value , class Algebra , class Operations >
class error_checker_explicit
{
public:
typedef Value value_type;
typedef Algebra algebra_type;
typedef Operations operations_type;
error_checker_explicit(
const value_type eps_abs = static_cast< value_type >( 1.0e-6 ) ,
const value_type eps_rel = static_cast< value_type >( 1.0e-6 ) ,
const value_type a_x = static_cast< value_type >( 1.0 ) ,
const value_type a_dxdt = static_cast< value_type >( 1.0 ) )
: m_eps_abs( eps_abs ) , m_eps_rel( eps_rel ) , m_a_x( a_x ) , m_a_dxdt( a_dxdt )
{ }
// template< class State1 , class State2 , class Err , class Time >
// value_type error( const State1 &x_old , const State2 &x , const Err &x_err , const Time &dt )
// {
// algebra.reduce3( x_err , x_old , x , typename operations_type::template rel_error< value_type >( m_eps_abs , m_eps_rel ) );
// }
template< class State1 , class State2 , class Deriv , class Err , class Time >
value_type error( const State1 &x_old , const State2 &x , const Deriv &dxdt_old , Err &x_err , const Time &dt )
{
// this overwrites x_err !
algebra.for_each3( x_err , x_old , dxdt_old ,
typename operations_type::template rel_error< value_type >( m_eps_abs , m_eps_rel , m_a_x , m_a_dxdt * detail::get_value( dt ) ) );
value_type res = algebra.reduce( x_err ,
typename operations_type::template maximum< value_type >() , static_cast< value_type >( 0.0 ) );
return res;
}
private:
value_type m_eps_abs;
value_type m_eps_rel;
value_type m_a_x;
value_type m_a_dxdt;
algebra_type algebra;
};
/*
* Error checker for controlled_error_stepper
*
* ToDo: rename to error_checker_explicit_new
*/
template< class Value , class Algebra , class Operations >
class error_checker_explicit_new
{
public:
typedef Value value_type;
typedef Algebra algebra_type;
typedef Operations operations_type;
error_checker_explicit_new(
algebra_type &algebra ,
const value_type eps_abs = static_cast< value_type >( 1.0e-6 ) ,
const value_type eps_rel = static_cast< value_type >( 1.0e-6 ) ,
const value_type a_x = static_cast< value_type >( 1.0 ) ,
const value_type a_dxdt = static_cast< value_type >( 1.0 ) )
: m_algebra( algebra ) , m_eps_abs( eps_abs ) , m_eps_rel( eps_rel ) , m_a_x( a_x ) , m_a_dxdt( a_dxdt )
{ }
template< class State1 , class State2 , class Deriv , class Err , class Time >
value_type error( const State1 &x_old , const State2 &x , const Deriv &dxdt_old , Err &x_err , const Time &dt )
{
// this overwrites x_err !
m_algebra.for_each3( x_err , x_old , dxdt_old ,
typename operations_type::template rel_error< value_type >( m_eps_abs , m_eps_rel , m_a_x , m_a_dxdt * detail::get_value( dt ) ) );
value_type res = m_algebra.reduce( x_err ,
typename operations_type::template maximum< value_type >() , static_cast< value_type >( 0.0 ) );
return res;
}
private:
algebra_type &m_algebra;
value_type m_eps_abs;
value_type m_eps_rel;
value_type m_a_x;
value_type m_a_dxdt;
};
}
}
}
#endif // BOOST_NUMERIC_ODEINT_STEPPER_ERROR_CHECKER_EXPLICIT_HPP_INCLUDED

119
libs/numeric/odeint/ideas/controller/boost/numeric/odeint/stepper/error_checker_l2_norm.hpp
View File

@@ -1,119 +0,0 @@
/*
[auto_generated]
boost/numeric/odeint/stepper/default_error_checker_l2_norm.hpp
[begin_description]
Default error checker for the use in the generic_controlled_steppers. Works with all error_steppers, explicit_error_steppers and
explicit_error_stepper_fsal;
[end_description]
Copyright 2009-2011 Karsten Ahnert
Copyright 2009-2011 Mario Mulansky
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or
copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#ifndef BOOST_NUMERIC_ODEINT_STEPPER_ERROR_CHECKER_L2_NORM_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_STEPPER_ERROR_CHECKER_L2_NORM_HPP_INCLUDED
#include <cmath>
#include <boost/numeric/odeint/util/ref_or_value_holder.hpp>
#include <boost/numeric/odeint/util/number_of_elements.hpp>
namespace boost {
namespace numeric {
namespace odeint {
/*
* Calculate the error in a generic way:
*
* for steppers: err = ||( | err[i] | / ( eps_abs + eps_rel * | x[i] | ) )||
* for explicit steppers: err = ||( | err[i] | ( eps_abs + eps_rel * ( a_x * | x[i] | + a_dxdt * | dxdt[i] * dt ] ) ) )||
* for explicit fsal steppers: err = ||( | err[i] | ( eps_abs + eps_rel * ( a_x * | x[i] | + a_dxdt * | dxdt[i] * dt ] ) ) )||
*/
template< class Value , class Algebra , class Operations , bool HoldsAlgebra = true >
class error_checker_l2_norm
{
public:
typedef Value value_type;
typedef Algebra algebra_type;
typedef Operations operations_type;
const static bool holds_algebra = HoldsAlgebra;
const static bool ref_algebra = !holds_algebra;
typedef ref_or_value_holder< algebra_type , ref_algebra > algebra_holder_type;
error_checker_l2_norm(
typename algebra_holder_type::constructor_type algebra ,
const value_type eps_abs = static_cast< value_type >( 1.0e-6 ) ,
const value_type eps_rel = static_cast< value_type >( 1.0e-6 ) ,
const value_type a_x = static_cast< value_type >( 1.0 ) ,
const value_type a_dxdt = static_cast< value_type >( 1.0 ))
: m_algebra( algebra ) ,
m_eps_abs( eps_abs ) , m_eps_rel( eps_rel ) , m_a_x( a_x ) , m_a_dxdt( a_dxdt )
{ }
// overload for steppers, x, x_old and x_err are available
template< class State1 , class State2 , class Err , class Time >
value_type error( const State1 &x_old , const State2 &x , const Err &x_err , const Time &dt )
{
using std::sqrt;
value_type result = m_algebra.get().reduce3( x_old , x , x_err ,
typename operations_type::template rel_error_l2< value_type >( m_eps_abs , m_eps_rel ) ,
static_cast< value_type >( 0.0 ) );
size_t num_elements = boost::numeric::odeint::number_of_elements( x_old );
if( num_elements != 0 ) return sqrt( result ) / static_cast< value_type >( num_elements );
return 0.0;
}
// overload for explicit steppers, x, x_old, dxdt_old and x_err are available
template< class State1 , class State2 , class Deriv , class Err , class Time >
value_type error( const State1 &x_old , const State2 &x , const Deriv &dxdt_old , Err &x_err , const Time &dt )
{
using std::sqrt;
value_type result = m_algebra.get().reduce4( x_old , x , dxdt_old , x_err ,
typename operations_type::template rel_error_l2_2< value_type >( m_eps_abs , m_eps_rel , m_a_x , m_a_dxdt * detail::get_value( dt ) ) ,
static_cast< value_type >( 0.0 ) );
size_t num_elements = boost::numeric::odeint::number_of_elements( x_old );
if( num_elements != 0 ) return sqrt( result ) / static_cast< value_type >( num_elements );
return 0.0;
}
// overload for explicit fsal steppers, x, x_old, dxdt, dxdt_old and x_err are available
template< class StateOld , class State , class DerivOld , class Deriv , class Err , class Time >
value_type error( const StateOld &x_old , const State &x , const DerivOld &dxdt_old , const Deriv &dxdt , const Err &x_err , const Time &dt )
{
using std::sqrt;
value_type result = m_algebra.get().reduce4( x_old , x , dxdt_old , x_err ,
typename operations_type::template rel_error_l2_2< value_type >( m_eps_abs , m_eps_rel , m_a_x , m_a_dxdt * detail::get_value( dt ) ) ,
static_cast< value_type >( 0.0) );
size_t num_elements = boost::numeric::odeint::number_of_elements( x_old );
if( num_elements != 0 ) return sqrt( result ) / static_cast< value_type >( num_elements );
return 0.0;
}
private:
algebra_holder_type m_algebra;
value_type m_eps_abs;
value_type m_eps_rel;
value_type m_a_x;
value_type m_a_dxdt;
};
}
}
}
#endif // BOOST_NUMERIC_ODEINT_STEPPER_ERROR_CHECKER_EXPLICIT_L2_NORM_HPP_INCLUDED

109
libs/numeric/odeint/ideas/controller/boost/numeric/odeint/stepper/error_checker_max_norm.hpp
View File

@@ -1,109 +0,0 @@
/*
[auto_generated]
boost/numeric/odeint/stepper/default_error_checker_max_norm.hpp
[begin_description]
Default error checker for the use in the generic_controlled_steppers. Works with all error_steppers, explicit_error_steppers and
explicit_error_stepper_fsal;
[end_description]
Copyright 2009-2011 Karsten Ahnert
Copyright 2009-2011 Mario Mulansky
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or
copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#ifndef BOOST_NUMERIC_ODEINT_STEPPER_ERROR_CHECKER_MAX_NORM_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_STEPPER_ERROR_CHECKER_MAX_NORM_HPP_INCLUDED
#include <boost/numeric/odeint/util/ref_or_value_holder.hpp>
namespace boost {
namespace numeric {
namespace odeint {
/*
* Calculate the error in a generic way:
*
* for steppers: err = max_i( | err[i] | / ( eps_abs + eps_rel * | x[i] | ) )
* for explicit steppers: err = max_i( | err[i] | ( eps_abs + eps_rel * ( a_x * | x[i] | + a_dxdt * | dxdt[i] * dt ] ) ) )
* for explicit fsal steppers: err = max_i( | err[i] | ( eps_abs + eps_rel * ( a_x * | x[i] | + a_dxdt * | dxdt[i] * dt ] ) ) )
*/
template< class Value , class Algebra , class Operations , bool HoldsAlgebra = true >
class error_checker_max_norm
{
public:
typedef Value value_type;
typedef Algebra algebra_type;
typedef Operations operations_type;
const static bool holds_algebra = HoldsAlgebra;
const static bool ref_algebra = !holds_algebra;
typedef ref_or_value_holder< algebra_type , ref_algebra > algebra_holder_type;
error_checker_max_norm(
typename algebra_holder_type::constructor_type algebra ,
const value_type eps_abs = static_cast< value_type >( 1.0e-6 ) ,
const value_type eps_rel = static_cast< value_type >( 1.0e-6 ) ,
const value_type a_x = static_cast< value_type >( 1.0 ) ,
const value_type a_dxdt = static_cast< value_type >( 1.0 ))
: m_algebra( algebra ) ,
m_eps_abs( eps_abs ) , m_eps_rel( eps_rel ) , m_a_x( a_x ) , m_a_dxdt( a_dxdt )
{ }
// overload for steppers, x, x_old and x_err are available
template< class State1 , class State2 , class Err , class Time >
value_type error( const State1 &x_old , const State2 &x , const Err &x_err , const Time &dt )
{
value_type result = m_algebra.get().reduce3( x_old , x , x_err ,
typename operations_type::template rel_error_max< value_type >( m_eps_abs , m_eps_rel ) ,
static_cast< value_type >( 0.0 ) );
return result;
}
// overload for explicit steppers, x, x_old, dxdt_old and x_err are available
template< class State1 , class State2 , class Deriv , class Err , class Time >
value_type error( const State1 &x_old , const State2 &x , const Deriv &dxdt_old , Err &x_err , const Time &dt )
{
value_type result = m_algebra.get().reduce4( x_old , x , dxdt_old , x_err ,
typename operations_type::template rel_error_max2< value_type >( m_eps_abs , m_eps_rel , m_a_x , m_a_dxdt * detail::get_value( dt ) ) ,
static_cast< value_type >( 0.0 ) );
return result;
}
// overload for explicit fsal steppers, x, x_old, dxdt, dxdt_old and x_err are available
template< class StateOld , class State , class DerivOld , class Deriv , class Err , class Time >
value_type error( const StateOld &x_old , const State &x , const DerivOld &dxdt_old , const Deriv &dxdt , const Err &x_err , const Time &dt )
{
value_type result = m_algebra.get().reduce4( x_old , x , dxdt_old , x_err ,
typename operations_type::template rel_error_max2< value_type >( m_eps_abs , m_eps_rel , m_a_x , m_a_dxdt * detail::get_value( dt ) ) ,
static_cast< value_type >( 0.0) );
return result;
}
private:
algebra_holder_type m_algebra;
value_type m_eps_abs;
value_type m_eps_rel;
value_type m_a_x;
value_type m_a_dxdt;
};
}
}
}
#endif // BOOST_NUMERIC_ODEINT_STEPPER_ERROR_CHECKER_EXPLICIT_MAX_NORM_HPP_INCLUDED

179
libs/numeric/odeint/ideas/controller/boost/numeric/odeint/stepper/generic_controlled_stepper.hpp
View File

@@ -1,179 +0,0 @@
/*
[auto_generated]
boost/numeric/odeint/stepper/generic_controlled_stepper.hpp
[begin_description]
Implementation of the generic_controlled_stepper for the error_stepper_tag. Specialized versions
for explicit_error_stepper and explicit_error_stepper_tag also exist.
[end_description]
Copyright 2009-2011 Karsten Ahnert
Copyright 2009-2011 Mario Mulansky
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or
copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#ifndef BOOST_NUMERIC_ODEINT_STEPPER_GENERIC_CONTROLLED_STEPPER_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_STEPPER_GENERIC_CONTROLLED_STEPPER_HPP_INCLUDED
#include <boost/numeric/odeint/stepper/stepper_categories.hpp>
#include <boost/numeric/odeint/util/is_resizeable.hpp>
#include <boost/numeric/odeint/util/state_wrapper.hpp>
#include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
#include <boost/numeric/odeint/stepper/generic_controlled_stepper_definition.hpp>
#include <boost/bind.hpp>
namespace boost {
namespace numeric {
namespace odeint {
template<
class ErrorStepper ,
class ErrorChecker ,
class Controller ,
class Resizer >
class generic_controlled_stepper< ErrorStepper , ErrorChecker , Controller , Resizer , error_stepper_tag >
{
public:
typedef ErrorStepper stepper_type;
typedef ErrorChecker error_checker_type;
typedef Controller controller_type;
typedef Resizer resizer_type;
typedef typename stepper_type::state_type state_type;
typedef typename stepper_type::value_type value_type;
typedef typename stepper_type::time_type time_type;
typedef typename stepper_type::order_type order_type;
typedef explicit_controlled_stepper_tag stepper_category;
typedef state_wrapper< state_type > wrapped_state_type;
generic_controlled_stepper(
const stepper_type &stepper = stepper_type() ,
const error_checker_type &error_checker = error_checker_type() ,
const controller_type &controller = controller_type() )
: m_stepper( stepper ) , m_error_checker( error_checker ) , m_controller( controller ) ,
m_xerr_resizer() , m_xnew_resizer() ,
m_xerr() , m_xnew() { }
/*
* Version 1 : try_step( sys , x , t , dt )
*
* The overloads are needed to solve the forwarding problem
*/
template< class System , class StateInOut >
controlled_step_result try_step( System system , StateInOut &x , time_type &t , time_type &dt )
{
m_xnew_resizer.adjust_size( x , boost::bind( &generic_controlled_stepper::template resize_m_xnew_impl< state_type > , boost::ref( *this ) , _1 ) );
controlled_step_result result = try_step( system , x , t , m_xnew.m_v , dt );
if( result == success )
{
boost::numeric::odeint::copy( m_xnew.m_v , x );
}
return result;
}
template< class System , class StateInOut >
controlled_step_result try_step( System system , const StateInOut &x , time_type &t , time_type &dt )
{
m_xnew_resizer.adjust_size( x , boost::bind( &generic_controlled_stepper::template resize_m_xnew_impl< state_type > , boost::ref( *this ) , _1 ) );
controlled_step_result result = try_step( system , x , t , m_xnew.m_v , dt );
if( result == success )
{
boost::numeric::odeint::copy( m_xnew.m_v , x );
}
return result;
}
/*
* Version 2 : try_step( sys , in , t , out , dt )
*
* this version does not solve the forwarding problem, boost.range can not be used
*/
template< class System , class StateIn , class StateOut >
controlled_step_result try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt )
{
// return try_step( system , in , m_dxdt.m_v , t , out , dt );
// ToDo : implement
m_xerr_resizer.adjust_size( in , boost::bind( &generic_controlled_stepper::template resize_m_xerr_impl< StateIn > , boost::ref( *this ) , _1 ) );
m_stepper.do_step( system , in , t , out , dt , m_xerr.m_v );
value_type error = m_error_checker.error( in , out , m_xerr.m_v , dt );
return m_controller( error , t , dt , m_stepper.stepper_order , m_stepper.error_order );
}
template< class StateType >
void adjust_size( const StateType &x )
{
resize_m_xerr_impl( x );
resize_m_xnew_impl( x );
m_stepper.adjust_size( x );
}
stepper_type& stepper( void ) { return m_stepper; }
const stepper_type stepper( void ) const { return m_stepper; }
private:
template< class StateIn >
bool resize_m_xerr_impl( const StateIn &x )
{
return adjust_size_by_resizeability( m_xerr , x , typename is_resizeable< wrapped_state_type >::type() );
}
template< class StateIn >
bool resize_m_xnew_impl( const StateIn &x )
{
return adjust_size_by_resizeability( m_xnew , x , typename is_resizeable< wrapped_state_type >::type() );
}
stepper_type m_stepper;
error_checker_type m_error_checker;
controller_type m_controller;
resizer_type m_xerr_resizer;
resizer_type m_xnew_resizer;
wrapped_state_type m_xerr;
wrapped_state_type m_xnew;
};
}
}
}
#endif // BOOST_NUMERIC_ODEINT_STEPPER_GENERIC_CONTROLLED_STEPPER_HPP_INCLUDED

43
libs/numeric/odeint/ideas/controller/boost/numeric/odeint/stepper/generic_controlled_stepper_definition.hpp
View File

@@ -1,43 +0,0 @@
/*
[auto_generated]
boost/numeric/odeint/stepper/generic_controlled_stepper_definition.hpp
[begin_description]
Definition of the generic controlled stepper consisting of an error stepper, an error checker,
and a controller.
[end_description]
Copyright 2009-2011 Karsten Ahnert
Copyright 2009-2011 Mario Mulansky
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or
copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#ifndef BOOST_NUMERIC_ODEINT_STEPPER_GENERIC_CONTROLLED_STEPPER_DEFINITION_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_STEPPER_GENERIC_CONTROLLED_STEPPER_DEFINITION_HPP_INCLUDED
namespace boost {
namespace numeric {
namespace odeint {
template<
class ErrorStepper ,
class ErrorChecker ,
class Controller ,
class Resizer ,
class ErrorStepperCategory = typename ErrorStepper::stepper_category >
class generic_controlled_stepper ;
}
}
}
#endif // BOOST_NUMERIC_ODEINT_STEPPER_GENERIC_CONTROLLED_STEPPER_DEFINITION_HPP_INCLUDED

207
libs/numeric/odeint/ideas/controller/boost/numeric/odeint/stepper/generic_controlled_stepper_explicit.hpp
View File

@@ -1,207 +0,0 @@
/*
[auto_generated]
boost/numeric/odeint/stepper/generic_controlled_stepper_explicit.hpp
[begin_description]
Specialization of the generic_controlled_stepper for the explicit_error_stepper_tag. This class is for
runge_kutta54_cash_karp or runge_kutta78_ehlberg.
[end_description]
Copyright 2009-2011 Karsten Ahnert
Copyright 2009-2011 Mario Mulansky
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or
copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#ifndef BOOST_NUMERIC_ODEINT_STEPPER_GENERIC_CONTROLLED_STEPPER_EXPLICIT_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_STEPPER_GENERIC_CONTROLLED_STEPPER_EXPLICIT_HPP_INCLUDED
#include <boost/bind.hpp>
#include <boost/numeric/odeint/util/is_resizeable.hpp>
#include <boost/numeric/odeint/util/state_wrapper.hpp>
#include <boost/numeric/odeint/stepper/stepper_categories.hpp>
#include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
#include <boost/numeric/odeint/stepper/generic_controlled_stepper_definition.hpp>
namespace boost {
namespace numeric {
namespace odeint {
template< class ErrorStepper , class ErrorChecker , class Controller , class Resizer >
class generic_controlled_stepper< ErrorStepper , ErrorChecker , Controller , Resizer , explicit_error_stepper_tag >
{
public:
typedef ErrorStepper stepper_type;
typedef ErrorChecker error_checker_type;
typedef Controller controller_type;
typedef Resizer resizer_type;
typedef typename stepper_type::state_type state_type;
typedef typename stepper_type::value_type value_type;
typedef typename stepper_type::deriv_type deriv_type;
typedef typename stepper_type::time_type time_type;
typedef typename stepper_type::order_type order_type;
typedef explicit_controlled_stepper_tag stepper_category;
typedef state_wrapper< state_type > wrapped_state_type;
typedef state_wrapper< deriv_type > wrapped_deriv_type;
generic_controlled_stepper(
const stepper_type &stepper = stepper_type() ,
const error_checker_type &error_checker = error_checker_type() ,
const controller_type &controller = controller_type() )
: m_stepper( stepper ) , m_error_checker( error_checker ) , m_controller( controller ) ,
m_dxdt_resizer() , m_xerr_resizer() , m_xnew_resizer() ,
m_dxdt() , m_xerr() , m_xnew() { }
/*
* Version 1 : try_step( sys , x , t , dt )
*
* The overloads are needed to solve the forwarding problem
*/
template< class System , class StateInOut >
controlled_step_result try_step( System system , StateInOut &x , time_type &t , time_type &dt )
{
return try_step_v1( system , x , t, dt );
}
template< class System , class StateInOut >
controlled_step_result try_step( System system , const StateInOut &x , time_type &t , time_type &dt )
{
return try_step_v1( system , x , t, dt );
}
/*
* Version 2 : try_step( sys , x , dxdt , t , dt )
*
* this version does not solve the forwarding problem, boost.range can not be used
*/
template< class System , class StateInOut , class DerivIn >
controlled_step_result try_step( System system , StateInOut &x , const DerivIn &dxdt , time_type &t , time_type &dt )
{
m_xnew_resizer.adjust_size( x , boost::bind( &generic_controlled_stepper::template resize_m_xnew_impl< StateInOut > , boost::ref( *this ) , _1 ) );
controlled_step_result res = try_step( system , x , dxdt , t , m_xnew.m_v , dt );
if( res == success )
{
boost::numeric::odeint::copy( m_xnew.m_v , x );
}
return res;
}
/*
* Version 3 : try_step( sys , in , t , out , dt )
*
* this version does not solve the forwarding problem, boost.range can not be used
*/
template< class System , class StateIn , class StateOut >
controlled_step_result try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt )
{
typename boost::unwrap_reference< System >::type &sys = system;
m_dxdt_resizer.adjust_size( in , boost::bind( &generic_controlled_stepper::template resize_m_dxdt_impl< StateIn > , boost::ref( *this ) , _1 ) );
sys( in , m_dxdt.m_v , t );
return try_step( system , in , m_dxdt.m_v , t , out , dt );
}
/*
* Version 4 : try_step( sys , in , dxdt , t , out , dt )
*
* this version does not solve the forwarding problem, boost.range can not be used
*/
template< class System , class StateIn , class DerivIn , class StateOut >
controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , time_type &dt )
{
m_xerr_resizer.adjust_size( in , boost::bind( &generic_controlled_stepper::template resize_m_xerr_impl< StateIn > , boost::ref( *this ) , _1 ) );
m_stepper.do_step( system , in , dxdt , t , out , dt , m_xerr.m_v );
value_type error = m_error_checker.error( in , out , dxdt , m_xerr.m_v , dt );
return m_controller( error , t , dt , m_stepper.stepper_order() , m_stepper.error_order() );
}
template< class StateType >
void adjust_size( const StateType &x )
{
resize_m_xerr_impl( x );
resize_m_dxdt_impl( x );
resize_m_xnew_impl( x );
m_stepper.adjust_size( x );
}
stepper_type& stepper( void ) { return m_stepper; }
const stepper_type stepper( void ) const { return m_stepper; }
private:
template< class StateIn >
bool resize_m_xerr_impl( const StateIn &x )
{
return adjust_size_by_resizeability( m_xerr , x , typename is_resizeable< wrapped_state_type >::type() );
}
template< class StateIn >
bool resize_m_dxdt_impl( const StateIn &x )
{
return adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable< wrapped_deriv_type >::type() );
}
template< class StateIn >
bool resize_m_xnew_impl( const StateIn &x )
{
return adjust_size_by_resizeability( m_xnew , x , typename is_resizeable< wrapped_state_type >::type() );
}
template< class System , class StateInOut >
controlled_step_result try_step_v1( System system , StateInOut &x , time_type &t , time_type &dt )
{
typename boost::unwrap_reference< System >::type &sys = system;
m_dxdt_resizer.adjust_size( x , boost::bind( &generic_controlled_stepper::template resize_m_dxdt_impl< StateInOut > , boost::ref( *this ) , _1 ) );
sys( x , m_dxdt.m_v ,t );
return try_step( system , x , m_dxdt.m_v , t , dt );
}
stepper_type m_stepper;
error_checker_type m_error_checker;
controller_type m_controller;
resizer_type m_dxdt_resizer;
resizer_type m_xerr_resizer;
resizer_type m_xnew_resizer;
wrapped_deriv_type m_dxdt;
wrapped_state_type m_xerr;
wrapped_state_type m_xnew;
};
}
}
}
#endif // BOOST_NUMERIC_ODEINT_STEPPER_GENERIC_CONTROLLED_STEPPER_EXPLICIT_HPP_INCLUDED

454
libs/numeric/odeint/ideas/controller/boost/numeric/odeint/stepper/generic_controlled_stepper_explicit_fsal.hpp
View File

@@ -1,454 +0,0 @@
/*
[auto_generated]
boost/numeric/odeint/stepper/generic_controlled_stepper_explicit_fsal.hpp
[begin_description]
Specialization of the generic_controlled_stepper for the explicit_error_stepper_fsal_tag. This class is for
runge_kutta_dopri5.
[end_description]
Copyright 2009-2011 Karsten Ahnert
Copyright 2009-2011 Mario Mulansky
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or
copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#ifndef BOOST_NUMERIC_ODEINT_STEPPER_GENERIC_CONTROLLED_STEPPER_EXPLICIT_FSAL_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_STEPPER_GENERIC_CONTROLLED_STEPPER_EXPLICIT_FSAL_HPP_INCLUDED
#include <boost/bind.hpp>
#include <boost/numeric/odeint/util/is_resizeable.hpp>
#include <boost/numeric/odeint/util/state_wrapper.hpp>
#include <boost/numeric/odeint/stepper/stepper_categories.hpp>
#include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
#include <boost/numeric/odeint/stepper/generic_controlled_stepper_definition.hpp>
namespace boost {
namespace numeric {
namespace odeint {
template< class ErrorStepper , class ErrorChecker , class Controller , class Resizer >
class generic_controlled_stepper< ErrorStepper , ErrorChecker , Controller , Resizer , explicit_error_stepper_fsal_tag >
{
public:
typedef ErrorStepper stepper_type;
typedef ErrorChecker error_checker_type;
typedef Controller controller_type;
typedef Resizer resizer_type;
typedef typename stepper_type::state_type state_type;
typedef typename stepper_type::value_type value_type;
typedef typename stepper_type::deriv_type deriv_type;
typedef typename stepper_type::time_type time_type;
typedef typename stepper_type::order_type order_type;
typedef explicit_controlled_stepper_fsal_tag stepper_category;
typedef state_wrapper< state_type > wrapped_state_type;
typedef state_wrapper< deriv_type > wrapped_deriv_type;
generic_controlled_stepper(
const stepper_type &stepper = stepper_type() ,
const error_checker_type &error_checker = error_checker_type() ,
const controller_type &controller = controller_type() )
: m_stepper( stepper ) , m_error_checker( error_checker ) , m_controller( controller ) ,
m_dxdt_resizer() , m_xerr_resizer() , m_xnew_resizer() , m_dxdt_new_resizer() ,
m_xerr() ,
m_first_call( true )
{ }
/*
* Version 1 : try_step( sys , x , t , dt )
*
* The two overloads are needed in order to solve the forwarding problem
*/
template< class System , class StateInOut >
controlled_step_result try_step( System system , StateInOut &x , time_type &t , time_type &dt )
{
return try_step_v1( system , x , t , dt );
}
template< class System , class StateInOut >
controlled_step_result try_step( System system , const StateInOut &x , time_type &t , time_type &dt )
{
return try_step_v1( system , x , t , dt );
}
/*
* Version 2 : try_step( sys , in , t , out , dt );
*
* This version does not solve the forwarding problem, boost::range can not be used.
*/
template< class System , class StateIn , class StateOut >
controlled_step_result try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt )
{
if( m_dxdt_resizer.adjust_size( in , boost::bind( &generic_controlled_stepper::template resize_m_dxdt_impl< StateIn > , boost::ref( *this ) , _1 ) ) || m_first_call )
{
typename boost::unwrap_reference< System >::type &sys = system;
sys( in , m_dxdt.m_v ,t );
m_first_call = false;
}
return try_step( system , in , m_dxdt.m_v , t , out , dt );
}
/*
* Version 3 : try_step( sys , x , dxdt , t , dt )
*
* This version does not solve the forwarding problem, boost::range can not be used.
*/
template< class System , class StateInOut , class DerivInOut >
controlled_step_result try_step( System system , StateInOut &x , DerivInOut &dxdt , time_type &t , time_type &dt )
{
m_xnew_resizer.adjust_size( x , boost::bind( &generic_controlled_stepper::template resize_m_xnew_impl< StateInOut > , boost::ref( *this ) , _1 ) );
m_dxdt_new_resizer.adjust_size( x , boost::bind( &generic_controlled_stepper::template resize_m_dxdt_new_impl< StateInOut > , boost::ref( *this ) , _1 ) );
controlled_step_result res = try_step( system , x , dxdt , t , m_xnew.m_v , m_dxdt_new.m_v , dt );
if( res == success )
{
boost::numeric::odeint::copy( m_xnew.m_v , x );
boost::numeric::odeint::copy( m_dxdt_new.m_v , dxdt );
}
return res;
}
/*
* Version 4 : try_step( sys , in , dxdt , t , out , dxdtout , dt )
*
* This version does not solve the forwarding problem, boost::range can not be used.
*/
template< class System , class StateIn , class DerivIn , class StateOut , class DerivOut >
controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt_in , time_type &t ,
StateOut &out , DerivOut &dxdt_out , time_type &dt )
{
m_xerr_resizer.adjust_size( in , boost::bind( &generic_controlled_stepper::template resize_m_xerr_impl< StateIn > , boost::ref( *this ) , _1 ) );
m_stepper.do_step( system , in , dxdt_in , t , out , dxdt_out , dt , m_xerr.m_v );
value_type error = m_error_checker.error( in , out , dxdt_in , m_xerr.m_v , dt );
return m_controller( error , t , dt , m_stepper.stepper_order() , m_stepper.error_order() );
}
template< class StateType >
void adjust_size( const StateType &x )
{
resize_m_xerr_impl( x );
resize_m_dxdt_impl( x );
resize_m_xnew_impl( x );
resize_m_dxdt_new_impl( x );
m_stepper.adjust_size( x );
}
stepper_type& stepper( void ) { return m_stepper; }
const stepper_type stepper( void ) const { return m_stepper; }
private:
template< class System , class StateInOut >
controlled_step_result try_step_v1( System system , StateInOut &x , time_type &t , time_type &dt )
{
if( m_dxdt_resizer.adjust_size( x , boost::bind( &generic_controlled_stepper::template resize_m_dxdt_impl< StateInOut > , boost::ref( *this ) , _1 ) ) || m_first_call )
{
typename boost::unwrap_reference< System >::type &sys = system;
sys( x , m_dxdt.m_v , t );
m_first_call = false;
}
return try_step( system , x , m_dxdt.m_v , t , dt );
}
template< class StateIn >
bool resize_m_xerr_impl( const StateIn &x )
{
return adjust_size_by_resizeability( m_xerr , x , typename is_resizeable< wrapped_state_type >::type() );
}
template< class StateIn >
bool resize_m_dxdt_impl( const StateIn &x )
{
return adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable< wrapped_deriv_type >::type() );
}
template< class StateIn >
bool resize_m_xnew_impl( const StateIn &x )
{
return adjust_size_by_resizeability( m_xnew , x , typename is_resizeable< wrapped_state_type >::type() );
}
template< class StateIn >
bool resize_m_dxdt_new_impl( const StateIn &x )
{
return adjust_size_by_resizeability( m_dxdt_new , x , typename is_resizeable< wrapped_deriv_type >::type() );
}
stepper_type m_stepper;
error_checker_type m_error_checker;
controller_type m_controller;
resizer_type m_dxdt_resizer;
resizer_type m_xerr_resizer;
resizer_type m_xnew_resizer;
resizer_type m_dxdt_new_resizer;
wrapped_deriv_type m_dxdt;
wrapped_deriv_type m_dxdt_new;
wrapped_state_type m_xerr;
wrapped_state_type m_xnew;
bool m_first_call;
};
//template< class ErrorStepper , class ErrorChecker , class Resizer >
//class controlled_runge_kutta< ErrorStepper , ErrorChecker , Resizer , explicit_error_stepper_fsal_tag >
//{
//
//public:
//
// typedef ErrorStepper stepper_type;
// typedef typename stepper_type::state_type state_type;
// typedef typename stepper_type::value_type value_type;
// typedef typename stepper_type::deriv_type deriv_type;
// typedef typename stepper_type::time_type time_type;
// typedef typename stepper_type::order_type order_type;
// typedef typename stepper_type::algebra_type algebra_type;
// typedef typename stepper_type::operations_type operations_type;
// typedef Resizer resizer_type;
// typedef ErrorChecker error_checker_type;
// typedef explicit_controlled_stepper_fsal_tag stepper_category;
// typedef typename stepper_type::wrapped_state_type wrapped_state_type;
// typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type;
//
// typedef controlled_runge_kutta< ErrorStepper , ErrorChecker , Resizer , explicit_error_stepper_tag > controlled_stepper_type;
//
// controlled_runge_kutta(
// const error_checker_type &error_checker = error_checker_type() ,
// const stepper_type &stepper = stepper_type()
// )
// : m_stepper( stepper ) , m_error_checker( error_checker ) ,
// m_first_call( true )
// { }
//
// /*
// * Version 1 : try_step( sys , x , t , dt )
// *
// * The two overloads are needed in order to solve the forwarding problem
// */
// template< class System , class StateInOut >
// controlled_step_result try_step( System system , StateInOut &x , time_type &t , time_type &dt )
// {
// return try_step_v1( system , x , t , dt );
// }
//
// template< class System , class StateInOut >
// controlled_step_result try_step( System system , const StateInOut &x , time_type &t , time_type &dt )
// {
// return try_step_v1( system , x , t , dt );
// }
//
//
//
// /*
// * Version 2 : try_step( sys , in , t , out , dt );
// *
// * This version does not solve the forwarding problem, boost::range can not be used.
// */
// template< class System , class StateIn , class StateOut >
// controlled_step_result try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt )
// {
// if( m_dxdt_resizer.adjust_size( in , boost::bind( &controlled_runge_kutta::template resize_m_dxdt_impl< StateIn > , boost::ref( *this ) , _1 ) ) || m_first_call )
// {
// typename boost::unwrap_reference< System >::type &sys = system;
// sys( in , m_dxdt.m_v ,t );
// m_first_call = false;
// }
// return try_step( system , in , m_dxdt.m_v , t , out , dt );
// }
//
//
// /*
// * Version 3 : try_step( sys , x , dxdt , t , dt )
// *
// * This version does not solve the forwarding problem, boost::range can not be used.
// */
// template< class System , class StateInOut , class DerivInOut >
// controlled_step_result try_step( System system , StateInOut &x , DerivInOut &dxdt , time_type &t , time_type &dt )
// {
// m_xnew_resizer.adjust_size( x , boost::bind( &controlled_runge_kutta::template resize_m_xnew_impl< StateInOut > , boost::ref( *this ) , _1 ) );
// m_dxdt_new_resizer.adjust_size( x , boost::bind( &controlled_runge_kutta::template resize_m_dxdt_new_impl< StateInOut > , boost::ref( *this ) , _1 ) );
// controlled_step_result res = try_step( system , x , dxdt , t , m_xnew.m_v , m_dxdtnew.m_v , dt );
// if( res == success )
// {
// boost::numeric::odeint::copy( m_xnew.m_v , x );
// boost::numeric::odeint::copy( m_dxdtnew.m_v , dxdt );
// }
// return res;
// }
//
//
// /*
// * Version 3 : try_step( sys , in , dxdt , t , out , dt )
// *
// * This version does not solve the forwarding problem, boost::range can not be used.
// */
// template< class System , class StateIn , class DerivIn , class StateOut , class DerivOut >
// controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt_in , time_type &t ,
// StateOut &out , DerivOut &dxdt_out , time_type &dt )
// {
// using std::max;
// using std::min;
// using std::pow;
//
// m_xerr_resizer.adjust_size( in , boost::bind( &controlled_runge_kutta::template resize_m_xerr_impl< StateIn > , boost::ref( *this ) , _1 ) );
//
// //fsal: m_stepper.get_dxdt( dxdt );
// //fsal: m_stepper.do_step( sys , x , dxdt , t , dt , m_x_err );
// m_stepper.do_step( system , in , dxdt_in , t , out , dxdt_out , dt , m_xerr.m_v );
//
// // this potentially overwrites m_x_err! (standard_error_checker does, at least)
// value_type max_rel_err = m_error_checker.error( m_stepper.algebra() , in , dxdt_in , m_xerr.m_v , 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 );
// return fail;
// }
// else
// {
// if( max_rel_err < 0.5 )
// {
// //error too small - increase dt and keep the evolution and limit scaling factor to 5.0
// t += dt;
// dt *= min( 0.9 * pow( max_rel_err , -1.0 / m_stepper.stepper_order() ) , 5.0 );
// return success;
// }
// else
// {
// t += dt;
// return success;
// }
// }
// }
//
//
//
//
// template< class StateType >
// void adjust_size( const StateType &x )
// {
// resize_m_xerr_impl( x );
// resize_m_dxdt_impl( x );
// resize_m_dxdt_new_impl( x );
// resize_m_xnew_impl( x );
// }
//
//
// stepper_type& stepper( void )
// {
// return m_stepper;
// }
//
// const stepper_type& stepper( void ) const
// {
// return m_stepper;
// }
//
//
//
//private:
//
//
// template< class StateIn >
// bool resize_m_xerr_impl( const StateIn &x )
// {
// return adjust_size_by_resizeability( m_xerr , x , typename is_resizeable<state_type>::type() );
// }
//
// template< class StateIn >
// bool resize_m_dxdt_impl( const StateIn &x )
// {
// return adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() );
// }
//
// template< class StateIn >
// bool resize_m_dxdt_new_impl( const StateIn &x )
// {
// return adjust_size_by_resizeability( m_dxdtnew , x , typename is_resizeable<deriv_type>::type() );
// }
//
// template< class StateIn >
// bool resize_m_xnew_impl( const StateIn &x )
// {
// return adjust_size_by_resizeability( m_xnew , x , typename is_resizeable<state_type>::type() );
// }
//
//
// template< class System , class StateInOut >
// controlled_step_result try_step_v1( System system , StateInOut &x , time_type &t , time_type &dt )
// {
// if( m_dxdt_resizer.adjust_size( x , boost::bind( &controlled_runge_kutta::template resize_m_dxdt_impl< StateInOut > , boost::ref( *this ) , _1 ) ) || m_first_call )
// {
// typename boost::unwrap_reference< System >::type &sys = system;
// sys( x , m_dxdt.m_v , t );
// m_first_call = false;
// }
// return try_step( system , x , m_dxdt.m_v , t , dt );
// }
//
//
// stepper_type m_stepper;
// error_checker_type m_error_checker;
//
// resizer_type m_dxdt_resizer;
// resizer_type m_xerr_resizer;
// resizer_type m_xnew_resizer;
// resizer_type m_dxdt_new_resizer;
//
// wrapped_deriv_type m_dxdt;
// wrapped_state_type m_xerr;
// wrapped_state_type m_xnew;
// wrapped_deriv_type m_dxdtnew;
// bool m_first_call;
//};
}
}
}
#endif // BOOST_NUMERIC_ODEINT_STEPPER_GENERIC_CONTROLLED_STEPPER_EXPLICIT_FSAL_HPP_INCLUDED

208
libs/numeric/odeint/ideas/controller/boost/numeric/odeint/stepper/generic_dense_output_stepper.hpp
View File

@@ -1,208 +0,0 @@
/*
[auto_generated]
boost/numeric/odeint/stepper/generic_dense_output_stepper.hpp
[begin_description]
Implementation of a generic dense output stepper
[end_description]
Copyright 2009-2011 Karsten Ahnert
Copyright 2009-2011 Mario Mulansky
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or
copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#ifndef BOOST_NUMERIC_ODEINT_STEPPER_GENERIC_DENSE_OUTPUT_STEPPER_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_STEPPER_GENERIC_DENSE_OUTPUT_STEPPER_HPP_INCLUDED
#include <boost/numeric/odeint/util/state_wrapper.hpp>
#include <boost/numeric/odeint/stepper/stepper_categories.hpp>
#include <boost/numeric/odeint/stepper/generic_dense_output_stepper_definition.hpp>
#include <boost/bind.hpp>
namespace boost {
namespace numeric {
namespace odeint {
template< class ControlledStepper >
class generic_dense_output_stepper< ControlledStepper , controlled_stepper_tag >
{
public:
void copy_variables( const generic_dense_output_stepper &stepper )
{
m_stepper = stepper.m_stepper;
m_x1 = stepper.m_x1;
m_x2 = stepper.m_x2;
if( stepper.m_current_state == ( & ( stepper.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 = stepper.m_t;
m_t_old = stepper.m_t_old;
m_dt = stepper.m_dt;
}
public:
typedef ControlledStepper controlled_stepper_type;
typedef typename controlled_stepper_type::stepper_type 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 typename stepper_type::resizer_type resizer_type;
typedef state_wrapper< state_type > wrapped_state_type;
typedef state_wrapper< deriv_type > wrapped_deriv_type;
typedef dense_output_stepper_tag stepper_category;
typedef generic_dense_output_stepper< ControlledStepper > dense_output_stepper_type;
generic_dense_output_stepper( const controlled_stepper_type &stepper = controlled_stepper_type() )
: m_stepper( stepper ) ,
m_x1() , m_x2() , m_current_state( &m_x1.m_v ) , m_old_state( &m_x2.m_v ) ,
m_t() , m_t_old() , m_dt()
{ }
generic_dense_output_stepper( const generic_dense_output_stepper &rb )
: m_current_state( &m_x1.m_v ) , m_old_state( &m_x2.m_v )
{ }
generic_dense_output_stepper& operator=( const generic_dense_output_stepper &rb )
{
copy_variables( rb );
return *this;
}
template< class StateType >
void initialize( const StateType &x0 , const time_type &t0 , const time_type &dt0 )
{
m_resizer.adjust_size( x0 , boost::bind( &dense_output_stepper_type::template resize_impl< StateType > , boost::ref( *this ) , _1 ) );
*m_current_state = x0;
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 = fail;
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( "generic_dense_output_stepper : too much iterations!");
}
while( res == fail );
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 );
}
template< class StateType >
void adjust_size( const StateType &x )
{
m_stepper.adjust_size( x );
resize_impl( x );
}
const state_type& current_state( void ) const
{
return *m_current_state;
}
const time_type& current_time( void ) const
{
return m_t;
}
const time_type& previous_state( void ) const
{
return *m_old_state;
}
const time_type& previous_time( void ) const
{
return m_t_old;
}
const time_type& current_time_step( void ) const
{
return m_dt;
}
private:
template< class StateIn >
bool resize_impl( const StateIn &x )
{
bool resized = false;
resized |= adjust_size_by_resizeability( m_x1 , x , typename is_resizeable<state_type>::type() );
resized |= adjust_size_by_resizeability( m_x2 , x , typename is_resizeable<state_type>::type() );
return resized;
}
controlled_stepper_type m_stepper;
resizer_type m_resizer;
wrapped_state_type m_x1 , m_x2;
state_type *m_current_state , *m_old_state;
time_type m_t , m_t_old , m_dt;
};
}
}
}
#endif // BOOST_NUMERIC_ODEINT_STEPPER_GENERIC_DENSE_OUTPUT_STEPPER_HPP_INCLUDED

39
libs/numeric/odeint/ideas/controller/boost/numeric/odeint/stepper/generic_dense_output_stepper_definition.hpp
View File

@@ -1,39 +0,0 @@
/*
[auto_generated]
boost/numeric/odeint/stepper/generic_dense_output_stepper_definition.hpp
[begin_description]
Definition of the generic dense output stepper encapsulating a stepper capable of calculating intermediate steps.
[end_description]
Copyright 2009-2011 Karsten Ahnert
Copyright 2009-2011 Mario Mulansky
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or
copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#ifndef BOOST_NUMERIC_ODEINT_STEPPER_GENERIC_DENSE_OUTPUT_STEPPER_DEFINITION_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_STEPPER_GENERIC_DENSE_OUTPUT_STEPPER_DEFINITION_HPP_INCLUDED
namespace boost {
namespace numeric {
namespace odeint {
template<
class Stepper ,
class StepperCategory = typename Stepper::stepper_category >
class generic_dense_output_stepper ;
}
}
}
#endif // BOOST_NUMERIC_ODEINT_STEPPER_GENERIC_DENSE_OUTPUT_STEPPER_DEFINITION_HPP_INCLUDED

55
libs/numeric/odeint/ideas/controller/boost/numeric/odeint/util/number_of_elements.hpp
View File

@@ -1,55 +0,0 @@
/*
[auto_generated]
boost/numeric/odeint/util/number_of_elements.hpp
[begin_description]
Generic number of elements (size) function.
[end_description]
Copyright 2009-2011 Karsten Ahnert
Copyright 2009-2011 Mario Mulansky
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or
copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#ifndef BOOST_NUMERIC_ODEINT_UTIL_NUMBER_OF_ELEMENTS_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_UTIL_NUMBER_OF_ELEMENTS_HPP_INCLUDED
namespace boost {
namespace numeric {
namespace odeint {
template< class T >
struct number_of_elements_impl
{
typedef typename T::size_type size_type;
typedef size_type result_type;
static result_type number_of_elements( const T &x )
{
return x.size();
}
};
template< class T >
typename number_of_elements_impl< T >::result_type number_of_elements( const T &x )
{
return number_of_elements_impl< T >::number_of_elements( x );
}
} // namespace odeint
} // namespace numeric
} // namespace boost
#endif // BOOST_NUMERIC_ODEINT_UTIL_NUMBER_OF_ELEMENTS_HPP_INCLUDED

82
libs/numeric/odeint/ideas/controller/boost/numeric/odeint/util/ref_or_value_holder.hpp
View File

@@ -1,82 +0,0 @@
/*
[auto_generated]
boost/numeric/odeint/util/ref_or_value_holder.hpp
[begin_description]
Util class for holding either a value or a reference.
[end_description]
Copyright 2009-2011 Karsten Ahnert
Copyright 2009-2011 Mario Mulansky
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or
copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#ifndef BOOST_NUMERIC_ODEINT_UTIL_REF_OR_VALUE_HOLDER_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_UTIL_REF_OR_VALUE_HOLDER_HPP_INCLUDED
namespace boost {
namespace numeric {
namespace odeint {
/*
* ToDo : find nice names for the types in ref_or_value_holder
*/
template< class T , bool RefOrVariable > struct ref_or_value_holder ;
// template< class T , bool RefOrVariable > struct cref_or_value_holder ;
template< class T >
struct ref_or_value_holder< T , true >
{
typedef T value_type;
typedef T& holder_type;
typedef T& constructor_type;
typedef T& reference_type;
holder_type m_t;
ref_or_value_holder( constructor_type t ) : m_t( t ) { }
reference_type get( void ) { return m_t; }
};
template< class T >
struct ref_or_value_holder< T , false >
{
typedef T value_type;
typedef T holder_type;
typedef const T& constructor_type;
typedef T& reference_type;
holder_type m_t;
ref_or_value_holder( constructor_type t ) : m_t( t ) { }
reference_type get( void ) { return m_t; }
};
//template< class T >
//struct cref_or_value_holder< T , true >
//{
// const T &m_t;
// cref_or_value_holder( const T &t ) : m_t( t ) { }
// const T& get( void ) { return m_t; }
//};
//
//template< class T >
//struct cref_or_value_holder< T , false >
//{
// const T m_t;
// cref_or_value_holder( const T &t ) : m_t( t ) { }
// const T& get( void ) { return m_t; }
//};
} // namespace odeint
} // namespace numeric
} // namespace boost
#endif // BOOST_NUMERIC_ODEINT_UTIL_REF_OR_VALUE_HOLDER_HPP_INCLUDED

142
libs/numeric/odeint/ideas/controller/controller.cpp
View File

@@ -1,142 +0,0 @@
#include <iostream>
#include <cmath>
#include <boost/timer.hpp>
#include <boost/numeric/odeint/algebra/range_algebra.hpp>
#include <boost/numeric/odeint/algebra/default_operations.hpp>
#include <boost/numeric/odeint/util/state_wrapper.hpp>
#include <boost/numeric/odeint/stepper/stepper_categories.hpp>
#include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
#include <boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp>
#include <boost/numeric/odeint/stepper/controlled_runge_kutta.hpp>
#include <boost/numeric/odeint/integrate/integrate_adaptive.hpp>
#include <boost/numeric/odeint/util/resizer.hpp>
#include <boost/numeric/odeint/stepper/error_checker_explicit.hpp>
#include <boost/numeric/odeint/stepper/generic_controlled_stepper.hpp>
#include <boost/numeric/odeint/stepper/generic_controlled_stepper_explicit.hpp>
#include <boost/numeric/odeint/stepper/controller/default_controller.hpp>
#include <boost/numeric/odeint/stepper/controller/pi_controller.hpp>
#define tab "\t"
namespace boost {
namespace numeric {
namespace odeint {
/*
* move this class into controller_runge_kutta.hpp
*/
template< class Stepper >
class controlled_runge_kutta_explicit :
public generic_controlled_stepper<
Stepper ,
error_checker_explicit_new< typename Stepper::value_type , typename Stepper::algebra_type , typename Stepper::operations_type > ,
default_controller ,
typename Stepper::resizer_type ,
typename Stepper::stepper_category >
{
public:
typedef generic_controlled_stepper<
Stepper ,
error_checker_explicit_new< typename Stepper::value_type , typename Stepper::algebra_type , typename Stepper::operations_type > ,
default_controller ,
typename Stepper::resizer_type ,
typename Stepper::stepper_category > base_type;
typedef error_checker_explicit_new< typename Stepper::value_type , typename Stepper::algebra_type , typename Stepper::operations_type > error_checker_type;
typedef typename Stepper::value_type value_type;
controlled_runge_kutta_explicit( const Stepper &stepper = Stepper() ,
const value_type eps_abs = static_cast< value_type >( 1.0e-6 ) ,
const value_type eps_rel = static_cast< value_type >( 1.0e-6 ) ,
const value_type a_x = static_cast< value_type >( 1.0 ) ,
const value_type a_dxdt = static_cast< value_type >( 1.0 ) )
: base_type( stepper ,
error_checker_type( base_type::stepper().algebra() , eps_abs , eps_rel , a_x , a_dxdt )
) { }
};
}
}
}
typedef boost::array< double , 3 > state_type;
void lorenz( const state_type &x , state_type &dxdt , double t )
{
const double sigma = 10.0;
const double R = 28.0;
const double b = 8.0 / 3.0;
dxdt[0] = sigma * ( x[1] - x[0] );
dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
dxdt[2] = -b * x[2] + x[0] * x[1];
}
int main( int argc , char **argv )
{
using namespace std;
using namespace boost::numeric::odeint;
generic_controlled_stepper<
runge_kutta_cash_karp54< state_type > ,
error_checker_explicit< double , range_algebra , default_operations > ,
default_controller ,
initially_resizer ,
explicit_error_stepper_tag > stepper1;
controlled_runge_kutta< runge_kutta_cash_karp54< state_type > > stepper2;
controlled_runge_kutta_explicit< runge_kutta_cash_karp54< state_type > > stepper3;
generic_controlled_stepper<
runge_kutta_cash_karp54< state_type > ,
error_checker_explicit< double , range_algebra , default_operations > ,
pi_controller< double > ,
initially_resizer ,
explicit_error_stepper_tag > stepper4;
state_type x1 = {{ 10.0 , 10.0 , 10.0 }};
state_type x2 = x1;
state_type x3 = x1;
state_type x4 = x1;
boost::timer timer;
const double t_max = 100000.0;
timer.restart();
size_t steps1 = integrate_adaptive( stepper1 , lorenz , x1 , 0.0 , t_max , 0.1 );
double t1 = timer.elapsed();
cout << steps1 << tab << t1 << tab << x1[0] << tab << x1[1] << tab << x1[2] << endl;
timer.restart();
size_t steps2 = integrate_adaptive( stepper2 , lorenz , x2 , 0.0 , t_max , 0.1 );
double t2 = timer.elapsed();
cout << steps2 << tab << t2 << tab << x2[0] << tab << x2[1] << tab << x2[2] << endl;
timer.restart();
size_t steps3 = integrate_adaptive( stepper3 , lorenz , x3 , 0.0 , t_max , 0.1 );
double t3 = timer.elapsed();
cout << steps3 << tab << t3 << tab << x3[0] << tab << x3[1] << tab << x3[2] << endl;
timer.restart();
size_t steps4 = integrate_adaptive( stepper4 , lorenz , x4 , 0.0 , t_max , 0.1 );
double t4 = timer.elapsed();
cout << steps4 << tab << t4 << tab << x4[0] << tab << x4[1] << tab << x4[2] << endl;
return 0;
}

168
libs/numeric/odeint/ideas/controller/dense_out_dopri5_controller.cpp
View File

@@ -1,168 +0,0 @@
#include <iostream>
#include <cmath>
#include <boost/timer.hpp>
#include <boost/numeric/odeint/algebra/range_algebra.hpp>
#include <boost/numeric/odeint/algebra/default_operations.hpp>
#include <boost/numeric/odeint/util/state_wrapper.hpp>
#include <boost/numeric/odeint/stepper/stepper_categories.hpp>
#include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
#include <boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp>
#include <boost/numeric/odeint/stepper/controlled_runge_kutta.hpp>
#include <boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp>
#include <boost/numeric/odeint/integrate/integrate_adaptive.hpp>
#include <boost/numeric/odeint/integrate/integrate_const.hpp>
#include <boost/numeric/odeint/util/resizer.hpp>
#include <boost/numeric/odeint/stepper/error_checker_explicit.hpp>
#include <boost/numeric/odeint/stepper/error_checker_max_norm.hpp>
#include <boost/numeric/odeint/stepper/generic_controlled_stepper_explicit_fsal.hpp>
#include <boost/numeric/odeint/stepper/controller/default_controller.hpp>
#include <boost/numeric/odeint/stepper/controller/pi_controller.hpp>
#define tab "\t"
namespace boost {
namespace numeric {
namespace odeint {
/*
* move this class into controller_runge_kutta.hpp
*/
template< class Stepper >
class controlled_runge_kutta_explicit :
public generic_controlled_stepper<
Stepper ,
error_checker_explicit_new< typename Stepper::value_type , typename Stepper::algebra_type , typename Stepper::operations_type > ,
default_controller ,
typename Stepper::resizer_type ,
typename Stepper::stepper_category >
{
public:
typedef generic_controlled_stepper<
Stepper ,
error_checker_explicit_new< typename Stepper::value_type , typename Stepper::algebra_type , typename Stepper::operations_type > ,
default_controller ,
typename Stepper::resizer_type ,
typename Stepper::stepper_category > base_type;
typedef error_checker_explicit_new< typename Stepper::value_type , typename Stepper::algebra_type , typename Stepper::operations_type > error_checker_type;
typedef typename Stepper::value_type value_type;
controlled_runge_kutta_explicit( const Stepper &stepper = Stepper() ,
const value_type eps_abs = static_cast< value_type >( 1.0e-6 ) ,
const value_type eps_rel = static_cast< value_type >( 1.0e-6 ) ,
const value_type a_x = static_cast< value_type >( 1.0 ) ,
const value_type a_dxdt = static_cast< value_type >( 1.0 ) )
: base_type( stepper ,
error_checker_type( base_type::stepper().algebra() , eps_abs , eps_rel , a_x , a_dxdt )
) { }
};
}
}
}
typedef boost::array< double , 3 > state_type;
void lorenz( const state_type &x , state_type &dxdt , double t )
{
const double sigma = 10.0;
const double R = 28.0;
const double b = 8.0 / 3.0;
dxdt[0] = sigma * ( x[1] - x[0] );
dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
dxdt[2] = -b * x[2] + x[0] * x[1];
}
int main( int argc , char **argv )
{
using namespace std;
using namespace boost::numeric::odeint;
dense_output_runge_kutta<
generic_controlled_stepper<
runge_kutta_dopri5< state_type > ,
error_checker_explicit< double , range_algebra , default_operations > ,
default_controller ,
initially_resizer ,
explicit_error_stepper_fsal_tag > > stepper1;
dense_output_runge_kutta< controlled_runge_kutta< runge_kutta_dopri5< state_type > > > stepper2;
dense_output_runge_kutta< controlled_runge_kutta_explicit< runge_kutta_dopri5< state_type > > > stepper3;
dense_output_runge_kutta<
generic_controlled_stepper<
runge_kutta_dopri5< state_type > ,
error_checker_explicit< double , range_algebra , default_operations > ,
pi_controller< double > ,
initially_resizer ,
explicit_error_stepper_fsal_tag > > stepper4;
typedef generic_controlled_stepper<
runge_kutta_dopri5< state_type > ,
error_checker_max_norm< double , range_algebra , default_operations > ,
default_controller ,
initially_resizer ,
explicit_error_stepper_fsal_tag > controller5_type;
typedef dense_output_runge_kutta< controller5_type > stepper5_type;
range_algebra al;
stepper5_type stepper5 = stepper5_type(
controller5_type( runge_kutta_dopri5< state_type >() ,
error_checker_max_norm< double , range_algebra , default_operations >( al ) ) );
state_type x1 = {{ 10.0 , 10.0 , 10.0 }};
state_type x2 = x1;
state_type x3 = x1;
state_type x4 = x1;
state_type x5 = x1;
boost::timer timer;
const double t_max = 1000000.0;
timer.restart();
size_t steps1 = integrate_const( stepper1 , lorenz , x1 , 0.0 , t_max , 10.0 );
double t1 = timer.elapsed();
cout << steps1 << tab << t1 << tab << x1[0] << tab << x1[1] << tab << x1[2] << endl;
timer.restart();
size_t steps2 = integrate_const( stepper2 , lorenz , x2 , 0.0 , t_max , 10.0 );
double t2 = timer.elapsed();
cout << steps2 << tab << t2 << tab << x2[0] << tab << x2[1] << tab << x2[2] << endl;
timer.restart();
size_t steps3 = integrate_const( stepper3 , lorenz , x3 , 0.0 , t_max , 10.0 );
double t3 = timer.elapsed();
cout << steps3 << tab << t3 << tab << x3[0] << tab << x3[1] << tab << x3[2] << endl;
timer.restart();
size_t steps4 = integrate_const( stepper4 , lorenz , x4 , 0.0 , t_max , 10.0 );
double t4 = timer.elapsed();
cout << steps4 << tab << t4 << tab << x4[0] << tab << x4[1] << tab << x4[2] << endl;
timer.restart();
size_t steps5 = integrate_const( stepper5 , lorenz , x5 , 0.0 , t_max , 10.0 );
double t5 = timer.elapsed();
cout << steps5 << tab << t5 << tab << x5[0] << tab << x5[1] << tab << x5[2] << endl;
return 0;
}

142
libs/numeric/odeint/ideas/controller/dopri5_controller.cpp
View File

@@ -1,142 +0,0 @@
#include <iostream>
#include <cmath>
#include <boost/timer.hpp>
#include <boost/numeric/odeint/algebra/range_algebra.hpp>
#include <boost/numeric/odeint/algebra/default_operations.hpp>
#include <boost/numeric/odeint/util/state_wrapper.hpp>
#include <boost/numeric/odeint/stepper/stepper_categories.hpp>
#include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
#include <boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp>
#include <boost/numeric/odeint/stepper/controlled_runge_kutta.hpp>
#include <boost/numeric/odeint/integrate/integrate_adaptive.hpp>
#include <boost/numeric/odeint/util/resizer.hpp>
#include <boost/numeric/odeint/stepper/error_checker_explicit.hpp>
#include <boost/numeric/odeint/stepper/generic_controlled_stepper_explicit_fsal.hpp>
#include <boost/numeric/odeint/stepper/controller/default_controller.hpp>
#include <boost/numeric/odeint/stepper/controller/pi_controller.hpp>
#define tab "\t"
namespace boost {
namespace numeric {
namespace odeint {
/*
* move this class into controller_runge_kutta.hpp
*/
template< class Stepper >
class controlled_runge_kutta_explicit :
public generic_controlled_stepper<
Stepper ,
error_checker_explicit_new< typename Stepper::value_type , typename Stepper::algebra_type , typename Stepper::operations_type > ,
default_controller ,
typename Stepper::resizer_type ,
typename Stepper::stepper_category >
{
public:
typedef generic_controlled_stepper<
Stepper ,
error_checker_explicit_new< typename Stepper::value_type , typename Stepper::algebra_type , typename Stepper::operations_type > ,
default_controller ,
typename Stepper::resizer_type ,
typename Stepper::stepper_category > base_type;
typedef error_checker_explicit_new< typename Stepper::value_type , typename Stepper::algebra_type , typename Stepper::operations_type > error_checker_type;
typedef typename Stepper::value_type value_type;
controlled_runge_kutta_explicit( const Stepper &stepper = Stepper() ,
const value_type eps_abs = static_cast< value_type >( 1.0e-6 ) ,
const value_type eps_rel = static_cast< value_type >( 1.0e-6 ) ,
const value_type a_x = static_cast< value_type >( 1.0 ) ,
const value_type a_dxdt = static_cast< value_type >( 1.0 ) )
: base_type( stepper ,
error_checker_type( base_type::stepper().algebra() , eps_abs , eps_rel , a_x , a_dxdt )
) { }
};
}
}
}
typedef boost::array< double , 3 > state_type;
void lorenz( const state_type &x , state_type &dxdt , double t )
{
const double sigma = 10.0;
const double R = 28.0;
const double b = 8.0 / 3.0;
dxdt[0] = sigma * ( x[1] - x[0] );
dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
dxdt[2] = -b * x[2] + x[0] * x[1];
}
int main( int argc , char **argv )
{
using namespace std;
using namespace boost::numeric::odeint;
generic_controlled_stepper<
runge_kutta_dopri5< state_type > ,
error_checker_explicit< double , range_algebra , default_operations > ,
default_controller ,
initially_resizer ,
explicit_error_stepper_fsal_tag > stepper1;
controlled_runge_kutta< runge_kutta_dopri5< state_type > > stepper2;
controlled_runge_kutta_explicit< runge_kutta_dopri5< state_type > > stepper3;
generic_controlled_stepper<
runge_kutta_dopri5< state_type > ,
error_checker_explicit< double , range_algebra , default_operations > ,
pi_controller< double > ,
initially_resizer ,
explicit_error_stepper_fsal_tag > stepper4;
state_type x1 = {{ 10.0 , 10.0 , 10.0 }};
state_type x2 = x1;
state_type x3 = x1;
state_type x4 = x1;
boost::timer timer;
const double t_max = 100000.0;
timer.restart();
size_t steps1 = integrate_adaptive( stepper1 , lorenz , x1 , 0.0 , t_max , 0.1 );
double t1 = timer.elapsed();
cout << steps1 << tab << t1 << tab << x1[0] << tab << x1[1] << tab << x1[2] << endl;
timer.restart();
size_t steps2 = integrate_adaptive( stepper2 , lorenz , x2 , 0.0 , t_max , 0.1 );
double t2 = timer.elapsed();
cout << steps2 << tab << t2 << tab << x2[0] << tab << x2[1] << tab << x2[2] << endl;
timer.restart();
size_t steps3 = integrate_adaptive( stepper3 , lorenz , x3 , 0.0 , t_max , 0.1 );
double t3 = timer.elapsed();
cout << steps3 << tab << t3 << tab << x3[0] << tab << x3[1] << tab << x3[2] << endl;
timer.restart();
size_t steps4 = integrate_adaptive( stepper4 , lorenz , x4 , 0.0 , t_max , 0.1 );
double t4 = timer.elapsed();
cout << steps4 << tab << t4 << tab << x4[0] << tab << x4[1] << tab << x4[2] << endl;
return 0;
}

27
libs/numeric/odeint/ideas/controller/libs/numeric/odeint/test/Jamfile
View File

@@ -1,27 +0,0 @@
# (C) Copyright 2010 : Karsten Ahnert, Mario Mulansky
# Distributed under the Boost Software License, Version 1.0.
# (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
# bring in rules for testing
import testing ;
use-project boost : $(BOOST_ROOT) ;
project
: requirements
<library>/boost/test//boost_unit_test_framework
<define>BOOST_ALL_NO_LIB=1
<include>../../../..
<link>static
# <cxxflags>-D_SCL_SECURE_NO_WARNINGS
;
test-suite "odeint"
:
[ run ref_or_value_holder.cpp ]
[ run error_checker_max_norm.cpp ]
[ run number_of_elements.cpp ]
[ run error_checker_l2_norm.cpp ]
: <testing.launcher>valgrind
;

110
libs/numeric/odeint/ideas/controller/libs/numeric/odeint/test/error_checker_l2_norm.cpp
View File

@@ -1,110 +0,0 @@
/*
* is_pair.cpp
*
* Created on: Feb 12, 2011
* Author: karsten
*/
#define BOOST_TEST_MODULE odeint_error_checker_max_norm
#include <utility>
#include <boost/test/unit_test.hpp>
#include <boost/static_assert.hpp>
#include <boost/numeric/odeint/algebra/range_algebra.hpp>
#include <boost/numeric/odeint/algebra/default_operations.hpp>
#include <boost/numeric/odeint/stepper/error_checker_l2_norm.hpp>
#include <boost/array.hpp>
using namespace boost::numeric::odeint;
#include <iostream>
using namespace std;
BOOST_AUTO_TEST_SUITE( error_checker_max_norm_test )
BOOST_AUTO_TEST_CASE( test_version1 )
{
// computes err = max( x_err / ( eps_abs + eps_rel * max( x , x_old ) ) );
range_algebra algebra;
error_checker_l2_norm< double , range_algebra , default_operations > error_checker( algebra , 1.0e-3 , 1.0e-4 );
double dt , err;
boost::array< double , 2 > x , x_old , x_err;
x[0] = 1.0 ; x[1] = 2.0;
x_old[0] = 0.5 ; x_old[1] = 1.5;
x_err[0] = 0.01 ; x_err[1] = 0.02;
// err1 = 0.01 / ( 0.001 + 0.0001 * max( 1.0 , 0.5 ) )
// err1 = 0.01 / ( 0.0011 )
// err1 = 9.090909090
// err2 = 0.02 / ( 0.001 + 0.0001 * max( 2.0 , 1 ) )
// err2 = 0.02 / ( 0.0012 )
// err2 = 16.6666666666666666
// err = sqrt( err1 * err1 + err2 + err2 )
err = error_checker.error( x_old , x , x_err , dt );
BOOST_CHECK_CLOSE( err , 9.49239703496 , 1.0e-10 );
x_err[0] *= -1.0;
err = error_checker.error( x_old , x , x_err , dt );
BOOST_CHECK_CLOSE( err , 9.49239703496 , 1.0e-10 );
x[0] *= -1.0;
err = error_checker.error( x_old , x , x_err , dt );
BOOST_CHECK_CLOSE( err , 9.49239703496 , 1.0e-10 );
x_old[0] *= -1.0;
err = error_checker.error( x_old , x , x_err , dt );
BOOST_CHECK_CLOSE( err , 9.49239703496 , 1.0e-10 );
}
BOOST_AUTO_TEST_CASE( test_version2 )
{
range_algebra algebra;
error_checker_l2_norm< double , range_algebra , default_operations > error_checker( algebra , 1.0e-3 , 1.0e-4 , 0.5 , 1.5 );
double dt = 0.02 , err;
boost::array< double , 2 > x , x_old , x_err , dxdt_old , dxdt;
x[0] = 1.0 ; x[1] = 2.0;
x_old[0] = 0.5 ; x_old[1] = 1.5;
x_err[0] = 0.01 ; x_err[1] = 0.02;
dxdt_old[0] = 2.0 ; dxdt_old[1] = 4.0;
dxdt[0] = 1000.0 , dxdt[1] = 255555.0;
// err1 = 0.01 / ( 0.001 + 0.0001 * ( 0.5 * 0.5 + 0.02 * 1.5 * 2.0 ) )
// err1 = 0.01 / ( 0.001 + 0.0001 * ( 0.25 + 0.06 ) )
// err1 = 0.01 / ( 0.001031
// err1 = 9.699321048
// err2 = 0.02 / ( 0.001 + 0.0001 * ( 0.5 * 1.5 + 0.02 * 1.5 * 4.0 ) )
// err2 = 0.02 / ( 0.001 + 0.0001 * ( 0.75 + 0.12 ) )
// err2 = 0.02 / ( 0.001087 )
// err2 = 18.399264029
err = error_checker.error( x_old , x , dxdt_old , x_err , dt );
BOOST_CHECK_CLOSE( err , 10.3996363591 , 1.0e-8 );
err = error_checker.error( x_old , x , dxdt_old , dxdt , x_err , dt );
BOOST_CHECK_CLOSE( err , 10.3996363591 , 1.0e-8 );
x_err[0] *= -1.0;
err = error_checker.error( x_old , x , dxdt_old , x_err , dt );
BOOST_CHECK_CLOSE( err , 10.3996363591 , 1.0e-8 );
x_old[0] *= -1.0;
err = error_checker.error( x_old , x , dxdt_old , x_err , dt );
BOOST_CHECK_CLOSE( err , 10.3996363591 , 1.0e-8 );
dxdt_old[0] *= 1.0;
err = error_checker.error( x_old , x , dxdt_old , x_err , dt );
BOOST_CHECK_CLOSE( err , 10.3996363591 , 1.0e-8 );
}
BOOST_AUTO_TEST_SUITE_END()

118
libs/numeric/odeint/ideas/controller/libs/numeric/odeint/test/error_checker_max_norm.cpp
View File

@@ -1,118 +0,0 @@
/*
* is_pair.cpp
*
* Created on: Feb 12, 2011
* Author: karsten
*/
#define BOOST_TEST_MODULE odeint_error_checker_max_norm
#include <utility>
#include <boost/test/unit_test.hpp>
#include <boost/static_assert.hpp>
#include <boost/numeric/odeint/algebra/range_algebra.hpp>
#include <boost/numeric/odeint/algebra/default_operations.hpp>
#include <boost/numeric/odeint/stepper/error_checker_max_norm.hpp>
#include <boost/array.hpp>
using namespace boost::numeric::odeint;
BOOST_AUTO_TEST_SUITE( error_checker_max_norm_test )
BOOST_AUTO_TEST_CASE( test_version1 )
{
// computes err = max( x_err / ( eps_abs + eps_rel * max( x , x_old ) ) );
range_algebra algebra;
error_checker_max_norm< double , range_algebra , default_operations > error_checker( algebra , 1.0e-3 , 1.0e-4 );
double dt , err;
boost::array< double , 2 > x , x_old , x_err;
x[0] = 1.0 ; x[1] = 2.0;
x_old[0] = 0.5 ; x_old[1] = 1.5;
x_err[0] = 0.01 ; x_err[1] = 0.02;
// err1 = 0.01 / ( 0.001 + 0.0001 * max( 1.0 , 0.5 ) )
// err1 = 0.01 / ( 0.0011 )
// err1 = 9.090909090
// err2 = 0.02 / ( 0.001 + 0.0001 * max( 2.0 , 1.5 ) )
// err2 = 0.02 / ( 0.0012 )
// err2 = 16.6666666666666666
// err = max( err1 , err2 )
err = error_checker.error( x_old , x , x_err , dt );
BOOST_CHECK_CLOSE( err , 16.66666666666 , 1.0e-10 );
x_err[1] = 0.00002;
err = error_checker.error( x_old , x , x_err , dt );
BOOST_CHECK_CLOSE( err , 9.0909090909091 , 1.0e-10 );
x_err[0] *= -1.0;
err = error_checker.error( x_old , x , x_err , dt );
BOOST_CHECK_CLOSE( err , 9.0909090909091 , 1.0e-10 );
x[0] *= -1.0;
err = error_checker.error( x_old , x , x_err , dt );
BOOST_CHECK_CLOSE( err , 9.0909090909091 , 1.0e-10 );
x_old[0] *= -1.0;
err = error_checker.error( x_old , x , x_err , dt );
BOOST_CHECK_CLOSE( err , 9.0909090909091 , 1.0e-10 );
}
BOOST_AUTO_TEST_CASE( test_version2 )
{
range_algebra algebra;
error_checker_max_norm< double , range_algebra , default_operations > error_checker( algebra , 1.0e-3 , 1.0e-4 , 0.5 , 1.5 );
double dt = 0.02 , err;
boost::array< double , 2 > x , x_old , x_err , dxdt_old , dxdt;
x[0] = 1.0 ; x[1] = 2.0;
x_old[0] = 0.5 ; x_old[1] = 1.5;
x_err[0] = 0.01 ; x_err[1] = 0.02;
dxdt_old[0] = 2.0 ; dxdt_old[1] = 4.0;
dxdt[0] = 1000.0 , dxdt[1] = 255555.0;
// err1 = 0.01 / ( 0.001 + 0.0001 * ( 0.5 * 0.5 + 0.02 * 1.5 * 2.0 ) )
// err1 = 0.01 / ( 0.001 + 0.0001 * ( 0.25 + 0.06 ) )
// err1 = 0.01 / ( 0.001031
// err1 = 9.699321048
// err2 = 0.02 / ( 0.001 + 0.0001 * ( 0.5 * 1.5 + 0.02 * 1.5 * 4.0 ) )
// err2 = 0.02 / ( 0.001 + 0.0001 * ( 0.75 + 0.12 ) )
// err2 = 0.02 / ( 0.001087 )
// err2 = 18.399264029
err = error_checker.error( x_old , x , dxdt_old , x_err , dt );
BOOST_CHECK_CLOSE( err , 18.399264029 , 1.0e-8 );
err = error_checker.error( x_old , x , dxdt_old , dxdt , x_err , dt );
// BOOST_CHECK_CLOSE( err , 18.399264029 , 1.0e-8 );
x_err[1] = 0.00002;
err = error_checker.error( x_old , x , dxdt_old , x_err , dt );
BOOST_CHECK_CLOSE( err , 9.699321048 , 1.0e-8 );
x_err[0] *= -1.0;
err = error_checker.error( x_old , x , dxdt_old , x_err , dt );
BOOST_CHECK_CLOSE( err , 9.699321048 , 1.0e-8 );
x_old[0] *= -1.0;
err = error_checker.error( x_old , x , dxdt_old , x_err , dt );
BOOST_CHECK_CLOSE( err , 9.699321048 , 1.0e-8 );
dxdt_old[0] *= 1.0;
err = error_checker.error( x_old , x , dxdt_old , x_err , dt );
BOOST_CHECK_CLOSE( err , 9.699321048 , 1.0e-8 );
}
BOOST_AUTO_TEST_SUITE_END()

57
libs/numeric/odeint/ideas/controller/libs/numeric/odeint/test/number_of_elements.cpp
View File

@@ -1,57 +0,0 @@
/*
* number_of_elements.cpp
*
* Created on: Feb 12, 2011
* Author: karsten
*/
#define BOOST_TEST_MODULE odeint_number_of_elements
#include <vector>
#include <list>
#include <boost/array.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/test/unit_test.hpp>
#include <boost/numeric/odeint/util/number_of_elements.hpp>
using namespace boost::numeric::odeint;
BOOST_AUTO_TEST_SUITE( number_of_elements_test )
BOOST_AUTO_TEST_CASE( test_number_of_elements_vector )
{
typedef std::vector< double > state_type;
state_type v( 10 );
BOOST_CHECK_EQUAL( size_t( 10 ) , number_of_elements( v ) );
}
BOOST_AUTO_TEST_CASE( test_number_of_elements_boost_array )
{
typedef boost::array< double , 5 > state_type;
state_type v;
BOOST_CHECK_EQUAL( size_t( 5 ) , number_of_elements( v ) );
}
BOOST_AUTO_TEST_CASE( test_number_of_elements_list )
{
typedef std::list< double > state_type;
state_type v( 10 );
BOOST_CHECK_EQUAL( size_t( 10 ) , number_of_elements( v ) );
}
BOOST_AUTO_TEST_CASE( test_number_of_elements_ublas_vector )
{
typedef boost::numeric::ublas::vector< double > state_type;
state_type v( 10 );
BOOST_CHECK_EQUAL( size_t( 10 ) , number_of_elements( v ) );
}
BOOST_AUTO_TEST_SUITE_END()

65
libs/numeric/odeint/ideas/controller/libs/numeric/odeint/test/ref_or_value_holder.cpp
View File

@@ -1,65 +0,0 @@
/*
* is_pair.cpp
*
* Created on: Feb 12, 2011
* Author: karsten
*/
#define BOOST_TEST_MODULE odeint_ref_or_value
#include <utility>
#include <boost/test/unit_test.hpp>
#include <boost/static_assert.hpp>
#include <boost/numeric/odeint/util/ref_or_value_holder.hpp>
using namespace boost::numeric::odeint;
BOOST_AUTO_TEST_SUITE( ref_or_value_holder_test )
BOOST_AUTO_TEST_CASE( test_value_holder )
{
int a = 1;
ref_or_value_holder< int , false > value( a );
value.get() = 2;
BOOST_CHECK_EQUAL( 1 , a );
BOOST_CHECK_EQUAL( 2 , value.get() );
}
BOOST_AUTO_TEST_CASE( test_ref_holder )
{
int a = 1;
ref_or_value_holder< int , true > value( a );
value.get() = 2;
BOOST_CHECK_EQUAL( 2 , a );
BOOST_CHECK_EQUAL( 2 , value.get() );
}
BOOST_AUTO_TEST_CASE( test_const_value_holder )
{
int a = 1;
ref_or_value_holder< const int , false > value( a );
BOOST_CHECK_EQUAL( 1 , a );
BOOST_CHECK_EQUAL( 1 , value.get() );
}
BOOST_AUTO_TEST_CASE( test_const_ref_holder )
{
int a = 1;
ref_or_value_holder< const int , true > value( a );
BOOST_CHECK_EQUAL( 1 , a );
BOOST_CHECK_EQUAL( 1 , value.get() );
}
BOOST_AUTO_TEST_SUITE_END()

223
libs/numeric/odeint/ideas/controller/rb_controller.cpp
View File

@@ -1,223 +0,0 @@
/*
* TODO:
*
* * implement rosenbrock_error_checker
* * implement rosenbrock_controller
* * implement generic_controlled_stepper
*/
#include <iostream>
#include <fstream>
#include <cmath>
#include <boost/timer.hpp>
#include <boost/numeric/odeint/stepper/rosenbrock4.hpp>
#include <boost/numeric/odeint/stepper/rosenbrock4_controller.hpp>
#include <boost/numeric/odeint/stepper/rosenbrock4_dense_output.hpp>
#include <boost/numeric/odeint/stepper/generation.hpp>
#include <boost/numeric/odeint/integrate/integrate_adaptive.hpp>
#include <boost/numeric/odeint/stepper/error_checker_explicit.hpp>
#include <boost/numeric/odeint/stepper/generic_controlled_stepper.hpp>
#include <boost/numeric/odeint/stepper/generic_controlled_stepper_explicit.hpp>
#include <boost/numeric/odeint/stepper/controller/default_controller.hpp>
#include <boost/numeric/odeint/stepper/controller/pi_controller.hpp>
#include <boost/numeric/odeint/stepper/controller/stiff_controller.hpp>
#define tab "\t"
namespace boost {
namespace numeric {
namespace odeint {
template< class Value , class Algebra , class Operations >
class error_checker
{
public:
typedef Value value_type;
typedef Algebra algebra_type;
typedef Operations operations_type;
error_checker(
const value_type eps_abs = static_cast< value_type >( 1.0e-6 ) ,
const value_type eps_rel = static_cast< value_type >( 1.0e-6 )
)
: m_eps_abs( eps_abs ) , m_eps_rel( eps_rel )
{ }
template< class State1 , class State2 , class Err , class Time >
value_type error( const State1 &x_old , const State2 &x , Err &x_err , const Time &dt )
{
// this overwrites x_err !
algebra.for_each3( x_err , x_old , x ,
typename operations_type::template default_rel_error< value_type >( m_eps_abs , m_eps_rel ) );
value_type res = algebra.reduce( x_err ,
typename operations_type::template maximum< value_type >() , static_cast< value_type >( 0.0 ) );
return res;
}
private:
value_type m_eps_abs;
value_type m_eps_rel;
algebra_type algebra;
};
template< class Value >
class rosenbrock_error_checker
{
public:
typedef Value value_type;
rosenbrock_error_checker(
const value_type eps_abs = static_cast< value_type >( 1.0e-6 ) ,
const value_type eps_rel = static_cast< value_type >( 1.0e-6 )
)
: m_eps_abs( eps_abs ) , m_eps_rel( eps_rel )
{ }
template< class State1 , class State2 , class Err , class Time >
value_type error( const State1 &x_old , const State2 &x , Err &x_err , const Time &dt )
{
value_type res = 0.0 , sk = 0.0;
for( size_t i=0 ; i<x_old.size() ; ++i )
{
sk = m_eps_abs + m_eps_rel * std::max( std::abs( x_old[i] ) , std::abs( x[i] ) );
res += x_err[i] * x_err[i] / sk / sk;
}
return sqrt( res / value_type( x_old.size() ) );
}
private:
value_type m_eps_abs;
value_type m_eps_rel;
};
/*
* ToDo: move this class into rosenbrock4_controller.hpp
*/
template< class Value >
class rosenbrock4_controller_tmp : public generic_controlled_stepper<
rosenbrock4< double > ,
rosenbrock_error_checker< double > ,
stiff_controller< double > ,
initially_resizer ,
error_stepper_tag >
{
};
}
}
}
typedef boost::numeric::odeint::rosenbrock4< double > stepper_type;
typedef stepper_type::state_type state_type;
typedef stepper_type::matrix_type matrix_type;
const double mu = 1000.0;
struct vdp_stiff
{
template< class State , class Deriv >
void operator()( const State &x , Deriv &dxdt , double t )
{
dxdt[0] = x[1];
dxdt[1] = -x[0] - mu * x[1] * (x[0]*x[0]-1.0);
}
};
struct vdp_stiff_jacobi
{
void operator()( const state_type &x , matrix_type &J , const double &t , state_type &dfdt )
{
J(0, 0) = 0.0;
J(0, 1) = 1.0;
J(1, 0) = -1.0 - 2.0*mu * x[0] * x[1];
J(1, 1) = -mu * ( x[0] * x[0] - 1.0);
dfdt[0] = 0.0;
dfdt[1] = 0.0;
}
};
struct streaming_observer
{
std::ostream& m_out;
streaming_observer( std::ostream &out ) : m_out( out ) { }
template< class State >
void operator()( const State &x , double t ) const
{
m_out << t;
for( size_t i=0 ; i<x.size() ; ++i ) m_out << "\t" << x[i];
m_out << "\n";
}
};
int main( int argc , char **argv )
{
using namespace std;
using namespace boost::numeric::odeint;
stepper_type stepper1;
typedef rosenbrock4_controller_tmp< double > stepper2_type;
state_type x1( 2 );
x1[ 0 ] = 1.0 ; x1[ 1 ] = 1.0;
state_type x2 = x1;
state_type x3 = x1;
state_type x4 = x1;
boost::timer timer;
const double t_max = 100000.0;
ofstream fout1( "vdp_rb_1.dat" );
timer.restart();
size_t steps1 = integrate_adaptive( make_controlled( 1.0e-6 , 1.0e-6 , stepper1 ) ,
make_pair( vdp_stiff() , vdp_stiff_jacobi() ) , x1 , 0.0 , t_max , 0.1 , streaming_observer( fout1 ) );
double t1 = timer.elapsed();
cout << steps1 << tab << t1 << tab << x1[0] << tab << x1[1] << endl;
ofstream fout2( "vdp_rb_2.dat" );
timer.restart();
size_t steps2 = integrate_adaptive( stepper2_type() ,
make_pair( vdp_stiff() , vdp_stiff_jacobi() ) , x2 , 0.0 , t_max , 0.1 , streaming_observer( fout2 ) );
double t2 = timer.elapsed();
cout << steps2 << tab << t2 << tab << x2[0] << tab << x2[1] << endl;
return 0;
}

77
libs/numeric/odeint/ideas/controller/rb_dense_output.cpp
View File

@@ -1,77 +0,0 @@
#include <iostream>
#include <cmath>
#include <boost/timer.hpp>
#include <boost/numeric/odeint/stepper/rosenbrock4_controller.hpp>
#include <boost/numeric/odeint/stepper/generic_dense_output_stepper.hpp>
#include <boost/numeric/odeint/stepper/generation.hpp>
#include <boost/numeric/odeint/integrate/integrate_adaptive.hpp>
#define tab "\t"
typedef boost::numeric::odeint::rosenbrock4< double > stepper_type;
typedef stepper_type::state_type state_type;
typedef stepper_type::matrix_type matrix_type;
const double mu = 1000.0;
struct vdp_stiff
{
template< class State , class Deriv >
void operator()( const State &x , Deriv &dxdt , double t )
{
dxdt[0] = x[1];
dxdt[1] = -x[0] - mu * x[1] * (x[0]*x[0]-1.0);
}
};
struct vdp_stiff_jacobi
{
void operator()( const state_type &x , matrix_type &J , const double &t , state_type &dfdt )
{
J(0, 0) = 0.0;
J(0, 1) = 1.0;
J(1, 0) = -1.0 - 2.0*mu * x[0] * x[1];
J(1, 1) = -mu * ( x[0] * x[0] - 1.0);
dfdt[0] = 0.0;
dfdt[1] = 0.0;
}
};
int main( int argc , char **argv )
{
using namespace std;
using namespace boost::numeric::odeint;
stepper_type stepper1;
state_type x1( 2 );
x1[ 0 ] = 1.0 ; x1[ 1 ] = 1.0;
state_type x2 = x1;
state_type x3 = x1;
state_type x4 = x1;
boost::timer timer;
const double t_max = 10000000.0;
timer.restart();
size_t steps1 = integrate_adaptive( make_dense_output( 1.0e-6 , 1.0e-6 , stepper1 ) ,
make_pair( vdp_stiff() , vdp_stiff_jacobi() ) , x1 , 0.0 , t_max , 0.1 );
double t1 = timer.elapsed();
cout << steps1 << tab << t1 << tab << x1[0] << tab << x1[1] << endl;
generic_dense_output_stepper< rosenbrock4_controller< rosenbrock4< double > > > stepper2;
timer.restart();
size_t steps2 = integrate_adaptive( stepper2 , make_pair( vdp_stiff() , vdp_stiff_jacobi() ) , x2 , 0.0 , t_max , 0.1 );
double t2 = timer.elapsed();
cout << steps2 << tab << t2 << tab << x2[0] << tab << x2[1] << endl;
return 0;
}

17
libs/numeric/odeint/ideas/fusion_runge_kutta/Jamfile
View File

@@ -1,17 +0,0 @@
# (C) Copyright 2010 : Karsten Ahnert, Mario Mulansky
# Distributed under the Boost Software License, Version 1.0.
# (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
import os ;
import modules ;
import path ;
project
: requirements
<define>BOOST_ALL_NO_LIB=1
<include>../../../../..
;
exe test_explicit_rk
: test_explicit_rk.cpp
;

197
libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_algebra.hpp
View File

@@ -1,197 +0,0 @@
/*
* fusion_algebra.hpp
*
* Created on: Apr 26, 2011
* Author: mario
*/
#ifndef FUSION_ALGEBRA_HPP_
#define FUSION_ALGEBRA_HPP_
#include <boost/array.hpp>
#include <iostream>
template< size_t n >
struct fusion_algebra
{
template< typename T , size_t dim >
inline static void foreach( boost::array< T , dim > &x_tmp , const boost::array< T , dim > &x ,
const boost::array< double , n > &a ,
const boost::array< T , dim > k_vector[n] , const double dt )
{
for( size_t i=0 ; i<dim ; ++i )
{
x_tmp[i] = x[i];// + a[0]*dt*k_vector[0][i];
for( size_t j = 0 ; j<n ; ++j )
x_tmp[i] += a[j]*dt*k_vector[j][i];
}
}
template< typename T , size_t dim >
inline static void foreach( boost::array< T , dim > &x_tmp ,
const boost::array< double , n > &a ,
const boost::array< T , dim > k_vector[n] , const double dt )
{
for( size_t i=0 ; i<dim ; ++i )
{
x_tmp[i] = a[0]*dt*k_vector[0][i];
for( size_t j = 1 ; j<n ; ++j )
x_tmp[i] += a[j]*dt*k_vector[j][i];
}
}
};
/** hand-wise implementation for performance improvement for n = 1..4 **/
/* !!!!!!! Actually, this is factor 3 slower with intel compiler, so we don'y use it !!!!!
* Update: It increases performance on msvc 9.0 by about 30%, so it is activated for MSVC
*/
//#ifdef BOOST_MSVC
template<>
struct fusion_algebra< 1 >
{
template< typename T , size_t dim >
inline static void foreach( boost::array< T , dim > &x_tmp , const boost::array< T , dim > &x ,
const boost::array< double , 1 > &a ,
const boost::array< T , dim > *k_vector , const double dt )
{
for( size_t i=0 ; i<dim ; ++i )
{
x_tmp[i] = x[i]
+ a[0]*dt*k_vector[0][i];
}
}
};
template<>
struct fusion_algebra< 2 >
{
template< typename T , size_t dim >
inline static void foreach( boost::array< T , dim > &x_tmp , const boost::array< T , dim > &x ,
const boost::array< double , 2 > &a ,
const boost::array< T , dim > *k_vector , const double dt )
{
for( size_t i=0 ; i<dim ; ++i )
{
x_tmp[i] = x[i]
+ a[0]*dt*k_vector[0][i]
+ a[1]*dt*k_vector[1][i];
}
}
};
template<>
struct fusion_algebra< 3 >
{
template< typename T , size_t dim >
inline static void foreach( boost::array< T , dim > &x_tmp , const boost::array< T , dim > &x ,
const boost::array< double , 3 > &a ,
const boost::array< T , dim > *k_vector , const double dt )
{
for( size_t i=0 ; i<dim ; ++i )
{
x_tmp[i] = x[i]
+ a[0]*dt*k_vector[0][i]
+ a[1]*dt*k_vector[1][i]
+ a[2]*dt*k_vector[2][i];
}
}
};
template<>
struct fusion_algebra< 4 >
{
template< typename T , size_t dim >
inline static void foreach( boost::array< T , dim > &x_tmp , const boost::array< T , dim > &x ,
const boost::array< double , 4 > &a ,
const boost::array< T , dim > *k_vector , const double dt )
{
for( size_t i=0 ; i<dim ; ++i )
{
x_tmp[i] = x[i]
+ a[0]*dt*k_vector[0][i]
+ a[1]*dt*k_vector[1][i]
+ a[2]*dt*k_vector[2][i]
+ a[3]*dt*k_vector[3][i];
}
}
};
template<>
struct fusion_algebra< 5 >
{
template< typename T , size_t dim >
inline static void foreach( boost::array< T , dim > &x_tmp , const boost::array< T , dim > &x ,
const boost::array< double , 5 > &a ,
const boost::array< T , dim > *k_vector , const double dt )
{
for( size_t i=0 ; i<dim ; ++i )
{
x_tmp[i] = x[i]
+ a[0]*dt*k_vector[0][i]
+ a[1]*dt*k_vector[1][i]
+ a[2]*dt*k_vector[2][i]
+ a[3]*dt*k_vector[3][i]
+ a[4]*dt*k_vector[4][i];
}
}
};
template<>
struct fusion_algebra< 6 >
{
template< typename T , size_t dim >
inline static void foreach( boost::array< T , dim > &x_tmp , const boost::array< T , dim > &x ,
const boost::array< double , 6 > &a ,
const boost::array< T , dim > *k_vector , const double dt )
{
for( size_t i=0 ; i<dim ; ++i )
{
x_tmp[i] = x[i]
+ a[0]*dt*k_vector[0][i]
+ a[1]*dt*k_vector[1][i]
+ a[2]*dt*k_vector[2][i]
+ a[3]*dt*k_vector[3][i]
+ a[4]*dt*k_vector[4][i]
+ a[5]*dt*k_vector[5][i];
}
}
template< typename T , size_t dim >
inline static void foreach(boost::array<T , dim> &x_tmp ,
const boost::array<double , 6> &a ,
const boost::array<T , dim> *k_vector , const double dt)
{
for (size_t i = 0 ; i < dim ; ++i)
{
x_tmp[i] = a[0] * dt * k_vector[0][i] + a[1] * dt * k_vector[1][i]
+ a[2] * dt * k_vector[2][i] + a[3] * dt * k_vector[3][i]
+ a[4] * dt * k_vector[4][i] + a[5] * dt * k_vector[5][i];
}
}
};
//#endif /* BOOST_MSVC */
#endif /* FUSION_ALGEBRA_HPP_ */

58
libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_explicit_error_rk.hpp
View File

@@ -1,58 +0,0 @@
/*
* fusion_explicit_error_rk.hpp
*
* Created on: Apr 29, 2011
* Author: mario
*/
#ifndef FUSION_EXPLICIT_ERROR_RK_HPP_
#define FUSION_EXPLICIT_ERROR_RK_HPP_
#include "fusion_explicit_rk_new.hpp"
#include "fusion_algebra.hpp"
namespace mpl = boost::mpl;
namespace fusion = boost::fusion;
using namespace std;
template< class StateType , size_t stage_count >
class explicit_error_rk : public explicit_rk< StateType , stage_count >
{
public:
typedef explicit_rk< StateType , stage_count > base;
typedef StateType state_type;
typedef typename base::stage_indices stage_indices;
typedef typename base::coef_a_type coef_a_type;
typedef typename base::coef_b_type coef_b_type;
typedef typename base::coef_c_type coef_c_type;
public:
explicit_error_rk( const coef_a_type &a ,
const coef_b_type &b ,
const coef_b_type &b2 ,
const coef_c_type &c )
: base( a , b , c ) , m_b2( b2 )
{ }
template< class System >
void inline do_step( System system , state_type &x , const double t , const double dt , state_type &x_err )
{
base::do_step( system , x , t , dt );
// compute error estimate
fusion_algebra< stage_count >::foreach( x_err , m_b2 , base::m_F , dt );
}
private:
const coef_b_type m_b2;
};
#endif /* FUSION_EXPLICIT_ERROR_RK_HPP_ */

211
libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_explicit_rk.hpp
View File

@@ -1,211 +0,0 @@
/*
* fusion_runge_kutta.hpp
*
* Created on: Apr 26, 2011
* Author: mario
*/
#ifndef FUSION_EXPLICIT_RK_HPP_
#define FUSION_EXPLICIT_RK_HPP_
#include <boost/mpl/vector.hpp>
#include <boost/mpl/push_back.hpp>
#include <boost/mpl/for_each.hpp>
#include <boost/mpl/range_c.hpp>
#include <boost/mpl/copy.hpp>
#include <boost/mpl/size_t.hpp>
#include <boost/fusion/container.hpp>
#include <boost/fusion/algorithm.hpp>
#include <boost/fusion/sequence.hpp>
#include <boost/fusion/adapted.hpp>
#include <algorithm>
#include <iostream>
#include <string>
#include <boost/array.hpp>
#include <typeinfo>
#include "fusion_algebra.hpp"
namespace mpl = boost::mpl;
namespace fusion = boost::fusion;
using namespace std;
struct intermediate_stage {};
struct last_stage {};
template< class T , class Constant >
struct array_wrapper
{
typedef typename boost::array< T , Constant::value > type;
};
template< class T , class Constant , class StageCategory >
struct stage_fusion_wrapper
{
typedef typename fusion::vector< size_t , T , boost::array< T , Constant::value > , StageCategory > type;
};
template< class StateType , size_t stage_count >
class explicit_rk
{
public:
typedef StateType state_type;
typedef mpl::range_c< size_t , 1 , stage_count > stage_indices;
typedef typename fusion::result_of::as_vector
<
typename mpl::copy
<
stage_indices ,
mpl::inserter
<
mpl::vector0< > ,
mpl::push_back< mpl::_1 , array_wrapper< double , mpl::_2 > >
>
>::type
>::type coef_a_type;
typedef boost::array< double , stage_count > coef_b_type;
typedef boost::array< double , stage_count > coef_c_type;
typedef typename fusion::result_of::as_vector
<
typename mpl::push_back
<
typename mpl::copy
<
stage_indices,
mpl::inserter
<
mpl::vector0<> ,
mpl::push_back< mpl::_1 , stage_fusion_wrapper< double , mpl::_2 , intermediate_stage > >
>
>::type ,
typename stage_fusion_wrapper< double , mpl::size_t< stage_count > , last_stage >::type
>::type
>::type stage_vector_base;
struct stage_vector : public stage_vector_base
{
struct do_insertion
{
stage_vector_base &m_base;
const coef_a_type &m_a;
const coef_c_type &m_c;
do_insertion( stage_vector_base &base , const coef_a_type &a , const coef_c_type &c )
: m_base( base ) , m_a( a ) , m_c( c ) { }
template< class Index >
void operator()( Index ) const
{
fusion::at_c< 0 >( fusion::at< Index >( m_base ) ) = Index::value;
fusion::at_c< 1 >( fusion::at< Index >( m_base ) ) = m_c[ Index::value ];
fusion::at_c< 2 >( fusion::at< Index >( m_base ) ) = fusion::at< Index >( m_a );
}
};
stage_vector( const coef_a_type &a , const coef_b_type &b , const coef_c_type &c )
{
typedef mpl::range_c< size_t , 0 , stage_count - 1 > indices;
mpl::for_each< indices >( do_insertion( *this , a , c ) );
fusion::at_c< 0 >( fusion::at_c< stage_count - 1 >( *this ) ) = stage_count - 1 ;
fusion::at_c< 1 >( fusion::at_c< stage_count - 1 >( *this ) ) = c[ stage_count - 1 ];
fusion::at_c< 2 >( fusion::at_c< stage_count - 1 >( *this ) ) = b;
}
};
template< class System >
struct calculate_stage
{
System &system;
state_type &x , &x_tmp;
state_type *k_vector;
const double t;
const double dt;
calculate_stage( System &_system , state_type &_x , state_type &_x_tmp , state_type *_k_vector ,
const double _t , const double _dt )
: system( _system ) , x( _x ) , x_tmp( _x_tmp ) , k_vector( _k_vector ) , t( _t ) , dt( _dt )
{}
template< typename T , size_t stage_number >
void operator()( fusion::vector< size_t , T , boost::array< T , stage_number > , intermediate_stage > const &stage ) const
//typename stage_fusion_wrapper< T , mpl::size_t< stage_number > , intermediate_stage >::type const &stage ) const
{
double c = fusion::at_c< 1 >( stage );
if( stage_number == 1 )
system( x , k_vector[stage_number-1] , t + c * dt );
else
system( x_tmp , k_vector[stage_number-1] , t + c * dt );
fusion_algebra<stage_number>::foreach( x_tmp , x , fusion::at_c< 2 >( stage ) , k_vector , dt);
}
template< typename T , size_t stage_number >
void operator()( fusion::vector< size_t , T , boost::array< T , stage_number > , last_stage > const &stage ) const
//void operator()( typename stage_fusion_wrapper< T , mpl::size_t< stage_number > , last_stage >::type const &stage ) const
{
double c = fusion::at_c< 1 >( stage );
if( stage_number == 1 )
system( x , k_vector[stage_number-1] , t + c * dt );
else
system( x_tmp , k_vector[stage_number-1] , t + c * dt );
fusion_algebra<stage_number>::foreach( x , x , fusion::at_c< 2 >( stage ) , k_vector , dt);
}
};
public:
explicit_rk( const coef_a_type &a ,
const coef_b_type &b ,
const coef_c_type &c )
: m_a( a ) , m_b( b ) , m_c( c ) ,
m_stages( a , b , c )
{ }
template< class System >
void do_step( System &system , state_type &x , double t , const double dt )
{
fusion::for_each( m_stages , calculate_stage< System >( system , x , m_x_tmp , m_k_vector , t , dt ) );
}
private:
const coef_a_type m_a;
const coef_b_type m_b;
const coef_c_type m_c;
const stage_vector m_stages;
state_type m_x_tmp;
state_type m_k_vector[stage_count];
};
#endif /* FUSION_EXPLICIT_RK_HPP_ */

212
libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_explicit_rk_new.hpp
View File

@@ -1,212 +0,0 @@
/*
* fusion_runge_kutta.hpp
*
* Created on: Apr 26, 2011
* Author: mario
*/
#ifndef FUSION_EXPLICIT_RK_HPP_
#define FUSION_EXPLICIT_RK_HPP_
#include <boost/mpl/vector.hpp>
#include <boost/mpl/push_back.hpp>
#include <boost/mpl/for_each.hpp>
#include <boost/mpl/range_c.hpp>
#include <boost/mpl/copy.hpp>
#include <boost/mpl/size_t.hpp>
#include <boost/fusion/container.hpp>
#include <boost/fusion/algorithm/iteration.hpp>
#include <boost/array.hpp>
#include "fusion_algebra.hpp"
//#include "fusion_foreach_performance.hpp"
namespace mpl = boost::mpl;
namespace fusion = boost::fusion;
using namespace std;
struct intermediate_stage {};
struct last_stage {};
template< class T , class Constant >
struct array_wrapper
{
typedef const typename boost::array< T , Constant::value > type;
};
template< class T , size_t i , class StageCategory >
struct stage
{
T c;
boost::array< T , i > a;
typedef StageCategory category;
};
template< class T , size_t i>
struct stage< T , i , last_stage >
{
T c;
boost::array< T , i > b;
typedef last_stage category;
};
template< class T , class Constant , class StageCategory >
struct stage_wrapper
{
typedef stage< T , Constant::value , StageCategory > type;
};
template< class StateType , size_t stage_count >
class explicit_rk
{
public:
typedef StateType state_type;
typedef mpl::range_c< size_t , 1 , stage_count > stage_indices;
typedef typename fusion::result_of::as_vector
<
typename mpl::copy
<
stage_indices ,
mpl::inserter
<
mpl::vector0< > ,
mpl::push_back< mpl::_1 , array_wrapper< double , mpl::_2 > >
>
>::type
>::type coef_a_type;
typedef boost::array< double , stage_count > coef_b_type;
typedef boost::array< double , stage_count > coef_c_type;
typedef typename fusion::result_of::as_vector
<
typename mpl::push_back
<
typename mpl::copy
<
stage_indices,
mpl::inserter
<
mpl::vector0<> ,
mpl::push_back< mpl::_1 , stage_wrapper< double , mpl::_2 , intermediate_stage > >
>
>::type ,
stage< double , stage_count , last_stage >
>::type
>::type stage_vector_base;
struct stage_vector : public stage_vector_base
{
struct do_insertion
{
stage_vector_base &m_base;
const coef_a_type &m_a;
const coef_c_type &m_c;
do_insertion( stage_vector_base &base , const coef_a_type &a , const coef_c_type &c )
: m_base( base ) , m_a( a ) , m_c( c ) { }
template< class Index >
void operator()( Index ) const
{
//fusion::at< Index >( m_base ) = stage< double , Index::value+1 , intermediate_stage >( m_c[ Index::value ] , fusion::at< Index >( m_a ) );
fusion::at< Index >( m_base ).c = m_c[ Index::value ];
fusion::at< Index >( m_base ).a = fusion::at< Index >( m_a );
}
};
stage_vector( const coef_a_type &a , const coef_b_type &b , const coef_c_type &c )
{
typedef mpl::range_c< size_t , 0 , stage_count - 1 > indices;
mpl::for_each< indices >( do_insertion( *this , a , c ) );
//fusion::at_c< 0 >( fusion::at_c< stage_count - 1 >( *this ) ) = stage_count - 1 ;
fusion::at_c< stage_count - 1 >( *this ).c = c[ stage_count - 1 ];
fusion::at_c< stage_count - 1 >( *this ).b = b;
}
};
template< class System >
struct calculate_stage
{
System &system;
state_type &x , &x_tmp;
state_type *F;
const double t;
const double dt;
calculate_stage( System &_system , state_type &_x , state_type &_x_tmp , state_type *_F ,
const double _t , const double _dt )
: system( _system ) , x( _x ) , x_tmp( _x_tmp ) , F( _F ) , t( _t ) , dt( _dt )
{}
template< typename T , size_t stage_number >
void inline operator()( stage< T , stage_number , intermediate_stage > const &stage ) const
//typename stage_fusion_wrapper< T , mpl::size_t< stage_number > , intermediate_stage >::type const &stage ) const
{
if( stage_number == 1 )
system( x , F[stage_number-1] , t + stage.c * dt );
else
system( x_tmp , F[stage_number-1] , t + stage.c * dt );
fusion_algebra<stage_number>::foreach( x_tmp , x , stage.a , F , dt);
}
template< typename T , size_t stage_number >
void inline operator()( stage< T , stage_number , last_stage > const &stage ) const
//void operator()( typename stage_fusion_wrapper< T , mpl::size_t< stage_number > , last_stage >::type const &stage ) const
{
if( stage_number == 1 )
system( x , F[stage_number-1] , t + stage.c * dt );
else
system( x_tmp , F[stage_number-1] , t + stage.c * dt );
fusion_algebra<stage_number>::foreach( x , x , stage.b , F , dt);
}
};
public:
explicit_rk( const coef_a_type &a ,
const coef_b_type &b ,
const coef_c_type &c )
: m_stages( a , b , c )
{ }
template< class System >
void inline do_step( System system , state_type &x , const double t , const double dt )
{
fusion::for_each( m_stages , calculate_stage< System >( system , x , m_x_tmp , m_F , t , dt ) );
}
private:
stage_vector m_stages;
state_type m_x_tmp;
protected:
state_type m_F[stage_count];
};
#endif /* FUSION_EXPLICIT_RK_HPP_ */

39
libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_foreach_performance.hpp
View File

@@ -1,39 +0,0 @@
#include <boost/fusion/sequence/intrinsic/begin.hpp>
#include <boost/fusion/sequence/intrinsic/end.hpp>
#include <boost/fusion/iterator/equal_to.hpp>
#include <boost/fusion/iterator/next.hpp>
#include <boost/fusion/iterator/deref.hpp>
#include <boost/fusion/iterator/distance.hpp>
#include <boost/mpl/bool.hpp>
#include <iostream>
namespace boost { namespace fusion {
namespace detail
{
template<>
struct for_each_unrolled<6>
{
template<typename I0, typename F>
static void inline call(I0 const& i0, F const& f)
{
f(*i0);
typedef typename result_of::next<I0>::type I1;
I1 i1(fusion::next(i0));
f(*i1);
typedef typename result_of::next<I1>::type I2;
I2 i2(fusion::next(i1));
f(*i2);
typedef typename result_of::next<I2>::type I3;
I3 i3(fusion::next(i2));
f(*i3);
typedef typename result_of::next<I3>::type I4;
I4 i4(fusion::next(i3));
f(*i4);
typedef typename result_of::next<I4>::type I5;
I5 i5(fusion::next(i4));
f(*i5);
}
};
}
} }

131
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/Jamfile
View File

@@ -1,131 +0,0 @@
# (C) Copyright 2010 : Karsten Ahnert, Mario Mulansky
# Distributed under the Boost Software License, Version 1.0.
# (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
import os ;
import modules ;
import path ;
project
: requirements
<define>BOOST_ALL_NO_LIB=1
<include>../../../../../..
: default-build release
;
lib libgsl : : <name>gsl ;
lib libgslcblas : : <name>gslcblas ;
#exe generic_rk4
# : generic_rk4.cpp
# : <toolset>intel:<cxxflags>-inline-forceinline
# ;
exe generic_rk4_lorenz
: generic_rk4_lorenz.cpp
: <toolset>intel:<cxxflags>-inline-forceinline
;
exe generic_odeint_rk4_lorenz
: generic_odeint_rk4_lorenz.cpp
: <toolset>intel:<cxxflags>-inline-forceinline
;
#exe odeint_rk4
# : odeint_rk4.cpp
# : <cxxflags>-Winline
# ;
exe odeint_rk4_lorenz
: odeint_rk4_lorenz.cpp
;
exe odeint_rk4_lorenz_def_alg
: odeint_rk4_lorenz_def_alg.cpp
;
#exe nr_rk4
# : nr_rk4.cpp
# ;
exe nr_rk4_lorenz
: nr_rk4_lorenz.cpp
;
#exe rt_generic_rk4
# : rt_generic_rk4.cpp
# ;
exe rt_generic_rk4_lorenz
: rt_generic_rk4_lorenz.cpp
;
#exe gsl_rk4
# : gsl_rk4.cpp libgsl libgslcblas
# ;
exe gsl_rk4_lorenz
: gsl_rk4_lorenz.cpp libgsl libgslcblas
;
exe generic_rk54ck
: generic_rk54ck.cpp
: <toolset>intel:<cxxflags>-inline-forceinline
;
exe generic_rk54ck_lorenz
: generic_rk54ck_lorenz.cpp
: <toolset>intel:<cxxflags>-inline-forceinline
;
exe odeint_rk54ck
: odeint_rk54ck.cpp
;
exe odeint_rk54ck_lorenz
: odeint_rk54ck_lorenz.cpp
;
exe odeint_rk54ck_def_alg
: odeint_rk54ck_def_alg.cpp
;
exe odeint_rk54ck_lorenz_def_alg
: odeint_rk54ck_lorenz_def_alg.cpp
;
exe nr_rk54ck
: nr_rk54ck.cpp
;
exe nr_rk54ck_lorenz
: nr_rk54ck_lorenz.cpp
;
exe gsl_rk54ck
: gsl_rk54ck.cpp libgsl libgslcblas
;
exe gsl_rk54ck_lorenz
: gsl_rk54ck_lorenz.cpp libgsl libgslcblas
;
exe generic_rk4_phase_lattice
: generic_rk4_phase_lattice.cpp
: <toolset>intel:<cxxflags>-inline-forceinline
;
exe odeint_rk4_phase_lattice
: odeint_rk4_phase_lattice.cpp
;
exe nr_rk4_phase_lattice
: nr_rk4_phase_lattice.cpp
;
exe rt_generic_rk4_phase_lattice
: rt_generic_rk4_phase_lattice.cpp
;

71
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/generic_odeint_rk4_lorenz.cpp
View File

@@ -1,71 +0,0 @@
#include <boost/array.hpp>
//#include <boost/numeric/odeint/stepper/explicit_generic_rk.hpp>
#include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
#include <boost/numeric/odeint/algebra/array_algebra.hpp>
#include "rk_performance_test_case.hpp"
#include "lorenz.hpp"
using namespace boost::numeric::odeint;
typedef boost::array< double , 3 > state_type;
/*typedef explicit_generic_rk< 4 , 4 , state_type , double , state_type , double , array_algebra > rk4_type;
typedef rk4_type::coef_a_type coef_a_type;
typedef rk4_type::coef_b_type coef_b_type;
typedef rk4_type::coef_c_type coef_c_type;
const boost::array< double , 1 > a1 = {{ 0.5 }};
const boost::array< double , 2 > a2 = {{ 0.0 , 0.5 }};
const boost::array< double , 3 > a3 = {{ 0.0 , 0.0 , 1.0 }};
const coef_a_type a = fusion::make_vector( a1 , a2 , a3 );
const coef_b_type b = {{ 1.0/6 , 1.0/3 , 1.0/3 , 1.0/6 }};
const coef_c_type c = {{ 0.0 , 0.5 , 0.5 , 1.0 }};
*/
typedef runge_kutta4< state_type , double , state_type , double , array_algebra > rk4_type;
class rk4_wrapper
{
public:
rk4_wrapper()
{ }
void reset_init_cond()
{
m_x[0] = 10.0 * rand() / RAND_MAX;
m_x[1] = 10.0 * rand() / RAND_MAX;
m_x[2] = 10.0 * rand() / RAND_MAX;
m_t = 0.0;
}
inline void do_step( const double dt )
{
m_stepper.do_step( lorenz(), m_x , m_t , dt );
}
double state( const size_t i ) const
{ return m_x[i]; }
private:
state_type m_x;
double m_t;
rk4_type m_stepper;
};
int main()
{
srand( 12312354 );
rk4_wrapper stepper;
run( stepper );
}

63
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/generic_rk4_lorenz.cpp
View File

@@ -1,63 +0,0 @@
#include <boost/array.hpp>
#include "../fusion_explicit_rk_new.hpp"
#include "rk_performance_test_case.hpp"
#include "lorenz.hpp"
typedef boost::array< double , 3 > state_type;
typedef explicit_rk< state_type , 4 > rk4_fusion_type;
typedef rk4_fusion_type::coef_a_type coef_a_type;
typedef rk4_fusion_type::coef_b_type coef_b_type;
typedef rk4_fusion_type::coef_c_type coef_c_type;
const boost::array< double , 1 > a1 = {{ 0.5 }};
const boost::array< double , 2 > a2 = {{ 0.0 , 0.5 }};
const boost::array< double , 3 > a3 = {{ 0.0 , 0.0 , 1.0 }};
const coef_a_type a = fusion::make_vector( a1 , a2 , a3 );
const coef_b_type b = {{ 1.0/6 , 1.0/3 , 1.0/3 , 1.0/6 }};
const coef_c_type c = {{ 0.0 , 0.5 , 0.5 , 1.0 }};
class fusion_wrapper
{
public:
fusion_wrapper() : m_stepper( a , b , c )
{ }
void reset_init_cond()
{
m_x[0] = 10.0 * rand() / RAND_MAX;
m_x[1] = 10.0 * rand() / RAND_MAX;
m_x[2] = 10.0 * rand() / RAND_MAX;
m_t = 0.0;
}
inline void do_step( const double dt )
{
m_stepper.do_step( lorenz(), m_x , m_t , dt );
}
double state( const size_t i ) const
{ return m_x[i]; }
private:
state_type m_x;
double m_t;
rk4_fusion_type m_stepper;
};
int main()
{
srand( 12312354 );
fusion_wrapper stepper;
run( stepper );
}

71
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/generic_rk4_phase_lattice.cpp
View File

@@ -1,71 +0,0 @@
/*
* generic_rk4_phase_lattice.cpp
*
* Created on: May 15, 2011
* Author: mario
*/
#include <boost/array.hpp>
#include "../fusion_explicit_rk_new.hpp"
#include "rk_performance_test_case.hpp"
#include "phase_lattice.hpp"
const size_t N = 1024;
typedef boost::array< double , N > state_type;
typedef explicit_rk< state_type , 4 > rk4_fusion_type;
typedef rk4_fusion_type::coef_a_type coef_a_type;
typedef rk4_fusion_type::coef_b_type coef_b_type;
typedef rk4_fusion_type::coef_c_type coef_c_type;
const boost::array< double , 1 > a1 = {{ 0.5 }};
const boost::array< double , 2 > a2 = {{ 0.0 , 0.5 }};
const boost::array< double , 3 > a3 = {{ 0.0 , 0.0 , 1.0 }};
const coef_a_type a = fusion::make_vector( a1 , a2 , a3 );
const coef_b_type b = {{ 1.0/6 , 1.0/3 , 1.0/3 , 1.0/6 }};
const coef_c_type c = {{ 0.0 , 0.5 , 0.5 , 1.0 }};
class fusion_wrapper
{
public:
fusion_wrapper() : m_stepper( a , b , c )
{ }
void reset_init_cond()
{
for( size_t i = 0 ; i<N ; ++i )
m_x[i] = 2.0*3.1415927*rand() / RAND_MAX;
m_t = 0.0;
}
inline void do_step( const double dt )
{
m_stepper.do_step( phase_lattice<N>(), m_x , m_t , dt );
}
double state( const size_t i ) const
{ return m_x[i]; }
private:
state_type m_x;
double m_t;
rk4_fusion_type m_stepper;
};
int main()
{
srand( 12312354 );
fusion_wrapper stepper;
run( stepper , 10000 , 1E-6 );
}

97
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/generic_rk54ck.cpp
View File

@@ -1,97 +0,0 @@
/*
* generic_rk78.cpp
*
* Created on: Apr 29, 2011
* Author: mario
*/
#include <iostream>
#include <fstream>
#include <boost/array.hpp>
#include <boost/accumulators/accumulators.hpp>
#include <boost/accumulators/statistics.hpp>
#include <boost/timer.hpp>
#include "../fusion_explicit_error_rk.hpp"
#include "lorenz.hpp"
#define tab "\t"
using namespace std;
using namespace boost::accumulators;
typedef accumulator_set<
double , stats< tag::mean , tag::variance >
> accumulator_type;
ostream& operator<<( ostream& out , accumulator_type &acc )
{
out << boost::accumulators::mean( acc ) << tab;
// out << boost::accumulators::variance( acc ) << tab;
return out;
}
typedef boost::timer timer_type;
typedef boost::array< double , 3 > state_type;
typedef explicit_error_rk< state_type , 6 > rk54ck_fusion_type;
//typedef explicit_rk< state_type , 6 > rk54ck_fusion_type;
const size_t loops = 20;
int main( int argc , char **argv )
{
typedef rk54ck_fusion_type::coef_a_type coef_a_type;
typedef rk54ck_fusion_type::coef_b_type coef_b_type;
typedef rk54ck_fusion_type::coef_c_type coef_c_type;
const boost::array< double , 1 > a1 = {{ 0.2 }};
const boost::array< double , 2 > a2 = {{ 3.0/40.0 , 9.0/40 }};
const boost::array< double , 3 > a3 = {{ 0.3 , -0.9 , 1.2 }};
const boost::array< double , 4 > a4 = {{ -11.0/54.0 , 2.5 , -70.0/27.0 , 35.0/27.0 }};
const boost::array< double , 5 > a5 = {{ 1631.0/55296.0 , 175.0/512.0 , 575.0/13824.0 , 44275.0/110592.0 , 253.0/4096.0 }};
const coef_a_type a = fusion::make_vector( a1 , a2 , a3 , a4 , a5 );
const coef_b_type b = {{ 37.0/378.0 , 0.0 , 250.0/621.0 , 125.0/594.0 , 0.0 , 512.0/1771.0 }};
const coef_b_type b2 = {{ b[0]-2825.0/27648.0 , b[1]-0.0 , b[2]-18575.0/48384.0 , b[3]-13525.0/55296.0 , b[4]-277.0/14336.0 , b[5]-0.25 }};
//const coef_b_type b = {{ 2825.0/27648 , 0.0 , 18575.0/48384 , 13525.0/55296 , 277.0/14336 , 1.0/4 }};
//const coef_b_type b2 = {{ b[0]-37.0/378 , b[1]-0.0 , b[2]-250.0/621 , b[3]-125.0/594 , b[4]-0.0 , b[5]-512.0/1771 }};
const coef_c_type c = {{ 0.0 , 0.2 , 0.3 , 0.6 , 1.0 , 7.0/8 }};
rk54ck_fusion_type rk54ck( a , b , b2 , c );
//rk54ck_fusion_type rk54ck( a , b , c );
const size_t num_of_steps = 20000000;
double dt = 1E-10;
accumulator_type acc;
timer_type timer;
srand( 12312354 );
for( size_t n=0 ; n<loops ; ++n )
{
state_type x = {{ 10.0 * rand()/RAND_MAX , 10.0 * rand()/RAND_MAX , 10.0 * rand()/RAND_MAX }};
state_type x_err;
double t = 0.0;
timer.restart();
for( size_t i=0 ; i<num_of_steps ; ++i, t+=dt )
{
rk54ck.do_step( lorenz() , x , t , dt , x_err );
if( i % 1000 == 0 ) // simulated stepsize control
dt += (dt*1E-3*rand())/RAND_MAX - dt*5E-4;
}
acc( timer.elapsed() );
clog.precision( 15 );
clog.width( 20 );
clog << acc << " " << x[0] << tab << " " << x_err[0] << endl;
}
cout << acc << endl;
return 0;
}

76
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/generic_rk54ck_lorenz.cpp
View File

@@ -1,76 +0,0 @@
/*
* generic_rk54ck_lorenz.cpp
*
* Created on: May 12, 2011
* Author: mario
*/
#include <boost/array.hpp>
#include "../fusion_explicit_error_rk.hpp"
#include "rk_performance_test_case.hpp"
#include "lorenz.hpp"
typedef boost::array< double , 3 > state_type;
typedef explicit_error_rk< state_type , 6 > rk54ck_fusion_type;
typedef rk54ck_fusion_type::coef_a_type coef_a_type;
typedef rk54ck_fusion_type::coef_b_type coef_b_type;
typedef rk54ck_fusion_type::coef_c_type coef_c_type;
const boost::array< double , 1 > a1 = {{ 0.2 }};
const boost::array< double , 2 > a2 = {{ 3.0/40.0 , 9.0/40 }};
const boost::array< double , 3 > a3 = {{ 0.3 , -0.9 , 1.2 }};
const boost::array< double , 4 > a4 = {{ -11.0/54.0 , 2.5 , -70.0/27.0 , 35.0/27.0 }};
const boost::array< double , 5 > a5 = {{ 1631.0/55296.0 , 175.0/512.0 , 575.0/13824.0 , 44275.0/110592.0 , 253.0/4096.0 }};
const coef_a_type a = fusion::make_vector( a1 , a2 , a3 , a4 , a5 );
const coef_b_type b = {{ 37.0/378.0 , 0.0 , 250.0/621.0 , 125.0/594.0 , 0.0 , 512.0/1771.0 }};
const coef_b_type b2 = {{ b[0]-2825.0/27648.0 , b[1]-0.0 , b[2]-18575.0/48384.0 , b[3]-13525.0/55296.0 , b[4]-277.0/14336.0 , b[5]-0.25 }};
const coef_c_type c = {{ 0.0 , 0.2 , 0.3 , 0.6 , 1.0 , 7.0/8 }};
class fusion_wrapper
{
public:
fusion_wrapper() : m_stepper( a , b , b2 , c )
{ }
void reset_init_cond()
{
m_x[0] = 10.0 * rand() / RAND_MAX;
m_x[1] = 10.0 * rand() / RAND_MAX;
m_x[2] = 10.0 * rand() / RAND_MAX;
m_t = 0.0;
}
inline void do_step( const double dt )
{
m_stepper.do_step( lorenz() , m_x , m_t , dt , m_x_err );
//m_t += dt;
}
double state( const size_t i ) const
{ return m_x[i]; }
private:
state_type m_x;
state_type m_x_err;
double m_t;
rk54ck_fusion_type m_stepper;
};
int main()
{
srand( 12312354 );
fusion_wrapper stepper;
run( stepper );
}

66
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/gsl_rk4_lorenz.cpp
View File

@@ -1,66 +0,0 @@
/*
* gsl_rk4_lorenz.cpp
*
* Created on: May 11, 2011
* Author: mario
*/
#include <gsl/gsl_odeiv.h>
#include "rk_performance_test_case.hpp"
#include "lorenz_gsl.hpp"
const size_t dim = 3;
class gsl_wrapper
{
public:
gsl_wrapper()
{
m_s = gsl_odeiv_step_alloc( gsl_odeiv_step_rk4 , dim);
m_sys.function = lorenz_gsl;
m_sys.jacobian = 0;
m_sys.dimension = dim;
m_sys.params = 0;
}
void reset_init_cond()
{
m_x[0] = 10.0 * rand() / RAND_MAX;
m_x[1] = 10.0 * rand() / RAND_MAX;
m_x[2] = 10.0 * rand() / RAND_MAX;
m_t = 0.0;
}
inline void do_step( const double dt )
{
gsl_odeiv_step_apply ( m_s , m_t , dt , m_x , m_x_err , 0 , 0 , &m_sys );
//m_t += dt;
}
double state( const size_t i ) const
{ return m_x[i]; }
~gsl_wrapper()
{
gsl_odeiv_step_free( m_s );
}
private:
double m_x[dim];
double m_x_err[dim];
double m_t;
gsl_odeiv_step *m_s;
gsl_odeiv_system m_sys;
};
int main()
{
gsl_wrapper stepper;
run( stepper , 20000000 / 3 , 1E-10 * 3);
}

77
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/gsl_rk54ck.cpp
View File

@@ -1,77 +0,0 @@
/*
* gsl_rk54ck.cpp
*
* Created on: Apr 29, 2011
* Author: mario
*/
#include <iostream>
#include <boost/accumulators/accumulators.hpp>
#include <boost/accumulators/statistics.hpp>
#include <boost/timer.hpp>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>
#include "lorenz_gsl.hpp"
#define tab "\t"
using namespace std;
using namespace boost::accumulators;
typedef accumulator_set<
double , stats< tag::mean , tag::variance >
> accumulator_type;
ostream& operator<<( ostream& out , accumulator_type &acc )
{
out << boost::accumulators::mean( acc ) << tab;
// out << boost::accumulators::variance( acc ) << tab;
return out;
}
typedef boost::timer timer_type;
const size_t loops = 20;
int main()
{
const size_t num_of_steps = 20000000;
double dt = 1E-10;
accumulator_type acc;
timer_type timer;
srand( 12312354 );
gsl_odeiv_step *s = gsl_odeiv_step_alloc( gsl_odeiv_step_rkck , 3);
gsl_odeiv_system sys = { lorenz_gsl , 0 , 3 , 0 };
for( size_t n=0 ; n<loops ; ++n )
{
double x[3] = { 10.0 * rand()/RAND_MAX ,
10.0 * rand()/RAND_MAX ,
10.0 * rand()/RAND_MAX };
double x_err[3];
double t = 0.0;
timer.restart();
for( size_t i=0 ; i<num_of_steps ; ++i, t+=dt )
{
gsl_odeiv_step_apply ( s , t , dt , x , x_err , 0 , 0 , &sys );
if( i % 1000 == 0 ) // simulated stepsize control
dt += (dt*1E-3*rand())/RAND_MAX - dt*5E-4;
}
acc( timer.elapsed() );
clog.precision( 15 );
clog.width( 20 );
clog << acc << " " << x[0] << tab << " " << x_err[0] << endl;
}
cout << acc << endl;
return 0;
}

66
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/gsl_rk54ck_lorenz.cpp
View File

@@ -1,66 +0,0 @@
/*
* gsl_rk54ck_lorenz.cpp
*
* Created on: May 12, 2011
* Author: mario
*/
#include <gsl/gsl_odeiv.h>
#include "rk_performance_test_case.hpp"
#include "lorenz_gsl.hpp"
const size_t dim = 3;
class gsl_wrapper
{
public:
gsl_wrapper()
{
m_s = gsl_odeiv_step_alloc( gsl_odeiv_step_rkck , dim);
m_sys.function = lorenz_gsl;
m_sys.jacobian = 0;
m_sys.dimension = dim;
m_sys.params = 0;
}
void reset_init_cond()
{
m_x[0] = 10.0 * rand() / RAND_MAX;
m_x[1] = 10.0 * rand() / RAND_MAX;
m_x[2] = 10.0 * rand() / RAND_MAX;
m_t = 0.0;
}
inline void do_step( const double dt )
{
gsl_odeiv_step_apply ( m_s , m_t , dt , m_x , m_x_err , 0 , 0 , &m_sys );
//m_t += dt;
}
double state( const size_t i ) const
{ return m_x[i]; }
~gsl_wrapper()
{
gsl_odeiv_step_free( m_s );
}
private:
double m_x[dim];
double m_x_err[dim];
double m_t;
gsl_odeiv_step *m_s;
gsl_odeiv_system m_sys;
};
int main()
{
gsl_wrapper stepper;
run( stepper );
}

62
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/java/rk4.java
View File

@@ -1,62 +0,0 @@
//package org.apache.commons.math.ode;
//package org.apache.commons.math.ode.nonstiff;
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math.ode.FirstOrderIntegrator;
import org.apache.commons.math.ode.nonstiff.ClassicalRungeKuttaIntegrator;
import org.apache.commons.math.ode.DerivativeException;
import org.apache.commons.math.ode.IntegratorException;
public class rk4 {
private static class Lorenz implements FirstOrderDifferentialEquations
{
private static double sigma = 10.0;
private static double R = 28.0;
private static double b = 8.0/3.0;
public int getDimension()
{
return 3;
}
public void computeDerivatives( double t , double[] y , double[] yDot )
{
yDot[0] = sigma * ( y[1] - y[0] );
yDot[1] = R * y[0] - y[1] - y[0] * y[2];
yDot[2] = y[0]*y[1] - b * y[2];
}
}
public static void main( String[] args )
{
final int loops = 20;
final double dt = 1E-10;
final int steps = 20000000;
double mean = 0.0;
FirstOrderIntegrator rk4 = new ClassicalRungeKuttaIntegrator( dt );
final FirstOrderDifferentialEquations ode = new Lorenz();
for( int n=0 ; n<loops+3 ; n++ )
{
double[] y = new double[] { 0.0, 1.0 , 1.0 }; // initial state
long start = System.currentTimeMillis();
try {
rk4.integrate(ode, 0.0, y, dt*steps , y);
} catch(DerivativeException de) {
System.out.println("wrong exception caught");
} catch(IntegratorException ie) {
}
System.out.println( y[0] + " " + y[1] + " " + y[2] );
long end = System.currentTimeMillis();
System.out.println( "Time: " + (end-start) + " ms" );
if( n > 2 )
mean += end - start;
}
System.out.println("Average run time: " + (mean/loops) + " ms" );
}
}

41
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/lorenz.hpp
View File

@@ -1,41 +0,0 @@
/*
* lorenz.hpp
*
* Created on: May 12, 2011
* Author: mario
*/
#ifndef LORENZ_HPP_
#define LORENZ_HPP_
#include <boost/array.hpp>
struct lorenz
{
template< class state_type >
void inline operator()( const state_type &x , state_type &dxdt , const double t ) const
{
const double sigma = 10.0;
const double R = 28.0;
const double b = 8.0 / 3.0;
dxdt[0] = sigma * ( x[1] - x[0] );
dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
dxdt[2] = x[0]*x[1] - b * x[2];
}
};
typedef boost::array< double , 3 > state_type;
inline void lorenz_func( const state_type &x , state_type &dxdt , const double t )
{
const double sigma = 10.0;
const double R = 28.0;
const double b = 8.0 / 3.0;
dxdt[0] = sigma * ( x[1] - x[0] );
dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
dxdt[2] = x[0]*x[1] - b * x[2];
}
#endif /* LORENZ_HPP_ */

25
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/lorenz_gsl.hpp
View File

@@ -1,25 +0,0 @@
/*
* lorenz_gsl.hpp
*
* Created on: May 12, 2011
* Author: mario
*/
#ifndef LORENZ_GSL_HPP_
#define LORENZ_GSL_HPP_
#include <gsl/gsl_matrix.h>
int lorenz_gsl( const double t , const double y[] , double f[] , void *params)
{
const double sigma = 10.0;
const double R = 28.0;
const double b = 8.0 / 3.0;
f[0] = sigma * ( y[1] - y[0] );
f[1] = R * y[0] - y[1] - y[0] * y[2];
f[2] = y[0]*y[1] - b * y[2];
return GSL_SUCCESS;
}
#endif

79
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/nr_rk4_lorenz.cpp
View File

@@ -1,79 +0,0 @@
/*
* nr_rk4_lorenz.cpp
*
* Created on: May 11, 2011
* Author: mario
*/
#include <boost/array.hpp>
#include "rk_performance_test_case.hpp"
#include "lorenz.hpp"
const size_t dim = 3;
typedef boost::array< double , dim > state_type;
template< class System , typename T , size_t dim >
void rk4_step( const System sys , boost::array< T , dim > &x , const double t , const double dt )
{ // fast rk4 implementation adapted from the book 'Numerical Recipes'
size_t i;
const double hh = dt*0.5;
const double h6 = dt/6.0;
const double th = t+hh;
boost::array< T , dim > dydx , dym , dyt , yt;
sys( x , dydx , t );
for( i=0 ; i<dim ; i++ )
yt[i] = x[i] + hh*dydx[i];
sys( yt , dyt , th );
for( i=0 ; i<dim ; i++ )
yt[i] = x[i] + hh*dyt[i];
sys( yt , dym , th );
for( i=0 ; i<dim ; i++ ) {
yt[i] = x[i] + dt*dym[i];
dym[i] += dyt[i];
}
sys( yt , dyt , t+dt );
for( i=0 ; i<dim ; i++ )
x[i] += h6*( dydx[i] + dyt[i] + 2.0*dym[i] );
}
class nr_wrapper
{
public:
void reset_init_cond()
{
m_x[0] = 10.0 * rand() / RAND_MAX;
m_x[1] = 10.0 * rand() / RAND_MAX;
m_x[2] = 10.0 * rand() / RAND_MAX;
m_t = 0.0;
}
inline void do_step( const double dt )
{
rk4_step( lorenz() , m_x , m_t , dt );
}
double state( const size_t i ) const
{ return m_x[i]; }
private:
state_type m_x;
double m_t;
};
int main()
{
nr_wrapper stepper;
run( stepper );
}

80
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/nr_rk4_phase_lattice.cpp
View File

@@ -1,80 +0,0 @@
/*
* nr_rk4_phase_lattice.cpp
*
* Created on: May 15, 2011
* Author: mario
*/
#include <boost/array.hpp>
#include "rk_performance_test_case.hpp"
#include "phase_lattice.hpp"
const size_t dim = 1024;
typedef boost::array< double , dim > state_type;
template< class System , typename T , size_t dim >
void rk4_step( const System sys , boost::array< T , dim > &x , const double t , const double dt )
{ // fast rk4 implementation adapted from the book 'Numerical Recipes'
size_t i;
const double hh = dt*0.5;
const double h6 = dt/6.0;
const double th = t+hh;
boost::array< T , dim > dydx , dym , dyt , yt;
sys( x , dydx , t );
for( i=0 ; i<dim ; i++ )
yt[i] = x[i] + hh*dydx[i];
sys( yt , dyt , th );
for( i=0 ; i<dim ; i++ )
yt[i] = x[i] + hh*dyt[i];
sys( yt , dym , th );
for( i=0 ; i<dim ; i++ ) {
yt[i] = x[i] + dt*dym[i];
dym[i] += dyt[i];
}
sys( yt , dyt , t+dt );
for( i=0 ; i<dim ; i++ )
x[i] += h6*( dydx[i] + dyt[i] + 2.0*dym[i] );
}
class nr_wrapper
{
public:
void reset_init_cond()
{
for( size_t i = 0 ; i<dim ; ++i )
m_x[i] = 2.0*3.1415927*rand() / RAND_MAX;
m_t = 0.0;
}
inline void do_step( const double dt )
{
rk4_step( phase_lattice<dim>() , m_x , m_t , dt );
}
double state( const size_t i ) const
{ return m_x[i]; }
private:
state_type m_x;
double m_t;
};
int main()
{
srand( 12312354 );
nr_wrapper stepper;
run( stepper , 10000 , 1E-6 );
}

121
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/nr_rk54ck.cpp
View File

@@ -1,121 +0,0 @@
/*
* nr_rk4.cpp
*
* Created on: Apr 28, 2011
* Author: mario
*/
#include <iostream>
#include <fstream>
#include <boost/array.hpp>
#include <boost/numeric/odeint/stepper/runge_kutta4_classic.hpp>
#include <boost/accumulators/accumulators.hpp>
#include <boost/accumulators/statistics.hpp>
#include <boost/timer.hpp>
#include "lorenz.hpp"
#define tab "\t"
using namespace std;
using namespace boost::accumulators;
typedef accumulator_set<
double , stats< tag::mean , tag::variance >
> accumulator_type;
ostream& operator<<( ostream& out , accumulator_type &acc )
{
out << boost::accumulators::mean( acc ) << tab;
// out << boost::accumulators::variance( acc ) << tab;
return out;
}
typedef boost::timer timer_type;
typedef boost::array< double , 3 > state_type;
template< class System , typename T , size_t dim >
void rk54ck_step( System sys , boost::array< T , dim > &x , const double t , const double dt , boost::array< T , dim > &xerr )
{ // fast rk54ck implementation adapted from the book 'Numerical Recipes'
size_t i;
static const double a2 = 0.2 , a3 = 0.3 , a4 = 0.6 , a5 = 1.0 , a6 = 0.875 ,
b21 = 0.2 , b31 = 3.0 / 40.0 , b32 = 9.0 / 40.0 , b41 = 0.3 , b42 =
-0.9 , b43 = 1.2 , b51 = -11.0 / 54.0 , b52 = 2.5 , b53 =
-70.0 / 27.0 , b54 = 35.0 / 27.0 , b61 = 1631.0 / 55296.0 ,
b62 = 175.0 / 512.0 , b63 = 575.0 / 13824.0 , b64 = 44275.0
/ 110592.0 , b65 = 253.0 / 4096.0 , c1 = 37.0 / 378.0 , c3 =
250.0 / 621.0 , c4 = 125.0 / 594.0 , c6 = 512.0 / 1771.0 ,
dc1 = c1 - 2825.0 / 27648.0 , dc3 = c3 - 18575.0 / 48384.0 , dc4 =
c4 - 13525.0 / 55296.0 , dc5 = -277.00 / 14336.0 , dc6 = c6
- 0.25;
const size_t n = dim;
boost::array< T , dim > dydx , ak2 , ak3 , ak4 , ak5 , ak6 , ytemp;
sys( x , dydx , t );
for (i=0;i<n;i++)
ytemp[i] = x[i] + b21 * dt * dydx[i];
sys( ytemp , ak2 , t+a2*dt );
for (i=0;i<n;i++)
ytemp[i] = x[i] + dt*(b31*dydx[i]+b32*ak2[i]);
sys( ytemp , ak3 , t+a3*dt );
for (i=0;i<n;i++)
ytemp[i] = x[i] + dt*(b41*dydx[i]+b42*ak2[i]+b43*ak3[i]);
sys( ytemp , ak4 , t+a4*dt );
for (i=0;i<n;i++)
ytemp[i] = x[i] + dt*(b51*dydx[i]+b52*ak2[i]+b53*ak3[i]+b54*ak4[i]);
sys( ytemp, ak5 , t+a5*dt );
for (i=0;i<n;i++)
ytemp[i] = x[i] + dt*(b61*dydx[i]+b62*ak2[i]+b63*ak3[i]+b64*ak4[i]+b65*ak5[i]);
sys( ytemp , ak6 , t+a6*dt );
for (i=0;i<n;i++)
x[i] += dt*(c1*dydx[i]+c3*ak3[i]+c4*ak4[i]+c6*ak6[i]);
for (i=0;i<n;i++)
xerr[i] = dt*(dc1*dydx[i]+dc3*ak3[i]+dc4*ak4[i]+dc5*ak5[i]+dc6*ak6[i]);
}
const size_t loops = 20;
int main( int argc , char **argv )
{
const size_t num_of_steps = 20000000;
double dt = 1E-10;
accumulator_type acc;
timer_type timer;
srand( 12312354 );
for( size_t n=0 ; n<loops ; ++n )
{
state_type x = {{ 10.0 * rand()/RAND_MAX ,
10.0 * rand()/RAND_MAX ,
10.0 * rand()/RAND_MAX }};
state_type x_err;
double t = 0.0;
timer.restart();
for( size_t i=0 ; i<num_of_steps ; ++i, t+=dt )
{
rk54ck_step( lorenz() , x , t , dt , x_err );
if( i % 1000 == 0 ) // simulated stepsize control
dt += (dt*1E-3*rand())/RAND_MAX - dt*5E-4;
}
acc( timer.elapsed() );
clog.precision( 15 );
clog.width( 20 );
clog << acc << " " << x[0] << tab << " " << x_err[0] << endl;
}
cout << acc << endl;
return 0;
}

96
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/nr_rk54ck_lorenz.cpp
View File

@@ -1,96 +0,0 @@
/*
* nr_rk54ck_lorenz.cpp
*
* Created on: May 12, 2011
* Author: mario
*/
#include <boost/array.hpp>
#include "rk_performance_test_case.hpp"
#include "lorenz.hpp"
const size_t dim = 3;
typedef boost::array< double , dim > state_type;
template< class System , typename T , size_t dim >
void rk54ck_step( System sys , boost::array< T , dim > &x , const double t , const double dt , boost::array< T , dim > &xerr )
{ // fast rk54ck implementation adapted from the book 'Numerical Recipes'
size_t i;
static const double a2 = 0.2 , a3 = 0.3 , a4 = 0.6 , a5 = 1.0 , a6 = 0.875 ,
b21 = 0.2 , b31 = 3.0 / 40.0 , b32 = 9.0 / 40.0 , b41 = 0.3 , b42 =
-0.9 , b43 = 1.2 , b51 = -11.0 / 54.0 , b52 = 2.5 , b53 =
-70.0 / 27.0 , b54 = 35.0 / 27.0 , b61 = 1631.0 / 55296.0 ,
b62 = 175.0 / 512.0 , b63 = 575.0 / 13824.0 , b64 = 44275.0
/ 110592.0 , b65 = 253.0 / 4096.0 , c1 = 37.0 / 378.0 , c3 =
250.0 / 621.0 , c4 = 125.0 / 594.0 , c6 = 512.0 / 1771.0 ,
dc1 = c1 - 2825.0 / 27648.0 , dc3 = c3 - 18575.0 / 48384.0 , dc4 =
c4 - 13525.0 / 55296.0 , dc5 = -277.00 / 14336.0 , dc6 = c6
- 0.25;
const size_t n = dim;
boost::array< T , dim > dydx , ak2 , ak3 , ak4 , ak5 , ak6 , ytemp;
sys( x , dydx , t );
for (i=0;i<n;i++)
ytemp[i] = x[i] + b21 * dt * dydx[i];
sys( ytemp , ak2 , t+a2*dt );
for (i=0;i<n;i++)
ytemp[i] = x[i] + dt*(b31*dydx[i]+b32*ak2[i]);
sys( ytemp , ak3 , t+a3*dt );
for (i=0;i<n;i++)
ytemp[i] = x[i] + dt*(b41*dydx[i]+b42*ak2[i]+b43*ak3[i]);
sys( ytemp , ak4 , t+a4*dt );
for (i=0;i<n;i++)
ytemp[i] = x[i] + dt*(b51*dydx[i]+b52*ak2[i]+b53*ak3[i]+b54*ak4[i]);
sys( ytemp, ak5 , t+a5*dt );
for (i=0;i<n;i++)
ytemp[i] = x[i] + dt*(b61*dydx[i]+b62*ak2[i]+b63*ak3[i]+b64*ak4[i]+b65*ak5[i]);
sys( ytemp , ak6 , t+a6*dt );
for (i=0;i<n;i++)
x[i] += dt*(c1*dydx[i]+c3*ak3[i]+c4*ak4[i]+c6*ak6[i]);
for (i=0;i<n;i++)
xerr[i] = dt*(dc1*dydx[i]+dc3*ak3[i]+dc4*ak4[i]+dc5*ak5[i]+dc6*ak6[i]);
}
class nr_wrapper
{
public:
void reset_init_cond()
{
m_x[0] = 10.0 * rand() / RAND_MAX;
m_x[1] = 10.0 * rand() / RAND_MAX;
m_x[2] = 10.0 * rand() / RAND_MAX;
m_t = 0.0;
}
inline void do_step( const double dt )
{
rk54ck_step( lorenz() , m_x , m_t , dt , m_x_err );
}
double state( const size_t i ) const
{ return m_x[i]; }
private:
state_type m_x;
state_type m_x_err;
double m_t;
};
int main()
{
nr_wrapper stepper;
run( stepper );
}

57
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4_lorenz.cpp
View File

@@ -1,57 +0,0 @@
/*
* odeint_rk4_lorenz.cpp
*
* Created on: May 11, 2011
* Author: mario
*/
#include <boost/array.hpp>
#include <boost/numeric/odeint/stepper/runge_kutta4_classic.hpp>
#include <boost/numeric/odeint/algebra/array_algebra.hpp>
#include "rk_performance_test_case.hpp"
#include "lorenz.hpp"
typedef boost::array< double , 3 > state_type;
typedef boost::numeric::odeint::runge_kutta4_classic< state_type , double , state_type , double ,
boost::numeric::odeint::array_algebra > rk4_odeint_type;
class odeint_wrapper
{
public:
void reset_init_cond()
{
m_x[0] = 10.0 * rand() / RAND_MAX;
m_x[1] = 10.0 * rand() / RAND_MAX;
m_x[2] = 10.0 * rand() / RAND_MAX;
m_t = 0.0;
}
inline void do_step( const double dt )
{
m_stepper.do_step( lorenz() , m_x , m_t , dt );
//m_t += dt;
}
double state( const size_t i ) const
{ return m_x[i]; }
private:
state_type m_x;
double m_t;
rk4_odeint_type m_stepper;
};
int main()
{
srand( 12312354 );
odeint_wrapper stepper;
run( stepper );
}

54
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4_lorenz_def_alg.cpp
View File

@@ -1,54 +0,0 @@
/*
* odeint_rk4_lorenz_def_alg.cpp
*
* Created on: Apr 28, 2011
* Author: mario
*/
#include <boost/array.hpp>
#include <boost/numeric/odeint/stepper/runge_kutta4_classic.hpp>
#include <boost/numeric/odeint/algebra/array_algebra.hpp>
#include "rk_performance_test_case.hpp"
#include "lorenz.hpp"
typedef boost::array< double , 3 > state_type;
typedef boost::numeric::odeint::runge_kutta4_classic< state_type > rk4_odeint_type;
class odeint_wrapper
{
public:
void reset_init_cond()
{
m_x[0] = 10.0 * rand() / RAND_MAX;
m_x[1] = 10.0 * rand() / RAND_MAX;
m_x[2] = 10.0 * rand() / RAND_MAX;
m_t = 0.0;
}
inline void do_step( const double dt )
{
m_stepper.do_step( lorenz() , m_x , m_t , dt );
//m_t += dt;
}
double state( const size_t i ) const
{ return m_x[i]; }
private:
state_type m_x;
double m_t;
rk4_odeint_type m_stepper;
};
int main()
{
odeint_wrapper stepper;
run( stepper );
}

60
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4_phase_lattice.cpp
View File

@@ -1,60 +0,0 @@
/*
* odeint_rk4_phase_lattice.cpp
*
* Created on: May 15, 2011
* Author: mario
*/
#include <cmath>
#include <boost/array.hpp>
#include <boost/numeric/odeint/stepper/runge_kutta4_classic.hpp>
#include <boost/numeric/odeint/algebra/array_algebra.hpp>
#include "rk_performance_test_case.hpp"
#include "phase_lattice.hpp"
const size_t N = 1024;
typedef boost::array< double , N > state_type;
typedef boost::numeric::odeint::runge_kutta4_classic< state_type , double , state_type , double ,
boost::numeric::odeint::array_algebra > rk4_odeint_type;
class odeint_wrapper
{
public:
void reset_init_cond()
{
for( size_t i = 0 ; i<N ; ++i )
m_x[i] = 2.0*3.1415927*rand() / RAND_MAX;
m_t = 0.0;
}
inline void do_step( const double dt )
{
m_stepper.do_step( phase_lattice<N>() , m_x , m_t , dt );
//m_t += dt;
}
double state( const size_t i ) const
{ return m_x[i]; }
private:
state_type m_x;
double m_t;
rk4_odeint_type m_stepper;
};
int main()
{
srand( 12312354 );
odeint_wrapper stepper;
run( stepper , 10000 , 1E-6 );
}

83
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk54ck.cpp
View File

@@ -1,83 +0,0 @@
/*
* odeint_rk54ck.cpp
*
* Created on: Apr 29, 2011
* Author: mario
*/
#include <iostream>
#include <fstream>
#include <boost/array.hpp>
#include <boost/numeric/odeint/stepper/runge_kutta_cash_karp54_classic.hpp>
#include <boost/numeric/odeint/algebra/array_algebra.hpp>
#include <boost/accumulators/accumulators.hpp>
#include <boost/accumulators/statistics.hpp>
#include <boost/timer.hpp>
#include "lorenz.hpp"
#define tab "\t"
using namespace std;
using namespace boost::accumulators;
typedef accumulator_set<
double , stats< tag::mean , tag::variance >
> accumulator_type;
ostream& operator<<( ostream& out , accumulator_type &acc )
{
out << boost::accumulators::mean( acc ) << tab;
// out << boost::accumulators::variance( acc ) << tab;
return out;
}
typedef boost::timer timer_type;
typedef boost::array< double , 3 > state_type;
//typedef boost::numeric::odeint::explicit_error_rk54_ck< state_type > rk54_ck_odeint_type;
typedef boost::numeric::odeint::runge_kutta_cash_karp54_classic< state_type , double , state_type , double ,
boost::numeric::odeint::array_algebra > rk54_ck_odeint_type;
const size_t loops = 20;
int main( int argc , char **argv )
{
rk54_ck_odeint_type rk54_ck_odeint;
const size_t num_of_steps = 20000000;
double dt = 1E-10;
accumulator_type acc;
timer_type timer;
srand( 12312354 );
for( size_t n=0 ; n<loops ; ++n )
{
state_type x = {{ 10.0 * rand()/RAND_MAX ,
10.0 * rand()/RAND_MAX ,
10.0 * rand()/RAND_MAX }};
state_type x_err;
double t = 0.0;
timer.restart();
for( size_t i=0 ; i<num_of_steps ; ++i, t+=dt )
{
rk54_ck_odeint.do_step( lorenz() , x , t , dt , x_err );
if( i % 1000 == 0 ) // simulated stepsize control
dt += (dt*1E-3*rand())/RAND_MAX - dt*5E-4;
}
acc( timer.elapsed() );
clog.precision( 15 );
clog.width( 20 );
clog << acc << " " << x[0] << tab << " " << x_err[0] << endl;
}
cout << acc << endl;
return 0;
}

81
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk54ck_def_alg.cpp
View File

@@ -1,81 +0,0 @@
/*
* odeint_rk54ck.cpp
*
* Created on: Apr 29, 2011
* Author: mario
*/
#include <iostream>
#include <fstream>
#include <boost/array.hpp>
#include <boost/numeric/odeint/stepper/runge_kutta_cash_karp54_classic.hpp>
#include <boost/numeric/odeint/algebra/array_algebra.hpp>
#include <boost/accumulators/accumulators.hpp>
#include <boost/accumulators/statistics.hpp>
#include <boost/timer.hpp>
#include "lorenz.hpp"
#define tab "\t"
using namespace std;
using namespace boost::accumulators;
typedef accumulator_set<
double , stats< tag::mean , tag::variance >
> accumulator_type;
ostream& operator<<( ostream& out , accumulator_type &acc )
{
out << boost::accumulators::mean( acc ) << tab;
// out << boost::accumulators::variance( acc ) << tab;
return out;
}
typedef boost::timer timer_type;
typedef boost::array< double , 3 > state_type;
typedef boost::numeric::odeint::runge_kutta_cash_karp54_classic< state_type > rk54_ck_odeint_type;
const size_t loops = 20;
int main( int argc , char **argv )
{
rk54_ck_odeint_type rk54_ck_odeint;
const size_t num_of_steps = 20000000;
double dt = 1E-10;
accumulator_type acc;
timer_type timer;
srand( 12312354 );
for( size_t n=0 ; n<loops ; ++n )
{
state_type x = {{ 10.0 * rand()/RAND_MAX ,
10.0 * rand()/RAND_MAX ,
10.0 * rand()/RAND_MAX }};
state_type x_err;
double t = 0.0;
timer.restart();
for( size_t i=0 ; i<num_of_steps ; ++i, t+=dt )
{
rk54_ck_odeint.do_step( lorenz() , x , t , dt , x_err );
if( i % 1000 == 0 ) // simulated stepsize control
dt += (dt*1E-3*rand())/RAND_MAX - dt*5E-4;
}
acc( timer.elapsed() );
clog.precision( 15 );
clog.width( 20 );
clog << acc << " " << x[0] << tab << " " << x_err[0] << endl;
}
cout << acc << endl;
return 0;
}

57
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk54ck_lorenz.cpp
View File

@@ -1,57 +0,0 @@
/*
* odeint_rk54ck_lorenz.cpp
*
* Created on: May 12, 2011
* Author: mario
*/
#include <boost/array.hpp>
#include <boost/numeric/odeint/stepper/runge_kutta_cash_karp54_classic.hpp>
#include <boost/numeric/odeint/algebra/array_algebra.hpp>
#include "rk_performance_test_case.hpp"
#include "lorenz.hpp"
typedef boost::array< double , 3 > state_type;
typedef boost::numeric::odeint::runge_kutta_cash_karp54_classic< state_type , double , state_type , double ,
boost::numeric::odeint::array_algebra > rk54_odeint_type;
class odeint_wrapper
{
public:
void reset_init_cond()
{
m_x[0] = 10.0 * rand() / RAND_MAX;
m_x[1] = 10.0 * rand() / RAND_MAX;
m_x[2] = 10.0 * rand() / RAND_MAX;
m_t = 0.0;
}
inline void do_step( const double dt )
{
m_stepper.do_step( lorenz() , m_x , m_t , dt , m_x_err );
//m_t += dt;
}
double state( const size_t i ) const
{ return m_x[i]; }
private:
state_type m_x;
state_type m_x_err;
double m_t;
rk54_odeint_type m_stepper;
};
int main()
{
srand( 12312354 );
odeint_wrapper stepper;
run( stepper );
}

56
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk54ck_lorenz_def_alg.cpp
View File

@@ -1,56 +0,0 @@
/*
* odeint_rk54ck_lorenz.cpp
*
* Created on: May 12, 2011
* Author: mario
*/
#include <boost/array.hpp>
#include <boost/numeric/odeint/stepper/runge_kutta_cash_karp54_classic.hpp>
#include <boost/numeric/odeint/algebra/array_algebra.hpp>
#include "rk_performance_test_case.hpp"
#include "lorenz.hpp"
typedef boost::array< double , 3 > state_type;
typedef boost::numeric::odeint::runge_kutta_cash_karp54_classic< state_type > rk54_odeint_type;
class odeint_wrapper
{
public:
void reset_init_cond()
{
m_x[0] = 10.0 * rand() / RAND_MAX;
m_x[1] = 10.0 * rand() / RAND_MAX;
m_x[2] = 10.0 * rand() / RAND_MAX;
m_t = 0.0;
}
inline void do_step( const double dt )
{
m_stepper.do_step( lorenz() , m_x , m_t , dt , m_x_err );
//m_t += dt;
}
double state( const size_t i ) const
{ return m_x[i]; }
private:
state_type m_x;
state_type m_x_err;
double m_t;
rk54_odeint_type m_stepper;
};
int main()
{
srand( 12312354 );
odeint_wrapper stepper;
run( stepper );
}

51
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/perf_analysis.py
View File

@@ -1,51 +0,0 @@
from pylab import *
rcParams.update( { 'xtick.labelsize': 20 , 'ytick.labelsize': 20 })
results_rk4 = array(
[[ 0.82 , 0.59 , 0.60 , 0.60 , 0.54 , 0.54 , 0.81 , 0.98 , 1.14 ] , #odeint
[ 0.64 , 0.56 , 0.54 , 0.60 , 0.77 , 0.65 , 0.85 , 0.98 , 1.09 ] , #generic
[ 0.96 , 0.64 , 0.62 , 0.60 , 0.80 , 0.47 , 1.05 , 0.63 , 0.76 ] , #nr
[ 1.52 , 1.61 , 1.35 , 1.77 , 1.07 , 1.06 , 1.07 , 1.28 , 1.62 ] , #gsl
[ 1.85 , 1.85 , 1.47 , 2.08 , 1.08 , 1.08 , 1.57 , 1.70 , 2.10 ]]) #rt gen
results_rk54ck = array(
[[ 1.34 , 1.14 , 1.38 , 1.00 , 0.95 , 1.53 , 1.91 ] , #odeint
[ 1.26 , 1.13 , 1.43 , 1.28 , 1.02 , 1.63 , 1.84 ] , #generic
[ 2.04 , 1.32 , 1.34 , 1.48 , 1.10 , 1.82 , 1.15 ] , #nr
[ 2.80 , 2.34 , 2.79 , 1.80 , 1.94 , 1.99 , 2.16 ]]) #gsl
means_rk4 = 100*ones( 5 )
error_rk4 = zeros( 5 )
for i in arange(1,5):
tmp = results_rk4[0] / results_rk4[i]
means_rk4[i] = 100*mean( tmp )
error_rk4[i] = 100*sqrt(var( tmp ))
means_rk54ck = 100*ones( 4 )
error_rk54ck = zeros( 4 )
for i in arange(1,4):
tmp = results_rk54ck[0] / results_rk54ck[i]
means_rk54ck[i] = 100*mean( tmp )
error_rk54ck[i] = 100*sqrt(var( tmp ))
bar_width = 0.6
figure(1)
title("Runge-Kutta 4" , fontsize=20)
bar( arange(5) , means_rk4 , bar_width , color='blue' , linewidth=4 , edgecolor='blue' , yerr = error_rk4 , ecolor='red') #, elinewidth=2, ecolor='red' )
xlim( -0.5 , 4.5+bar_width )
xticks( arange(5)+bar_width/2 , ('odeint' , 'generic' , 'NR' , 'GSL' , 'rt gen' ) )
ylabel('Performance in %' , fontsize=20)
figure(2)
title("Runge-Kutta 5(4) Cash-Karp" , fontsize=20)
bar( arange(4) , means_rk54ck , bar_width , color='blue' , linewidth=4 , edgecolor='blue' , yerr = error_rk54ck , ecolor='red' )
xlim( -0.5 , 3.5+bar_width )
xticks( arange(4)+bar_width/2 , ('odeint' , 'generic' , 'NR' , 'GSL' ) )
ylabel('Performance in %' , fontsize=20)
show()

39
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/perf_tests.py
View File

@@ -1,39 +0,0 @@
from os import popen
from os import system
from os.path import isfile
#toolset = "gcc-4.4"
#toolset = "intel"
#toolset = "msvc"
toolset = "msvc-9.0express"
#bin_path = "bin/gcc-4.4/release/"
#bin_path = "bin/intel-linux/release/"
bin_path = "bin\\msvc-9.0express\\release\\threading-multi\\"
extension = ".exe"
#extension = ""
bins = [ "odeint_rk4_lorenz" , "odeint_rk4_lorenz_def_alg" , "generic_odeint_rk4_lorenz" , "generic_rk4_lorenz" , "nr_rk4_lorenz" , "gsl_rk4_lorenz" , "rt_generic_rk4_lorenz" ,
"odeint_rk54ck_lorenz" , "odeint_rk54ck_lorenz_def_alg" , "generic_rk54ck_lorenz" , "nr_rk54ck_lorenz" , "gsl_rk54ck_lorenz" ]
results = []
print "Performance tests for " , bin_path
print
for bin in bins:
system( "bjam toolset=" + toolset + " -a " + bin );
if isfile( bin_path + bin + extension):
print "Running" , bin
res = popen( bin_path+bin+extension ).read()
print bin , res
results.append( res )
else:
print "no executable found:" , bin_path + bin + extension
results.append( " -- " )
print "Results from" , bin_path
print
for i in range(len(bins)):
print bins[i] , results[i]

33
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/phase_lattice.hpp
View File

@@ -1,33 +0,0 @@
/*
*/
#include <cmath>
#include <boost/array.hpp>
template< size_t N >
struct phase_lattice
{
typedef double value_type;
typedef boost::array< value_type , N > state_type;
value_type m_epsilon;
state_type m_omega;
phase_lattice() : m_epsilon( 6.0/(N*N) ) // should be < 8/N^2 to see phase locking
{
for( size_t i=1 ; i<N-1 ; ++i )
m_omega[i] = m_epsilon*(N-i);
}
void inline operator()( const state_type &x , state_type &dxdt , const double t ) const
{
dxdt[0] = m_omega[0] + sin( x[1] - x[0] );
for( size_t i=1 ; i<N-1 ; ++i )
dxdt[i] = m_omega[i] + sin( x[i] - x[i-1] ) + sin( x[i+1] - x[i] );
dxdt[N-1] = m_omega[N-1] + sin( x[N-1] - x[N-2] );
}
};

41
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/results.dat
View File

@@ -1,41 +0,0 @@
Results for Runge-Kutta 4 with different implementations and compilers
| odeint | def alg | generic | nr | gsl | rt gen |
-------------------------------------------------------------------------------------
gcc 4.3.2 | 0.82 | 0.84 | 0.64 | 0.96 | 1.52 | 1.85 | Core2Quad Q9550 @ 2.83 GHz
gcc 4.4.1 | 0.59 | 0.59 | 0.56 | 0.64 | 1.61 | 1.85 | PhenomII X4 945 @ 3 GHz
gcc 4.4.1 | 0.60 | 0.60 | 0.54 | 0.62 | 1.35 | 1.47 | Core i7 930 @ 2.80 GHz
gcc 4.4.1 | 0.60 | 0.59 | 0.60 | 0.60 | 1.77 | 2.08 | Opteron 2224 @ 3 GHz
gcc 4.4.1 | 0.731 | 0.7315 | 0.729 | 0.7785 | 1.7755 | 2.8295 | Core 2 Duo P8400 @ 2.26GHz
gcc 4.4.1 | 0.7125 | 0.7095 | 0.7085 | 0.7555 | 2.0065 | 2.915 | Core 2 Quad Q6600 @ 2.40GHz
gcc 4.5.0 | 0.54 | 0.92 | 0.51 | 0.47 | 1.13 | 1.16 | Core i7 870 @ 2.93 GHz
gcc 4.6.0 | 0.54 | 0.54 | 0.65 | 0.47 | 1.06 | 1.08 | Core i7 870 @ 2.93 GHz
icc 12.0.2 | 0.90 | 0.81 | 0.85 | 1.05 | 1.07 | 1.57 | Core i7 870 @ 2.93 GHz
icc 11.1 | 1.11 | 0.98 | 0.98 | 0.63 | 1.28 | 1.70 | Xeon X5650 @ 2.67 GHz
msvc 9.0 | 6.91 | | 7.28 | 5.59 | ---- | 13.2 | Via Nano @ 1.60 GHz
msvc 9.0 | 1.957 | 2.2625 | 2.2953 | 1.2531 | ---- | 4.35155 | Core 2 Quad Q6600 @ 2.40GHz (32bit)
icc 11.1 | 1.24 | 1.14 | 1.09 | 0.76 | 1.62 | 2.10 | PhenomII X4 945 @ 3 GHz
Results for Runge-Kutta 54 Cash Karp (including simulated stepsize control)
| odeint | def alg | generic | nr | gsl |
-------------------------------------------------------------------------
gcc 4.3.2 | 1.34 | 1.42 | 1.26 | 2.04 | 2.80 | Core2Quad Q9550 @ 2.83 GHz
gcc 4.4.1 | 1.06 | 1.05 | 1.15 | | 2.63 | PhenomII X4 945 @ 3 GHz
gcc 4.4.1 | 1.14 | 1.14 | 1.13 | 1.32 | 2.34 | Core i7 930 @ 2.80 GHz
gcc 4.4.1 | 1.38 | 1.38 | 1.43 | 1.34 | 2.79 | Opteron 2224 @ 3 GHz
gcc 4.4.1 | 1.412 | 1.4135 | 1.4265 | 1.571 | 2.955 | Core 2 Duo P8400 @ 2.26GHz
gcc 4.4.1 | 1.37 | 1.387 | 1.424 | 2.287 | 4.3235 | Core 2 Quad Q6600 @ 2.40GHz
gcc 4.5.0 | 1.00 | 1.13 | 1.28 | 1.48 | 1.80 | Core i7 870 @ 2.93 GHz
gcc 4.6.0 | 0.95 | 0.95 | 1.02 | 1.10 | 1.94 | Core i7 870 @ 2.93 GHz
icc 12.0.2 | 1.53 | 1.82 | 1.63 | 1.82 | 1.99 | Core i7 870 @ 2.93 GHz
icc 11.1 | 1.91 | 1.93 | 1.84 | 1.15 | 2.16 | Xeon X5650 @ 2.67 GHz
gcc 4.4.1 | | | | | | Core2Quad Q6600 @ 2.40GHz
msvc 9.0 | 12.8 | | 15.4 | | ---- | Via Nano @ 1.60 GHz
msvc 9.0 | 4.21015 | 4.5484 | 4.79375 | 2.8133 | ---- | Core 2 Quad Q6600 @ 2.40GHz (32bit)
icc 11.1 | 1.80 | 2.44 | 2.39 | | 2.69 | PhenomII X4 945 @ 3 GHz
NOTE: additional Intel Compiler option: -inline-forceinline
see Jamfile

56
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/rk_performance_test_case.hpp
View File

@@ -1,56 +0,0 @@
/*
* rk_performance_test_case.hpp
*
* Created on: May 11, 2011
* Author: mario
*/
#include <iostream>
#include <boost/accumulators/accumulators.hpp>
#include <boost/accumulators/statistics.hpp>
#include <boost/timer.hpp>
#define tab "\t"
using namespace std;
using namespace boost::accumulators;
typedef accumulator_set<
double , stats< tag::mean , tag::variance >
> accumulator_type;
ostream& operator<<( ostream& out , accumulator_type &acc )
{
out << boost::accumulators::mean( acc ) << tab;
// out << boost::accumulators::variance( acc ) << tab;
return out;
}
typedef boost::timer timer_type;
template< class Stepper >
void run( Stepper &stepper , const size_t num_of_steps = 20000000 , const double dt = 1E-10 )
{
const size_t loops = 20;
accumulator_type acc;
timer_type timer;
srand( 12312354 );
for( size_t n=0 ; n<loops ; ++n )
{
stepper.reset_init_cond( );
timer.restart();
for( size_t i = 0 ; i < num_of_steps ; ++i )
stepper.do_step( dt );
acc(timer.elapsed());
clog.precision(3);
clog.width(5);
clog << acc << " " << stepper.state(0) << endl;
}
cout << acc << endl;
}

73
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/rt_generic_rk4_lorenz.cpp
View File

@@ -1,73 +0,0 @@
/*
* rt_generic_rk4_lorenz.cpp
*
* Created on: May 11, 2011
* Author: mario
*/
#include <boost/array.hpp>
#include "../rt_explicit_rk.hpp"
#include "rk_performance_test_case.hpp"
#include "lorenz.hpp"
typedef boost::array< double , 3 > state_type;
typedef rt_explicit_rk< state_type > rk_stepper_type;
const size_t stage_count = 4;
class rt_generic_wrapper
{
public:
rt_generic_wrapper() : m_stepper( stage_count )
{
rk_stepper_type::coeff_a_type a( stage_count-1 );
a[0].resize(1); a[0][0] = 0.5;
a[1].resize(2); a[1][0] = 0.0; a[1][1] = 0.5;
a[2].resize(3); a[2][0] = 0.0; a[2][1] = 0.0; a[2][2] = 1.0;
rk_stepper_type::coeff_b_type b( stage_count );
b[0] = 1.0/6; b[1] = 1.0/3; b[2] = 1.0/3; b[3] = 1.0/6;
rk_stepper_type::coeff_c_type c( stage_count );
c[0] = 0.0; c[1] = 0.5; c[2] = 0.5; c[3] = 1.0;
m_stepper.set_params( a , b , c );
}
void reset_init_cond()
{
m_x[0] = 10.0 * rand() / RAND_MAX;
m_x[1] = 10.0 * rand() / RAND_MAX;
m_x[2] = 10.0 * rand() / RAND_MAX;
m_t = 0.0;
}
inline void do_step( const double dt )
{
m_stepper.do_step( lorenz() , m_x , m_t , dt );
//m_t += dt;
}
double state( const size_t i ) const
{ return m_x[i]; }
private:
state_type m_x;
double m_t;
rk_stepper_type m_stepper;
};
int main()
{
rt_generic_wrapper stepper;
run( stepper );
}

76
libs/numeric/odeint/ideas/fusion_runge_kutta/performance/rt_generic_rk4_phase_lattice.cpp
View File

@@ -1,76 +0,0 @@
/*
* rt_generic_rk4_phase_lattice.cpp
*
* Created on: May 15, 2011
* Author: mario
*/
#include <boost/array.hpp>
#include "../rt_explicit_rk.hpp"
#include "rk_performance_test_case.hpp"
#include "phase_lattice.hpp"
const size_t N = 1024;
typedef boost::array< double , N > state_type;
typedef rt_explicit_rk< state_type > rk_stepper_type;
const size_t stage_count = 4;
class rt_generic_wrapper
{
public:
rt_generic_wrapper() : m_stepper( stage_count )
{
rk_stepper_type::coeff_a_type a( stage_count-1 );
a[0].resize(1); a[0][0] = 0.5;
a[1].resize(2); a[1][0] = 0.0; a[1][1] = 0.5;
a[2].resize(3); a[2][0] = 0.0; a[2][1] = 0.0; a[2][2] = 1.0;
rk_stepper_type::coeff_b_type b( stage_count );
b[0] = 1.0/6; b[1] = 1.0/3; b[2] = 1.0/3; b[3] = 1.0/6;
rk_stepper_type::coeff_c_type c( stage_count );
c[0] = 0.0; c[1] = 0.5; c[2] = 0.5; c[3] = 1.0;
m_stepper.set_params( a , b , c );
}
void reset_init_cond()
{
for( size_t i = 0 ; i<N ; ++i )
m_x[i] = 2.0*3.1415927*rand() / RAND_MAX;
m_t = 0.0;
}
inline void do_step( const double dt )
{
m_stepper.do_step( phase_lattice<N>() , m_x , m_t , dt );
//m_t += dt;
}
double state( const size_t i ) const
{ return m_x[i]; }
private:
state_type m_x;
double m_t;
rk_stepper_type m_stepper;
};
int main()
{
srand( 12312354 );
rt_generic_wrapper stepper;
run( stepper , 10000 , 1E-6 );
}

Some files were not shown because too many files have changed in this diff Show More