Code

The three-body problem in thirty lines of Python, no libraries, and it draws its own picture. Copy it, run it, change the numbers.

The program below is the same arithmetic as everything else on this site: work out the pull on each body, kick the speeds by half a step, drift the positions a whole step, kick again. It starts the three bodies on the figure-eight orbit, runs one full period, prints the track as characters, and then prints the energy at the start and at the end so you can see for yourself that the stepping did not leak.

Download three_body.py On GitHub

#!/usr/bin/env python3
"""The whole three-body problem, in thirty lines, with no libraries.

Run it:      python3 three_body.py
It prints the figure-eight orbit as characters, and the energy, which should not move.

The physics is the line every page on this site is about: each body is pulled toward
each other body, by their weights divided by the square of the distance. The stepping
is leapfrog — half a kick, a full drift, half a kick — which is what keeps the energy
from walking away.  https://nanobotco.github.io/three-body/
"""

G = 1.0
M = [1.0, 1.0, 1.0]                                    # three equal weights
R = [[0.97000436, -0.24308753],                        # the starting places
     [-0.97000436, 0.24308753], [0.0, 0.0]]
P = [[0.466203685, 0.43236573]] * 2 + [[-0.93240737, -0.86473146]]    # the starting speeds
T, H = 6.32591398, 0.0002                              # one full period, and the step


def pull(r):
    """The acceleration of each body: the line the whole problem is written in."""
    a = [[0.0, 0.0] for _ in r]
    for i in range(len(r)):
        for j in range(len(r)):
            if i != j:
                dx, dy = r[j][0] - r[i][0], r[j][1] - r[i][1]
                d3 = (dx * dx + dy * dy) ** 1.5
                a[i][0] += G * M[j] * dx / d3
                a[i][1] += G * M[j] * dy / d3
    return a


def energy(r, p):
    e = sum(0.5 * M[i] * (p[i][0] ** 2 + p[i][1] ** 2) for i in range(3))
    for i in range(3):
        for j in range(i + 1, 3):
            dx, dy = r[j][0] - r[i][0], r[j][1] - r[i][1]
            e -= G * M[i] * M[j] / (dx * dx + dy * dy) ** 0.5
    return e


def step(r, p, h):
    """Half a kick, a whole drift, half a kick."""
    a = pull(r)
    for i in range(3):
        for k in (0, 1):
            p[i][k] += 0.5 * h * a[i][k]
            r[i][k] += h * p[i][k]
    a = pull(r)
    for i in range(3):
        for k in (0, 1):
            p[i][k] += 0.5 * h * a[i][k]


if __name__ == "__main__":
    r = [q[:] for q in R]
    p = [q[:] for q in P]
    e0 = energy(r, p)
    grid = [[" "] * 78 for _ in range(15)]
    marks = ".oO"
    for n in range(int(T / H)):
        step(r, p, H)
        for i, (x, y) in enumerate(r):
            col, rowi = int(38 + x * 36), int(7 - y * 19)
            if 0 <= col < 78 and 0 <= rowi < 15:
                grid[rowi][col] = marks[i]
    print("\n".join("".join(row) for row in grid))
    print(f"one period of the figure eight, {int(T / H):,} steps of {H}")
    print(f"energy at the start {e0:.9f}   at the end {energy(r, p):.9f}")
    print(f"the three bodies are back within "
          f"{max(abs(r[i][k] - R[i][k]) for i in range(3) for k in (0, 1)):.2e} of where they set off")

What it prints

       ................                              OOOOOOOOOOOOOOOO         
   .....              .....                      OOOOO              OOOOO     
 ooo                      ....                OOOO                      OOO   
oo                           ...            OOO                           OO  
o                              ...        OOO                              OO 
o                                ...    OOO                                 O 
o                                  ...OOO                                   O 
o                                  ooo...                                   O 
o                                ooo    ...                                 O 
o                              ooo        ...                              OO 
oo                           ooo            ...                           OO  
 ooo                      oooo                ....                      .OO   
   ooooo              ooooo                      .....              .....     
       oooooooooooooooo                              ................         
                                                                              
one period of the figure eight, 31,629 steps of 0.0002
energy at the start -1.287141992   at the end -1.287141992
the three bodies are back within 1.06e-04 of where they set off

Things to change

The step, H. Make it ten times bigger and watch the energy at the end move. Make it a hundred times bigger and the orbit falls apart. That is the whole numbers chapter in one edit.
The starting speeds, P. Nudge the last digit and run it again. Nothing happens at first, and by the end of the period the picture is different — that is the chaos chapter, at home, in one edit.
The weights, M. Make one of them 1.2. The figure eight needs all three equal, so it comes apart; how long it takes to is a fair question and the program answers it.
Swap the stepper. Replace the three lines in step with position first, speed second — the obvious way — and watch the energy climb instead of hold.

The rest of the toolchain

The site is built by the same kind of code, a little longer:

tools/physics.pythe equations, the leapfrog, a fourth-order stepper for when the answer has to be good, and the batch version that moves forty thousand trajectories at once tools/find_orbits.pythe scan over starting speeds and the Gauss-Newton solver that walks each near-miss in until it closes tools/draw.pythe photographs that are not photographs: every raster on the site is a trajectory laid down as light tools/svg.pythe diagrams, computed while the page is written, which is why the numbers in the captions match the pictures above them tools/anim.jsthe demos: the same steppers in the browser, including the fourth-order one and the step chosen from the closest pair tools/build.pythis page and the rest of them

All of it is on GitHub, MIT for the code and CC BY 4.0 for the words, pictures and data. It needs Python 3, numpy and pillow; the program above needs nothing at all.