The world’s leading publication for data science, AI, and ML professionals.

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

Use Python SymPy to turn Math and Physics into Programming

This is the third 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. For a video overview of the articles, see this talk from PyData 2023. To play with the engine yourself, try this live demo.

This article is about getting the computer to do math for us. Not just numerical calculations, but also algebraic and calculus-related manipulations. This is useful because:

  • I – and perhaps you – don’t remember all my high school and college math.
  • In school, our teachers only gave us problems we could solve. I want to solve problems I can’t solve.

Recall that in Part 1 and Part 2 of this series we developed the top-level of a perfect physics engine. It uses expressions such as 175–9*sqrt(3)/2 to represent time, position, and velocity exactly, that is, with no approximations. We applied the engine to Newton’s Cradle, a tennis ball and basketball drop, and a billiards break. The simulations resulted in surprising-to-me discoveries about physics and, even, philosophy.

Here in Part 3, we’ll use the computer to derive the needed expressions from the basic laws of physics. We will do this work in five steps:

  1. Model two moving circles. Use SymPy to find exactly when two moving circles will "just touch". Understand why SymPy returns multiple answers.
  2. Model a circle and a wall. Find when they will just touch. Add a "secret" sixth equation about angles (but don’t use angles).
  3. Try to model circle-to-circle collisions with conservation of momentum and energy. You can’t. Add a "secret" fourth equation about angles (but don’t use angles).
  4. Assume walls have infinite mass. A circle bouncing off a wall violates the conservation of momentum. Fix this with limits.
  5. Two objects can be touching and yet moving away from each other. Use limits to find these. Remove them from the list of colliding objects.

We start with two circles.

Step 1: Model two moving circles. Use SymPy to find exactly when two moving circles will "just touch". Understand why SymPy returns multiple answers.

We’re going to model the world as circles moving on an infinite two-dimensional plane. Circles move at constant velocity until they bump. Circles can bump into each other or into walls. Walls are (infinitely long) lines.

The first question we need to answer is when (if ever) will two moving circles just start to touch. To answer this question, we’ll engineer a set of equations. We’ll then have the computer solve the equations and give us the answer.

To engineer our first equation, think about two circles, a and b. Each has a position (x, y), a velocity (vx, vy), and a radius, r. Can you predict a circle’s future position (x’, y’) for any time t?

Our first equation says that circle a‘s x position at time t is its current x position plus its vx velocity multiplied by time t. In math notation:

Equation 2 will tell us circle a‘s y position in the future. Equations 3 and 4 will tell us b‘s x and y position in the future.

As a first step toward turning math into programming, we’ll write these four equations with Python’s SymPy package. SymPy is a computer algebra system, similar to Mathematica and Maple. Unlike those systems, however, SymPy is free and available inside Python.

Aside #1: Don’t confuse SymPy (a Python package for symbolic mathematics) with SimPy (a Python package for simulation).

Aside #2: Computer algebra systems represent expressions as nested computer objects. This "symbolic" representation makes some things easy for the computer – for example, working with fractions or taking derivatives of other expressions. Even with a symbolic representation, some tasks remain difficult for the computer – for example, taking integrals or solving systems of nonlinear equations. This tracks with what is easy and hard for people. I like using computer algebra systems because they are better than I at most tasks, easy or hard.

Here are our four equations in Python & SymPy:

from sympy import symbols, Eq

# define symbols
t = symbols("t")
a_x, a_y, a_vx, a_vy, a_r, aprime_x, aprime_y = symbols("a_x, a_y, a_vx, a_vy, a_r, a'_x, a'_y")
b_x, b_y, b_vx, b_vy, b_r, bprime_x, bprime_y = symbols("b_x, b_y, b_vx, b_vy, b_r, b'_x, b'_y")

# define equations
eq1 = Eq(aprime_x, a_x + a_vx * t)
eq2 = Eq(aprime_y, a_y + a_vy * t)
eq3 = Eq(bprime_x, b_x + b_vx * t)
eq4 = Eq(bprime_y, b_y + b_vy * t)p

