📚 Perfect, Infinite-Precision, Game Physics in Python (Part 4)

Improving speed and features, perfectly and imperfectly

This is the fourth of four articles showing you how to program a perfect physics engine in Python. It is a tiny step in my grand ambition to turn all of physics, mathematics, and even philosophy into programming. All code is available on GitHub. Part 1 introduced the high-level engine and applied it to physics worlds such as Newton’s Cradle and the tennis ball & basketball drop. Part 2 applied the engine to the billiard’s break, uncovering unexpected non-determinism. In Part 3, we created the low-level engine by engineering equations then having the SymPy package do the hard math.

This article is about the limitations of our perfect physics engine.

It is slow. We will see how to make it a bit faster, but not fast enough to make it practical. Unless … we replace its perfect calculations with numerical approximations. We can do this, but if we do, we would be removing the engine’s reason for being.

We’ll also look at its many other limitations — for example, being only 2D, not having gravity, etc. For each limitation, we’ll discuss a possible extension to fix it. Some extensions would be easy. Others would be very difficult.

Speed

When I started this physics engine, I hoped that all the expressions would simplify to something like 8*sqrt(3)/3. My hopes were dashed when I looked into the expressions. Consider three circles inscribed by a triangle:

We’ll look at the expression for one ball’s x position after some number of collision events:

After 10 collision events, the expression looks good. After 40, we see a square root inside a square root, which seems ominous. By 50 events, the expression explodes.

The next plot shows the length of the expression vs. the number of collision events. Notice that the y-axis is logarithmic, so we’re seeing at least an exponential explosion.

I implemented the following speed ups. They make the engine several times faster. Note, however, that a “several times faster” exponential, is still exponentially slow.

• Use floats when safe: Where safe, evaluate an expression as an (approximate) float or complex. For example, to check if two objects are moving toward each other, we find their instantaneous relative speed (see Part 3). We need to know if speed>0. I assume that if float(speed)<-0.00001, then speed is not positive. Similarly, if the (approximate) imaginary part of a speed is more than 0.00001 from 0, I assume the speed is complex. This optimization helps. However, some values will be too close to zero to use it and will still need to be evaluated symbolically.
• Recycle some next collisions: In the world below, A and B will collide next. A bit later, C will collide with the wall. After A and B collide, we need to compute their next collisions with every object in the scene. We don’t, however, need to recompute the collisions of C with the wall, because it was not involved in A and B’s collision. In many situations, avoiding this recomputation provides a huge speed up (from quadradic to linear in the number of objects). We can also avoid a bit of computation by skipping pairs for which both pairs are stationary.

• Use some parallelism: Several parts of the engine’s work can be done in parallel. Specifically, finding the time until the next collision for a list of object pairs. Also, after saving all collision events, we can render the video frames independently of each other. I used mapreduce1, a multithreading/multiprocessing library that I originally wrote for genomics processing. Surprisingly, running SymPy functions in multiple processes gave me no net gain. I don’t know the cause, but suspect SymPy’s automatic caching.

So, is a perfect infinite-precision physics engine doomed to be impractical? Yes! However, imperfect and finite-precision can be practical.

Numerical simulation: I didn’t implement this. However, if you are interested, here are some steps for converting to a numerical simulator:

• Generate numeric code (with efficient common subexpressions). For example, here we generate numeric code for the time until the next circle-circle collision. It generates 35 expressions.

