第一段引用上面的摘要
在进行大规模粒子模拟时，例如模拟行星轨道，存储所有时间步的位置和速度数据会消耗大量的内存。为了解决这个问题，我们需要找到一种方法，只保存每隔 N 个时间步的数据，从而在保证一定精度的前提下，显著降低内存占用。
原代码存在的问题在于，保存数据的时间点与数据更新的时间点不一致，导致保存的是未更新的数据。具体来说，在 counter % intervals == 0 条件成立时，t 的值已经比 counter 大 1，而 pos[:, t] 还没有被更新，因此保存的是零值。
以下是几种修改方案，可以解决这个问题：
方案一：保存 pos[:, counter]
最直接的解决方法是将 saved_pos.append(pos[:,t].copy()) 修改为 saved_pos.append(pos[:, counter].copy())。因为 counter 记录的是已经完成计算的时间步，所以 pos[:, counter] 包含的是正确的数据。
# Save the pos and vel here
if counter % intervals == 0:
saved_pos.append(pos[:, counter].copy())
saved_vel.append(vel[:, counter].copy())方案二：调整判断条件
另一种方法是调整判断条件，使得在 t 的值为 intervals 的整数倍时保存数据。由于 t 从 1 开始，因此判断条件应该改为 if t % intervals == 1:。
# Save the pos and vel here
if t % intervals == 1:
saved_pos.append(pos[:, t-1].copy())
saved_vel.append(vel[:, t-1].copy())需要注意的是，由于t是从1开始，所以这里保存的是pos[:, t-1]。
方案三：调整 t 的递增位置
还可以将 t += 1 移动到保存数据之后，这样可以保证在保存数据时，pos[:, t] 包含的是正确的数据。
# Save the pos and vel here
if counter % intervals == 0:
saved_pos.append(pos[:,t].copy())
saved_vel.append(vel[:,t].copy())
t += 1完整代码示例（基于方案一）
以下是修改后的完整代码示例，基于方案一：
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# Constants
M_Sun = 1.989e30 #Solar Mass
G = 6.67430e-11 # m^3 kg^(-1) s^(-2)
yr = 365 * 24 * 60 * 60 #1 year in seconds
# Number of particles
num_particles = 8
# Initial conditions for the particles (m and m/s)
initial_pos = np.array([
[57.9e9, 0, 0], #Mercury
[108.2e9, 0, 0], #Venus
[149.6e9, 0, 0], #Earth
[228e9, 0, 0], #Mars
[778.5e9, 0, 0], #Jupiter
[1432e9, 0, 0], #Saturn
[2867e9, 0, 0], #Uranus
[4515e9, 0, 0] #Neptune
])
initial_vel = np.array([
[0, 47400, 0],
[0, 35000, 0],
[0, 29800, 0],
[0, 24100, 0],
[0, 13100, 0],
[0, 9700, 0],
[0, 6800, 0],
[0, 5400, 0]
])
# Steps
t_end = 0.004 * yr #Total time of integration
dt_constant = 0.1
intervals = 10000 #Number of outputs of pos and vel to be saved
# Arrays to store pos and vel
pos = np.zeros((num_particles, int(t_end), 3))
vel = np.zeros((num_particles, int(t_end), 3))
# Leapfrog Integration (2nd Order)
pos[:, 0] = initial_pos
vel[:, 0] = initial_vel
saved_pos = []
saved_vel = []
t = 1
counter = 0
while t < int(t_end):
r = np.linalg.norm(pos[:, t - 1], axis=1)
acc = -G * M_Sun / r[:, np.newaxis]**3 * pos[:, t - 1] #np.newaxis for broadcasting with pos[:, i-1]
# Calculate the time step for the current particle
current_dt = dt_constant * np.sqrt(np.linalg.norm(pos[:, t - 1], axis=1)**3 / (G * M_Sun))
min_dt = np.min(current_dt) # Use the minimum time step for all particles
half_vel = vel[:, t - 1] + 0.5 * acc * min_dt
pos[:, t] = pos[:, t - 1] + half_vel * min_dt
# Recalculate acceleration with the new position
r = np.linalg.norm(pos[:, t], axis=1)
acc = -G * M_Sun / r[:, np.newaxis]**3 * pos[:, t] #np.newaxis for broadcasting with pos[:, i-1]
vel[:, t] = half_vel + 0.5 * acc * min_dt
t += 1
counter += 1
# Save the pos and vel here
if counter % intervals == 0:
saved_pos.append(pos[:,counter].copy())
saved_vel.append(vel[:,counter].copy())
saved_pos = np.array(saved_pos)
saved_vel = np.array(saved_vel)
# Orbit Plot
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(0, 0, 0, color='yellow', marker='o', s=50, label='Sun')
for particle in range(num_particles):
x_particle = pos[particle, :, 0]
y_particle = pos[particle, :, 1]
z_particle = pos[particle, :, 2]
ax.plot(x_particle, y_particle, z_particle, label=f'Particle {particle + 1} Orbit (km)')
ax.set_xlabel('X (km)')
ax.set_ylabel('Y (km)')
ax.set_zlabel('Z (km)')
ax.legend(loc='upper right', bbox_to_anchor=(1.1, 1.1))
ax.set_title('Orbits of Planets around Sun (km)')
plt.show()注意事项
- 时间步长选择： 时间步长 dt_constant 的选择会影响模拟的精度和稳定性。需要根据具体问题进行调整。
- 数据保存频率： intervals 的值越大，保存的数据越少，内存占用越低，但可能会损失部分精度。需要在内存占用和精度之间进行权衡。
- 内存管理： 即使只保存部分时间步的数据，大规模模拟仍然可能占用大量内存。可以考虑使用更高效的数据结构或算法，例如使用稀疏矩阵存储数据，或者使用流式计算，将数据写入文件而不是全部加载到内存中。
总结
通过修改保存数据的时间点或调整判断条件，可以有效地解决粒子模拟中内存占用过大的问题。选择合适的方案，并根据具体问题进行调整，可以提高模拟效率，并获得更准确的结果。同时，需要注意时间步长选择和数据保存频率，以在精度和内存占用之间取得平衡。