We still need to engineer the notion of two circles "just touching". We do this with a fifth equation. It says that two circles just touch when the distance between their centers equals the sum of the radii. (We avoid a sqrt function, by squaring both sides.) In math notation:

In Python and SymPy:

eq5 = Eq((aprime_x - bprime_x) ** 2 + (aprime_y - bprime_y) ** 2,
         (a_r + b_r)**2)

If we give the first four equations information about two circles and a time, they will tell us the positions of the circles at that time. Equation five can then tell us of the two circles are just touching at that time.

In other words, given a time, these equations can answer the question "Do these two circles just touch at time t?" But that’s not what we want. We want to turn that inside out and answer the question "At what time will these two circles just touch?"

To answer the question of interest, Mathematics tells us to just "solve" the equations. Sadly, to paraphrase Barbie, "Math (that solves equations containing squares) is hard." Happily, the math is not hard for SymPy:

from sympy import nonlinsolve
from perfectphysics import savecc_all_solutions = nonlinsolve([eq1, eq2, eq3, eq4, eq5], t, aprime_x, aprime_y, bprime_x, bprime_y)
cc_time_solutions = [t for t, aprime_x, ap_y, bp_x, bp_y in cc_all_solutions]
save(cc_time_solutions, "cc_time_solutions.sympy")
cc_time_solutions[0]

In about 30 seconds, SymPy finds the two general, symbolic solutions. Here is just bit of the first solution, cc_time_solutions[0] in math notation:

This looks complicated to us. It is not, however, complicated for the computer.

Why are there two solutions? Evaluating them can help us understand. Consider this situation: two circles of radius 1 with the blue ball moving at speed 2:

from perfect_physics import Circle, plot
a = Circle(x=0,y=0,r=1,vx=2,vy=0)
b = Circle(x=4,y=0,r=1,vx=0,vy=0)
plot([a, b], colors=["tab:blue","tab:red"], xlim=(-2,8), ylim=(-2,2), figsize=(4,2), font_scale=1)

We use the SymPy .subs (substitute) method to find the two times when the circles a and b will just touch. We also plot the circles at those times.

from perfect_physics import load
cc_time_solutions = load("cc_time_solutions.sympy")
times = [time_solution.subs([("a_x", a.x), ("a_y", a.y), ("a_r", a.r), ("a_vx", a.vx), ("a_vy", a.vy),
                             ("b_x", b.x), ("b_y", b.y), ("b_r", b.r), ("b_vx", b.vx), ("b_vy", b.vy)])
         for time_solution in cc_time_solutions]
print(times)
for time in times:
    plot([a.tick_clone(time), b.tick_clone(time)], clock=time,
    colors=["tab:blue","tab:red"], xlim=(-2,8), ylim=(-2,2), figsize=(4,2), font_scale=1)

Python prints [3, 1] and plots:

In other words, at time 3 (ignoring collisions), the blue circle will just be touching the red. Likewise, at time 1. If we want the time until the next collision, we just take the smaller time, right?

Sadly, the smaller time might not be right either. Consider this situation, in which the circles are moving away from each other. When will they just touch?

a = Circle(x=2,y=2,r=1,vx=1,vy=1)
b = Circle(x=0,y=0,r=1,vx=0,vy=0)
plot([a, b], colors=["tab:blue","tab:red"], xlim=(-2,4), ylim=(-2,4), figsize=(3,3), font_scale=1)
times = [time_solution.subs([("a_x", a.x), ("a_y", a.y), ("a_r", a.r), ("a_vx", a.vx), ("a_vy", a.vy),
                             ("b_x", b.x), ("b_y", b.y), ("b_r", b.r), ("b_vx", b.vx), ("b_vy", b.vy)])
         for time_solution in cc_time_solutions]
print(times)

