Why is my simple RK4 failing to capture adequate spring deformation?

Question

I am setting up a very simple physical model for two colliding particles modeled as elastic spring's obeying Hooke's law. I have set up a numerical model for the collision and while the end stages of the model are correct something is going wrong during the collision interaction. I was hoping someone could point out where I am going wrong, I imagine it is in my implementation of the Runge-Kutta method to solve the equations of motion, but I will walk through my whole set-up.

My work is just a recreation of the solution to the question posted here, but I provide all the necessary information in this post, so there is no reason to flip back and forth. Ultimately I can't get my simulated results to match the results of the person who posted the answer:

Energy transfer during perfectly elastic collision between bodies of same mass

Finding the Equations of Motion

key variables

• $m_1, m_2, x_1, x_2, v_1, v_2$ mass position and velocity of particles one and two respectively
• $\Delta x_1$ width of particle 1
• $\Delta x_2$ width of particle 2
• $\Delta x_0$ equilibrium length of particles (same for both)
• $m_1,m_2$ masses of particles one and two respectively
• $k$ spring constant (same for both)

During collision we have: $$ F_{12}=-k(\Delta x_1 - \Delta x_0) $$ $$ F_{21} = k(\Delta x_2 - \Delta x_0) $$

Because both particles have the same spring constant newtons third law guarantees that the deformation, now just called $\Delta x$ is equivalent : $$ \Delta x _1 = \Delta x_2 \equiv \Delta x $$

Finally, we relate this deformation to the particles position as follows: $$ \Delta x = \frac{\Delta x_1}{2} + \frac{\Delta x_2}{2}=x_2-x_1 $$

Now we can write our equations of motion! $$ a_1 = \frac{k}{m_1}(\Delta x_2 - \Delta x_0) = \frac{k}{m_1}(x_2-x_1-\Delta x_0) $$ $$ a_2 = -\frac{k}{m_2}(\Delta x_2 - \Delta x_0) = -\frac{k}{m_2}(x_2-x_1-\Delta x_0) $$

Now for setting up the RK4 to solve these equations of motion

Numerical Solution to Equations of Motion

I set up a multi-variable RK4 of the form $$ X' = F(X) $$

Where $X=[x_1,x_2,v_1,v_2]^T$, $X'=[\frac{dx_1}{dt},\frac{dx_2}{dt},\frac{dv_1}{dt},\frac{dv_2}{dt}]^T$, and F is a vector of function linking the two $F=[f_1(x_1),f_2(x_2),f_3(v_1),f_4(v_2)]$, more specifically:

$$ f_1(x_1)=v_1 $$ $$ f_2(x_2)=v_2 $$ $$ f_3(v_1) = \frac{k}{m_1}(x_2-x_1-\Delta x_0) \text{(or 0 if }x_2-x_1>\Delta x_0) $$ $$ f_4(v_2) =-\frac{k}{m_2}(x_2-x_1-\Delta x_0) \text{(or 0 if }x_2-x_1>\Delta x_0) $$

Now to actually run the RK4 I use the following initial conditions: $x_1=-1 m $, $x_2=1 m$, $m_1=1 kg$, $m_2=1 kg$, $k=5 \frac{N}{m}$, $v_1=3 m/s$, $v_2= 0m/s$.

I set h = 0.1, and $t_{end}=3.0 s$, the number of steps equal to $N=t_{end}/h$, $X_1$ equal to the relevant initial conditions and then set up the following loop:

for i in 2:N $$ k_1 = F(X_{i-1}) $$ $$ k_2 = F(X_{i-1}+\frac{h}{2}k_1) $$ $$ k_3 = F(X_{i-1}+\frac{h}{2}k_2) $$ $$ k_4 = F(X_{i-1}+ h*k_3) $$ $$ X_i = X_{i-1}+\frac{h}{6}(k_1+2k_2+2k_3+k_4) $$ end

My end results recover complete momentum transfer but the springs end up compressing almost entirely with $x_2-x_1-\Delta x_0$ getting very close to zero. These are summarized in the gif below:

