我正在嘗試使用 Velocity-Verlet 積分器在勢內積分恒星的軌道。然而，我目前的演算法每個軌道大約需要 25 秒，我的目標是整合大約 1000 個軌道，這需要大約 7 個小時。因此我希望有人可以幫助我優化演算法：
Velocity-Verlet 積分器定義為（參見verlet_integration()實作）：
我想使用定義為的簡單對數勢 (SLP) 來整合軌道：
該值v_0始終等于，1并且該值q可以是0.7或0.9（請參閱SLP()）。
通過使用該電位，可以計算出該積分器所需的加速度（參見apply_forces()）
在此之后，我選擇x = 0作為所有軌道的起始值，以及能量E = 0。vx使用可以計算的起始值（請參閱calc_vx()）
為了獲得足夠準確的軌道，時間步長需要1E-4更小或更小。我需要稍后分析這些軌道，因為它們需要足夠長，因此我在 和 之間進行t=0積分t=200。
我需要計算（y，vy）的整個允許相空間的值。允許的相空間是該方法calc_vx()不會產生負數平方根的地方。因此，需要整合大量軌道。我希望能夠計算出至少 1000 個值。但是 10,000 肯定會更理想。雖然也許這要求太多了。如果您對性能改進有任何想法，請告訴我。使用另一種語言不是一種選擇，所以請不要建議。
可以在這里看到這些軌道的外觀示例：
運行代碼所需的一切都應該在下面找到，以及運行它的起始值。
更新：我已經實施了 mivkov 的建議，這將時間減少到 9 秒，快了 3 秒，但它仍然很慢。任何其他建議仍然歡迎
import numpy as np
def calc_vx(y, vy, q):
"""
Calculate starting value of x velocity
"""
vx2 = -np.log((y / q) ** 2) - vy ** 2
return np.sqrt(vx2)
def apply_forces(x, y, q):
"""
Apply forces to determine the accelerations
"""
Fx = -x / (y ** 2 / q ** 2 x ** 2)
Fy = -y / (q ** 2 * x ** 2 y ** 2)
return Fx, Fy
def verlet_integration(start, dt, steps, q):
# initialise an array and set the first value to the starting value
vals = np.zeros((steps, *np.shape(start)))
vals[0] = start
# run through all elements and apply the integrator to each value
for i in range(steps - 1):
x_vec, v_vec, a_vec = vals[i]
new_x_vec = x_vec dt * (v_vec 0.5 * a_vec * dt)
new_a_vec = apply_forces(*new_x_vec, q)
new_v_vec = v_vec 0.5 * (a_vec new_a_vec) * dt
vals[i 1] = new_x_vec, new_v_vec, new_a_vec
# I return vals.T so i can use the following to get arrays for the position, velocity and acceleration
# ((x, vx, ax), (y, vy, ay)) = verlet_integration_vec( ... )
return vals.T
def integration(y, vy, dt, t0, t1, q):
# calculate the starting values
vx = calc_vx(y, vy, q)
ax, ay = apply_forces(0, y, q)
start = [(0, y), (vx, vy), (ax, ay)]
steps = round((t1 - t0) / dt) # bereken het aantal benodigde stappen
e = verlet_integration(start, dt, steps, q) # integreer
return e
((x, vx, ax), (y, vy, ay)) = integration(0.1, 0.2, 1E-4, 0, 100, 0.7)
uj5u.com熱心網友回復：
我在這里看到兩件事可以提供幫助。
(1) 考慮使用庫中的 ODE 求解器，而不是您手動撰寫的 Verlet 方法。Runge-Kutta 方法及其變體（例如 Runge-Kutta-Fehlberg）廣泛適用，我會先嘗試。庫方法很可能會用 C 代碼撰寫，因此速度更快，而且有人已經解決了這些錯誤。這兩個方面都將對此問題有用。
(2) 如果您需要自己撰寫 ODE 求解器，在 Python 中加快實作的一個想法是對軌跡進行矢量化，以便您可以使用 Numpy 陣列操作。也就是說，創建一個代表所有 1000 或 10,000 條軌跡的陣列，然后一次將所有軌跡推進一步。我假設您可以以矩陣形式重寫運動方程。
順便說一句，我的建議是保留您已經制定的解決方案，因為它似乎正在作業，并使用想法（1）或（2）開始一個單獨的實作，并使用您的原始實作來驗證結果。
uj5u.com熱心網友回復：
首先要做的是使用分析器來分析代碼中的熱點。你可以用它cProfile.run來做到這一點。一份基本報告顯示，幾乎所有時間都花在了以下函式上：x_func和v_func（apply_forces呼叫了 200_000 次）：
ncalls tottime percall cumtime percall filename:lineno(function)
200000 0.496 0.000 0.496 0.000 apply_forces
199999 0.279 0.000 2.833 0.000 func_verlet
199999 0.713 0.000 0.713 0.000 x_func
199999 0.780 0.000 0.780 0.000 v_func
1 0.000 0.000 0.000 0.000 SLP
199999 0.222 0.000 0.718 0.000 a_func
1 0.446 0.446 3.279 3.279 verlet_integration
[...]
快速分析x_func并v_func顯示您正在微型陣列上呼叫 Numpy 函式。問題是Numpy 沒有針對處理這么小的陣列進行優化（Numpy 對輸入進行許多檢查，與微型陣列的計算時間相比，這些檢查是昂貴的）。
解決這個問題的主要方法是使用向量化的 Numpy 呼叫在更大的陣列上運行。問題是這個解決方案需要完全重新設計您的代碼，以便擺脫for i in range(steps - 1)問題根源所在的回圈。
另一種解決方案是使用 Numba 以避免 Numpy 開銷。這個解決方案要簡單得多（盡管如果目標是學習 Numpy，你的老師可能不會想到這一點）。這是一個例子：
import numpy as np
import numba as nb
@nb.njit
def SLP(x, y, v0, q):
return 0.5 * v0 ** 2 * np.log(x ** 2 (y / q) ** 2)
@nb.njit
def calc_vx(x, y, vy, q):
"""
Calculate starting value of x velocity
"""
vx2 = -2 * SLP(x, y, 1, q) - vy ** 2
return np.sqrt(vx2)
@nb.njit
def apply_forces(x, y, v0, q):
"""
Apply forces to determine the accelerations
"""
Fx = -(v0 ** 2 * x) / (y ** 2 / q ** 2 x ** 2)
Fy = -(v0 ** 2 * y) / (q ** 2 * x ** 2 y ** 2)
return np.array([Fx, Fy])
@nb.njit
def x_func(x_vec, v_vec, a_vec, dt):
return x_vec dt * (v_vec 0.5 * a_vec * dt)
@nb.njit
def v_func(v_vec, a_vec, new_a_vec, dt):
return v_vec 0.5 * (a_vec new_a_vec) * dt
@nb.njit
def a_func(x_vec, dt, q):
x, y = x_vec
return apply_forces(x, y, 1, q)
# The parameter is a signature of the function that provides the input type to
# Numba so it can eagerly compile the function and all the dependent functions.
# Please read the Numba documentation for more information about this.
@nb.njit('(float64[:], float64[:], float64[:], float64, float64)')
def func_verlet(x_vec, v_vec, a_vec, dt, q):
# calculate the new position, velocity and acceleration
new_x_vec = x_func(x_vec, v_vec, a_vec, dt)
new_a_vec = a_func(new_x_vec, dt, q)
new_v_vec = v_func(v_vec, a_vec, new_a_vec, dt)
out = np.empty((len(new_x_vec), 3))
out[:,0] = new_x_vec
out[:,1] = new_v_vec
out[:,2] = new_a_vec
return out
def verlet_integration(start, f, dt, steps, q):
# initialise an array and set the first value to the starting value
vals = np.zeros((steps, *np.shape(start)))
vals[0] = start
# run through all elements and apply the integrator to each value
for i in range(steps - 1):
vals[i 1] = f(*vals[i].T, dt, q)
# I return vals.T so i can use the following to get arrays for the position, velocity and acceleration
# ((x, y), (vx, vy), (ax, ay)) = verlet_integration_vec( ... )
return vals.T
def integration(y, vy, dt, t0, t1, q):
# calculate the starting values
x = 0
vx = calc_vx(x, y, vy, q)
ax, ay = apply_forces(x, y, 1, q)
start = [(x, vx, ax), (y, vy, ay)]
steps = round((t1 - t0) / dt) # calculate the number of necessary steps
return verlet_integration(start, func_verlet, dt, steps, q)
((x_, y_), (vx_, vy_), (ax_, ay_)) = integration(0.1, 0.2, 1E-4, 0, 20, 0.7)
請注意，某些功能已被修改，以避免使用 Numba 不支持的功能/運算子。例如，*不支持打開 Numpy 陣列或串列的一元運算子，但無論如何它通常效率很低。
生成的代碼快 5 倍。
由于 for 回圈的開銷（包括運算子和函式呼叫），分析器現在顯示該函式verlet_integration負責執行的主要部分。*由于 lambda，這部分不能輕易移植到 Numba。我認為如果你成功地重新設計這部分以避免 lambda 和展開，它可以加快兩倍。事實上，對包含 2 個專案的陣列進行操作是非常低效的，即使使用 Numba 也是如此。在標量上操作會使代碼的可讀性降低一些，但速度會更快（當然無論有沒有 Numba）。我想代碼可以再快幾倍。
更新：使用更新的代碼，Numba 可以提供更好的幫助，因為現在主要的性能瓶頸已得到修復。這是新的 Numba 版本：
import numpy as np
import numba as nb
@nb.njit
def calc_vx(y, vy, q):
vx2 = -np.log((y / q) ** 2) - vy ** 2
return np.sqrt(vx2)
@nb.njit
def apply_forces(x, y, q):
Fx = -x / (y ** 2 / q ** 2 x ** 2)
Fy = -y / (q ** 2 * x ** 2 y ** 2)
return np.array([Fx, Fy])
@nb.njit('(float64[:,:], float64, int_, float64)')
def verlet_integration(start, dt, steps, q):
vals = np.zeros((steps, 3, 2))
vals[0] = start
for i in range(steps - 1):
x_vec, v_vec, a_vec = vals[i]
new_x_vec = x_vec dt * (v_vec 0.5 * a_vec * dt)
x, y = new_x_vec
new_a_vec = apply_forces(x, y, q)
new_v_vec = v_vec 0.5 * (a_vec new_a_vec) * dt
vals[i 1, 0] = new_x_vec
vals[i 1, 1] = new_v_vec
vals[i 1, 2] = new_a_vec
return vals.T
def integration(y, vy, dt, t0, t1, q):
vx = calc_vx(y, vy, q)
ax, ay = apply_forces(0, y, q)
start = [(0, y), (vx, vy), (ax, ay)]
steps = round((t1 - t0) / dt)
e = verlet_integration(np.array(start), dt, steps, q)
return e
((x, vx, ax), (y, vy, ay)) = integration(0.1, 0.2, 1E-4, 0, 100, 0.7) # 9.7
這比問題的更新代碼快 36 倍。在我的機器上只需要 0.27 秒，而初始代碼需要 9.7 秒。
uj5u.com熱心網友回復：
在我的基準測驗中，行內apply_forces將 @Jér?meRichard 的解決方案加速了約 2 倍。
import numpy as np
import numba as nb #tested with numba 0.55.2
@nb.njit
def calc_vx(y, vy, q):
vx2 = -np.log((y / q) ** 2) - vy ** 2
return np.sqrt(vx2)
@nb.njit
def verlet_integration(start, dt, steps, q):
vals = np.zeros((steps, 3, 2))
vals[0] = start
for i in range(steps - 1):
x_vec, v_vec, a_vec = vals[i]
x, y = x_vec dt * (v_vec 0.5 * a_vec * dt)
ax, ay = -x / (y ** 2 / q ** 2 x ** 2), -y / (q ** 2 * x ** 2 y ** 2)
vx, vy = v_vec[0] 0.5 * (a_vec[0] ax) * dt, v_vec[1] 0.5 * (a_vec[1] ay) * dt
vals[i 1, 0] = x, y
vals[i 1, 1] = vx, vy
vals[i 1, 2] = ax, ay
return vals.T
@nb.njit
def integration(y, vy, dt, t0, t1, q):
vx = calc_vx(y, vy, q)
ax, ay = 0, -y / y**2
start = np.array([[0, y], [vx, vy], [ax, ay]])
steps = round((t1 - t0) / dt)
return verlet_integration(start, dt, steps, q)
((x, vx, ax), (y, vy, ay)) = integration(0.1, 0.2, 1E-4, 0, 100, 0.7)
轉載請註明出處，本文鏈接：https://www.uj5u.com/ruanti/485474.html