Python answers [-2 + sqrt(2), -2 - sqrt(2)]. (Note that it answers with infinite precision rather than something approximate such as [-0.59, -3.41]). These times are negative, meaning the circles will just touch in the past. We don’t care about past collisions, so we ignore negative times.

To summarize: The answer will give us two times and we keep the smallest, non-negative, time, right? Nope. We must also watch out for two more complications. First, when the circles are stationary or moving together at the same speed, the times will be [nan, nan], that is, floating point "not a number".

Second, when the circles are moving such that their paths never cross, the times will be something like [2 + sqrt(2)*I, 2 - sqrt(2)*I], containing imaginary numbers.

In the end, our procedure is to use SymPy to get two times and then:

  • filter out answers containing imaginary numbers
  • filter out nan answers
  • filter out negative answers
  • keep the smallest remaining answer (if any)

Will this result be the infinite-precision time of the next collision we desire? Almost. When it says that they will just touch in time zero, the circles could be moving away from each other. That is a situation that doesn’t count as a collision. We’ll look at this problem in Step 5.

Note that if we can handle two circles, we can handle any number of circles. We just need to look at all pairs of circles and return the smallest time (if any).

I think what we’ve done in this step is very cool. We just had to engineer five simple equations about circles moving forward in time and checking a distance. The Python SymPy package then did the hard math for us. It created a function that tells us exactly when (if ever) the circles will touch.

In Step 3, we’ll see how touching changes the speed of each circle. Next, however, let’s work with SymPy to find when a circle will run into a wall.

Step 2: Model a circle and a wall. Find when they will just touch. Add a "secret" sixth equation about angles (but don’t use angles).

Suppose we have circle a and wall w. As before, the circle has a position (x, y), a velocity (vx, vy), and a radius, r. The wall is an (infinitely long) line. We specify the wall with any two distinct points along it, (x0, y0) and (x1, y1).

To find when the circle and wall will just touch, we again model the circle’s position (x’, y’) at time t.

from sympy import symbols, Eq
t = symbols("t")
a_x, a_y, a_vx, a_vy, a_r, aprime_x, aprime_y = symbols("a_x, a_y, a_vx, a_vy, a_r, a'_x, a'_y")
eq1 = Eq(aprime_x, a_x + a_vx * t)
eq2 = Eq(aprime_y, a_y + a_vy * t)

Next, we model the place on the wall where the circle and wall will touch. Call this point (x2, y2). It will be some proportion p between wall-defining points (x0, y0) and (x1, y1). For example, if p is 1/2, then (x2, y2) is the midpoint of the two wall-defining points. Also, by letting p range less than zero and above 1, we can move (x2, y2) to any place along the wall:

p = symbols("p")
x0, y0, x1, y1, x2, y2 = symbols("x_0, y_0, x_1, y_1, x_2, y_2")
eq3 = Eq(x2, x0 + (x1-x0) * p)
eq4 = Eq(y2, y0 + (y1-y0) * p)

The circle will touch the wall when the distance between its center point (x’, y’) and the place on the wall it touches (x2, y2) is equal to the circle’s radius. In other words:

eq5 = Eq((x2-aprime_x)**2 + (y2-aprime_y)**2, a_r**2)

When we ask SymPy to solve and apply these five equations, however, it gives us strange solutions like this:

SymPy thinks this solves our equations because the circle’s center is one radius away from some point on the wall (the two red dots). This is not, however, the solution we want because the circle overlaps the wall.

A sixth equation can fix the problem. It will require that angle α between the wall and that white line be 90° (see figure above). We can use slopes to say this without explicit angles or trigonometric functions. The slope of the wall is (y1-y0)/(x1-x0). The slope of the figure’s white line is (aprime_y-y2)/(aprime_x-x2). Two lines intersect at 90° when the slope of one is the negative reciprocal of the other. In SymPy, we say:

slope_0 = (y1-y0)/(x1-x0)
slope_1 = (aprime_y-y2)/(aprime_x-x2)
eq6 = Eq(slope_0, -1/slope_1)

