# How can I numerically integrate the Kepler problem?

I tried to solve a simple Kepler problem numerically.

I have discrete time steps, a starting position $$(x_0,y_0)$$ and starting velocity $$(u_0, v_0)$$.

I used this iteration by calculating the forces and applying Newton's law:

$$F = k/r^2 = k/(x^2+y^2)$$ $$F_x = F \cdot x/\sqrt{x^2+y^2}$$ $$F_y = F \cdot y/\sqrt{x^2+y^2}$$

$$u_{i+1} = u_{i}+ F_x/m \cdot \Delta t$$ $$v_{i+1} = v_{i}+ F_y/m \cdot \Delta t$$

$$x_{i+1} = x_{i}+ (u_{i+1}+u_i)/2 \cdot \Delta t$$ $$y_{i+1} = y_{i}+ (v_{i+1}+v_i)/2 \cdot \Delta t$$

I used the mean velocity between the last and actual time to calculate the new (x, y).

Naively, I thought, that by making the timestep small enough this would give good results. But I observed, that this method is completely bad, because the total energy tends to increase ant what I get is not what I expect: What method for integrating such problem is the standard? My idea is obviously complete garbage because it is the numerically most unstable. When I calculate the total energy it is not constant at all. Even when I make the timestep incredible small, I get totally dissatisfying results.

I read the comment of Daniel Shapero and had a look on his website:

Symplectic integrators

I did a run of Python code using what is referred to as "semi_explicit_euler" for a problem with two fixed masses and one charge coming from the left:

Masses are at x,y=0,5 and 0,-5, initial position: -10,-1, initial speed: 1, 0.

The result is as follows:  You see that the moving mass stops at the right and gets back. But shouldn't the way back be the same as the way as to the point where it stops? The motion must be symmetric to time reversal. Interestingly this "discrepancy" is stable and can not be made smaller by increasing the number of steps. The same picture is obtained when the number of steps is only 25% of that: The code is as follows:

import numpy as np
from numpy import pi as π
import tqdm
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection

q_0 = np.array([-10.0, -1])
p_0 = np.array([1.0, 0])
p_1 = np.array([3.0, 0])
p_2 = np.array([5.0, 0])

r1 = np.array([0, 5])
r2 = np.array([0, -5])

final_time = 80.0
num_steps = 2000000*4
dt = final_time / num_steps

def gravitational_force(q, r1, r2):

return -4 * π ** 2 * (q-r1) / np.sqrt(np.dot((q-r1), (q-r1))) ** 3  -4 * π ** 2 * (q-r2) / np.sqrt(np.dot((q-r2), (q-r2))) ** 3
#return  -4 * π ** 2 * (q-r2) / np.sqrt(np.dot((q-r2), (q-r2))) ** 3

def semi_explicit_euler(q, p, dt, num_steps, r1, r2, force, progressbar=True):
qs = np.zeros((num_steps + 1,) + q.shape)
ps = np.zeros((num_steps + 1,) + p.shape)

qs = q
ps = p

iterator = tqdm.trange(num_steps) if progressbar else range(num_steps)
for t in iterator:
qs[t + 1] = qs[t] + dt * ps[t]
ps[t + 1] = ps[t] + dt * force(qs[t + 1], r1, r2)

return qs, ps

def plot_trajectory(q, start_width=1.0, end_width=3.0, **kwargs):
points = q.reshape(-1, 1, 2)
segments = np.concatenate([points[:-1], points[1:]], axis=1)
widths = np.linspace(start_width, end_width, len(points))
return LineCollection(segments, linewidths=widths, **kwargs)

def energies(qs, ps, r1, r2):
kinetic = 0.5 * np.sum(ps ** 2, axis=1)
potential = 0
potential = -4 * π ** 2 / np.sqrt(np.sum((qs-r1) ** 2, axis=1))
potential = potential - 4 * π ** 2 / np.sqrt(np.sum((qs-r2) ** 2, axis=1))
return kinetic + potential

q_se0, p_se0 = semi_explicit_euler(q_0, p_0, dt, num_steps, r1, r2, gravitational_force)
#q_se1, p_se1 = semi_explicit_euler(q_0, p_1, dt, num_steps, r1, r2, gravitational_force)
#q_se2, p_se2 = semi_explicit_euler(q_0, p_2, dt, num_steps, r1, r2, gravitational_force)

fig, ax = plt.subplots()
ax.set_aspect("equal")
ax.set_xlim((-30, +30))
ax.set_ylim((-30, +30))
#ax.get_xaxis().set_visible(False)
#ax.get_yaxis().set_visible(False)
ax.grid(True)

ax.legend(loc="upper right")

fig, ax = plt.subplots()
ts = np.linspace(0.0, final_time, num_steps + 1)
ax.plot(ts, energies(q_se0, p_se0, r1, r2), color="tab:green", label="v=1")
#ax.plot(ts, energies(q_se1, p_se1, r1, r2), color="tab:blue",  label="v=3")
#ax.plot(ts, energies(q_se2, p_se2, r1, r2), color="tab:red",   label="v=5")
ax.set_xlabel("time (steps)")
ax.set_ylabel("Energy")
ax.legend()

plt.show()


I would have thought that errors gets smaller when increasing number of steps. But this isn't the case - why?

This second question has been solved: The mass doesn't stop, as pointed out in a comment below.

I'm going to assume for the moment that your code is correctly implemented and that this problem isn't a bug.

