两体问题中的Python Euler方法实现不起作用 [英] Python Euler Method implementation in Two Body Problem not working
问题描述
我目前正在尝试解决一个两体问题，然后我可以升级到更多的行星，但是它不起作用。它向我输出不可能的职位。有人知道是什么原因吗？
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坐标：
coordinate 1（41.14712335398 1485，-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)
但是它输出以下内容：
Sun
0.0，0.0，0.0
Sun 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更新更改为辛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 Euler方法实现不起作用的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
