图片-1000字
您代码中的直接错误是
您计算出的力向错误，应该是错误的方向，
rx = b[n].x-b[x].x
或者稍后需要删除减号。您在单个坐标中的计算会引起复制粘贴错误，如下所示
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
其中，在
y仍然配合你的使用
vx。中间舍入到整数值没有意义，只会稍微降低精度。
我更改了代码，以将numpy数组用作向量，将加速度计算与Euler更新分离，在数值模拟过程中删除了无意义的舍入为整数值，删除了未使用的变量和字段，删除了用于力/加速度计算的中间变量，直接更新加速度字段，更改循环以使用时间来通知已过去一年（或10年）的时间（您的代码以0.1天的增量迭代了100年，这是预期的吗？），…并将金星添加到主体和添加的代码以生成图像，结果请参见上文。
这种螺旋对于欧拉方法是典型的。您可以通过将Euler更新更改为辛欧拉更新来轻松地改进该模式，这意味着首先更新速度并使用新速度计算位置。在其他所有条件相同的情况下，这给出了图像
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()
