Python 欧拉方法在二体问题中的实施不起作用
Python Euler Method implementation in Two Body Problem not working
我目前正在尝试解决二体问题，然后我可以升级到更多行星，但它不起作用。它正在输出我不可能的位置。有谁知道是什么原因造成的？
这是我使用的代码：
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.147123353981485, -2812171.2728945166)
coordinate 2 (150013715707.77917, 2374319765.821534)
但它输出如下：
Sun
0.0, 0.0, 0.0
Earth
149600000000.0, 0.0, 0.0
注意：问题可能出在 for 循环或数组总和的舍入中，但我不确定。
我认为你问题的核心是你没有把它当做一个状态引擎。
假设"Bodies"是一个完全不可改变的值，它决定了系统在某个时间点的状态：
bodies_at_time_0 = ((sun, position, velocity, mass), (earth, position, velocity, mass))
你会得到下一个状态：
bodies_at_time_1 = apply_euler_method_for_one_tick( bodies_at_time_0 )
因此你的 "Bodies" 一次完全固定，下一次你计算一个全新的 "Bodies"。在计算中，您总是使用输入中的数据，这就是它们现在所在的位置。你正在做的是移动一些东西，然后根据错误的数字计算将其他东西移动到哪里（因为你已经移动了其他东西）。
一旦您确定您的函数使用输入状态，returns 输出状态，您就可以更轻松地分解它：
# advance all bodies one time interval, using their frozen state
def compute(bodies):
new_bodies = []
for body in bodies:
new_bodies.append(compute_one_body(body, bodies))
return new_bodies
# figure out where one body will move to, return its new state
def compute_one_body(start, bodies):
end = math stuff using the fixed state in bodies
return end
# MAIN
bodies = initial_state
for timepoint in whatever:
bodies = compute(bodies)
我喜欢将元组用于此类事情，以避免意外更改其他范围内的列表（因为列表是可变的）。
图片 - 1000 字
您的代码中的直接错误是
你计算的力方向错误，应该是rx = b[n].x-b[x].x等等，或者你需要稍后去掉负号。
您在单坐标中的计算会引起 copy-paste 错误，如
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 数组作为向量，将加速度计算与欧拉更新分开，在数值模拟期间删除了 non-sensical 舍入为整数值，删除了未使用的变量和字段，删除了中间变量对于直接更新加速度场的 force/acceleration 计算，将循环更改为使用时间来通知一年（或 10 年）已经过去（您的代码以 0.1 天的增量迭代 100 多年，这是故意的吗？） , ... 并将 Venus 添加到 body 中并添加代码以生成图像，结果见上文。
这种螺旋是欧拉方法的典型特征。您可以通过将欧拉更新更改为辛欧拉更新来轻松改进该模式，这意味着首先更新速度并使用新速度计算位置。在其他一切相同的情况下，这给出了图像
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()
我目前正在尝试解决二体问题，然后我可以升级到更多行星，但它不起作用。它正在输出我不可能的位置。有谁知道是什么原因造成的？
这是我使用的代码：
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.147123353981485, -2812171.2728945166)
coordinate 2 (150013715707.77917, 2374319765.821534)
但它输出如下：
Sun 0.0, 0.0, 0.0
Earth 149600000000.0, 0.0, 0.0
注意：问题可能出在 for 循环或数组总和的舍入中，但我不确定。
我认为你问题的核心是你没有把它当做一个状态引擎。
假设"Bodies"是一个完全不可改变的值，它决定了系统在某个时间点的状态：
bodies_at_time_0 = ((sun, position, velocity, mass), (earth, position, velocity, mass))
你会得到下一个状态：
bodies_at_time_1 = apply_euler_method_for_one_tick( bodies_at_time_0 )
因此你的 "Bodies" 一次完全固定，下一次你计算一个全新的 "Bodies"。在计算中，您总是使用输入中的数据，这就是它们现在所在的位置。你正在做的是移动一些东西，然后根据错误的数字计算将其他东西移动到哪里（因为你已经移动了其他东西）。
一旦您确定您的函数使用输入状态，returns 输出状态，您就可以更轻松地分解它：
# advance all bodies one time interval, using their frozen state
def compute(bodies):
new_bodies = []
for body in bodies:
new_bodies.append(compute_one_body(body, bodies))
return new_bodies
# figure out where one body will move to, return its new state
def compute_one_body(start, bodies):
end = math stuff using the fixed state in bodies
return end
# MAIN
bodies = initial_state
for timepoint in whatever:
bodies = compute(bodies)
我喜欢将元组用于此类事情，以避免意外更改其他范围内的列表（因为列表是可变的）。
图片 - 1000 字
您的代码中的直接错误是
你计算的力方向错误，应该是
rx = b[n].x-b[x].x等等，或者你需要稍后去掉负号。您在单坐标中的计算会引起 copy-paste 错误，如
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 数组作为向量，将加速度计算与欧拉更新分开，在数值模拟期间删除了 non-sensical 舍入为整数值，删除了未使用的变量和字段，删除了中间变量对于直接更新加速度场的 force/acceleration 计算，将循环更改为使用时间来通知一年（或 10 年）已经过去（您的代码以 0.1 天的增量迭代 100 多年，这是故意的吗？） , ... 并将 Venus 添加到 body 中并添加代码以生成图像，结果见上文。
这种螺旋是欧拉方法的典型特征。您可以通过将欧拉更新更改为辛欧拉更新来轻松改进该模式，这意味着首先更新速度并使用新速度计算位置。在其他一切相同的情况下，这给出了图像
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()