An adaptive step size solver for an ODE

Question

I am trying solve an ordinary differential equation numerically, $\,dy/dt=10e^{-(t-2)^2/2(0.075)^2}-0.6y\,\,$ with an initial value and initial step size between $t=0\, and\, t=4$. In my code I implement an embedded Runge-Kutta formulation: Calculating two tentative estimates using RK2 and RK3 formulae with the same step size and the difference between them is the error for RK2. Then I compare this difference (will simply called 'error') with a value (an acceptable error) which actually is not a set number but is obtained through $10^{-3}*\lvert y_{current}\rvert$ or $10^{-6}$ whichever is larger ($y_{current}$ is the present estimate) If the error is greater than the acceptable error I update step size with $h_{new}=h_{old}(\frac{acceptable\, error}{error})^{\frac{1}{3}} $ and if the error is less than or equal to the acceptable error, the latest step size $h_{old}\,$is used to estimate the function value and the independent variable t is increased by that amount. The issue starts after that: I want to increase the step size and do not want to proceed with $h_{old}\,$ for the next set of RK calculations after establishing $y$. What are convenient ways of increasing my step size? Would using the formula $h_{new}=h_{old}(\frac{acceptable\, error}{error})^{\frac{1}{3}} $in reverse help?

Answer 1

Assume that you compute the step update resulting in $v_2,v_3,err$ so that $y_{k+1}=y_k+hv_2$ is the second order RK2 step, $\tilde y_{k+1}=y_k+hv_3$ is the third order RK3 step, and $y_{k+1}-\tilde y_{k+1}=h(v_2-v_3)=h\cdot err$ is the estimated error of the second order step.

Then $err=v_2-v_3$ is the error per unit step, or the error density. For consistent results this is the quantity to compare against the error tolerance. The optimal step size for the next step in the second-order method or a repetition of the current step is then $$ h_{opt}=\left(\frac{\max(atol,|y_k|\cdot rtol)}{|err|}\right)^{\frac12}\cdot h. $$ This is a best guess, one may reduce this slightly to increase the probability that the actual error growth remains below the level set by the tolerances.

One can use the 3rd order step with the step sizes for the 2nd order method, this will in general reduce the actual error to the level $(tol)^{\frac32}$. In the case of stiff equations there may be stability issues with this approach.

One can try to guess an optimal step size for the 3rd order method using $h\cdot err$ as proxy of its unit step error. This gives the formula $$ h_{opt}=\left(\frac{\max(atol,|y_k|\cdot rtol)}{|h\cdot err|}\right)^{\frac13}\cdot h. $$ This seems to be the formula that you are using. To repeat, this only works when using the RK3 step for the next value.

There is nothing to do in the case that a step is accepted, as if the error is smaller than the tolerance the first factor in the formula is larger than 1, so the step size increases.

There is of course a lot of uncertainty in this extrapolation approach, the error is not entirely correct for the method, the step size is only optimal for the last step. This might lead to chaotic oscillations in the sequence of the step sizes. Some kind of running average on the step sizes or the error quotients, or some other kind of dampening can improve the integration performance.