I am currently working on a rigid body simulator. Before getting proper collision and rotations working, I decided to simulate the collision of spheres as axis-aligned objects.

The simulation works by, each frame, calculating all forces on the sphere (just gravity for now), calculating the sphere's acceleration, applying the acceleration to change the velocity, and then applying the velocity to change the position. Both applications take into account the varying time to compute each frame to produce a real-time simulation.

After each sphere has moved, it is then checked for a collision with all other spheres based on the idea that if two spheres are colliding then the sum of their radii is greater than the distance between their centers. When two spheres collide, vector projection (dot products) is used to separate the velocity vectors in the direction between them. Then, using conservation of momentum and velocity in one dimension along this direction, the velocities of both spheres after the collision are calculated. The equations used are as follows, found with the Maxima computer algebra system:

Maxima command:

], [l_1, l_2]);


$$ {l_1}=\frac{2 {m_2} {v_2}+\left( {m_1}-{m_2}\right) {v_1}}{{m_2}+{m_1}},{l_2}=\frac{\left( {m_2}-{m_1}\right) {v_2}+2 {m_1} {v_1}}{{m_2}+{m_1}} $$

The balls are simulated inside an invisible box. Collisions with this box are simulated by changing the sign on the component of the velocity perpendicular to the wall when a ball collides. This should not violate conservation of energy as that component must be squared to find the magnitude of the velocity, making its sign irrelevant.

The issue is that, after the simulation runs for a few minutes, all of the balls have a very clearly diminished velocity, implying that kinetic energy is being destroyed.

To attempt to fix this, I tried the following:

  • I wrote a program to verify that kinetic energy is conserved in 3 dimensions. There is no loss of kinetic energy with the above formulas with my dot product system in the 1 billion tested collisions.
  • I temporarily disabled gravity (and gave randomized initial velocities). In this scenario, the issue goes away. I don't know why.
  • I temporarily disabled collisions between balls (balls still collide with walls). This also "fixes" the issue. I don't know why.
  • I considered that floating point and other errors inherent in the approximate nature of the simulation may cause this issue. This, however, does not make sense as floating point errors should round up and down in approximately equal proportions, and thus not cause issues like this, and should only be relevant around a bifurcation point in the state space, or when performing very precise calculations. It seems to me that neither scenario is the case in this situation. Additionally, this does not explain why disabling gravity fixes the issue.

Here is a video of the simulation:



closed as off-topic by Bill N, Sebastian Riese, Buzz, ZeroTheHero, Kyle Kanos Dec 12 '18 at 11:08

  • This question does not appear to be about physics within the scope defined in the help center.
If this question can be reworded to fit the rules in the help center, please edit the question.

  • $\begingroup$ Just a wild guess: what happens if more than two balls collide simultaneously? Can you log the total energy for each frame and look for drops? $\endgroup$ – Jasper Dec 11 '18 at 20:33
  • 4
    $\begingroup$ I ran across a similar problem once, and in my case the issue was energy drift, which is the non-conservation of energy in step-wise integrators. The integration scheme you describe in you second paragraph is known to violate energy conservation. One way to reduce (eliminate?) the problem is to use a symplectic integrator, as mentioned in the Wikipedia article above. There is a simple symplectic integrator that involve simply changing the order of operations and/or changing a few signs. It's been a while, details are forgotten. :-( $\endgroup$ – garyp Dec 11 '18 at 21:07
  • 1
    $\begingroup$ @garyp You mean taking the "obvious" Euler integrator $x(t+1)=x(t)+h\cdot dx(t), dx(t+1)=dx(t)+h\cdot f(t)$ and flipping the order of evaluation to $dx(t+1)=dx(t)+h\cdot f(t), x(t+1)=x(t)+h\cdot dx(t+1)$ ? $\endgroup$ – Anders Sandberg Dec 11 '18 at 23:09
  • 5
    $\begingroup$ I'm voting to close this question as off-topic because this is about debugging software. $\endgroup$ – Bill N Dec 12 '18 at 0:20
  • 1
    $\begingroup$ @BillN Although software is a component in this, physics is also a component. People here are more likely to know at least some programming than people on more programming focused stack exchanges are likely to know physics because 1. this is stack exchange -- a lot of people know programming, and 2. programming simulations is something useful to professional physicists in order to perform experiments that are impractical to perform in real life. In other words, a good physicist can at least muddle through programming even though programming is perfectly reasonable without physics knowledge. $\endgroup$ – john01dav Dec 12 '18 at 1:14