BS stepper: correct exponents for optimal h

Becoming suspicious by the difference of the exponents used for computing the new step size in the BS and BS denseout stepper (see 0f943fb) I checked again the Hairer book and I'm now convinced there was a mistake in our implementation and both steppers should use 1/(2*k+1) as exponent. The background is the this exponent represents the order of the error of the k-th iteration, and this order is always 2k+1, independent of the interval_sequence. This error is computed from the difference of the k-th and k-1 - th iteration, which have the orders 2k+2 2k respectively, which means the computed error has order 2k+1.
headmyshoulder · Dec 28, 2015 · cc9b196 · cc9b196
1 parent ae53eb1
commit cc9b196
Show file tree

Hide file tree

Showing 2 changed files with 3 additions and 5 deletions.
diff --git a/include/boost/numeric/odeint/stepper/bulirsch_stoer.hpp b/include/boost/numeric/odeint/stepper/bulirsch_stoer.hpp
@@ -113,7 +113,7 @@ class bulirsch_stoer {
             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 , 1.0 / static_cast< value_type > ( 2*i+1 ) );
+            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] );

diff --git a/include/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp b/include/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp
@@ -128,13 +128,11 @@ class bulirsch_stoer_dense_out {
             if( i == 0 )
             {
                 m_cost[i] = m_interval_sequence[i];
-                m_facmin_table[0] = 1; // never to be used
-            }
-            else
+            } else
             {
-                m_facmin_table[i] = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , static_cast<value_type>(1)/(m_interval_sequence[i-1]) );
                 m_cost[i] = m_cost[i-1] + m_interval_sequence[i];
             }
+            m_facmin_table[i] = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , static_cast< value_type >(1) / static_cast< value_type >( 2*i+1 ) );
             m_coeff[i].resize(i);
             for( size_t k = 0 ; k < i ; ++k  )
             {