Believe it or not, gradual increase of the energy is the expected behavior of most simple explicit time integration schemes. If you were to hold the final time fixed, halve the timestep, and double the number of steps, you would get a better approximation to the trajectory -- in other words, if you're using a convergent discretization scheme, the numerical trajectory will approach the (energy-conservative) true trajectory in the limit of small timesteps. That doesn't mean that any individual numerical trajectory with finite timestep will be energy-conservative, especially for very long integration times. Higher-order explicit Runge-Kutta methods can give a vast improvement but they still suffer from the same problem.

There are, however, certain time integration schemes that do a better job solving Hamiltonian systems; these are called symplectic methods. While there is no scheme that guarantees exact conservation of an arbitrary Hamiltonian $$H$$, what makes symplectic methods so convenient is that they sample exactly from the trajectory of a slightly perturbed Hamiltonian $$H' = H + \delta t\cdot\delta H$$. This implies certain really nice properties that we expect from the theory of Hamiltonian systems, i.e. phase space volume preservation, bounded orbits, all equilibria are centers or saddles, etc.

The simplest symplectic method is the semi-implicit Euler scheme, which conveniently is very close to the one you already wrote down: \begin{align} \mathbf{u}_{i + 1} & = \mathbf{u}_i + m^{-1}\delta t\cdot\mathbf{F}(\mathbf{x}_i) \\ \mathbf{x}_{i + 1} & = \mathbf{x}_i + \delta t\cdot\mathbf{u}_{i + 1} \end{align} I wrote a bit more about this on my website, including why symplectic integrators can exist at all, how to come up with a 2nd-order scheme, and an example from molecular dynamics.

• will take a while for reading your website, but it looks very informative. Dec 1, 2022 at 22:48
• Sympletic methods are ideal for systems with low (or no) "friction". With higher friction it is a harder call. Dec 2, 2022 at 8:42
• @Daniel Shapero: I used your code (see your website) and updated my question. Can you have a look on it? The basic question is: Why does decreasing stepsize not lead to a converging result? both green trajectories are equally wrong, although there is a factor 4 in number of steps... Dec 3, 2022 at 10:09
• The question behind the last comment has been solved. Seems that the curve is true and what I identified as a stop isn't one. Dec 3, 2022 at 11:31

For a properly formulated second-order two-stage method that is suitable for this field, take the midpoint method. Set $$\mu=-k/m$$ to stay compatible with other texts.

def acc(x): r3 = sum(x**2)**1.5; return (-mu/r3)*x

def RK2step(x,v,dt):
k1v = acc(x)               ## k1x = v
k2v = acc(x+0.5*v*dt)      ## k2x = v + 0.5*k1v*dt
return x+(v+0.5*k1v*dt)*dt, v+k2v*dt


Then the simulation of a circular orbit with $$\mu=1$$, $$x_0=(0,1)$$, $$v_0=(1,0)$$ gives for moderately large step sizes over about 10 revolutions and with some other step sizes and 60 revolutions As one can see, the range of step sizes from where the solution explodes to where it is circular is rather small.

With higher order methods this visualization exercise gets more boring and misleading.

• With a step size that ensures a dense sampling of the circle for plotting as polygon, the numerical error is so small that clear deviations become visible only after some very long time.

• Conversely, choosing a step size that shows visible deviations from the circle in a short time, it has also to be so large that there will be only some 10 or less points per revolution. But that can only be properly appreciated when using some cubic or higher order interpolation on the segments. Otherwise, for step sizes in-between that still give points on the circle for a long time, the line segments dip into the circle so that in the overlay of a larger number of revolutions only an annulus can be seen, and not the expected sharp circle.

• Thanks, have to read it carefully. Maybe you can have a look on my updated question? Dec 3, 2022 at 10:11
• This would be better as a new topic, the addendum now dominates the original question. I see nothiing wrong in code and results. Note that the second fly-by is much closer to the fixed body than the first, so the trajectory is not actually symmetric. Also, if you zoom in to the return point, you should see that it is not a cusp but a fold, so the return path is different from the arriving path from the start. Also, what is the x coordinate of the last computed point, does it even end up close to the start? Dec 3, 2022 at 10:31
• The loop ends close to the start for end time $105$, and it is not periodic. // In some sense the result is quite ridiculous, you give a step size of $h=10^{-6}$, the energy should be generally constant up to trends of $O(h^2)$ with bounded oscillations of size $O(h)$. The plot, ignoring the spikes, is more compatible with a step size $10^{-4}$. Dec 3, 2022 at 11:36
• To maintain the level of ridicule, I tested with a variable step Heun (using a 3rd order error estimator and applying the 3rd order update) and unit-step error target $\tau=10^{-6}$. This results in about 74000 steps, average size 1e-3 to 1e-2, with spikes down to 1e-5. The energy slowly rises in this interval with a difference of 3e-8, without spikes. The actual plots of course look identical to yours. Dec 3, 2022 at 11:39
• The second Hairer book, Hairer, Wanner: "Solving ODE II: stiff and DAE" is a good general resource. For simple step size control with low-order methods like the Fehlberg23 method, see personal.math.ubc.ca/~feldman/math/vble.pdf Dec 4, 2022 at 11:56

There are better approximations to the velocity update that you can do than your naive implementation, which is called the Euler method and belongs to a family of numerical algorithms called the Runge–Kutta algorithms. The higher-order Runge-Kutta methods do more intelligent weighting of the velocities at various points along your $$dt$$.

See e.g. Runge-Kutta Method.

I'm rather certain that each step along the orbit will overshoot in the tangential direction, and that's why you see this characteristic "loosening" of the orbit.

• I was going to add that for the Kepler problem, you need a symplectic Runge Kutta method which preserves the Hamiltonian, but nevertheless regular RK4 would be a vast improvement over Euler's. Dec 1, 2022 at 19:39