Runge-Kutta Method for Coulomb's Trouble

Question

我试图用Euler方法和Runge-Kutta四阶方法来模拟二维两电荷情况。使用这两种方法，我得到了相对预期的答案。但我从未尝试过在不同的初始条件下使用RK4方法。你知道吗

从一个正的初始速度vx0和一个负的初始x位置x0开始，代码似乎工作得很好。但是当我翻转它们，使vx0为负，x0为正，我得到不同的答案，当它们应该是对称的。你知道吗

我确定我的欧拉方法对初始条件的两种变化都有效，以确认问题出在我的RK4函数上。最初我很难将RK4方法应用于二维运动，所以出现这种错误并不奇怪。下面的图片可能会让问题更清楚一些。你知道吗

这里是我计算常数的地方：

# Initial conditions
time[0] = t = t0
vx1[0] = vx = vx0
vy1[0] = vy = vy0
x1[0] = x = x0
y1[0] = y = y0

for i in range(1, n + 1):
    # Calculate our constants
    k1vx = step * ax(x, y, q1, q2, m)
    k1vy = step * ay(x, y, q1, q2, m)
    k1x = step * vx
    k1y = step * vy
    k2vx = step * ax(x + 0.5 * k1x, y + 0.5 * k1y, q1, q2, m)
    k2vy = step * ay(x + 0.5 * k1x, y + 0.5 * k1y, q1, q2, m)
    k2x = step * (vx + (k1vx / 2))
    k2y = step * (vy + (k1vy / 2))
    k3vx = step * ax(x + 0.5 * k2x, y + 0.5 * k2y, q1, q2, m)
    k3vy = step * ay(x + 0.5 * k2x, y + 0.5 * k2y, q1, q2, m)
    k3x = step * (vx + (k2vx / 2))
    k3y = step * (vy + (k2vy / 2))
    k4vx = step * ax(x + k3x, y + k3y, q1, q2, m)
    k4vy = step * ay(x + k3x, y + k3y, q1, q2, m)
    k4x = step * (vx + k3vx)
    k4y = step * (vy + k3vy)

    # Update the values based on our calculated constants
    vx1[i] = vx = vx + (k1vx + k2vx + k2vx + k3vx + k3vx + k4vx) / 6
    vy1[i] = vy = vy + (k1vx + k2vy + k2vy + k3vy + k3vy + k4vy) / 6
    x1[i] = x = x + ((k1x + 2 * k2x + 2 * k3x + k4x) / 6)
    y1[i] = y = y + ((k1y + 2 * k2y + 2 * k3y + k4y) / 6)

    # Update the time
    time[i] = t = t0 + i * step

下面是我在前面的代码中用于ax和ay的函数

def accel_rk4_x(x, y, q1, q2, m):
    const = (q1 * q2) / (4 * math.pi * 8.854e-12 * m)
    return const * (x / ((x ** 2 + y ** 2) ** 1.5))


def accel_rk4_y(x, y, q1, q2, m):
    const = (q1 * q2) / (4 * math.pi * 8.854e-12 * m)
    return const * (y / ((x ** 2 + y ** 2) ** 1.5))

我很感激你能帮我解决这个问题！我可能只需要一双眼睛。你知道吗