We ask SymPy to solve this for us like this:

from sympy import nonlinsolve
from perfect_physics import save
cw_all_solutions = nonlinsolve([eq1, eq2, eq3, eq4, eq5, eq6], t, p, aprime_x, aprime_y, x2, y2)
cw_time_solutions = [t for t, p, aprime_x, aprime_y, x2, y2 in cw_all_solutions]
save(cw_time_solutions, "cw_time_solutions.sympy")

After working for 7 minutes (better it than me!), it returns two solutions. Here is the first solution in math notation (more readable than last time):

When we apply these solutions to the situation above, Python returns time spans [2, 3]:

_Aside: The solution expressions look like something from the quadradic equation. That suggests that a good applied mathematician could have figured this out manually. Good for them. I, however, am happy to solve the problem via computer._

We’ll apply the same rules as in Step 1, filtering out negative, complex, nan, and keeping the smallest time (if any). We may run across one new issue when a circle grazes a wall:

This world will return [nan, zoo]. where zoo is the SymPy symbol for complex infinity. We’ll filter it out because we don’t care about grazing collisions.

Let’s add this result to the result from Step 1. We now can find the exact, infinite-precision, time of the next "just touch" between a circle and either a circle or a wall.

Next, let’s figure out how a circle-circle touch can change the speed and direction of the circles.

Step 3: Try to model circle-to-circle collisions with conservation of momentum and energy. You can’t. Add a "secret" fourth equation about angles (but don’t use angles).

The first time I tried to model the collision between two circles, I thought I would need just two laws: conservation of momentum and conservation of energy. I was wrong. Let’s see why.

Conservation of energy says that the energy before and after a collision should be the same. For our worlds, the only type of energy is kinetic, so E= mv²/2, where m is mass and v is velocity. We’ll assume collisions are elastic, that is, no energy is converted to heat, sound, etc. In math notation, this becomes:

where the mass of circle a is aₘ, its velocity before the collision is (aᵥₓ, aᵥᵧ), and its velocity after the collision is (âᵥₓ, âᵥᵧ). Likewise for circle b.

In SymPy, we write:

from sympy import symbols, Eq
a_x, a_y, a_vx, a_vy, a_r, a_m, ahat_vx, ahat_vy = symbols("a_x, a_y, a_vx, a_vy, a_r, a_m, ahat_vx, ahat_vy")
b_x, b_y, b_vx, b_vy, b_r, b_m, bhat_vx, bhat_vy = symbols("b_x, b_y, b_vx, b_vy, b_r, b_m, bhat_vx, bhat_vy")
energy_before = a_m * (a_vx**2 + a_vy**2) / 2 + b_m * (b_vx**2 + b_vy**2) / 2
energy_after = a_m * (ahat_vx**2 + ahat_vy**2) / 2 + b_m * (bhat_vx**2 + bhat_vy**2) / 2
eq1 = Eq(energy_before, energy_after)

Conservation of momentum gives us two more equations, one for the x dimension and one for the y dimension. Momentum is mass times velocity, so:

eq2 = Eq(a_m * a_vx + b_m * b_vx, a_m * ahat_vx + b_m * bhat_vx)
eq3 = Eq(a_m * a_vy + b_m * b_vy, a_m * ahat_vy + b_m * bhat_vy)

We now have three equations and want SymPy to give us solutions for four unknowns, namely, the after-collision velocities: âᵥₓ, âᵥᵧ, b̂ᵥₓ, b̂ᵥᵧ. That will not work. We need as many equations as unknowns.

The first time I noticed this missing equation, it stumped me. I couldn’t think of any other applicable laws of physics. I eventually figured out the fourth missing "law":

When circles collide, their speed 90° from the direction of impact doesn’t change.

Consider this situation:

from perfect_physics import Circle, plot
a = Circle(x=0, y=0, r=1, vx=1, vy=2, m=1)
b = Circle(x=2, y=0, r=1, vx=0, vy=0, m=1)
plot([a, b], colors=["tab:red", "tab:blue"], xlim=(-2, 4), ylim=(-2, 2), figsize=(4, 2), font_scale=1)

