3
$\begingroup$

I am implementing adaptive step-size Runge–Kutta methods using embedded pairs (orders ( p ) and ( p+m )). The standard structure is:

  • Compute two solutions: $$ y^{[p]} \text{ and } y^{[p+m]} $$

  • Estimate local error: $$ \text{err} = \left| y^{[p+m]} - y^{[p]} \right|$$

  • Compute scaling factor:
    $$ s = \beta \left( \frac{\text{tol}}{\text{err}} \right)^{1/p} $$

  • Update step:
    $$ h_{\text{new}} = s \cdot h $$

My question is about the accept/reject criterion:

Two alternatives I found:

  1. Compare using ( s ):
if (s >= 1) accept step
else reject
  1. Compare using the error directly:
if (err <= tol) accept step
else reject

From what I understand, these should be mathematically equivalent if ( s ) is defined as above.

However:

  • In Shampine-style implementations, it seems the comparison is done directly on the error.
  • In some lecture notes and sources, the acceptance is expressed using: $$ s \ge 1 $$

My concern (practical behavior):

I tested both approaches on the double pendulum problem, counting rejected steps. Using:

  • $$ \text{tol} = 10^{-5} $$
  • $$ h_{\max} = 5 \times 10^{-3} $$

I observed:

  • Using ( s $\ge$ 1 ):

    • ~183,000 rejected steps
    • ~6,000 accepted steps

  • Using ( $ \text{err} \le \text{tol} $):

    • < 100 rejected steps
    • Same number of accepted points (~6,000)

Both solutions remain within the requested tolerance band.

Questions:

  1. Are these two acceptance criteria strictly equivalent in practice, or can numerical/implementation details make them behave differently?
  2. Is one approach considered more robust or standard in modern implementations?
  3. Are there known situations (e.g., stiff/chaotic systems) where comparing via ( s ) is preferable?

If anyone can clarify whether there is a theoretical or practical reason to prefer one over the other, I would really appreciate it.

Thanks!

Improve this question
$\endgroup$
1
  • 1
    $\begingroup$ The two methods are only identical if $\beta=1$. $\endgroup$ Commented Apr 3 at 21:08

1 Answer 1

Reset to default
5
$\begingroup$

The standard approach in most ODE solver packages is the following. First compute the error using a norm weighted by tolerances: $$ err = \sqrt{\frac{1}{d} \sum_{i=1}^{d} \left( \frac{y_i^{[p+m]} - y_i^{[p]}}{atol_i + rtol_i |y_i^{[p]}|} \right)^2} $$

Here, $atol$ is an absolute tolerance and $rtol$ is a relative tolerance. To be consistent with your criteria 2, you can take $atol_i = tol, rtol_i = 0$.

If $err \leq 1$, the step is accepted, and if $err > 1$, the step is rejected. In either case, the time step is updated by

$$ h_{new} = \beta \, h \, err^{-1/(\min(p+m, p) + 1)} $$

Note the extra +1 in the exponent that the formula in your question was missing. Typically there's some bounds placed on the growth or shrinking of the step too. I recommend reading pages 167-168 of https://doi.org/10.1007/978-3-540-78862-1

So to answer your questions

Are these two acceptance criteria strictly equivalent in practice

Since $\beta < 1$, they're not quite equivalent. Criteria 1 seems dubious if not outright misguided because you cannot shrink the time step without rejecting the step. With criteria 2 you can. In the test results you report, this likely explains why criteria 1 was significantly worse.

Is one approach considered more robust or standard in modern implementations

Criteria 1 is not robust, 2 is somewhat robust, and what I provided above most robust.

Are there known situations (e.g., stiff/chaotic systems) where comparing via ( s ) is preferable?

I struggle to see any situation where criteria 1 is preferable.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.