from sympy import cse, numbered_symbols
from perfect_physics import load
cc_time_solutions = load("cc_time_solutions.sympy")
cse(cc_time_solutions, symbols=numbered_symbols("temp"))
# Outputs ...
([(temp0, a_vx**2),
  (temp1, a_vy**2),
  (temp2, b_vx**2),
  (temp3, b_vy**2),
  (temp4, a_vx*b_vx),
...

• Combine the x position of every circle into a single x NumPy vector. Likewise, y, vx, vy, m, r. For example, if we have six circles, and one is moving with a vx of 1 and the others are stationary, we would say: vx = np.array([1, 0, 0, 0, 0, 0]). Use NumPy matrix and vector operators to automatically compute in parallel. Here we compute the first five subexpressions above for all 36 pairs of circles in our Newton’s Cradle. Instead of 5 x 6 x 6 calculations, we need just three fast NumPy steps:

import numpy as np

vx = np.array([1, 0, 0, 0, 0, 0])
vy = np.array([0, 0, 0, 0, 0, 0])
temp0 = vx**2 # this covers temp2
temp1 = vy**2 # this covers temp3
temp4 = np.outer(vx, vx)
temp4
# Outputs:
array([[1, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0]])

If we do all the steps this way, we’ll end up with two 6 x 6 matrices that tell us the times to collisions between any two circles. It may seem inefficient to compute both the time when circle A collides with B and B with A (and A with A), but with NumPy this is efficient.

Beyond being slow, our simulations are limited in other ways. Let’s look at these limitations and consider the practicality of improvements.

Limitations and Extensions

Here are some of the engine’s limitations, roughly in the order that I think they could/should be fixed.

Explicit Units

All physics engines should give units — meters, feet, meters per second, miles per hour, etc. The engine currently works for any consistent set of units, for example, SI. It does not, however, make these units explicit. Adopting SI and making these units explicit should be done and should be easy to do.

More interesting and fun would be to use SymPy’s unit library. I believe this would let the engine use any consistent unit system and easily convert between unit systems.

Simple Inelastic Collisions

Our simulations look unrealistic because the circles never slow down. A simple fix: after each collision reduce (kinetic) energy by a given percentage. This is done by reducing vx and vy proportionally to each other.

Downward Constant Gravity

The current simulator doesn’t know about gravity. This means, for example, that the tennis ball and basketball “drop” (Part 1) doesn’t actually drop.

Without gravity, circles move in straight lines between collisions. If we added downward constant gravity, the circles would instead move in parabolas. This would add more non-linearities to the equations we had SymPy solve. As far as I know, SymPy should be able to handle these.

Particles

Some physics papers discuss collisions between particles. I had hoped that the current engine could model a particle as a circle with zero radius. Sadly, it failed because colliding particles end up at exactly the same position but without information about which “side” of the particle is touching the other particle. I’m optimistic that the engine could be extended to fix this.

3D

It would be fun to extend the engine to work with spheres in 3D. Walls would then be infinite planes defined by three non-colinear points. I think our equations (Part 3) could straightforwardly be extended to 3D and that SymPy could again solve them for us.

To implement this extension, I’d personally need to brush up on my geometry of planes. Recall, however, that for 2D we used slopes and unit vectors rather than trigonometric functions. I think that gives us a smooth path to 3D. Also, instead of introducing a new variable z to our existing x and y variables, we should instead replace all three with a single variable position, a 3-component tuple. (Likewise for velocity.)

Rendering a 3D scene would be more involved than our current 2D plotting, but I presume packages exist for that.

Unbound Variables

We currently allow the properties of a circle, say vx, to be an expression, such as 3*sqrt(3)/7. We don’t, however, allow it to just stay vx, that is to remain unbound. This stops us from solving physics problems for which concrete values are not given or needed.

I think we could create a version of the engine that could handle unbound variables. The engine, however, would need to report its answer in cases. For example, if the left-most ball is going a_vx, what will be the final b_vx of the right-most ball?

If you answered b_vx=a_vx then you remember Newton’s Cradle. A better answer, however, is: b_vx=a_vx if a_vx > 0; b_vx=0, otherwise. In other words, we must consider the case in which a_vx is zero or negative and the left most ball never hits another ball.

Angular Momentum

I’d love to create a perfect physics engine that understands something about spinning objects and angular momentum. I would start with a new engine that only knew about angular momentum and nothing about the linear momentum of the current engine. Doing both linear and angular perfectly seems very hard.

Other Shapes

Why just circles; why not squares, etc.? Why just infinite walls; why not finite walls?

To do these perfectly would be very hard. Squares spin with angular momentum. Finite walls can be missed or hit the narrow way. Both squares and finite walls introduce convex corners.

More Limitations

Some limitations seem so hard that I can’t imagine extensions to fix them perfectly:

• Realistic friction and inelastic collisions
• Realistic collisions via springs or other non-rigid bodies.
• Mutual gravity (the three-body problem, known to be very, very hard to do perfectly)
• Speed of sound

Summing Up Part 4

So, we made our engine faster by adding numeric approximations where safe and by recycling collision information. Even with this, however, it is not fast enough to be practical for most uses. We also enumerated its limitations that could be overcome with extensions. The extensions of most interest to me are adding:

• Units
• Downward constant gravity
• 3D spheres and infinite planes as walls

Thank you for joining me on this journey to create a perfect physics engine. I encourage you to implement your own engine or to extend this one. You can download the code from CarlKCarlK/perfect-physics (github.com). Let me know if there is interest and I’ll create a nicer installer.

Follow Carl M. Kadie — Medium for notifications of new articles. On YouTube, I have older (approximate) physics simulations and some trying-to-be-humorous nature videos.

Perfect, Infinite-Precision, Game Physics in Python (Part 4) was originally published in Towards Data Science on Medium, where people are continuing the conversation by highlighting and responding to this story.