picture - 1000 words
The direct errors in your code are
You compute the force in the wrong direction, it should be rx = b[n].x-b[x].x etc. or you need to remove the minus sign some lines later.
Your computation in single coordinates invites copy-paste errors as in
x = int(Bodies[n].x) + int(Bodies[n].vx) * dt
y = int(Bodies[n].y) + int(Bodies[n].vx) * dt
z = int(Bodies[n].z) + int(Bodies[n].vz) * dt
where in the y coordinate you still use vx. The intermediate rounding to integer values makes no sense, it only reduces the accuracy somewhat.
I changed your code to use numpy arrays as vectors, separated the acceleration computation from the Euler updates, removed the non-sensical rounding to integer values during the numerical simulation, removed unused variables and fields, removed intermediate variables for the force/acceleration computations to directly update the acceleration field, changed the loop to use the time to notice when a year (or 10) has passed (your code iterates over 100 years in 0.1day increments, was that intended?), ... and added Venus to the bodies and added code to produce images, result see above.
This spiraling is typical for the Euler method. You can easily improve that pattern by changing the Euler update to the symplectic Euler update, which means to update the velocity first and compute the position with the new velocity. This gives, with everything else the same, the image
day = 60*60*24
# Constants
G = 6.67408e-11
au = 1.496e11
class CelBody(object):
# Constants of nature
# Universal constant of gravitation
def __init__(self, id, name, x0, v0, mass, color, lw):
# Name of the body (string)
self.id = id
self.name = name
# Mass of the body (kg)
self.M = mass
# Initial position of the body (au)
self.x0 = np.asarray(x0, dtype=float)
# Position (au). Set to initial value.
self.x = self.x0.copy()
# Initial velocity of the body (au/s)
self.v0 = np.asarray(v0, dtype=float)
# Velocity (au/s). Set to initial value.
self.v = self.v0.copy()
self.a = np.zeros([3], dtype=float)
self.color = color
self.lw = lw
# All Celestial Bodies
t = 0
dt = 0.1*day
Bodies = [
CelBody(0, 'Sun', [0, 0, 0], [0, 0, 0], 1.989e30, 'yellow', 10),
CelBody(1, 'Earth', [-1*au, 0, 0], [0, 29783, 0], 5.9742e24, 'blue', 3),
CelBody(2, 'Venus', [0, 0.723 * au, 0], [ 35020, 0, 0], 4.8685e24, 'red', 2),
]
paths = [ [ b.x[:2].copy() ] for b in Bodies]
# loop over ten astronomical years
v = 0
while t < 10*365.242*day:
# compute forces/accelerations
for body in Bodies:
body.a *= 0
for other in Bodies:
# no force on itself
if (body == other): continue # jump to next loop
rx = body.x - other.x
r3 = sum(rx**2)**1.5
body.a += -G*other.M*rx/r3
for n, planet in enumerate(Bodies):
# use the symplectic Euler method for better conservation of the constants of motion
planet.v += planet.a*dt
planet.x += planet.v*dt
paths[n].append( planet.x[:2].copy() )
#print("%10s x:%53s v:%53s"%(planet.name,planet.x, planet.v))
if t > v:
print("t=%f"%t)
for b in Bodies: print("%10s %s"%(b.name,b.x))
v += 30.5*day
t += dt
plt.figure(figsize=(8,8))
for n, planet in enumerate(Bodies):
px, py=np.array(paths[n]).T;
plt.plot(px, py, color=planet.color, lw=planet.lw)
plt.show()