The "direction of impact" is the direction from one center to the other, here horizontal. Circle a’s speed 90° from horizontal is aᵥᵧ which is 2. So, we expect the after-collision speed in this direction, âᵥᵧ to be 2 also.

But what if the circle centers don’t lie horizontally or vertically? Consider this situation:

from sympy import sqrt
from perfect_physics import Circle, plot
a = Circle(x=0, y=0, r=1, vx=1, vy=0, m=1)
b = Circle(x=sqrt(2), y=sqrt(2), r=1, vx=0, vy=0, m=1)
plot([a, b], colors=["tab:red", "tab:blue"], xlim=(-2, 4), ylim=(-2, 3), figsize=(4, 2), font_scale=1)

We could find the speed 90° from the direction of impact with the sine and cosine function. However, doing it with unit vectors lets us avoid introducing those functions. We want a vector <ux,uv> that is one unit long and that runs 90° from the line between the circle centers. We can get the length to be one by dividing by d, the distance between the centers. We can get the direction we want with the same -1/slope we used in Step 2. Specifically, this says the velocity 90° from the direction of impact doesn’t change:

from sympy import sqrt
d = sqrt((b_x-a_x)**2 + (b_y-a_y)**2)
ux = -(b_y-a_y)/d
uy = (b_x-a_x)/d
eq4 = Eq((b_vx - a_vx) * ux + (b_vy - a_vy) * uy, (bhat_vx-ahat_vx) * ux + (bhat_vy-ahat_vy) * uy)

We have SymPy solve our equations:

from sympy import nonlinsolve
from perfect_physics import save
cc_velocity_solutions = nonlinsolve([eq1, eq2, eq3, eq4], ahat_vx, ahat_vy, bhat_vx, bhat_vy)
cc_velocity_solutions = list(cc_velocity_solutions)
save(cc_velocity_solutions[1], "cc_velocity_solution.sympy")

After 2 or 3 minutes, it returns two solutions. Surprisingly, the first solution says that the after velocities equal the before velocities. In other words, it is the solution for when the circles miss each other.

The second solution, the one we care about, is long. I won’t show it here.

Let’s put the solution to work on the two examples from above. In the first example, the velocities before were a.vx=1, a.vy=2, b.vx=0, b.vy=0. After they become a.vx=0, a.vy=2, b.vx=1, b.vy=0.

In the second example, the velocities before were a.vx=1, a.vy=0, b.vx=0, b.vy=0. After they become a.vx=1/2, a.vy=-1/2, b.vx=1/2, b.vy=1/2.

Cool! So, we can now predict how circles will bounce off each other. We’ll see later that this solution works for circles of different masses and different radii.

We next try to apply these conservation laws to circles bouncing off walls. Surprisingly, the bounce seems to violate "conservation of momentum".

Step 4: Assume walls have infinite mass. A circle bouncing off a wall violates the conservation of momentum. Fix this with limits.

For a ball bouncing off a wall, we could just cheat. The answer, which you may remember from school, is that the "angle of incidence is equal to the angle of reflection":

But how can this be right? The circle’s y-momentum (that is, mass times y-velocity) before the bounce is -1. After the bounce it is 1. That violates conservation of momentum. Including the wall’s momentum doesn’t help because its mass is infinite and its y-velocity is 0 giving us an undefined momentum.

Maybe we can just model walls as big circles and use the solution from Step 3. Here we give the big circle a radius of 20 and a mass of 1000.

a = Circle(x=0, y=1, r=1, vx=1, vy=-1, m=1)
b = Circle(x=0, y=-20, r=20, vx=0, vy=0, m=1000)

This almost works. The after velocities are a.vx=1, a.vy=999/1001, b.vx=0, b.vy=-2/1001. What if we stick infinity in as the big circle’s mass and/or radius. Sadly, the result is undefined values.

