Existing implementations

Implementations

general (classical):

• error
  • val = || err / sk ||
  • with || x || is some norm
    • || x || = sqrt( 1 / n sum_i x_i^2 )
    • || x || = max_i |x_i|
  • sk is the scale function
    • sk = abs_tol + rel_tol | xold_i |
    • sk = abs_tol + rel_tol max( | xold_i | , | x_i | )
    • sk = abs_tol + rel_tol ( a_x | xold_i | + a_dxdt | dxdtold_i | )
• acceptance
  • val < 1 : accept
  • val > 1 : reject
• new step size
  • ...

controlled_runge_kutta:

• error:
  • val = max_i( | err[i] | / ( eps_abs + eps_rel * ( a_x | xold[i] | + a_dxdt | dxdtold[i] | ) )
• acceptance
  • val > 1 reject
  • val < 1 accept
• new step size
  • val > 1 : dt = dt * max( 0.9 * pow( val , -1 / ( error_order - 1 ) ) , 0.2 );
  • val < 0.5 : dt = dt * min( 0.9 * pow( val , -1 / stepper_order ) , 5 );
  • else : dt = dt
• stepper requirements
  • error stepper, out of place, explicit (dxdt enters do_step)
    • m_stepper.do_step( system , in , dxdt , t , out , dt , m_xerr.m_v );
  • error stepper fsal, output of place, explicit (if step is reject, new derivative needs to be rejected)
    • m_stepper.do_step( system , in , dxdt_in , t , out , dxdt_out , dt , m_xerr.m_v );

rosenbrock4:

• error:
  • val = sqrt( 1 / n * sum_i( err[i]^2 / ( eps_abs + eps_rel * max( | x[i] | , | xold[i] | ) )^2 )
• acceptance
  • val > 1 reject*
  • val < 1 accept
• new step size
  • safe = 0.9 , fac1 = 5 , fac2 = 1 / 6
  • fac = max( fac2 , min( fac1 , pow( val , 0.25 ) / safe ) )
  • val > 1 : complicated, see code below
  • val < 1 : dt = dt / fac
• stepper requirements:
  • error stepper, out of places
    • m_stepper.do_step( sys , x , t , xout , dt , m_xerr.m_v );

numercal recipes c++ (dopri5): (check code)

numerical recipes c++ (dopri853): (check code from webnote 20)

numerical recipes c (rk54ck): (check code nr_c)

numerical recipes c (rk54ck): (check code nr_c++)

gsl:

hairer page 168:

• error
  • val = sqrt( 1 / n * sum_i( err[i]^2 / ( eps_abs + eps_rel * max( | x[i] | , | xold[i] | ) )^2 )
• acceptance:
  • val > 1 reject
  • val < 1 accept
• new step size
  • dt = dt min( facmax , max( facmin , fac * (1/val)^(1/(q+1)) ) )
    • q = min( error_order , stepper_order )
    • fac, safety factor, 0.8, 0.9, 0.25^(1/(q+1)), 0.38^(1/(q+1))
    • facmax, maximal step increase (1.5-5)
      • 1 after an rejected step
    • facmin, minimal step decrease (0.2-0.6)

hairer page

Code

controlled_runge_kutta code:

`val = max_i( | err[i] | / ( eps_abs + eps_rel * ( a_x | xold[i] | + a_dxdt | dxdtold[i] | ) )

if( val > 1.0 )
{
    dt = dt * max( 0.9 * pow( val , -1.0 / ( m_stepper.error_order() - 1.0 ) ) , 0.2 );
    reject
}
else if( val < 0.5 )
{
    dt = dt * min( 0.9 * pow( val , -1.0 / m_stepper.stepper_order() ) , 5.0 );
    success;
}
else ( if( val > 0.5 && val < 1.0 )
{
    success;
}`

rosenbrock4 code:

`err = sqrt( sum_i( err[i] * err[i] / ( eps_abs + eps_rel * max( | x[i] | , | xold[i] | ) ) / n )

safe = 0.9 , fac1 = 5.0 , fac2 = 1.0 / 6.0
fac = max( fac2 , min( fac1 , pow( err , 0.25 ) / safe ) )
dt_new = dt / fac;
if ( err <= 1.0 )
{
    if( m_first_step )
    {
        m_first_step = false;
    }
    else
    {
        value_type fac_pred = ( m_dt_old / dt ) * pow( err * err / 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 , err );
    if( m_last_rejected )
        dt_new = ( dt >= 0.0 ? std::min( dt_new , dt ) : std::max( dt_new , dt ) );
    dt = dt_new;
    success;
}
else
{
    dt = dt_new;
    fail;
}`