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

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 :)

• I have faced the same issue myself. Note that when the spring engages it changes the system, and the RK4 assumes there is no change in the system during each time step. To model the correctly you need to create a time step at the moment the physics change. Commented Jun 27, 2022 at 20:19
• That's a good point! I'll start the RK4 at the time of collision and see if I can't recover the elongation that Kvante Kaffe (person who answered linked question) was able to capture in his model. However, can I not capture this transition by making my derivative functions appropriately defined piecewise? Give a small enough h it would be similar to just starting a new RK4 at the transition point, no? Commented Jun 27, 2022 at 20:43
• I started the simulation right at initiation of contact, and got effectively the same elongations :/ Commented Jun 27, 2022 at 20:53

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.

• Wow, thanks rob! I'll do that pen and paper exercise now, but one more quick question. I was checking my expectations using an energy argument as follows. The initial Energy in the system is 9J. The KE is at a min when compression is at a max, when the velocity curves intersect (1.5 ms). At that point in time the KE is 2.25. Meaning that the remaining 2.25J are stored in the springs. Thus 2.25J = E_spring1 +E_spring2 = 1/2k(x)^2+1/2k(x)^2 which tells us maximum compression should be .67, and maximum elongation should be ~.33, so dose this mean that simple NM can't capture this interaction? Commented Jun 28, 2022 at 17:23
• I found that the error rapidly decreases as the particles become more rigid, so maybe this phase space of particle collision is just super sensitive and hard to capture with numerical analysis! I updated this relationship between error and k value in my answer Commented Jun 28, 2022 at 17:43
• Your energy-conservation approach is more robust than my impulse approach (but you should get similar answers from both). However, you have an arithmetic error. You're correct that, at maximum compression, the system has only one-quarter of its initial kinetic energy. However you have typoed that one-quarter (instead of three-quarters) of the initial energy should be stored in the springs. This factor of $\sqrt3$ in length gives an unphysical squashing that I don't see numerically with $k=5$ N/m, but I do see with smaller $k$.
– rob
Commented Jun 28, 2022 at 18:02
• Like you, I got more reasonable results with stiffer $k$; I put those on my plot but I apparently left them out of the text.
– rob
Commented Jun 28, 2022 at 18:03
• Oh right, the initial energy is 4.5 J. Well thank you for working through this with me. It feels good to have arrived at some sort of consensus. I was going a bit mad trying to recreate the answer from that one post :) Commented Jun 28, 2022 at 18:10