How can I numerically integrate the Kepler problem?

Question

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[0] = q
    ps[0] = 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.add_collection(plot_trajectory(q_se0, color="tab:green", label="v=1"))
#ax.add_collection(plot_trajectory(q_se1, color="tab:blue", label="v=3"))
#ax.add_collection(plot_trajectory(q_se2, color="tab:red", label="v=5"))
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.

Answer 1

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.

Answer 2

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.

Answer 3

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.