Adaptive Step Size in RK45 for Second-Order ODE

Question

My question pertains to the answer given in this post, but I am implementing the RK45 or RKF45 algorithm. Following the explanation in that post for a second order ODE, as an example I am looking at a pendulum whose equation of motion is $$\ddot{\theta} = -\Omega_0^2 \sin\theta$$ Setting $\Omega_0$ (or converting to dimensionless variables) for simplicity, I separate this into two equation $$\dot{v} = -\sin\theta$$ $$\dot{\theta} = v$$ I change this into a vector relationship, $$ \begin{pmatrix} \dot{v} \\ \dot{\theta} \end{pmatrix} = \begin{pmatrix} -\sin\theta(t) \\ v(t) \end{pmatrix} $$ Let's define the vector $\vec{y} = (v,\theta)^T$, then the left hand side is $\dot{\vec{y}}$ and the right side $\vec{f}(\vec{y}) = \vec{f}(\theta,v)$. The initial condition is given by $\vec{y}(t_0)$. Note there is no explicit time dependence, therefore the RK45 equations become \begin{align*} \vec{k}_1 &= \delta t\,\vec{f}(\vec{y}(t_n)) \\[4pt] \vec{k}_2 &= \delta t\,\vec{f}\left(\vec{y}(t_n) + \tfrac{1}{4}\vec{k}_1\right) \\[4pt] \vec{k}_3 &= \delta t\,\vec{f}\left(\vec{y}(t_n) + \tfrac{3}{32}\vec{k}_1 + \tfrac{9}{32}\vec{k}_2\right) \\[4pt] \vec{k}_4 &= \delta t\,\vec{f}\left(\vec{y}(t_n) + \tfrac{1932}{2197}\vec{k}_1 - \tfrac{7200}{2197} \vec{k}_2 + \tfrac{7296}{2197}\vec{k}_3\right) \\[4pt] \vec{k}_5 &= \delta t\,\vec{f}\left(\vec{y}(t_n) + \tfrac{439}{216}\vec{k}_1 - 8\vec{k}_2 + \tfrac{3680}{513}\vec{k}_3 - \tfrac{845}{4104}\vec{k}_4\right) \\[4pt] \vec{k}_6 &= \delta t\,\vec{f}\left(\vec{y}(t_n) - \tfrac{8}{27}\vec{k}_1 + 2\vec{k}_2 - \tfrac{3544}{2565}\vec{k}_3 + \tfrac{1859}{4104}\vec{k}_4 - \tfrac{11}{40}\vec{k}_5\right) \end{align*} where $\vec{y}(t_n) = (v(t_n),\theta(t_n))^T$. The new values are given to fourth order by $$\vec{y}^{(4)}(t_{n+1}) = \vec{y}(t_n) + \left(\tfrac{25}{216} \vec{k}_1 + \tfrac{1408}{2565}\vec{k}_3 + \tfrac{2197}{4101}\vec{k}_4 - \tfrac{1}{5} \vec{k}_5\right)$$ or to fifth order by \begin{align*} \vec{y}^{(5)}(t_{n+1}) &= \vec{y}(t_n) + \left(\tfrac{16}{135}\vec{k}_1 + \tfrac{6656}{12825}\vec{k}_3 + \tfrac{28561}{56430}\vec{k}_4 - \tfrac{9}{50} \vec{k}_5 + \tfrac{2}{55} \vec{k}_6\right) \end{align*} Here is where my question is. In computing the new step size, $\delta t \to s \delta t$, and the formula generalizes to $$s = \left(\frac{\epsilon\,\delta t_{old}}{2\left|\vec{y}^{(5)}(t_{n+1}) - \vec{y}^{(4)}(t_{n+1})\right|}\right)^{1/4}$$ I assume then that I can just treat the magnitude as the vector norm, but I am not certain this is correct. I have written a code which works (produces reasonable results with a fixed step size) but the adaptive step size seems to always want to rapidly decrease to zero, even if my error tolerance is relatively large, say $10^{-4}$. I'm not sure if this is a bug in my code or if I am misunderstanding how to apply the adaptive step size algorithm. Any help appreciated.

Answer 1

See Butcher tableau for RKF45 to find that there is a typo in one denominator for $y^{(4)}$. Using $4101$ reproduced the observed behavior, changing to $4104$ as in the other fractions and the source gave regular behavior, the step-size is corrected once for the first step to $h=0.55109$ for tolerance $10^{-4}$ (with my implementation details, about $0.3$ to $0.6$ should remain true) and remains that way during the full integration.