When compared to the expected results posted at the top of the question we see that my model is tacking too long to transfer momentum. The particles compress too much and x1 travels too far. A simple check using conservation of energy show this is not realistic given the starting conditions. Why is my simulation failing to recover the results that the answer in the above question recovered? I have tried running this for up to h =0.000001 and the results are very stable, but seemingly wrong...i get almost complete compression every time. Below are my elongation results for different timesteps:

I fell like I'm missing something simple in my model...can you see what it is?

EDIT: Rob provided an excellent answer, and I think that the post in question must have just used a different k-value to create their results?

I was skeptical of my results because an energy argument shows that we should expect much less deformation from the spring. Thee energy argument is as follows:

• The Initially the total energy is just equal to $KE=1/2m_1v_1^2=4.5 J$
• Maximum compression of the springs occurs when spring energy is at a maximum and KE is at a minimum. KE is at a minimum when the velocity curves intersect and we have $KE_{min}= 1/2m_1v_1^2+1/2m_1v_2^2=2.25 J$
• Thus we can get the expected compression of each particle with $2.25 J = Energy_{spring1}+Energy_{spring2}=k\Delta x^2 \rightarrow \Delta x = 0.67$ and thus elongation $\approx .33$

This is the theoretical compression of the spring. I can extract the experimental compression of the spring from my simulations. Calculating the percent error for different values of k shows that our model is actually pretty good for rigid particles, it just fails in the phase space of highly deform able particles. This is summarized in the figure below:

TLDR; If you want to model particle collisions using simple numerical analysis,you'll be fine as long as your particles are rigid :)

Answer 1

It's not your Runge-Kutta implementation: I reproduce your result. A quick-n-dirty Python implementation, using RK45 from scipy.integrate with its default tolerances and adaptive stepsize:

def dxdt(t, y, k=5, dx=1, m=(1,1,)):
    x1,x2,v1,v2 = y
    m1, m2 = m
    if abs(x2-x1) > dx:
        a1, a2 = 0,0
    else:
        a1 = +k*(x2-x1-dx) / m1
        a2 = -k*(x2-x1-dx) / m2
    return v1,v2,a1,a2

trange = (0,2)
result = scipy.integrate.solve_ivp(dxdt, trange, (-1,1,3,0), dense_output=True)

plot(result.t, result.y[:2].T, '.')
t_dense = linspace(*trange)
plot(t_dense, result.sol(t_dense)[:2].T, '-')

Note that taking a zillion tiny steps is not the issue: the algorithm is only explicitly evaluated at about a dozen points, and the rest are filled in reasonably well with interpolation. An eighth-order RK method and a backwards-differentiation implicit method give comparable results, though they (a) adaptively compute a bunch more steps around the discontinuity in the force, and (b) do a little better at having the first object actually come to rest afterwards.

I think the difference in results is a physics problem. For the first object to stop, it must undergo an impulse

$$ \vec J = \vec p_\text{final} - \vec p_\text{initial} = \int \vec F\,\mathrm dt $$

Within a factor of a few, we'll have

$$ J_1 = -m_1 v_1 \sim k\,\Delta x\,\Delta t $$

where the length scale $\Delta x$ is like the size of the particle, and the time scale $\Delta t$ is the duration of the collision. The distance traveled from the start of the collision to the end will go like

$$ v_1\Delta t \sim \frac{k\,\Delta x}{m_1} (\Delta t)^2 $$

In your plots (and in my duplication of them) the duration $\Delta t$ of the interaction is nearly a second long. In the results you are comparing to, the duration of the collision is only about a half-second long. Furthermore, your model predicts a smaller compression scale $\Delta x$ than the other model. Since you are in principle removing the same initial momentum, I think that at least one of the two models is effectively using a different value for the spring constant $k$.

A useful pen-and-paper exercise would be to make an educated guess for the average force you expect. For a bell-curve-ish shape like your parabolic force curve, a good order-of-magnitude approximation is a triangle force, $F_\text{average} \approx \frac12 k\Delta x_\text{max}$. Then you can figure out right away how long the duration must last so that the impulse $F\,\Delta t$ can remove the initial momentum.