How to simulate a universe (of dots on a screen)
Netflix premiered a new sci-fi show this year, 3 Body Problem, which brought attention to a whole host of nerdy Wikipedia articles. [No spoilers!] One of those topics was the three-body problem. The TL;DR version is that if you have a gravitational system of three similar-size bodies, e.g. stars, the result is almost always chaos. Regardless of initial conditions, the bodies resist settling into a neat and predictable orbit, like the moon around the earth, or the earth around the sun.
I’d like to simulate a three-body system and visualize its chaotic motion. By defining the fundamental forces (gravity) and everything in it (the stars), we can essentially create our own miniature universe. Of course it’s a massive simplification of the universe we currently occupy, and just talking about it this way opens up a Costco-size can of philosophical worms. But that’s part of the fun.
1. Newton’s Law of Universal Gravitation.
To build a gravitational system we’ll need a basic understanding of gravity. Let’s start with something we can notice intuitively:
- Velocity is the rate of change of position.
- Acceleration is the rate of change of velocity.
It follows that if we know an object’s acceleration then we can calculate its velocity and, in turn, its position. At least assuming we know its initial state.
Okay, but how would we know acceleration?
Isaac Newton realized that the force of gravity affects bodies in space just as it affects apples falling from trees. They only seem like different processes because the scale is so different. He was able to reason that the gravitational force between two objects depends on (1) their mass and (2) the distance between them. More specifically, it’s inversely proportional to the square of the distance between them.
For example, if you tripled the distance between earth and the moon, the earth’s pull on the moon would be 1/9 as strong. And it’s true in the other direction as well. The moon’s pull on earth would be cut in nine.
Newton’s Law of Universal Gravitation can be expressed this way:
Where m1 and m2 are the objects’ masses, d is the distance between their centers, and G is the gravitational constant.
2. Newton’s Second Law.
It turns out that Newton was a really smart guy and he also discovered that force is equal to mass multiplied by acceleration. As we said above, if we know acceleration then we’ll also know velocity and position. This equation includes acceleration so we must be making progress!
Now look at what happens when we apply Newton’s Second Law to the previous formula for gravitational force:
In this case m1 is the mass of a given object and m2 is the other mass applying a force to it. m1 cancels and we’re left with a formula for acceleration due to gravity. Notably, the mass of the object in question makes no difference. It’s accelerated just the same. This helps explain the lesson we’re taught as kids: bowling balls and feathers fall at the same rate (ignoring air resistance).
Of course that skips over a few details, not to mention Einstein’s superseding theory of general relativity. My explanation probably wouldn’t cut it in a classroom.
But the math does closely approximate the universe around us. We can use it to predict, with surprisingly high accuracy, the position of objects in our solar system. Or to predict more tangible things here on earth, like how far a baseball will travel if it’s hit with X velocity at Y angle.
Let’s use the formula to model our own three-body system.
3. Build a simulation.
I want to write the code with minimal imports—nothing beyond basic math functions and Matplotlib to visualize. Since we’ll be working with vectors, it would be easier and more efficient to use a library like NumPy. But I like the idea of breaking down the math as much as possible for this post. Python’s standard library has a built-in option to display graphics: turtle. If this were a real-life demonstration it would be great choice. Instead I’ll turn the output into a GIF and display it on this page.
Start by creating a Star class, of which we’ll make three instances. Recall that we can calculate a body’s acceleration and predict how it will move, but we also need to know where it was and how fast it was moving in the previous moment. So pass an initial state during initialization—position, velocity, and mass (and color for plotting).
The values I set are entirely arbitrary. I played around with the model and found a system that interested me. Freedom to create is one of the best parts of simulating your own universe!
from math import sqrt, sin, cos, atan2
class Star:
def __init__(self, initial_position, initial_velocity, mass, color):
self.x = initial_position[0]
self.y = initial_position[1]
self.velocity_x = initial_velocity[0]
self.velocity_y = initial_velocity[1]
self.mass = mass
self.color = color
self.acceleration_x = 0
self.acceleration_y = 0
[...]
john = Star(initial_position=(-10, 60),
initial_velocity=(0, -1.5),
mass=200,
color="#5F72C9")
paul = Star(initial_position=(10, 60),
initial_velocity=(0, 1.5),
mass=100,
color="#5FC97B")
ringo = Star(initial_position=(0, 85),
initial_velocity=(-3, 0),
mass=2,
color="#DB5C82")
Next we’ll calculate how each Star is affected by the other two. The model’s general structure should look like the code below. It’s important to update every Star’s acceleration before moving any of them.
star_list = [john, paul, ringo]
NUMBER_OF_FRAMES = 275
for n in range(NUMBER_OF_FRAMES):
for star in star_list:
star.update_acceleration()
for star in star_list:
star.update_position()
I’m getting ahead of myself because we haven’t yet defined update_acceleration() and update_position() methods.
Updating acceleration is where we implement Newton’s Law of Universal Gravitation. With our Star instances packaged into a star_list, we can iterate through the list and look at interactions between each pair of celestial bodies. For example, John experiences two forces—one from Paul and one from Ringo. We calculate the two acceleration vectors and add them together (example below).
Remember, at this point we’re only updating acceleration_x and acceleration_y attributes. We do this for all three bodies before changing positions.
class Star:
[...]
def update_acceleration(self):
x_components = []
y_components = []
for star in star_list:
if star != self:
distance_x = star.x - self.x
distance_y = star.y - self.y
distance = sqrt(distance_x**2 + distance_y**2)
acceleration = G_CONSTANT * star.mass / distance**2
angle = atan2(distance_y, distance_x)
x_components.append(acceleration * cos(angle))
y_components.append(acceleration * sin(angle))
self.acceleration_x = sum(x_components)
self.acceleration_y = sum(y_components)
[...]
G_CONSTANT = 1.0
Now we know what every Star is going to do at the next moment in time, so we call update_position().
This method is much simpler than the previous one. It’s the same bit of intuition that got us started: acceleration is the rate of change of velocity and, in turn, velocity is the rate of change of position.
class Star:
[...]
def update_position(self):
self.velocity_x += self.acceleration_x
self.velocity_y += self.acceleration_y
self.x += self.velocity_x
self.y += self.velocity_y
That’s all we need to model the system. Let’s work on visualizing the result.
Each moment (frame) of the simulation is plotted and saved in a frames folder. Once the script finishes I’ll open the images in GIMP and save them as an animated GIF.
Let’s include a “progress bar” at the top of the figure. I think it makes it a little easier to see what’s happening when the GIF loops. Set clip_on=False to allow Matplotlib to draw outside axes limits.
Each Star’s mass is represented by its size on the plot. Although since the smallest mass is only 1% of the largest, we should tweak the scaling to improve visibility.
import matplotlib.pyplot as plt
plt.style.use("orbit.mplstyle")
NUMBER_OF_FRAMES = 275
for n in range(NUMBER_OF_FRAMES):
for star in star_list:
star.update_acceleration()
for star in star_list:
star.update_position()
fig, ax = plt.subplots()
ax.plot([-100, -100 + 200 * (n + 1) / NUMBER_OF_FRAMES], [103, 103], clip_on=False)
for star in star_list:
ax.scatter([star.x], [star.y], s=star.mass * 0.7 + 15, c=star.color, lw=0.5, ec="black")
ax.set(xticks=range(-100, 120, 20),
yticks=range(-100, 120, 20),
xlim=(-101, 101),
ylim=(-101, 101))
plt.savefig(f"frames/{n}.png")
4. The output.
Seeing this in action puts a big smile on my face. I think it’s really cool that we can create such complex behavior from 50-ish lines of Python. And to connect it back to Newton’s brilliant discoveries 300 years ago.
The two larger bodies behave like a binary star system. They almost look like a pair of dancers twirling across the dance floor. Meanwhile the small mass makes wide, chaotic swoops in and out of the action, before being summarily ejected from the system.
And the simulation isn’t limited to three stars or 275 frames. We could write a loop that instantiates a dozen bodies and appends each to star_list. We’re only limited by CPU power and patience.
Now that we’ve created a miniature “universe” with functional gravity, can you imagine adding more forces and additional complexity until it seems as real as our own universe? What if our reality is just someone else’s code running on a much more powerful computer? I know people like to speculate about life being a simulation. It’s a fun idea and I’d like to know the answer, but I don’t think it would change my behavior. Our little bubble is worth caring about regardless of how it came to exist.
Download the Matplotlib style.
Full code:
from math import sqrt, sin, cos, atan2
import matplotlib.pyplot as plt
class Star:
def __init__(self, initial_position, initial_velocity, mass, color):
self.x = initial_position[0]
self.y = initial_position[1]
self.velocity_x = initial_velocity[0]
self.velocity_y = initial_velocity[1]
self.mass = mass
self.color = color
self.acceleration_x = 0
self.acceleration_y = 0
def update_acceleration(self):
x_components = []
y_components = []
for star in star_list:
if star != self:
distance_x = star.x - self.x
distance_y = star.y - self.y
distance = sqrt(distance_x**2 + distance_y**2)
acceleration = G_CONSTANT * star.mass / distance**2
angle = atan2(distance_y, distance_x)
x_components.append(acceleration * cos(angle))
y_components.append(acceleration * sin(angle))
self.acceleration_x = sum(x_components)
self.acceleration_y = sum(y_components)
def update_position(self):
self.velocity_x += self.acceleration_x
self.velocity_y += self.acceleration_y
self.x += self.velocity_x
self.y += self.velocity_y
john = Star(initial_position=(-10, 60),
initial_velocity=(0, -1.5),
mass=200,
color="#5F72C9")
paul = Star(initial_position=(10, 60),
initial_velocity=(0, 1.5),
mass=100,
color="#5FC97B")
ringo = Star(initial_position=(0, 85),
initial_velocity=(-3, 0),
mass=2,
color="#DB5C82")
star_list = [john, paul, ringo]
G_CONSTANT = 1.0
NUMBER_OF_FRAMES = 275
plt.style.use("orbit.mplstyle")
for n in range(NUMBER_OF_FRAMES):
print(n)
for star in star_list:
star.update_acceleration()
for star in star_list:
star.update_position()
fig, ax = plt.subplots()
ax.plot([-100, -100 + 200 * (n + 1) / NUMBER_OF_FRAMES], [103, 103], clip_on=False)
for star in star_list:
ax.scatter([star.x], [star.y], s=star.mass * 0.7 + 15, c=star.color, lw=0.5, ec="black")
ax.set(xticks=range(-100, 120, 20),
yticks=range(-100, 120, 20),
xlim=(-101, 101),
ylim=(-101, 101))
plt.savefig(f"frames/{n}.png")