My problem in brief: in some situations, the Runge Kutta 4th order method (RK4) doesn't seem to give 4th order improvement when using a smaller time step. I wonder how this worse-than-expected result can be explained.

Consider two identical massive steel spherical balls with radius $R$. The values of the relevant properties of the steel (density, Young's modulus and Poisson ratio) are given. We study what happens during an elastic, central collision of the balls when they come in contact in free space with an initial relative speed $v_0$.

During the collision, each ball is deformed slightly. The deformation $x$ is the difference between $2\cdot R$ and the momentary distance between the centres of the balls. As a result of the deformation, the balls exert a kind of restoring 'spring force' $F$ on each other.
Heinrich Hertz derived a relation of the form $F=-C*x^{3/2}$ between $F$ and $x$. The constant $C$ depends on the material properties of the steel, the ball radius $R$ and the initial relative speed $v_0$.
Hertz also derived an analytic expression for the total time duration $\tau$ of the collision, with $x$ varying from $x=0$ to a maximum $x=x_{max}$ (speed $v=0$) and back to $x=0$.
An expression for $x_{max}$ follows from conservation of mechanical energy.

We can interpret Newton's second law $F=m\cdot a$ as a second order differential equation.
So, using centre-of-mass coordinates, $\frac{\text{d}^2 x}{\text{d}t^2}=-C/\mu*x^{3/2}$ where $\mu=\frac{m_1\cdot m_2}{m_1+m_2}$ is the reduced mass of the two balls.

The usual RK4 iteration scheme for such a second order differential equation with a time step dt is $\text{d}x_1=v\cdot \text{d}t\qquad\text{d}v_1=-C/\mu\cdot x^{3/2}\cdot\text{d}t$






$x=x+\text{d}x\qquad v=v+\text{d}v\qquad t=t+\text{d}t$

The iteration is ended when $x$ is very nearly the analytical $x_{max}$, using a greatly reduced time step dt near the end. Then $v$ turns out to be very nearly 0.

Using RK4, I calculate the time duration $T$ from $x=0$ to $x_{max}$.
I calculate the same duration $T=\tau/2$ using Hertz' analytical expression.
The difference between these values of $T$ is called $\Delta T$.
Using RK4, I expect that $\Delta T$ becomes $10^4$ times smaller when I repeat the RK4 iteration with a time step dt which is 10 times smaller than the previous dt. Instead, $\Delta T$ turns out to be $10^{2.500}$ times smaller when dt is 10 times smaller. So RK4 only seems to achieve an improvement of order $10^{2.500}$, not $10^4$.

This does not seem to be due to a coding error in RK4. The same code produces an improvement with a factor of $10^4$ when I only change the restoring force from $F=-C\cdot x^{3/2}$ to $F=-C\cdot x$ or $F=-C\cdot x^2$.
All RK4 calculations are done in more than 50 correct significant digits.

My questions: Is there a mathematical or computational explanation why the improvement of RK4 with $F=-C\cdot x^{3/2}$ is of order $10^{2.500}$ instead of $10^4$ when reducing the time step?
Under which conditions can I expect 4th order improvement when reducing the time step?

I realise there are other possibilities for improvement, such as using an adaptive time step or higher order Runge Kutta. My curiosity concerns the 'bare' RK4 method, though.

  • 5
    $\begingroup$ $x^{3/2}$ has unbounded derivatives near $x=0$, which the other cases that you tried do not. RK error bounds are based on Taylor series, so these bounds likely are invalid in the presence of unbounded derivatives. $\endgroup$
    – whpowell96
    Commented Sep 6, 2023 at 20:57


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.