Dormand–Prince 5(4): How to update the stepsize and make accept/reject decision?

Question

https://en.wikipedia.org/wiki/Dormand–Prince_method

I want to implement the Dormand-Prince 4(5) version to solve Initial Value problems. Using regular notation I have $A$ matrix and the $c,b,\hat{b}$ vectors.

I also have the error tolerence of $ETOL$ absolute error tolerence $ATOL$ and relative error tolerance $RTOL$

At each time step I compute $x_n$ by use of $b$ and $\hat{x}_n$ by us of $\hat{b}$ and estimate the error by as the difference $e_n=x_n-\hat{x}_n$

The way I understand this method, when the error is approximated then I have to do the following to things:

• Decide whether to accept the step just computed or not
• update the stepsize

Which exact decision rules/algorithm should I use to do the those steps?

the book I am using (Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations by Petzold and Ascher) does not explain this part.

Side Note:

And I assume there is a way to update the stepsize regarding what kind of RUnge Kutta method you use.

Answer 1

The basic idea is that

1. You use the estimated error given to you (cheaply) by the embedded methods;
2. You use a metric to define acceptance using a user-defined relative and absolute tolerance;
3. Based on the order properties of the code and this metric, a new step size is computed;
4. You avoid large differences in step sizes from one step to the next by limiting the maximum increase and maximum decrease.

From "Solving Ordinary Differential Equations, part I" by Hairer, Norsett and Wanner :

For two approximations to the solution $y_{1i}$ (using the lower order) and $\hat{y}_{1i}$ (using the higher order), we want the error to satisfy (componentwise) $$|y_{1i} - \hat{y}_{1i}| \leq sc_i|, \quad sc_i = Atol_i + \max\left(|y_{0i}|,|y_{1i}|\right)\cdot Rtol_i.$$

As a measure of the error, one takes $$err = \sqrt{\frac{1}{n}\sum_{i=1}^n\left(\frac{y_{1i}-\hat{y}_{1i}}{sc_i}\right)^{2}}$$ and this value is compared to one. From the expected error behaviour $err \propto C\cdot h^{q+1}$ the optimal step size is obtained as $$h_{\text{opt}} = h\cdot\left(\frac{1}{err}\right)^{\frac{1}{q+1}}$$ where $h$ is the previous step size.

Typically, one would add a kind of safety factor so that the probability of having selected a good step size is increased. Typically this safety factor $fac$ is taken $fac = 0.8$ or some power of the order like $fac = \left(0.25\right)^{\frac{1}{q+1}}$ so that $$h_{\text{opt}} = fac\cdot h\cdot\left(\frac{1}{err}\right)^{\frac{1}{q+1}}$$ .

As a second measure of safety, one wants to prohibit the stepsize to vary too much from one step to the next. Typically, one limits the new stepsize to be less than a factor five of the previous stepsize. And in some codes, a stepsize increase is forbidden if the previous step was rejected (so if in a previous step, $h$ was decreased, it can only increase again if a step with that step size is accepted).

The initial step size for the algorithm can be set by the user or it could be estimated using a procedure based on the order of the method. Hindmarsch and Shampine have developed methods to do so.