如何優化python中的時間步進算法？

Question

我有一段代碼，它爲所謂Lorenz95 model（由Ed Lorenz於1995年發明）增加了一個時間步長。它通常作爲40-變量模型實現，並顯示混沌行爲。我已經編寫了時間步的算法如下：

class Lorenz: 
    '''Lorenz-95 equation''' 

    global F, dt, SIZE 
    F = 8 
    dt = 0.01 
    SIZE = 40 

    def __init__(self): 
     self.x = [random.random() for i in range(SIZE)] 

    def euler(self): 
     '''Euler time stepping''' 
     newvals = [0]*SIZE 
     for i in range(SIZE-1): 
      newvals[i] = self.x[i] + dt * (self.x[i-1] * (self.x[i+1] - self.x[i-2]) - self.x[i] + F) 
     newvals[SIZE-1] = self.x[SIZE-1] + dt * (self.x[SIZE-2] * (self.x[0] - self.x[SIZE-3]) - self.x[SIZE-1] + F) 
     self.x = newvals 

此功能歐拉不慢，但不幸的是，我的代碼需要做出一個非常大量的調用它。有沒有一種方法可以編寫時間步進以使其運行更快？

非常感謝。

Answer 1

至少有兩種可能的優化：以更智能的方式工作（算法改進）並加快工作速度。

• 在算法側，您使用的Euler method，這是一階方法（因此全局誤差正比於步長大小），並且具有一個短小穩定區域。也就是說，效率不高。

• 另一方面，如果您使用標準的CPython實現，這種代碼將會非常緩慢。爲了解決這個問題，你可以試着在PyPy下運行它。它的即時編譯器可以使數字代碼運行速度提高100倍。您也可以編寫自定義的C或Cython擴展。

但有一個更好的辦法。求解常微分方程組是非常普遍的，因此scipy是Python中的核心科學庫之一，它包裝了快速的，經過測試的Fortran庫來解決它們。通過使用scipy，您可以同時獲得算法改進（因爲集成商將擁有更高的訂單）和快速實施。

解決洛倫茨95模型的一組的擾動的初始條件是這樣的：

import numpy as np 


def lorenz95(x, t): 
    return np.roll(x, 1) * (np.roll(x, -1) - np.roll(x, 2)) - x + F 

if __name__ == '__main__': 
    import matplotlib.pyplot as plt 
    from scipy.integrate import odeint 
    SIZE = 40 
    F = 8 
    t = np.linspace(0, 10, 1001) 
    x0 = np.random.random(SIZE) 
    for perturbation in 0.1 * np.random.randn(5): 
     x0i = x0.copy() 
     x0i[0] += perturbation 
     x = odeint(lorenz95, x0i, t) 
     plt.plot(t, x[:, 0]) 
    plt.show() 

和輸出（設定np.random.seed(7)，你可以是不同的）是很好的混亂。初始條件下的小擾動（僅在他的座標之一中）產生了非常不同的解決方案：

但是，它真的比歐拉時間步進快嗎？對於dt = 0.01它看起來幾乎快了三倍，但解決方案除了在最開始時不匹配。

如果dt減小，由歐拉方法提供的解決方案變得越來越類似於odeint解決方案，但它需要更長的時間。請注意較小的dt，後來的歐拉解決方案是如何鬆動odeint解決方案的。最精確的歐拉解決方案花費了600倍的時間來計算t = 6時的解決方案，而不是t = 10時的解決方案。請參閱完整腳本here。

最後，這個系統是如此不穩定，我甚至不認爲odeint解決方案在所有繪圖時間都是準確的。