# Integration error in high velocity

I've implemented a simple simulation of two planets (simple 2D disks really) in which the only force is gravity and there is also collision detection/response (collisions are completely elastic). I can launch one planet into orbit of the other just fine.

The collision detection code though does not work so well. I noticed that when one planet hits the other in a free fall it speeds backward and goes much higher than its original position. Some poking around convinced me that the simplistic Euler integration is causing the error.

Consider this case. One object has a mass of 1kg and the other has a mass equal to earth. Say the object is 10 meters above ground. Assume that our dt (delta t) is 1 second. The object goes to the height of 9 meters at the end of the first iteration, 7 at the end of the second, 4 at the end of the third and 0 at the end of the fourth iteration.

At this points it hits the ground and bounces back with the speed of 10 meters per second. The problem is with dt=1, on the first iteration it bounces back to a height of 10. It takes several more steps to make the object change its course.

So my question is, what integration method can I use which fixes this problem. Should I split dt to smaller pieces when velocity is high? Or should I use another method altogether? What method do you suggest?

EDIT: You can see the source code here at github:https://github.com/elektito/diskworld/

• Please try symplectic euler integration. You will also likely need to use sphere to capsule collision detection in order to compute a time of impact for high speed collisions. If you find that high speed planets zipping around other planets at close proximity (but not colliding) to be a source of error, then you will need an integration scheme that takes into account the non-constant acceleration. There are some nice resources online for RK4 integration. May 5, 2014 at 19:27
• Hmmm. I just tried RK4 and turned out replacing my simple integration method with something else is not so easy, because integration code and collision detection code are so much mingled together. It seems I'll have to separate those somehow. I gave up for now! I'll be very grateful if you have any advice on this, too. May 5, 2014 at 19:41
• I should have mentioned that high speed planets zipping around others do not seem to be the source of the problem. The simulation seems to be working just fine when no collisions occur. Collision detection by itself seems to be working perfectly too. It's the combination that causes the problem, or so it seems. May 5, 2014 at 19:46
• If the problem is with integration error, then all this planet and orbit stuff obscures the question. Could you post a minimal example where the problem stil occurs?
– Anko
Jul 5, 2014 at 8:31

## 2 Answers

Say the object is 10 meters above ground. Assume that our dt (delta t) is 1 second. The object goes to the height of 9 meters at the end of the first iteration

Here lies your problem. It is true that the velocity at the end of the first iteration is $$\1 m.s^{-1}\$$. However during that time the object has not travelled $$\1 m\$$.

In fact, since the acceleration is constant, the average object velocity is simply $$\0.5 m.s^{-1}\$$ and thus the object has only travelled $$\0.5 m\$$. The next frames will be:

$$\begin{array}{c|c|c|l} \text{time} & \;\; v \;\; & v(avg) & \text{height} \\ \hline 0 & 0 & - & \;\;10.0 \\ 1 & 1 & 0.5 & \;\;9.5 \\ 2 & 2 & 1.5 & \;\;8.0 \\ 3 & 3 & 2.5 & \;\;5.5 \\ 4 & 4 & 3.5 & \;\;2.0 \\ 5 & 5 & 4.5 & -2.5 \text{(collision happened)} \end{array}$$

So, to know the new object position, use the average velocity instead of the new velocity. It will greatly improve your accuracy (in fact, in the case of constant acceleration, it will even give you exact results).

At this points it hits the ground and bounces back with the speed of 10 meters per second.

This is incorrect, too. Where does the value 10 come from? With your integration method, the velocity should be 4.

• The 10m/s likely comes from an incorrect collision response. A fairly common collision response is to make everything "rubbery" and use hooke's law to separate collisions (F = kx, where x is the penetration vector, and k is a spring constant). The problem is, this is no longer constant acceleration, and treating it as if it were usually ends up adding energy to the system. Aug 5, 2014 at 1:13

Just guessing from the description of your issue, since I don't have any actual code, I would also look in to the following if you haven't already:

What happens if you use a smaller time step, say 1/30 sec?

Check that you apply the time step variable (dt) in the appropriate parameters of your simulation calculation. Since you mentioned that the "combination causes the problem" and the physics are mangled with the collisions it might be necessary to apply the time step to both physics and collision detection/response code.

Make sure that the objects don't overlap during the collision. So the collision detection code will not run multiple times that may add the collision response force multiple times.

Make sure that physics calculations occur in fixed time step because if you use variable time step the integration produces small errors proportional to each time step that sum up during multiple frames.

Try to smooth the results of the physics calculations by averaging the results of the last two or three frames. That way you will also average the small errors and the end result will be perceived as more accurate.

• Let's see.Shortening the time step as expected alleviates the problem but it is still there. I do believe that dt is applied appropriately, that objects don't overlap. I fixed the time-step and it doesn't seem to help. And finally I'm not sure how to smoothe the results when there is a collision. After all the collision system is applying a sudden impulse to the object. Changing that won't be correct. How do you propose I do that? May 6, 2014 at 3:33
• I've added a link to my source code by the way. It uses python and pygame. Right clicking an empty space moves the smaller disk there. Dragging the smaller disk throws it. Mouse wheel zooms. Dragging empty space pans the camera. "p" pauses/un-pauses. May 6, 2014 at 3:35