提问人：Vaeu 提问时间：4/30/2022 更新时间：4/30/2022 访问量：521
使用 Velocity Verlet 模拟 3D 三体运动
Using Velocity Verlet to simulate 3-body motion in 3D
问：
正如标题所述，我想使用速度 Verlet 算法来模拟 3D 物体的轨迹。代码基于此。我的问题是输出的轨迹似乎是以直线形式传播的小振荡，而我应该得到一些混乱的东西，比如随机波浪线。我不确定问题出在我实现的速度 Verlet 算法上，还是在定义 3 个物体之间的引力或其他任何东西的数学上。
import numpy as np
import matplotlib.pyplot as plt
import numpy.linalg as lag
from mpl_toolkits.mplot3d import Axes3D
plt.style.use('dark_background')
# masses of stars
m1 = 10e30
m2 = 20e30
m3 = 30e30
# starting coordinates for suns
r1_initial = np.array([-10e10, 10e10, -11e10])
v1_initial = np.array([3e10, 0, 0])
r2_initial = np.array([0, 0, 0])
v2_initial = np.array([0, 3e10, 0])
r3_initial = np.array([10e10, 10e10, 12e10])
v3_initial = np.array([-3e10, 0, 0])
def force_grav_3bodies(r1,r2,r3,m1,m2,m3,G):
force_A = -G*m1*m2*(r1-r2)/lag.norm(r1-r2)**3 - G*m1*m3*(r1-r3)/lag.norm(r1-r3)**3
force_B = -G*m2*m1*(r2-r3)/lag.norm(r2-r3)**3 - G*m2*m3*(r2-r1)/lag.norm(r2-r1)**3
force_C = -G*m3*m1*(r3-r1)/lag.norm(r3-r1)**3 - G*m3*m2*(r3-r2)/lag.norm(r3-r2)**3
return force_A, force_B, force_C
def accel_3sun(r1,r2,r3,m1,m2,m3,G):
fg3bA, fg3bB, fg3bC = force_grav_3bodies(r1,r2,r3,m1,m2,m3,G)
accel_sunA, accel_sunB, accel_sunC = fg3bA/m1, fg3bB/m2, fg3bC/m3
return accel_sunA, accel_sunB, accel_sunC
dt = 0.1
Tmax = 2000
steps = int(Tmax/dt)
G = 6.67e-11
r1 = np.zeros([steps,3])
v1 = np.zeros([steps,3])
a1 = np.zeros([steps,3])
r2 = np.zeros([steps,3])
v2 = np.zeros([steps,3])
a2 = np.zeros([steps,3])
r3 = np.zeros([steps,3])
v3 = np.zeros([steps,3])
a3 = np.zeros([steps,3])
# initial conditions
r1[0], r2[0], r3[0] = r1_initial, r2_initial, r3_initial
v1[0], v2[0], v3[0] = v1_initial, v2_initial, v3_initial
a1[0], a2[0], a3[0] = accel_3sun(r1[0],r2[0],r3[0],m1,m2,m3,G)
# velocity verlet algorithm
for k in range(1,steps-1):
r1[k] = r1[k-1] + v1[k-1]*dt + 0.5*a1[k-1]*dt*dt
r2[k] = r2[k-1] + v2[k-1]*dt + 0.5*a2[k-1]*dt*dt
r3[k] = r3[k-1] + v3[k-1]*dt + 0.5*a3[k-1]*dt*dt
a1[k], a2[k], a3[k] = accel_3sun(r1[k],r2[k],r3[k],m1,m2,m3,G)
v1[k] = v1[k-1] + 0.5*(a1[k-1] + a1[k])*dt
v2[k] = v2[k-1] + 0.5*(a2[k-1] + a2[k])*dt
v3[k] = v3[k-1] + 0.5*(a3[k-1] + a3[k])*dt
fig = plt.figure(figsize=(8, 8))
ax = fig.gca(projection='3d')
plt.gca().patch.set_facecolor('black')
plt.plot([i[0] for i in r1], [j[1] for j in r1], [k[2] for k in r1] , '^', color='r', lw = 0.5, markersize = 0.1, alpha=0.5)
plt.plot([i[0] for i in r2], [j[1] for j in r2], [k[2] for k in r2] , '^', color='y', lw = 0.5, markersize = 0.1, alpha=0.5)
plt.plot([i[0] for i in r3], [j[1] for j in r3], [k[2] for k in r3] , '^', color='c', lw = 0.5, markersize = 0.1, alpha=0.5)
ax.scatter(r1[steps-2][0],r1[steps-2][1],r1[steps-2][2],color="r",marker="^",s=100)
ax.scatter(r2[steps-2][0],r2[steps-2][1],r2[steps-2][2],color="y",marker="^",s=100)
ax.scatter(r3[steps-2][0],r3[steps-2][1],r3[steps-2][2],color="c",marker="^",s=100)
ax.scatter(r1[0][0],r1[0][1],r1[0][2],color="r",marker="o",s=100,label="A")
ax.scatter(r2[0][0],r2[0][1],r2[0][2],color="y",marker="o",s=100,label="B")
ax.scatter(r3[0][0],r3[0][1],r3[0][2],color="c",marker="o",s=100,label="C")
plt.legend()
plt.axis('on')
ax.w_xaxis.set_pane_color((0.0, 0.0, 0.0, 1.0)), ax.w_yaxis.set_pane_color((0.0, 0.0, 0.0, 1.0)), ax.w_zaxis.set_pane_color((0.0, 0.0, 0.0, 1.0))
plt.show()
答：
0赞
Smoren
12/1/2023
#1
以下是 2D 中 Lennard Jones 势的 Verlet 算法实现示例。从 2 维到 3 维没有问题。
import numpy as np
def lennard_jones_potential(r, sigma: float, eps: float) -> np.ndarray:
"""
Lennard-Jones potential
"""
return 4 * eps * np.power(sigma, 12) / np.power(r, 12) - 4 * eps * np.power(sigma, 6) / np.power(r, 6)
def lennard_jones_force(r, sigma: float, eps: float) -> np.ndarray:
"""
Lennard-Jones force (Derivative of LJ potential)
"""
return 24 * eps * np.power(sigma, 6) / np.power(r, 7) - 48 * eps * np.power(sigma, 12) / np.power(r, 13)
def get_kinetic_energy(m, v):
"""
Get kinetic energy of particle
"""
return m * (v ** 2) / 2
def get_verlet_next_x(x, v, f, m, dt):
"""
Get next position of particle
"""
a = f / m
return x + v * dt + a * (dt ** 2) / 2
def get_verlet_next_f(r, rx, sigma: float, eps: float):
"""
Get next LJ interaction force of particle
"""
return np.sum(lennard_jones_force(r, sigma, eps) * rx / r)
def get_verlet_next_u(rs, sigma: float, eps: float):
"""
Get next LJ interaction potential of particle
"""
return np.sum(lennard_jones_potential(rs, sigma, eps)) / 2
def get_verlet_next_v(v, f, f_next, m, dt):
"""
Get next velocity of particle using Verlet algorithm
"""
a, a_next = f / m, f_next / m
return v + ((a + a_next) * dt) / 2
def get_verlet_next(p: np.ndarray, rx: np.ndarray, ry: np.ndarray, dt: float, sigma: float, eps: float):
"""
Get next state values of particle using Verlet algorithm
(force, velocity, potential, kinetic energy)
"""
r = np.sqrt(np.power(rx, 2) + np.power(ry, 2))
fx_next = get_verlet_next_f(r, rx, sigma, eps)
fy_next = get_verlet_next_f(r, ry, sigma, eps)
vx_next = get_verlet_next_v(p[P_VX], p[P_FX], fx_next, p[P_M], dt)
vy_next = get_verlet_next_v(p[P_VY], p[P_FY], fy_next, p[P_M], dt)
v_next = np.sqrt(np.power(vx_next, 2) + np.power(vy_next, 2))
u_next = get_verlet_next_u(r, sigma, eps)
ek_next = get_kinetic_energy(p[P_M], v_next)
return fx_next, fy_next, vx_next, vy_next, v_next, u_next, ek_next
# Indexes of particle properties
P_ID, P_M, P_X, P_Y, P_VX, P_VY, P_FX, P_FY = 0, 1, 2, 3, 4, 5, 6, 7
P_X_NEXT, P_Y_NEXT, P_VX_NEXT, P_VY_NEXT, P_FX_NEXT, P_FY_NEXT = 8, 9, 10, 11, 12, 13
# Indexes of result properties
R_ID, R_M, R_X, R_Y, R_VX, R_VY, R_FX, R_FY, R_U, R_EK, R_ET, R_T = 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11
def verlet_algo(particles: np.ndarray, simulation_time: float, dt: float, sigma: float, eps: float) -> np.ndarray:
"""
Verlet algorithm for N particles in 2D space
"""
time = np.arange(0, simulation_time + dt, dt)
result = np.zeros(shape=(particles.shape[0], time.shape[0], 12))
for i in range(particles.shape[0]):
result[i, :, R_ID] = np.repeat(particles[i, P_ID], result.shape[1])
result[i, :, R_M] = np.repeat(particles[i, P_M], result.shape[1])
result[:, 0, R_X] = particles[:, P_X]
result[:, 0, R_Y] = particles[:, P_Y]
result[:, 0, R_VX] = particles[:, P_VX]
result[:, 0, R_VY] = particles[:, P_VY]
result[:, 0, R_FX] = particles[:, P_FX]
result[:, 0, R_FY] = particles[:, P_FY]
result[:, 0, R_U] = 0
result[:, 0, R_EK] = 0
result[:, 0, R_ET] = 0
result[:, :, R_T] = time
# Calculate the initial forces and energies
for i in range(particles.shape[0]):
p = particles[i]
p_id = p[P_ID].astype(np.int64)
neighbours = particles[particles[:, P_ID] != p[P_ID]]
rx = neighbours[:, P_X] - p[P_X]
ry = neighbours[:, P_Y] - p[P_Y]
fx, fy, vx, vy, v, u, ek = get_verlet_next(p, rx, ry, dt, sigma, eps)
particles[i, P_FX] = fx
particles[i, P_FY] = fy
result[p_id, 0, R_FX] = fx
result[p_id, 0, R_FY] = fy
result[p_id, 0, R_U] = u
result[p_id, 0, R_EK] = ek
result[p_id, 0, R_ET] = ek + u
for curr in range(len(time) - 1):
next = curr + 1
# Get following coordinates of all particles:
particles[:, P_X_NEXT] = get_verlet_next_x(particles[:, P_X], particles[:, P_VX], particles[:, P_FX],
particles[:, P_M], dt)
particles[:, P_Y_NEXT] = get_verlet_next_x(particles[:, P_Y], particles[:, P_VY], particles[:, P_FY],
particles[:, P_M], dt)
for i in range(particles.shape[0]):
p = particles[i]
p_id = p[P_ID].astype(np.int64)
neighbours = particles[particles[:, P_ID] != p[P_ID]]
# Get distances between current particle and all other particles:
rx = neighbours[:, P_X_NEXT] - p[P_X_NEXT]
ry = neighbours[:, P_Y_NEXT] - p[P_Y_NEXT]
# Get next force, velocity and energies for current particle using verlet algorithm:
fx_next, fy_next, vx_next, vy_next, v_next, u_next, ek_next = get_verlet_next(p, rx, ry, dt, sigma, eps)
# Save next values of current particle:
p[P_VX_NEXT] = vx_next
p[P_VY_NEXT] = vy_next
p[P_FX_NEXT] = fx_next
p[P_FY_NEXT] = fy_next
# Save next values of current particle to result array:
result[p_id, next, R_X] = p[P_X_NEXT]
result[p_id, next, R_Y] = p[P_Y_NEXT]
result[p_id, next, R_VX] = vx_next
result[p_id, next, R_VY] = vy_next
result[p_id, next, R_FX] = fx_next
result[p_id, next, R_FY] = fy_next
result[p_id, next, R_U] = u_next
result[p_id, next, R_EK] = ek_next
result[p_id, next, R_ET] = ek_next + u_next
# Apply next values to all particles:
particles[:, P_X] = particles[:, P_X_NEXT]
particles[:, P_Y] = particles[:, P_Y_NEXT]
particles[:, P_VX] = particles[:, P_VX_NEXT]
particles[:, P_VY] = particles[:, P_VY_NEXT]
particles[:, P_FX] = particles[:, P_FX_NEXT]
particles[:, P_FY] = particles[:, P_FY_NEXT]
return result
particles = np.array([
#id,m, x, y, vx,vy,fx,fy,ẋ, ẏ, vẋ,vẏ,fẋ,fẏ (ẋ == x_next, ...)
[0, 2, 11, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[1, 1, 10, 13, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[2, 3, 12, 13, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
], dtype=np.float64)
result = verlet_algo(particles, simulation_time=50, dt=0.002, sigma=1, eps=0.7, animate=True)
能量图：
您可以在我自己的 GitHub 存储库上看到更多详细信息。
评论
3赞
Adrian Mole
12/1/2023
对特定产品/资源的过度宣传可能会被社区视为垃圾邮件。看看帮助中心，特别是用户应该有什么样的行为？最后一节：避免公开的自我推销。您可能还对如何不成为垃圾邮件发送者和如何在 Stack Overflow 上做广告感兴趣？。
0赞
Lutz Lehmann
12/2/2023
#2
基本示例轨道或其变体是半径为 1 的圆，中心质量位于中心。如果选择时间尺度使轨道速度为 ，则方程中的常数为 。让我们将其分发为 ，因此也 .1mu = GMq'' = -mu*q/|q|^3mu=1M=1G=1
现在改成链接柱的刻度，半径约为 10，总质量假定为 、 、 时间刻度的中心质量。在新的位置和质量中，方程显示为M=60G=9.8T
T^2*(q/10)'' = -(G/9.8)*(M/60)*(q/10)/|q/10|^3
q'' = -G*M*q/|q|^3 * 50/(3*9.8)/T^2
固定几乎消除了最后一个因素，因此新配置的轨道速度约为 ，这是链接页面上速度的比例。那里的轨迹不是圆形的，身体“相互坠落”。T=1.3R/T=10/1.3=7.7
您希望将系统缩放到太阳系的尺寸，并使用自然重力常数。所以半径变为 或 .总质量为30个太阳质量，时间单位变为一些未指定的。与上面相同的舞蹈将重力公式更改为R=100e9 m = 100e6 km2/3 AUM=60e30 kgT
T^2*(q/1e11)'' = -(G/6.67e-11)*(M/6e31)*(q/1e11)/|q/1e11|^3
q'' = -G*M*q/|q|^3 * 1e13/(6*6.67)/T^2
要将最后一个因数设置为 1，需要 大约，它给出的轨道速度为 或大约 。因此，要获得与网页上类似的图，您需要在这个范围内的速度，例如，通过将那里的速度缩放一个系数，使 的速度分量变为 。你的初始速度要高得多，所以物体会非常迅速地远离彼此，而不会对彼此产生重大影响。T=5e5R/T=10e10/5e5=2e5200 km/s2e5/8=25000375000
上述速度给出了合理但不是非常有趣的结果。将速度分量减小到尺寸可在图片中提供更多交互3e4
评论