太阳系8行星模拟：多列表可读存储及输出异常求助

问题核心 #

你的代码存在三个关键问题导致输出数据量不足，且不符合t, x, y, z, Vx, Vy, Vz的格式要求：

1. 固定大小的pos/vel数组不匹配变步长模拟逻辑
2. 保存条件的间隔设置过大，且索引引用错误
3. 未记录时间数据，且未按要求合并位置与速度列

错误点拆解与修正 #

1. 动态存储模拟数据 #

原代码用int(t_end)作为pos/vel数组的步数维度，这是错误的——因为采用变步长（每次取min_dt），实际模拟步数不等于t_end的整数。改用列表动态存储每一步的时间、位置、速度：

# 替换原固定数组为动态列表
pos_history = []
vel_history = []
time_history = []

# 初始化第一步
pos_history.append(initial_pos.copy())
vel_history.append(initial_vel.copy())
time_history.append(0.0)

2. 合理设置保存间隔与索引 #

原intervals=100000远大于模拟总步数（t_end≈126144秒，实际步数约为数千次），导致仅触发一次保存。建议改为时间间隔而非步数间隔（变步长下更合理），同时用当前模拟时间判断是否保存：

# 改为时间间隔（比如每0.0001年保存一次）
save_interval = 0.0001 * yr
next_save_time = save_interval

在循环中累计模拟时间，达到间隔时保存数据：

# 累计模拟时间
current_time += min_dt

# 达到保存间隔时存储数据
if current_time >= next_save_time:
    pos_output.append(pos_history[-1].copy())
    vel_output.append(vel_history[-1].copy())
    time_output.append(current_time)
    next_save_time += save_interval

3. 合并数据为要求格式 #

将每个行星的时间、位置、速度整合成一行，最终生成(保存次数×8)行的数组，每列对应t, x, y, z, Vx, Vy, Vz：

# 整合所有数据
output_data = []
for t, pos_batch, vel_batch in zip(time_output, pos_output, vel_output):
    for p, v in zip(pos_batch, vel_batch):
        output_data.append([t, p[0], p[1], p[2], v[0], v[1], v[2]])

# 转换为numpy数组并保存
output_data = np.array(output_data)
np.savetxt('planet_data.txt', output_data, fmt='%.6e', delimiter='\t', header='t\tx\ty\tz\tVx\tVy\tVz')

完整修正代码 #

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]
])

# Simulation parameters
t_end = 0.004 * yr  # Total time of integration
dt_constant = 0.1
save_interval = 0.0001 * yr  # 时间间隔，每0.0001年保存一次

# 动态存储每一步的模拟数据
pos_history = []
vel_history = []
current_time = 0.0

# 初始化第一步
pos_history.append(initial_pos.copy())
vel_history.append(initial_vel.copy())

# 用于保存输出的列表
pos_output = []
vel_output = []
time_output = []
next_save_time = save_interval

# Leapfrog Integration (2nd Order)
while current_time < t_end:
    # 获取当前位置和速度
    current_pos = pos_history[-1]
    current_vel = vel_history[-1]

    # 计算加速度
    r = np.linalg.norm(current_pos, axis=1)
    acc = -G * M_Sun / r[:, np.newaxis]**3 * current_pos

    # 计算变步长
    current_dt = dt_constant * np.sqrt(np.linalg.norm(current_pos, axis=1)**3 / (G * M_Sun))
    min_dt = np.min(current_dt)

    # Leapfrog步骤
    half_vel = current_vel + 0.5 * acc * min_dt
    new_pos = current_pos + half_vel * min_dt

    # 重新计算加速度
    r_new = np.linalg.norm(new_pos, axis=1)
    acc_new = -G * M_Sun / r_new[:, np.newaxis]**3 * new_pos
    new_vel = half_vel + 0.5 * acc_new * min_dt

    # 更新历史数据和时间
    pos_history.append(new_pos)
    vel_history.append(new_vel)
    current_time += min_dt

    # 达到保存间隔时存储数据
    if current_time >= next_save_time:
        pos_output.append(new_pos.copy())
        vel_output.append(new_vel.copy())
        time_output.append(current_time)
        next_save_time += save_interval

# 整合为要求的格式：t, x, y, z, Vx, Vy, Vz
output_data = []
for t, pos_batch, vel_batch in zip(time_output, pos_output, vel_output):
    for p, v in zip(pos_batch, vel_batch):
        output_data.append([t, p[0], p[1], p[2], v[0], v[1], v[2]])

# 转换为numpy数组并保存
output_data = np.array(output_data)
np.savetxt('planet_data.txt', output_data, fmt='%.6e', delimiter='\t', header='t\tx\ty\tz\tVx\tVy\tVz')

# 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] for pos in pos_history]
    y_particle = [pos[particle, 1] for pos in pos_history]
    z_particle = [pos[particle, 2] for pos in pos_history]
    ax.plot(x_particle, y_particle, z_particle, label=f'Planet {particle + 1} Orbit (m)')

ax.set_xlabel('X (m)')
ax.set_ylabel('Y (m)')
ax.set_zlabel('Z (m)')
ax.legend(loc='upper right', bbox_to_anchor=(1.1, 1.1))
ax.set_title('Orbits of Planets around Sun')
plt.show()

效果说明 #

• 现在会按照save_interval设置的时间间隔保存数据，总保存次数为t_end / save_interval（约40次），每个时间点保存8个行星的数据，最终输出文件有40×8=320行。
• 文件格式完全符合t, x, y, z, Vx, Vy, Vz的要求，可直接用Excel或Python读取分析。
• 动态存储解决了变步长下的数组大小不匹配问题，避免了内存浪费或越界错误。

内容的提问来源于stack exchange，提问作者bluebee09r