以可讀格式將多個列表保存到一個文件

我現在花了很多時間處理你的代碼，所以我希望你能感謝我的輸入。也許我也誤解了一切，請告訴我。

編輯：在評論中提供了信息后，我更改了原來的答案，以更適合您的問題。

回答您的問題

因此，首先回答你的問題：我認為主要問題是你混淆了時間步長和時間：t是你當前時間步長的索引。然而，通過提示t < int(t_end)，您假設每個時間步長都有一秒的長度（由于您的自適應時間步長，它顯然沒有）。因此，您需要將while條件更改為小于t_end的實際時間（而不是索引）。我將引入一個時間數組來跟蹤時間，這意味著我們可以將其更改為：

while time[-1] <= t_end:

注意：但請記住，這將在時間過去后邁出一步t_end。

由于自適應時間步長，我們不知道我們最終會有多少時間步長。這意味著我們不能像您（以及我在之前的回答中）那樣初始化pos和vel arrays，而是需要將每個新結果添加到arrray中。這意味著我們只從兩個arrays的形狀(8, 1, 3)（包含初始條件）開始。有趣的是，你使用了一個自適應的時間步長和一個像蛙跳一樣的辛積分器。正因為如此，你的時間步長并不總是一樣的。通過初始化時間數組來跟蹤當前時間步長當前所處的時間可能是一種很好的做法，該時間數組將是1d，當然從0開始。

# Arrays to store pos and vel
pos = np.zeros((num_particles, 1, 3)) # shape: (#particles, timesteps taken + 1, #coordinates)
vel = np.zeros((num_particles, 1, 3)) # shape: (#particles, timesteps taken + 1, #coordinates)

# Initial conditions
time = np.zeros(1) # time array
pos[:, 0] = initial_pos
vel[:, 0] = initial_vel

然后在循環內部為每個時間步長添加一個相同形狀的數組。也意味著我們不再需要變量t。相反，軸1上的索引-1可以訪問先前時間步長的數量。我將為每個時間步長引入當前位置、速度和時間的新內部變量：

_pos_t = pos[:, -1] + half_vel * min_dt # shape: (#particles, #coordinates)
_vel_t = half_vel + 0.5 * acc * min_dt # shape: (#particles, #coordinates)
_time_t = time[-1] + min_dt # shape: timesteps taken + 1

并在從（8，3）變為（8，1，3）后將它們連接到完整的數組：

pos = np.concatenate((pos, _pos_t[:, np.newaxis]), axis=1)
vel = np.concatenate((vel, _vel_t[:, np.newaxis]), axis=1)
time = np.append(time, _time_t)

每隔一段時間存儲到新列表

現在我明白你在if counter % intervals == 0:后面的進球了。但你又一次混淆了時間步和時間。counter只是第二個索引，在我看來這是很難理解的。但是intervals你是如何描述它的，是一個時間間隔，所以以單位秒為單位。然而，用實際時間替換counter也是不正確的，因為我們不知道時間是否真的達到了間隔的倍數。我認為一個簡單的解決方案是引入一個我們已經達到的新變量'next_interval_time(you can also changeintervalsif you prefer) that keeps track which multiple ofintervals。如果時間經過一個間隔，請將當前位置和速度附加到輸出列表中，并跟蹤時間：

next_interval_time = intervals
while ... :
    # Check if the elapsed time has surpassed the next interval
    if _time_t >= next_interval_time:
        pos_output.append(_pos_t.copy())
        vel_output.append(_vel_t.copy())
        t_output.append(_time_t.copy())
        next_interval_time += intervals

循環結束后，我會將其轉換為numpy數組并轉置以獲得相同的軸順序：

pos_output = np.array(pos_output).transpose((1, 0, 2)) # shape: (#objects, #stored timesteps, #coordinates)
vel_output = np.array(vel_output).transpose((1, 0, 2)) # shape: (#objects, #stored timesteps, #coordinates)

正在將arrays存儲到文件

