Numerical issues when implementing physics

Contents

  • 1. The engine
  • 2. Numerical problems
    • 2.1. Invisible contact
    • 2.2. Breaking contact

Written by Bert Peers, feedback very welcome on bpeers@acm.org ! Thanks & enjoy..
Last update : 31/05/2000 18:44

1. The Engine

This article discusses a few issues that can arise when implementing a physics engine. Ofcourse, there are numerous ways to implement such an engine, but this article talks specifically about the classic Baraff-like setup (see for instance Siggraph Course Notes 97).

In this setup, collision detection is first used to compute the time at which a penetration will occur when integrating the current velocities and orientations. The scene is then integrated up to the point of non-penetrating contact, and impulses are applied to remove any negative (ie, penetrating) velocities at all contact points. Next, the engine computes what forces and torques are acting on the contact points, verifying if any of these exercise a negative acceleration (negative again as in resulting in a penetrating velocity when integrated). If so, the resulting linear complimentarity problem is formulated and solved, resulting in a vector of contact forces that have the net effect of neutralising any negative acceleration.

Although this is a very brief description, it's probably clear that there is a lot of code involved, with quite a few numerical techniques. Still, it may come as a surprise that even with all this sophistication, and with tested and verified code routines, it is still possible for such an engine to freeze up or go crazy on the simplest of scenes (such as illustrated below), due to numerical precision problems. This text demonstrates two of these problems.

Stuck on a simple case

2. Numerical problems

The problems discussed here arise from the fact that the code must always use a finite approximation to the continuous functions commonly used in physics papers. For instance, it is easy for a paper to specify that the object positions can be integrated up to the point of contact, but in practice the collision detection will use e.g. binary rootfinding to find a time t at which point two objects are 'getting close'. For a realtime engine, one commonly choses an epsilon defining 'close' which will result in a guaranteed upperbound on the time necessary for finding the contact time. However, even in an offline simulation there is a limit to how small epsilon can be, simply because of the limited precision of a float or double (except if one uses arbitrary precision packages). Thus, the arguments below are anything but 'theoretical worst-cases' that can be solved by just using 'more cycles' : the problems can arise for arbitrarily small definitions of 'close', 'timesteps' or 'touching'.

2.1. Invisible contact

A well-known problem with non-analytical collision detection is 'tunneling' : when we integrate a certain timestep, verifying that no collision occurs at time t+1, it is very well possible that there was a collision at time t+0.5, but we just missed it because at time t+1 the collision is already over; an example is a bullet flying at a wall with such high velocity that is in fact already at the other side of the wall at t+1.

A more subtle variation of this problem is the following. Suppose that we have integrated up to time t, and have detected a point penetrating a surface at time t+epsilon. For reasons of accuracy, epsilon is the smallest timestep that we can take and thus we have to take measures at time t to prevent the penetration from occuring. Normally, we'd apply an impulse to the would-be-penetrating point, therefore altering it's course so it no longer penetrates at t+epsilon. However, to apply this impulse we must know what the point is that will be penetrating 'in a moment'. Commonly, there is a collision detection routine that will find all the points (closest features) within a certain tolerance distance delta. The idea then would be to look at the situation at time t, ask the CD routine to list all points that are already 'dangerously' close, and apply an impulse where necessary. Hopefully, this list will include the point that is about to penetrate.

The achilles'heel here is the word 'hopefully' : there is no guarantee that the point p which will penetrate a surface S at time t+epsilon is also within a distance delta from S at time t ! If that sounded confusing, read it again, because it's the whole foundation of the problem. In a perfect world, we would just keep refining epsilon until we're at a moment in time where p is exactly touching S, or at the very least is within a distance delta of it. But, as explained above, epsilon is really just a given constant, and if our collision detection routine does not report any points (or reports a set of points that does not enlist p) at time t, then p will not receive an impulse, and thus will still be penetrating at time t+epsilon. The whole system is now stuck at time t : it cannot integrate any further in time because it finds a penetration, but it is also unable to resolve the penetration with an impulse. This problem can be considered identical to the bullet-through-wall example given above : the 'bullet' p flies through the 'safety wall' which has thickness delta, covering S.


Figure 1 : A colliding point goes undetected

