summaryrefslogtreecommitdiff
path: root/boost/numeric/odeint/stepper/bulirsch_stoer.hpp
diff options
context:
space:
mode:
Diffstat (limited to 'boost/numeric/odeint/stepper/bulirsch_stoer.hpp')
-rw-r--r--boost/numeric/odeint/stepper/bulirsch_stoer.hpp23
1 files changed, 12 insertions, 11 deletions
diff --git a/boost/numeric/odeint/stepper/bulirsch_stoer.hpp b/boost/numeric/odeint/stepper/bulirsch_stoer.hpp
index 4b908333ba..02c37492b9 100644
--- a/boost/numeric/odeint/stepper/bulirsch_stoer.hpp
+++ b/boost/numeric/odeint/stepper/bulirsch_stoer.hpp
@@ -98,6 +98,7 @@ public:
m_interval_sequence( m_k_max+1 ) ,
m_coeff( m_k_max+1 ) ,
m_cost( m_k_max+1 ) ,
+ m_facmin_table( 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 )
{
@@ -112,21 +113,14 @@ public:
else
m_cost[i] = m_cost[i-1] + m_interval_sequence[i];
m_coeff[i].resize(i);
+ m_facmin_table[i] = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , static_cast< value_type >(1) / static_cast< value_type >( 2*i+1 ) );
for( size_t k = 0 ; k < i ; ++k )
{
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
}
-
- // crude estimate of optimal order
-
- m_current_k_opt = 4;
- /* no calculation because log10 might not exist for value_type!
- const value_type logfact( -log10( max BOOST_PREVENT_MACRO_SUBSTITUTION( eps_rel , static_cast< value_type >(1.0E-12) ) ) * 0.6 + 0.5 );
- m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>( 1 ) , min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>( m_k_max-1 ) , logfact ));
- */
}
-
+ reset();
}
@@ -282,7 +276,7 @@ public:
{
m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<int>(m_k_max-1) , static_cast<int>(m_current_k_opt)+1 );
new_h = h_opt[k];
- new_h *= m_cost[m_current_k_opt]/m_cost[k];
+ new_h *= 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;
@@ -344,6 +338,12 @@ public:
{
m_first = true;
m_last_step_rejected = false;
+ // crude estimate of optimal order
+ m_current_k_opt = 4;
+ /* no calculation because log10 might not exist for value_type!
+ const value_type logfact( -log10( max BOOST_PREVENT_MACRO_SUBSTITUTION( eps_rel , static_cast< value_type >(1.0E-12) ) ) * 0.6 + 0.5 );
+ m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>( 1 ) , min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>( m_k_max-1 ) , logfact ));
+ */
}
@@ -417,7 +417,7 @@ private:
BOOST_USING_STD_MAX();
using std::pow;
value_type expo( 1.0/(2*k+1) );
- value_type facmin = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , expo );
+ value_type facmin = m_facmin_table[k];
value_type fac;
if (error == 0.0)
fac=1.0/facmin;
@@ -506,6 +506,7 @@ private:
int_vector m_interval_sequence; // stores the successive interval counts
value_matrix m_coeff;
int_vector m_cost; // costs for interval count
+ value_vector m_facmin_table; // for precomputed facmin to save pow calls
state_table_type m_table; // sequence of states for extrapolation