Put masses 3, 4, and 5 at the vertices opposite the matching sides of a 3–4–5 triangle. Set every velocity to zero. Burrau's Pythagorean problem turns that schoolroom triangle into repeated near-collisions, temporary binaries, and an eventual escape. The classical benchmark uses singular point masses; this interactive experiment uses a disclosed Plummer softening length of 10⁻⁴ so its float64 trajectories can continue through the closest passages.
Nothing in the playback below is a prescribed orbit. It contains twenty complete trajectories read from AMBER RunResults.agent_vars: five choices for the first mass crossed with four horizontal offsets. Each one continues until a bound pair and an outward, positive-energy escaper remain separated for ten time units. Move either slider and the page loads another recorded AMBER run. No force calculation occurs in the browser.
A tiny nudge changes who escapes
In the softened control, masses 3 and 4 form the final binary while mass 5 escapes at t=184. Shift mass 3 left by only 0.0005 length units and mass 3 escapes at t=131.6 instead. Shift it right by 0.001 and the system keeps scattering until t=464.6 before mass 3 escapes. The early motion can look nearly identical while later close encounters completely reorder both the ending and how long it takes to arrive.
| m₁ | Δx₁ | Stable at | |ΔE/E| | Outcome |
|---|---|---|---|---|
| 3.0 | 0 | 184.0 | 3.45×10⁻⁵ | 3+4 binary; mass 5 escapes |
| 3.0 | −0.0005 | 131.6 | 1.24×10⁻⁵ | 4+5 binary; mass 3 escapes |
| 3.0 | +0.001 | 464.6 | 3.09×10⁻⁴ | 4+5 binary; mass 3 escapes |
| 2.95 | −0.0005 | 84.0 | 1.22×10⁻⁵ | 2.95+5 binary; mass 4 escapes |
These are softened sensitivity experiments, not the arbitrary-precision reference solution. Boekholt and Portegies Zwart used the unsoftened problem to test N-body reliability. They needed a tolerance of 10⁻¹⁴ and 88-bit arithmetic to converge the first three decimal places; a double-precision Hermite run could conserve energy near 10⁻⁸ and still take a different path after the final encounter. The historical calculation likewise required explicit close-encounter regularization. Here the 10⁻⁴ softening is the regularizer, and its effect on a chaotic outcome is part of what the experiment exposes.
Exact initial state
# Exact baseline; controls select finite grids for m1 and dx
mass = np.array([m1, 4.0, 5.0]) # m1 ∈ {2.9, 2.95, 3, 3.05, 3.1}
x = np.array([1.0 + dx, -2.0, 1.0]) # dx ∈ {-.001, -.0005, 0, .001}
y = np.array([ 3.0, -1.0, -1.0])
vx = np.zeros(3)
vy = np.zeros(3)
# Translate each choice to its center-of-mass frame.
x -= np.sum(mass*x)/np.sum(mass)
y -= np.sum(mass*y)/np.sum(mass)
# baseline pair distances: 5, 4, 3
# each mass equals the opposite side length
# Plummer softening: epsilon = 1e-4
# point-mass energy limit: -769/60 = -12.816666...At the center slider settings, the separations are exactly 5 between masses 3 and 4, 4 between masses 3 and 5, and 3 between masses 4 and 5. With zero initial velocities, the angular momentum is zero; the point-mass energy limit is −769/60, while the softened value differs below the displayed precision. The other nineteen settings change only the selected mass and offset, then translate the system back to its center-of-mass frame.
What AMBER actually does
AMBER is the agent engine, recorder, parameter layer, and execution lane; it is not a built-in celestial-mechanics package. The softened force law, adaptive step-doubling integrator, and stability test are explicit model code. At each AMBER step, the model advances 0.02 time units, commits all four state columns together, and records positions, velocities, energy, and closest separation.
def derivatives(state, masses):
x, y, vx, vy = np.split(state, 4)
ax, ay = accelerations(x, y, masses) # Plummer-softened gravity, ε=1e-4
return np.r_[vx, vy, ax, ay]
def rk4_step(state, h, masses):
k1 = derivatives(state, masses)
k2 = derivatives(state + .5*h*k1, masses)
k3 = derivatives(state + .5*h*k2, masses)
k4 = derivatives(state + h*k3, masses)
return state + h*(k1 + 2*k2 + 2*k3 + k4)/6
def advance_interval(state, masses, interval, rtol, h):
remaining = interval
while remaining > 1e-14:
h = min(h, remaining)
whole = rk4_step(state, h, masses)
half = rk4_step(rk4_step(state, h/2, masses), h/2, masses)
scale = .05*rtol + rtol*np.maximum(abs(state), abs(half))
error = np.max(abs(half - whole)/(15*scale))
if error <= 1:
state = half + (half - whole)/15
remaining -= h
h = min(.01, h*(2.5 if error == 0 else min(2.5, .9*error**-.2)))
else:
h *= max(.1, .9*error**-.2)
return state, h
class PythagoreanTriple(am.Model):
params = {'rtol': (float, 1e-9),
'mass3': (float, 3.0),
'dx': (float, 0.0)}
agent_reporters = ['x', 'y', 'vx', 'vy']
model_reporters = {
'energy': energy_reporter,
'min_separation': separation_reporter,
}
record_initial = True
def setup(self):
self.next_substep = .01
x, y, mass = initial_conditions(self.p.mass3, self.p.dx)
self.add_agents(3, x=x, y=y,
vx=np.zeros(3), vy=np.zeros(3),
mass=mass)
def step_vectorized(self):
x, _ = self.agents.borrow('x'); y, _ = self.agents.borrow('y')
vx, _ = self.agents.borrow('vx'); vy, _ = self.agents.borrow('vy')
mass, _ = self.agents.borrow('mass')
state, self.next_substep = advance_interval(
np.r_[x, y, vx, vy], masses=mass, interval=.02,
rtol=self.p.rtol, h=self.next_substep)
x, y, vx, vy = np.split(state, 4)
self.agents.commit(x=x, y=y, vx=vx, vy=vy)
results = PythagoreanTriple({'steps': 40000, 'rtol': 1e-9,
'mass3': 3.0, 'dx': 0.0001}) \
.cpu(mode='vectorized').run()The stopping test
# Evaluated on recorded AMBER positions and velocities.
stable_now = (
pair_energy < 0 # a bound binary
and escaper_energy > 0 # unbound from that binary
and outward_speed > 0
and escaper_distance >= 12
and escaper_distance >= 8*binary_separation
)
streak = streak + 0.1 if stable_now and same_pair else 0
if streak >= 10:
stable_at = time # playback endpoint
The complete generator runs all twenty grid points against ambr==0.4.4, tests stability from recorded positions and velocities, and retains synchronized playback frames for comparison. Sampling is dense through t=80 and coarser during long excursions. The browser interpolates those immutable coordinates and draws their trails; it contains no force calculation. The payload checksum is written into the generated data itself.
What is—and is not—being claimed
- All twenty selectable trajectories were genuinely executed and recorded by AMBER 0.4.4.
- The model uses G=1 and a disclosed Plummer softening length of 10⁻⁴; it is not the singular point-mass reference problem.
- “Stable” requires one bound pair plus a positive-energy outward escaper at least 12 units and eight binary separations away, continuously for ten time units.
- The browser performs smooth interpolation only; it never advances the dynamical state.
- The float64 grid is a sensitivity explorer, not an arbitrary-precision converged reference solution.
Sources: Szebehely and Peters' 1967 complete numerical solution, the modern N-body reliability study, and the Scholarpedia overview of the general problem.