A possible solution would be to keep increasing delta, say by a factor 2, until p is detected. This is easier said than done : how, in general, do you define the condition that must be met to stop increasing delta and be happy with the list of collision points that comes out of the CD routine ? In the example, we know that we should increase delta until p is a part of the collision list L. It is not sufficient to just increase delta until L is not empty. If there is geometry elsewhere in the scene that is making contact without penetration, those contact points will be part of L even if delta is zero. It is also not sufficient to increase delta until at least one feature pair is added to L which has a negative relative velocity, since there may not be such a pair (see 2.2). A possible criterium is to increase delta until at least one pair is added which has a relative velocity which is not zero (ie outside a small tolerance of zero), thus ignoring any touching but non-moving contact points.

An alternative solution would be to avoid using a delta to compute (or more accurately, guess) the list of will-be-penetrating points. This can be done by integrating to time t+epsilon, computing the list of penetrating points L', and transforming these points back to time t to yield L. This guarantees that a point p that will be penetrating is known at this time t. This is not so easy in practice : depending on the integration scheme used it may not be trivial at all to reverse the integration. On the other hand, there is nothing in the impulse formulas that actually requires the points that are being pushed to be on the object, and if a simple reverse-integration is used (such as Euler), it is not a problem if a point in L' is reversed to a point p' which is not really p at time t. Note though that using two different integrators (one for going forward, another simpler one for going backwards) no longer guarantees that applying an impulse to p' prevents the 'real' p from penetrating S at t+epsilon -- so you might still be stuck.

2.2. No impulses

A second problem that can occur, somewhat related to 2.1, is that even when the closest features are known, and these features are about to penetrate, none of them have a negative velocity right now (at time t). More precisely, it is possible that two objects are close enough together that they have a closest feature pair, and these features are safely moving away from eachother, but nevertheless, at time t+epsilon the two objects are penetrating. The figure illustrates how this can occur : the object in the lower right is not moving at all, while the object in the upperleft is slightly translating up and to the right, slowly rotating. As a result, the marked point is moving away from the stationary object's vertex, but the object is penetrating at time t+epsilon because it's edge has rotated into the stationary object.


Figure 2 : Rotating edge penetrates with a point that is not the closest right now

Again, this is a problem of precision. If the two closest points on the two objects are moving away from each other at time t, but at time t+epsilon the two objects are penetrating, then clearly there must a be a certain point in time t', t < t' < t+epsilon where the closest features are at exactly distance zero, with a negative velocity. Because we cannot step to time t', we're once more stuck with a situation were the 'present' does not hold a 'dangerous' condition, even though the 'near future' is in an invalid state (penetrating). This situation is not as theoretical as it may seem. In particular, when two objects of approximately the same mass and inertia collide, it is common (as observed in reality too) for them to align themselves surface (or edge) against surface, with low tangential velocities.

This problem is unrelated to contact forces : if the problem would be that the point which is currently breaking apart is in fact accelerating back (ie negatively) to the surface it's breaking apart from, then penetration would only occur at a later point in time, as the velocity must first become negative again; only then will the position start moving back to the surface, taking another amount of time to penetrate. For the same reason, it is not possible to apply a force (such as a penalty) to resolve the collision, because that would require integration of the velocities to have a result, but we're at a situation where an integration timestep epsilon already causes problems, immediately.

The simplest hack seems to be the only possible solution here : if a list of contact points L is computed at time t, and a penetration is detected at time t+epsilon, but none of the points in L have a negative relative velocity, then simply apply a small impulse anyway that'll push the points in L a bit apart even further, hopefully breaking the immanent penetration in the process. This can be iterated until either the penetration is indeed averted, or 'normal' impulses are possible (ie, points appear in L that have truely negative velocities).

One caveat here is that it may be a good idea to apply these impulses only selectively. In particular, L will also contain points that are caused by contact between two objects in rest (eg a cube resting on a larger cube). Applying the impulse would destabilise this rest-configuration, pushing the cubes apart. It's a better idea to only apply the impulse-penalty to those feature pairs that are actually going to penetrate at time t+epsilon. A simple way to do this is to actually integrate to t+epsilon, remember which object pairs (A, B) are penetrating, and apply the impulse to only those points p in L related to A and B.

Homepage