A three-body simulation can produce a plot that looks “shifted”: the paths seem displaced, misaligned, or unlike the expected motion. That visual symptom does not identify a single cause. It can come from how the initial state is built, how the solver output is interpreted, or how the result is plotted.
This guide gives a small example and a diagnostic process. It uses SciPy’s solve_ivp interface, but the checks also apply to other numerical integrators. The example is illustrative, not a claim that it reproduces the reported issue. Check the documentation for the versions of Python and SciPy you use; Python’s documentation alone does not define SciPy solver behaviour.
First clarify what “shifted” means
Look at the plot and describe the mismatch precisely. Is every body displaced by nearly the same vector? Does the motion start at the wrong position? Does the path have the right shape but appear at the wrong time? Or does it diverge progressively from the expected trajectory?
- Same displacement for every body: check the coordinate origin, initial positions and any translation applied while plotting.
- Wrong starting points: check the order and units of the initial-state values.
- Correct-looking path, wrong time: check the time array and whether the plot uses the solver’s returned times.
- Trajectory changes shape or diverges: check the equations, units, integration interval and numerical accuracy.
A plot can be internally consistent while differing from a reference plot that uses a different origin, coordinate frame, initial condition or time window. Compare those choices before concluding that the solver is wrong.
A compact example to inspect
The following skeleton uses a two-dimensional state. Each body has position (x, y) and velocity (vx, vy). The state is packed as six positions followed by six velocities. The gravitational constant and masses must use units consistent with the positions and velocities.
import numpy as np
from scipy.integrate import solve_ivp
# Illustrative values only. Choose physically consistent units and conditions.
masses = np.array([1.0, 1.0, 1.0])
G = 1.0
# State layout:
# [x1, y1, x2, y2, x3, y3, vx1, vy1, vx2, vy2, vx3, vy3]
y0 = np.array([
-1.0, 0.0,
1.0, 0.0,
0.0, 1.0,
0.0, -0.2,
0.0, 0.2,
0.1, 0.0,
], dtype=float)
def derivative(t, state):
positions = state[:6].reshape(3, 2)
velocities = state[6:].reshape(3, 2)
accelerations = np.zeros_like(positions)
for i in range(3):
for j in range(3):
if i == j:
continue
delta = positions[j] - positions[i]
distance = np.linalg.norm(delta)
# This simple model is singular at a collision.
if distance == 0:
raise ValueError("Two bodies occupy the same position")
accelerations[i] += (
G * masses[j] * delta / distance**3
)
return np.concatenate((velocities.ravel(), accelerations.ravel()))
solution = solve_ivp(
derivative,
t_span=(0.0, 10.0),
y0=y0,
method="RK45",
dense_output=False,
)
if not solution.success:
raise RuntimeError(solution.message)
positions = solution.y[:6].T.reshape(-1, 3, 2)
times = solution.t
for body in range(3):
x = positions[:, body, 0]
y = positions[:, body, 1]
# Plot x against y, using your preferred plotting library.
The crucial convention is that the derivative returns values in the same order as the state: position derivatives first (the velocities), then velocity derivatives (the accelerations). If your own state uses another layout, update both the unpacking and the returned derivative to match it exactly.
Likely causes of a shifted-looking result
1. Initial conditions are not packed or unpacked consistently
A common source of confusion is switching between per-body records such as [x, y, vx, vy] and grouped arrays such as [all positions, all velocities]. Both layouts can work, but mixing them silently gives the solver a different system from the one you intended.
Write down the state layout in one comment and use it in the initial-state construction, derivative function and plotting code. Check a few entries by hand: confirm that the first position values really belong to body 1 and that the corresponding velocity values are read from the intended indices.
2. The plot applies an origin or coordinate change
Plotting code may subtract a centre of mass, shift coordinates to a chosen body, or display a transformed frame. Those operations can make a correct trajectory look displaced relative to raw coordinates. Search for expressions such as x - x.mean(), position - reference_position, or a centre-of-mass calculation.
As a diagnostic, plot the raw positions directly, with no translation or rotation. If that changes the apparent offset, inspect the transformation and label which coordinate frame each plot uses. Translating all bodies by the same constant changes the displayed origin, not their relative positions.
3. Solver output is indexed incorrectly
With solve_ivp, the returned solution array stores state components by row and output time by column. In the example, solution.y[:6].T turns the first six state components into one row per time sample. Omitting or misplacing the transpose can lead to plotting the wrong dimension or pairing values incorrectly.
Use solution.t for the returned time points. Do not assume that the solver returns samples at every point in a separate plotting array unless you explicitly requested or constructed those times. If you use t_eval, check that the requested times are ordered and lie within the integration interval.
4. Units, constants or equations do not match the reference
Newtonian gravitational acceleration depends on the masses, gravitational constant and separation. Mixing units—for example, positions in one unit system and velocities or the gravitational constant in another—changes the trajectory. A sign error in the displacement vector or an incorrect exponent in the distance term can also alter the motion.
Compare the equations and units independently of the integrator. In particular, verify that acceleration points toward the other body and that the distance factor is the one intended by your model. Do not change solver tolerances to compensate for an equation or unit mismatch.
5. The numerical result is being compared too strictly
Numerical trajectories are approximations, and a three-body system can be sensitive to small changes in initial conditions. A small difference can grow over time, especially in a close encounter. First compare the early part of the integration and use the same initial conditions, units, frame and time interval as the reference.
Only after those checks should you investigate numerical settings. Review the current SciPy documentation for the available solve_ivp options and their meaning. Do not assume that changing the method or tightening a tolerance will fix a plotting or model error.
Safe diagnostic path
- Record the state convention. List the indices for each position and velocity. Confirm that
y0follows that exact order. - Check the initial sample. Compare the first plotted position of each body with the initial positions. If they differ, inspect reshaping, slicing and any coordinate transformation.
- Check array dimensions. Print or inspect the shapes of
solution.t,solution.yand the reshaped positions. There should be one position sample per returned time for each body. - Plot raw coordinates. Temporarily remove centring, normalisation and animation offsets. Plot each body’s raw
xagainst rawy. - Check solver completion. Inspect
solution.successandsolution.message. A completed integration does not prove that the equations are physically correct, but a failed one should not be treated as a complete trajectory. - Check physical consistency. Verify masses, units, gravitational constant, signs and initial velocities. Compare like with like when using a reference trajectory.
- Investigate accuracy last. If the model and plotting checks pass but the trajectory still differs, consult the SciPy documentation for the installed version and review its solver settings.
How to verify a correction
Verification should test more than whether the picture looks familiar:
- Confirm that the first plotted sample matches the specified initial positions.
- Confirm that the time values used for the plot come from the same integration output as the positions.
- Plot raw coordinates before applying any reference-frame transformation.
- Check that the result changes in the expected way when you deliberately alter one initial position or velocity.
- If using a physical reference, match its masses, units, initial conditions, coordinate frame and integration interval before comparing paths.
These checks help isolate data-layout and plotting errors. They do not, by themselves, prove that an implementation is an accurate physical model.
Failure cases to keep in mind
- Close approaches: the force law in the example is singular when two positions coincide. A collision or near-collision needs careful modelling and numerical treatment; do not silently replace zero distance with an arbitrary value.
- Long integrations: a trajectory may become sensitive to small numerical differences. A visual mismatch late in the run is not enough to identify a bug.
- Invalid output assumptions: if the solver fails or stops before the requested end time, the returned data may cover only part of the intended interval. Check the completion status and actual returned times.
Common follow-up questions
Should I switch away from RK45?
Not as the first step. First verify the state layout, equations, units and plotting. If those are correct, consult the current SciPy documentation and choose a method and settings suitable for the problem and accuracy needs. The right choice depends on the model and integration, so a method change is not a universal fix.
Why does the reference plot have a different centre?
It may use a different coordinate origin or reference frame. Compare the raw initial positions and then check whether either plot subtracts a reference point or centre of mass. A uniform translation can move the whole picture without changing relative separations between bodies.
References
- Source question about a three-body RK45 plot that appears shifted
- Official Python documentation (consult the documentation for your Python version)
- SciPy documentation for
solve_ivp(consult the documentation matching your installed SciPy version)
Dont SPAM