1
$\begingroup$

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.

$\endgroup$
0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thanks that did it! $\endgroup$ – Kai Mar 21 '18 at 17:34

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.