最好將numpy arrays存儲為txt（精度較低，效率較低），但最好使用numpy存儲arrays的方法：

np.save('array.npy', array)

可以加載的

np.load('array.npy')

最后得到完整的數組。列表也是如此。在那里，您可以使用Python泡菜方法。但是，您將所有列表轉換為arrays。

我建議的代碼

我解決了您的問題，并實施了我提出的所有更改。我改變了t_end =165*yr，這是一個海王星年，所以應該給你看粒子8的一個軌道。我還改為intervals=1000000，因為我發現時間步長對于間隔來說太大了（這意味著我們存儲了每個時間步長）。下面是完整的代碼：

# -- setup as defined by you --
# Steps
t_end = 165*yr # Total time of integration
dt_constant = 0.1
intervals = 1000000 # seconds after which to store pos & vel
next_interval_time = intervals

# Arrays to store pos and vel
pos = np.zeros((num_particles, 1, 3)) # shape: (#particles, timesteps taken, #coordinates)
vel = np.zeros((num_particles, 1, 3)) # shape: (#particles, timesteps taken, #coordinates)

# Initial conditions
time = np.zeros(1) # time array
pos[:, 0] = initial_pos
vel[:, 0] = initial_vel
pos_output = []
vel_output = []
t_output = []

while time[-1] <= t_end: # NOTE: We end up doing one more timestep after t_end
    r = np.linalg.norm(pos[:, -1], axis=1)
    acc = -G * M_Sun / r[:, np.newaxis]**3 * pos[:, -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[:, -1], axis=1)**3 / (G * M_Sun))
    min_dt = np.min(current_dt)  # Use the minimum time step for all particles

    # leap frog integration
    # calculate velocity at half timestep
    half_vel = vel[:, -1] + 0.5 * acc * min_dt
    # current position only used in this timestep
    _pos_t = pos[:, -1] + half_vel * min_dt # shape: (#particles, #coordinates)

    # Recalculate acceleration with the new position
    r = np.linalg.norm(_pos_t, axis=1)
    # Acceleration at timestep
    acc = -G * M_Sun / r[:, np.newaxis]**3 * _pos_t #np.newaxis for broadcasting with pos[:, i-1]
    # current velocity only used in this timestep
    _vel_t = half_vel + 0.5 * acc * min_dt # shape: (#particles, #coordinates)
    
    # time at timestep t
    _time_t = time[-1] + min_dt

   # Check if the elapsed time has surpassed the next interval
    if _time_t >= next_interval_time:
        pos_output.append(_pos_t.copy())
        vel_output.append(_vel_t.copy())
        t_output.append(_time_t.copy())
        next_interval_time += intervals

    # add axis at position 1 to allow concatenation
    pos = np.concatenate((pos, _pos_t[:, np.newaxis]), axis=1)
    vel = np.concatenate((vel, _vel_t[:, np.newaxis]), axis=1)
    time = np.append(time, _time_t)

    # show current status by printing timestep number (-1 because initial conditions)
    print(f'timestep: {time.size -1} [progress: {_time_t/t_end*100:.3f}%]')

pos_output = np.array(pos_output).transpose((1, 0, 2)) # shape: (#objects, #stored timesteps, #coordinates)
np.save('pos_output.npy', pos_output)
vel_output = np.array(vel_output).transpose((1, 0, 2)) # shape: (#objects, #stored timesteps, #coordinates)
np.save('vel_output.npy', vel_output)
t_output = np.array(t_output)
np.save('t_output.npy', t_output)

# Orbit Plot
# --- your plotting routine ---

Result

我最終得到了預期的情節：

輸出arrays具有形狀(8, 5203, 3)。結束時間為165*yr = 5.203E9，這意味著我們預計有5203個輸出(end_time/interval=5203.44）。我們可以稍后加載arrays，如前所述。

Improvements

也許我想到了一些需要改進的地方：

• 考慮使用無量綱量進行模擬
• 考慮使用極坐標來最小化數值誤差

Cheers!