Пытаетесь сделать простой двухмерный симулятор орбиты с двумя телами, но все продолжает вертеться в космос?

Question

Я делаю небольшой двухмерный симулятор орбиты с двумя телами. С исходными позициями (0,0) и (1,0) я бы ожидал, что два тела будут перемещаться взад и вперед между двумя позициями, но вместо этого они взлетают в космос.

import numpy

import pandas

import matplotlib.pyplot as plt

import math

# variables right here
m1 = 1
m2 = 1
x1 = 0
y1 = 0
x2 = 1
y2 = 0
vx1 = 0
vy1 = 0
vx2 = 0
vy2 = 0
t = 0
dt = 0.01
G = 10
ax1 = 0
ay1 = 0
ax2 = 0
ay2 = 0
r2 = 0
#a = Gm/r^2
#a = Gm/((x2-x1)^2 + (y2-y1)^2)
#tan theta = (y2-y1)/(x2-x1)
#theta = arctan(dy/dx)
#ax = a cos arctan dy/dx
#ay = a sin arctan dy/dx

for t in range(0, 100):
    # just the r^2 calculation so I don't have to add it a bunch of times
    r2 = (x2 - x1)**2 + (y2 - y1)**2
    # quick calculation of theta, just for debugging
    theta = math.atan((y2 - y1) / (x2 - x1))
    # calculate the accelerations
    if x1 > x2:
        ax1 = -1 * math.cos(numpy.arctan2(y2 - y1, x2 - x1)) * G * m2 / r2
    else:
        ax1 = 1 * math.cos(numpy.arctan2(y2 - y1, x2 - x1)) * G * m2 / r2
    if  y1 > y2:
        ay1 = -1 * math.sin(numpy.arctan2(y2 - y1, x2 - x1)) * G * m2 / r2
    else:
        ay1 = 1 * math.sin(numpy.arctan2(y2 - y1, x2 - x1)) * G * m2 / r2
    if x1 > x2:
        ax2 = 1 * math.cos(numpy.arctan2(y2 - y1, x2 - x1)) * G * m1 / r2
    else:
        ax2 = -1 * math.cos(numpy.arctan2(y2 - y1, x2 - x1)) * G * m1 / r2
    if y1 > y2:
        ay2 = 1 * math.sin(numpy.arctan2(y2 - y1, x2 - x1)) * G * m1 / r2
    else:
        ay2 = -1 * math.sin(numpy.arctan2(y2 - y1, x2 - x1)) * G * m1 / r2
    # now change the velocities
    vx1 = vx1 + (ax1 * dt)
    vy1 = vy1 + (ay1 * dt)
    vx2 = vx2 + (ax2 * dt)
    vy2 = vy2 + (ay2 * dt)
    # now the positions
    x1 = x1 + (vx1 * dt)
    y1 = y1 + (vy1 * dt)
    x2 = x2 + (vx2 * dt)
    y2 = y2 + (vy2 * dt)
    # change t
    t = t + dt
    print x1, x2    

Я бы ожидал, что x1, x2 почти «подпрыгнут» между 0 и 1, но вместо этого я получаю ускорение x1 до бесконечности и ускорение x2 до отрицательной бесконечности.