Let’s see if we can start with the conservation laws and discover the right answer. We’ll assume that the wall has finite mass and speed 0. We put in conservation of energy (eq1), conservation of momentum (eq2, eq3), and that the relative speed of the circle along the wall shouldn’t change.

from sympy import symbols, Eq, sqrt
a_x, a_y, a_vx, a_vy, a_r, a_m, ahat_vx, ahat_vy = symbols("a_x, a_y, a_vx, a_vy, a_r, a_m, ahat_vx, ahat_vy")
w_x0, w_y0, w_x1, w_y1, w_m, what_vx, what_vy = symbols("w_x0, w_y0, w_x1, w_y1, w_m, what_vx, what_vy")
energy_before = a_m * (a_vx**2 + a_vy**2) / 2 + w_m * 0
energy_after = a_m * (ahat_vx**2 + ahat_vy**2) / 2 + w_m * (what_vx**2 + what_vy**2) / 2
eq1 = Eq(energy_before, energy_after)
eq2 = Eq(a_m * a_vx + w_m * 0, a_m * ahat_vx + w_m * what_vx)
eq3 = Eq(a_m * a_vy + w_m * 0, a_m * ahat_vy + w_m * what_vy)
d = sqrt((w_x1-w_x0)**2 + (w_y1-w_y0)**2)
ux = (w_x1-w_x0)/d
uy = (w_y1-w_y0)/d
eq4 = Eq(a_vx * ux + a_vy * uy, (ahat_vx-what_vx) * ux + (ahat_vy-what_vy) * uy)

We again apply nonlinsolve and get one result where they don’t hit (ignored) and a second in which they do. Our next step is new. We use SymPy to take the limit as the wall’s mass goes to infinity.

from sympy import nonlinsolve, limit, oo, Tuple
from perfect_physics import save
no_hit, cw_velocity_solution = nonlinsolve(
    [eq1, eq2, eq3, eq4], ahat_vx, ahat_vy, what_vx, what_vy
)
a_vx_limit, a_vy_limit, w_vx_limit, w_vy_limit = [
    limit(velocity, w_m, oo) for velocity in cw_velocity_solution
]
assert w_vx_limit == 0 and w_vy_limit == 0
cw_velocity_limits = Tuple(a_vx_limit, a_vy_limit)
save(cw_velocity_limits, "cw_velocity_limits.sympy")

The result is the after vx and vy for the circle and the wall. As hoped, the wall’s velocities are zero. We save the formulas for the circle’s after velocities.

When we apply the formula, we get the expect result of "angle of incidence is equal to the angle of reflection":

and

So, are we now ready to solve physics homework problems, create animations, and make games? Not quite. I ignored one issue: two objects touching does not always mean they are colliding. We’ll look at that next.

Step 5: Two objects can be touching and yet moving away from each other. Use limits to find these. Remove them from the list of colliding objects.

How do we detect if two circles are moving toward each other? I originally thought we could just check their relative speed. The relative speed, however, of the red circle is the same in these two situations. In the first they are moving toward each other and in the second they are not.

A better approach is to see if the centers are getting closer together. But when? They may start by getting closer but later move away from each other. Might we ask if they are getting closer at time zero? That doesn’t work because at any moment in time, their distance is unchanging.

The solution is to again use limits. We’ll test if in the instant after time zero, they are getting closer.

from sympy import symbols, limit, sqrt
from perfect_physics import save
a_x, a_y, a_vx, a_vy, b_x, b_y, b_vx, b_vy, t = symbols("a_x, a_y, a_vx, a_vy, b_x, b_y, b_vx, b_vy, t")
d0 = sqrt((a_x - b_x) ** 2 + (a_y - b_y) ** 2)
d1 = sqrt((a_x + t * a_vx - (b_x + t * b_vx)) ** 2 + (a_y + t * a_vy - (b_y + t * b_vy)) ** 2)
speed_toward = (d0 - d1) / t
instant_speed = limit(speed_toward, t, 0)
save(instant_speed, "instant_speed.sympy")

