I'm developing my own generic Runge-Kutta solver, and I'm currently implementing the adaptive step-size routine. I say generic because I want to be able to test different RK implementations by only passing the solver a Butcher tableau, such as the following.

\begin{array}{c|cccc} c_{1}&a_{11}&a_{12}&\dots &a_{1s}\\ c_{2}&a_{21}&a_{22}&\dots &a_{2s}\\ \vdots &\vdots &\vdots &\ddots &\vdots \\c_{s}&a_{s1}&a_{s2}&\dots &a_{ss}\\ \hline &b_{1}&b_{2}&\dots &b_{s}\\ &b_{1}^*&b_{2}^*&\dots &b_{s}^*\\ \end{array}

Currently, I'm only focused on adaptive-explicit methods. As far as I understand it, a non-adaptive method would only have only the first row of $b$, and an explicit method would be zero on the diagonal and upper triangle.

Essentially what I'm doing within a single iteration of my solver is passing a timestep span [$t_0, t_1$], and initial guess of the step ($h$), and an initial value of $y_0$, which could be a vector, and the function $\frac{dy}{dt} = f(t, y)$.

I calculate $y(t_0 + h)$ and $y^*(t_0 + h)$ using coefficients of the Butcher tableau, calculate the element-wise error between them, and then calculate a scaling factor of the step using the maximum absolute error. The scaling factor is used to determine if $h$ is too small, a good size, or too big, and is adjusted accordingly. Currently, I'm calculating the scaling factor like this (described in this youtube video):

$$ s = \frac{\epsilon \cdot h}{2 \cdot (t_1 - t_0) \cdot max(|y-y^*|)}$$

where $t_0$ and $t_1$ are the beginning and end of the time step of interest, and epsilon is some error parameter set previously (e.g. 1e-8).

The resulting value of $h$ would be scaled using $s$: $$h = sh$$

If s is larger than 2, then I assume the error between $y(t_0 + h)$ and $y^*(t_0 + h)$ is small and I can continue with the next step. If s is less than 1, then I adjust h without moving forward and try the same timestep again. If $s$ is between 1 and 2, then I move on as well, but I don't adjust h because it's good enough.

What as not readily obvious to me at the start, but now becomes so as the $tspan$ increases ($t_1-t_0$ goes up), is that the acceptable h is kept really far down. This also seems to be greatly influenced by the RK method used.

For example, the following system takes about 2 seconds on my laptop with 100 pts between [0, $4\pi$]:

\begin{array} {rl} y_1'(t) &= y_2 \\ y_2'(t) &= -y_1 \\ y_1(0) &= 0 \\ y_2(0) &= 0.1 \\ \end{array}

Using 100 points

Decreasing the number of sampling points to 10 within the same range ([0, $4\pi$]) now takes 16 seconds because I think the value of s is too small, due to the $t_1-t_0$ term.

It may be that I'm only interested in a very few number of points (or even the final value), but this is causing the solution time to increase, even for this relatively small problem.

Using 10 points

Can someone point me to a reference that goes into more detail about this scaling ($s$) parameter? Or is there another way to solve this approach this problem?

  • 3
    $\begingroup$ For the theory, i suggest the classic ODE text by Hairer and Wanner. For the practical implementation details, have a look at the code KPP, downloadable from here. The most relevant bits for you are the Runge-Kutta integrators that you can find under, e.g., /int/kpp_dvode.f (disclaimer: the code was developed by my former PhD adviser). $\endgroup$ – GoHokies Nov 8 '16 at 19:41
  • $\begingroup$ I have implemented a generic integrator using a specified Butcher Tableau, just as you are aiming to do. I use the following formula typically: $h = \gamma h (\frac{\tau}{e})^{1/p}$, where $\tau$ is the tolerance, $ 0 \le \gamma \le 1$ (0.9 in my code), $e$ is the error you compute between solutions, and $p$ is the higher order of accuracy of the two schemes. There are PID control based approaches to refining step size that you could investigate. $\endgroup$ – spektr Nov 8 '16 at 22:32
  • 2
    $\begingroup$ The approach you're describing is still pretty standard, but it has been superseded by much better control algorithms. These have yet to penetrate the wider community and production codes, but if you want to do it right, read Gustaf Soderlind's papers on step control. $\endgroup$ – David Ketcheson Nov 9 '16 at 5:54
  • $\begingroup$ I mean this, although there is even more recent work. $\endgroup$ – David Ketcheson Nov 9 '16 at 6:00
  • $\begingroup$ @DavidKetcheson thanks for the links - very interesting! $\endgroup$ – GoHokies Nov 9 '16 at 13:13

Your Answer

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

Browse other questions tagged or ask your own question.