Recently I found a beautiful, although unphysical behavior when I was implementing the dynamics of elastic billiard balls ( hard spheres in physics lingo). And the bug was not at all where I expected it to be. But since while debugging I did not find any resource helping me out with this, and I think there is a good general lesson to be learned, I thought I write it up.
Pretty images
Okay. Mindset we want to simulate balls flying around and bouncing off of things, with the total kinetic energy being conserved. There is nothing else acting on the balls. No gravity, no nada.
So. Have a look at the animation below. What do you see?
Well. This is how it should look like
Did you notice something? Well besides the boring flying straight and occasional bouncing.
You may not, it takes a bit of patience and keen eyes. Both of which I lost during the debugging and looking at many of those animations to try and understand the problem.
Hence the impatient route. Let’s just trace out all the trajectories. Note, this only is visually manageable for a few balls for a few time increments, learned this the hard way.
Trace plot of the first animation
Trace plot of the second animation
You see it now don’t ya? :)
In each trace plot each ball has its own curve in the x-y plane. The color of the curve indicates the time the ball was at that location, starting off blue turning green and finally yellow. Each change in trajectory of a ball is denoted with a small black dot.
Now we can easily see that in the first trace plot that there are a lot of dots of two curves very nearby. Those are orbiting balls! Side note: THIS is how you keep sane analyzing those trajectories.
Equipped with trace plots let’s run the simulation for longer, just because we can and it neatly underlines the difference between the buggy implementation and the correct (boring) one.
We actually see something that resembles coils in the buggy trace plot. Pretty, no?
Event Driven Molecular Dynamics
The above thing is also referred to as Event Driven Molecular Dynamics (e.g. here). Juuust in case you want to research yourself a bit. Helps greatly to know the name of the things when researching.
When you have hard spheres without any external field, you are in the lucky position that you can just solve the equations of motion and jump the entire ball setup forward to the next event, instead of having to increment time by a small constant amount, hoping your time increment is small enough to not screw you over with numerical errors. The latter is the norm for molecular dynamics.
So event driven molecular dynamics works because hard spheres are totally unrealistic :P. I say unrealistic because normally you would expect that some effect applies between the balls which depends on their distance in some continuous fashion, like gravity or electro dynamics or what have you. Anyway.
The bug
The bug snuck in when implementing the estimate when the next ball-ball collision will be. It did because I copied the solution I found in a book, a paper and also on this page. Maybe I wasn’t paying close enough attention or something, but it still happened, so it may happen to you ^^. It was not, as I thought for too long, the numerical side, e.g. rounding errors.
This caused my bug
Δ t is our time increment, Δ v is the difference in velocity between our balls of interest, Δ r the difference in their position and d = (Δ v *Δr)^2 - (Δ v * Δv ) (Δr * Δr - 4 * sigma^2), for a prettier version see the same page again immediately below Δ t definition I copied. sigma is the radius of each ball.
Do you see the bug? If so I congratulate you. I didn’t. So let’s step through it.
The Δ t is derived from a single equation.
So a collision occurs once two balls are twice their radius, sigma, apart. The distance, Δ r between their positions is
using balls k and l in two dimensions, that’s where the 1 and 2 indices come from. Both positions are a function of time like
So now all we have to do is go backwards and plug x(t) into our Δ r equation and solve it for Δ t. This coincidentally is a quadratic equation in Δ t
with
Whose solution is ALMOST what had above. To prevent annoying scrolling, here it is again
it should be
with a +- sign and conditioned on the result being greater equal zero. By discarding the minus square root of d we ignore a second solution, which leads to some balls orbiting each other!!! So the original solution looks reasonable, but is technically not correct. Which was important here.
So what’s the general lesson I mentioned in the beginning?
Do use your intuition to understand solutions in papers and text books, but if you are implementing them try to understand the derivation of those solutions step by step. Otherwise you may end up orbiting bugs for some time. *ba dump tsss*