In the first example above, applying the resultant formula returns 1. In the second example, it returns -1. Both results are as expected.

To see if a circle is moving toward a wall, we find (x₂, y₂), the point on the wall where the circle’s center might hit. That requires four equations that can be solved with SymPy linsolve. We then use the instant_speed formula that we just found to see how fast the circle is moving towards or away from that point:

from sympy import symbols, linsolve, Eq, simplify
from perfect_physics import load, save
a_x, a_y, a_vx, a_vy, t = symbols("a_x, a_y, a_vx, a_vy, t")
x0, y0, x1, y1, x2, y2, p = symbols("x_0, y_0, x_1, y_1, x_2, y_2, p")
eq1 = Eq(x2, a_x + a_vx * t)
eq2 = Eq(y2, a_y + a_vy * t)
eq3 = Eq(x2, x0 + (x1-x0) * p)
eq4 = Eq(y2, y0 + (y1-y0) * p)
x2_formula, y2_formula, _, _ = list(linsolve([eq1, eq2, eq3, eq4], (x2, y2, p, t)))[0]
instant_speed = load("instant_speed.sympy")
instant_speed_wall = simplify(instant_speed.subs({"b_x": x2_formula, "b_y": y2_formula, "b_vx": 0, "b_vy": 0}))
save(instant_speed_wall, "instant_speed_wall.sympy")

For example, for this situation we find the circle is moving toward the wall at speed sqrt(5).

from sympy import sqrt
from perfect_physics import Circle, Wall, plot, load
a = Circle(x=-sqrt(2) / 2, y=sqrt(2) / 2, r=1, vx=2, vy=1)
w = Wall(x0=-1, y0=-1, x1=1, y1=1)
plot(circle_list=[a], wall_list=[w], xlim=(-3, 3), ylim=(-2, 3), figsize=(4, 2), font_scale=1)
instant_speed_wall = load("instant_speed_wall.sympy")
instant_speed_wall.subs({"a_x": a.x, "a_y": a.y, "a_vx": a.vx, "a_vy": a.vy,
                         "x_0": w.x0, "y_0": w.y0, "x_1": w.x1, "y_1": w.y1})

We now have everything we need to build a physics engine with infinite precision.

Summing Up Part 3

In this part, we’ve seen how to use the Python SymPy package to find the low-level expressions needed to create a perfect physics engine for our 2-D worlds of circles and wall. We found the expressions for the time when two objects will just touch (if they ever do). When they do touch, we found the expressions for their new velocities.

By going through the details, we saw:

  • Engineering equations is relatively easy but manually solving systems of nonlinear equations would be hard. Happily, the computer does this hard work for us.
  • Conservation of energy and momentum is not enough to predict post-collision velocities. We need a fourth equation saying the velocity 90° from the collision doesn’t change.
  • Our equations often give strange answers. We saw answers of nan, negative time spans, complex numbers, and complex infinity. We interpret these strange answers and then handle them. For example, a negative time span means a collision happened in the past and, so, can be ignored.
  • Twice we needed 90° angles. In both cases, we avoid introducing sine, cosine, etc. by using slopes and unit vectors.
  • Once we needed an infinitesimal time and once we needed an infinite mass. In both cases, the SymPy limit function lets us use the extreme time and mass we need.

We have now created the low-level expressions that provide the foundation we need for a high-level layer. We saw the high-level layer in Part 1 and Part 2. When we put both layers together, we get a perfect physics engine. It will, however, be painfully slow. In Part 4, we’ll speed up the engine a bit (but not enough) and discuss other limitations.

If you have ideas for simulations that you’d like me to run, please send them to me. They may become the basis of a part 5.

You can download this 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 the next parts. All parts will be free to read. Finally, on YouTube, I have older (approximate) physics simulations and some trying-to-be-humorous nature videos.


Related Articles