问题描述
我目前正在尝试解决两体问题，然后我可以升级到更多行星，但它不起作用.它正在输出我不可能的位置.有谁知道是什么原因造成的?
I am currently trying to get a two body problem to work, that I can then upgrade to more planets, but it is not working. It is outputting me impossible positions. Does anyone know what is causing that?
这是我使用的代码:
day = 60*60*24
# Constants
G = 6.67408e-11
dt = 0.1*day
au = 1.496e11
t = 0
class CelBody:
def __init__(self, id, name, x0, y0, z0, vx0, vy0, vz0, mass, vector, ax0, ay0, az0, totalforcex, totalforcey, totalforcez):
self.ax0 = ax0
self.ay0 = ay0
self.az0 = az0
self.ax = self.ax0
self.ay = self.ay0
self.az = self.az0
# Constants of nature
# Universal constant of gravitation
self.G = 6.67408e-11
# Name of the body (string)
self.id = id
self.name = name
# Initial position of the body (au)
self.x0 = x0
self.y0 = y0
self.z0 = z0
# Position (au). Set to initial value.
self.x = self.x0
self.y = self.y0
self.z = self.z0
# Initial velocity of the body (au/s)
self.vx0 = vx0
self.vy0 = vy0
self.vz0 = vz0
# Velocity (au/s). Set to initial value.
self.vx = self.vx0
self.vy = self.vy0
self.vz = self.vz0
# Mass of the body (kg)
self.M = mass
# Short name
self.vector = vector
self.totalforcex = totalforcex
self.totalforcey = totalforcey
self.totalforcez = totalforcez
# All Celestial Bodies
forcex = 0
forcey = 0
forcez = 0
Bodies = [
CelBody(0, 'Sun', 1, 1, 1, 0, 0, 0, 1.989e30, 'sun', 0, 0, 0, 0, 0, 0),
CelBody(1, 'Mercury', 1*au, 1, 1, 0, 29780, 0, 3.3e23, 'earth', 0, 0, 0, 0, 0, 0),
]
leftover_bin = []
templistx = []
templisty = []
templistz = []
for v in range(365242):
for n in range(len(Bodies)):
#Need to initialize the bodies
planetinit = Bodies[n]
for x in range(len(Bodies)):
# Temporary lists and initial conditions
planet = Bodies[x]
if (planet == planetinit):
pass
else:
rx = Bodies[x].x - Bodies[n].x
ry = Bodies[x].y - Bodies[n].y
rz = Bodies[x].z - Bodies[n].z
r3 = (rx**2+ry**2+rz**2)**1.5
gravconst = G*Bodies[n].M*Bodies[x].M
fx = -gravconst*rx/r3
fy = -gravconst*ry/r3
fz = -gravconst*rz/r3
# Make a temporary list of the total forces and then add them to get the resulting force
templistx.append(fx)
templisty.append(fy)
templistz.append(fz)
forcex = sum(templistx)
forcey = sum(templisty)
forcez = sum(templistz)
templistx.clear()
templisty.clear()
templistz.clear()
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
Bodies[n].x = x
Bodies[n].y = y
Bodies[n].z = z
vx = int(Bodies[n].vx) + forcex/int(Bodies[n].M)*dt
vy = int(Bodies[n].vy) + forcey/int(Bodies[n].M)*dt
vz = int(Bodies[n].vz) + forcez/int(Bodies[n].M)*dt
Bodies[n].vx = vx
Bodies[n].vy = vy
Bodies[n].vz = vz
t += dt
print(Bodies[0].name)
print(Bodies[0].x)
print(Bodies[0].y)
print(Bodies[0].z)
print(Bodies[1].name)
print(Bodies[1].x)
print(Bodies[1].y)
print(Bodies[1].z)
它应该输出类似于这里的坐标，然后还有一个 z 坐标:坐标1 (41.147123353981485, -2812171.2728945166)坐标2 (150013715707.77917, 2374319765.821534)
It should output something like the coordinates here, but then also a z coordinate:coordinate 1 (41.147123353981485, -2812171.2728945166)coordinate 2 (150013715707.77917, 2374319765.821534)
但它输出以下内容:
太阳0.0, 0.0, 0.0
地球149600000000.0, 0.0, 0.0
Earth 149600000000.0, 0.0, 0.0
注意:问题可能出在 for 循环或数组总和的舍入中，但我不确定.
Note: The problem is probably in the for loops or in the rounding of the sum of the arrays but I am not sure.
推荐答案
图片 - 1000字
代码中的直接错误是
你计算的力方向错误，应该是
rx = b[n].xb[x].x等，或者你需要去掉一些行的减号稍后.
You compute the force in the wrong direction, it should be
rx = b[n].x-b[x].xetc. 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
在 y 坐标中的哪个位置仍然使用 vx.中间舍入到整数值没有意义，它只会在一定程度上降低准确性.
where in the y coordinate you still use vx. The intermediate rounding to integer values makes no sense, it only reduces the accuracy somewhat.
我将您的代码更改为使用 numpy 数组作为向量，将加速度计算与 Euler 更新分开，在数值模拟期间删除了无意义的整数值舍入，删除了未使用的变量和字段，删除了力的中间变量/加速计算以直接更新加速度字段，更改循环以使用时间来注意一年(或 10)已过去(您的代码以 0.1 天的增量迭代 100 多年，是有意的吗?)，...并添加维纳斯到身体并添加代码以生成图像，结果见上.
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.
这种螺旋是欧拉方法的典型特征.您可以通过将 Euler 更新更改为辛 Euler 更新来轻松改进该模式，这意味着首先更新速度并使用新速度计算位置.这给出了，在其他一切都相同的情况下，图像
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()
这篇关于两体问题中的 Python 欧拉方法实现不起作用的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！