# Why does a two-body simulation result in no change of the y-component?

I've been attempting to create a model of a heliocentric orbit based on Newton's law of gravitation:

$$\frac{d^2 \vec r}{dt^2} = -\frac{GM}{|r|^2} \hat r$$

This is what I have so far:

import numpy as np
import matplotlib.pyplot as plt

# Heliocentric orbit
u0 = np.array([0.0, 1.515e11, 0.0])

# Approximately 31 km/s
# only x-component because the tangential
# velocity is parallel to x-axis at t=0
v0 = np.array([-3.1e4, 0, 0.0])

N = 5000
u = np.zeros((N, len(u0)))
v = np.zeros((N, len(v0)))
u[0] = u0
v[0] = v0
dt = 0.001
M = 2e30
G = 6.67e-11

# Euler's method
for i in range(N - 1):
u[i + 1] = u[i] + v[i] * dt
v[i + 1] = v[i] -(G * M) / (np.linalg.norm(u[i]) ** 2) * dt

x = u[:, 0]
y = u[:, 1]
z = u[:, 2]

fig = plt.figure()
ax.plot(x, y, z)
# show Sun
ax.scatter(0, 0, 0, c="black")
plt.show()


What I don't understand is why it results in this plot:

This orbit has movement along the $$x$$ and $$z$$ axes, but if I run y[-1] - y[0], it returns 0.0. Since every component contributes to the position vector in the equation of motion, why wouldn't there be movement along one axis?

• Where is your implementation of $\hat r$? Note that dt and especially its square is numerically/floating-point insignificant relative to the size of u0[1]. You are integrating a realistic orbit for 5 seconds in millisecond steps. Note also the scales on the axes, if you set the aspect ratios to 1, you would get a very short nearly straight line segment. Oct 30, 2023 at 6:35
• For RK4 you could take time steps of one day, for Euler a time step of about an hour should be sufficient. Try the suggested changes back-to-front, for an integration interval of some years the orbit without correcting the force should look interesting. Oct 30, 2023 at 6:35
• @LutzLehmann Yes, indeed that fixes it. Would you mind copying your comment as an answer so I can accept it? Oct 30, 2023 at 6:42

The time issue is that you use constants in a meter-second-kilogram system. This means that you integrate for 5 seconds in steps of milliseconds. The step update in the $$y$$ will be below the machine precision relative to the existing value, so that indeed the value remains the same.
The formula issue is that you lost the unit vector in radial direction. This is implicitly replaced with the vector $$(1,1,1)^T$$ per the broadcasting rules of numpy.
In the correct system there are two vectors as input and all changes are via linear combinations of these. This means that the trajectory and the velocity remain in the plane defined by the initial vectors. Here you have 3 vectors as input, so that the trajectory moves in all directions, with the influence of the initial vectors diminishing a little in every step. In the end the planet moves away in the diagonal $$(1,1,1)^T$$ direction. In the transition phase this could possibly give interesting solution pictures.