2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ I kindly ask people not to use use Scientific papers when answering since I have very limited experience in reading these papers/articles and I have already looked at some paper to find answer to this problem $\endgroup$ – k.dkhk May 2 at 7:45
  • 1
    $\begingroup$ There are simple methods, but there are also a lot of options as this is still an active area of research. It's hard to know which method you are looking for. What is your goal? $\endgroup$ – David Ketcheson May 2 at 12:21
  • $\begingroup$ @DavidKetcheson My goal is to replicate the results from Matllabs ode45 and SciPy's rk45. I get that there might be more than one option, in that case I am interested in the most common way for dopri54 (what is the consensus). I could have asked the question. If you needed to implement this dopir5, how would YOU have decided the step sizes $\endgroup$ – k.dkhk May 3 at 17:30
3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Do note that this is close to but not the same as the algorithm used in dopri5. $\endgroup$ – Chris Rackauckas May 3 at 7:17
  • $\begingroup$ Thank you very much. I will implement this method and comment if I need more. This was more or less what I was looking for! $\endgroup$ – k.dkhk May 3 at 17:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.