Volume preservation in phase space in leapfrog algorithm

Based on the information provided in these lecture notes, the leapfrog integration algorithm should preserve area (for single-body systems) in phase space. I implemented the algorithm for a Coulomb force and while the energy is conserved as it should be, the area of the parallelogram defined by four different initial conditions doesn't stay constant as I would expect.

DT = 0.1
T = 20

def inverse_square_law(x, k=1):
    return k / (x ** 2)

def KE(v, m=1):
    return 0.5 * m * v**2

def PE(x, k=1):
    return k / x

def area(a, b):
    return np.abs(a[0] * b[1] - a[1] * b[0])

x0 = np.array([2, 2, 1, 1])
v0 = np.array([1, 0, 1, 0])

x = x0
v = v0

times = []
ke_list = []
pe_list = []
total_energy_list = []
all_positions = [[] for _ in range(len(x))]
all_velocities = [[] for _ in range(len(v))]
areas = []

v = v + 0.5 * DT * inverse_square_law(x)

for t in np.arange(0, T, DT):
    x = x + DT * v
    a = inverse_square_law(x)
    v = v + DT * a

    v_full = v - 0.5 * DT * a
    ke = KE(v_full[0])
    pe = PE(x[0])
    total_energy = ke + pe

    times.append(t)
    ke_list.append(ke)
    pe_list.append(pe)
    total_energy_list.append(total_energy)

    for i in range(len(x)):
        all_positions[i].append(x[i])
        all_velocities[i].append(v_full[i])

    p0 = np.array([x[0], v_full[0]])
    p1 = np.array([x[1], v_full[1]])
    p2 = np.array([x[2], v_full[2]])
    p3 = np.array([x[3], v_full[3]])

    area_val = area(p2 - p0, p1 - p0)
    areas.append(area_val)

Why is that? I should also mention that if I run the algorithm on a simple harmonic oscillator instead, the area is conserved. The graph below illustrates how the area changes over time. Interestingly, the plot I get using Euler's method looks exactly the same.

In addition to the approximate nature of the method, a supplemental source of errors in your simulation is the fact that the shape of the initial rectangular area as it evolves in time is no longer a quadrilateral, as the points on the initial straight borders will no longer be aligned, so the figure transforms into a "quadrilateral with curved boundaries".

Here's a modification of your program showing how the area evolves, by taking multiple points along the borders.

import numpy as np
import matplotlib.pyplot as plt

DT = 0.01
T = 10

def inverse_square_law(x, k=1):
    return k / (x ** 2)

def KE(v, m=1):
    return 0.5 * m * v**2

def PE(x, k=1):
    return k / x

def area(a, b):
    return np.abs(a[0] * b[1] - a[1] * b[0])

x0 = np.concatenate((
    np.linspace(2, 2, num=10),
    np.linspace(2, 1, num=10),
    np.linspace(1, 1, num=10),
    np.linspace(1, 2, num=10)))
v0 = np.concatenate((
    np.linspace(1, 0, num=10),
    np.linspace(0, 0, num=10),
    np.linspace(0, 1, num=10),
    np.linspace(1, 1, num=10)))

x = x0
v = v0

v = v + 0.5 * DT * inverse_square_law(x)

for t in np.arange(0, T, DT):
    x = x + DT * v
    a = inverse_square_law(x)
    v = v + DT * a

    v_full = v - 0.5 * DT * a

# plt.plot(x0, v0, label = "Initial")
plt.plot(x, v, label = "at T=10")
plt.show()

The integration was done until T = 10 in the image above.

An interesting visualisation of the space phase figure as it evolves in time and the conservation of its area, for a pendulum can be found on this Ohio State "About physics 5300" github page by Richard J. Furnstahl.

Only a (theoretical) elementary area dA may be a rectangle that transforms into a quadrilateral. Of course, the smaller the area, the closer we are to the elementary area model (or the elementary area with straight borders becomes a better approximation).

To verify if that is significant, I did some experimentation with your code; I used to represent the relative error the quantity abs(area_val/area0 - 1) where area0 is the initial value of area_val. The same with energy, abs(total_energy/energy0 - 1).

Then I started reducing the size of the zone used to compute the (initial) area, from 1 x 1 by an order of magnitude on both dimensions: 0.1 x 0.1, 0.01 x 0.01, ...
At 1e-5 x 1e-5 the errors in energy and area were of the same approximate magnitude. Reducing the initial area even further, reduced the magnitude of the error in area well below that of energy.

Here's the code and result:

import numpy as np
import matplotlib.pyplot as plt

DT = 0.01
T = 20

def inverse_square_law(x, k=1):
    return k / (x ** 2)

def KE(v, m=1):
    return 0.5 * m * v**2

def PE(x, k=1):
    return k / x

def area(a, b):
    return np.abs(a[0] * b[1] - a[1] * b[0])

# x0 = np.array([2, 2, 1, 1])
# v0 = np.array([1, 0, 1, 0])
# x0 = np.array([1.1, 1.1, 1, 1])
# v0 = np.array([0.1, 0, 0.1, 0])
x0 = np.array([1.00001, 1.00001, 1, 1])
v0 = np.array([0.00001, 0, 0.00001, 0])
# x0 = np.array([1.000001, 1.000001, 1, 1])
# v0 = np.array([0.000001, 0, 0.000001, 0])

x = x0
v = v0

times = []
ke_list = []
pe_list = []
energy_err_list = []
all_positions = [[] for _ in range(len(x))]
all_velocities = [[] for _ in range(len(v))]
area_err_list = []

v = v + 0.5 * DT * inverse_square_law(x)
area0 = -1
energy0 = 0

for t in np.arange(0, T, DT):
    x = x + DT * v
    a = inverse_square_law(x)
    v = v + DT * a

    v_full = v - 0.5 * DT * a
    ke = KE(v_full[0])
    pe = PE(x[0])
    total_energy = ke + pe
    if area0 < 0:
        energy0 = total_energy
        energy_err = 0
    else:
        energy_err = abs(total_energy / energy0 - 1)

    times.append(t)
    ke_list.append(ke)
    pe_list.append(pe)
    energy_err_list.append(energy_err)

    for i in range(len(x)):
        all_positions[i].append(x[i])
        all_velocities[i].append(v_full[i])

    p0 = np.array([x[0], v_full[0]])
    p1 = np.array([x[1], v_full[1]])
    p2 = np.array([x[2], v_full[2]])
    p3 = np.array([x[3], v_full[3]])

    area_val = area(p2 - p0, p1 - p0)
    if area0 < 0:
        area0 = area_val
        area_err = 0
    else:
        area_err = abs(area_val / area0 - 1)
    area_err_list.append(area_err)

plt.plot(times, energy_err_list, label = "Energy Error")
plt.plot(times, area_err_list, label = "Area Error")
plt.legend(loc="lower right")
plt.show()