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 = fig.add_subplot(projection="3d")
ax.plot(x, y, z)
# show Sun
ax.scatter(0, 0, 0, c="black")

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

Orbit 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?

  • 1
    $\begingroup$ 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. $\endgroup$ Oct 30, 2023 at 6:35
  • 1
    $\begingroup$ 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. $\endgroup$ Oct 30, 2023 at 6:35
  • $\begingroup$ @LutzLehmann Yes, indeed that fixes it. Would you mind copying your comment as an answer so I can accept it? $\endgroup$
    – JS4137
    Oct 30, 2023 at 6:42

1 Answer 1


You have two separate issues, a formula error and the time. Rectifying the time issue might have highlighted the formula issue.

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.

With RK4 one gets satisfying results with time steps of 5 days, if Earth had a larger excentricity (but then we would not be here to discuss this) the time step might need to be reduced to one day. For the lower-order Euler method, I estimate that a time step of half a day down to an hour should be sufficient.

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.