向量化微分方程組的正向歐拉法

Question

我對一階微分方程組的x（t）數值求解。該系統是：

DX/DT = Y
DY/DT = -x - 一個* Y（X^2 + Y^2 -1）

我已經實現向前歐拉方法來解決這個問題如下：

def forward_euler(): 
    h = 0.01 
    num_steps = 10000 

    x = np.zeros([num_steps + 1, 2]) # steps, number of solutions 
    y = np.zeros([num_steps + 1, 2]) 
    a = 1. 

    x[0, 0] = 10. # initial condition 1st solution 
    y[0, 0] = 5. 

    x[0, 1] = 0. # initial condition 2nd solution 
    y[0, 1] = 0.0000000001 

    for step in xrange(num_steps): 
     x[step + 1] = x[step] + h * y[step] 
     y[step + 1] = y[step] + h * (-x[step] - a * y[step] * (x[step] ** 2 + y[step] ** 2 - 1)) 

    return x, y 

現在我想進一步向量化的代碼，並保持在同一陣列中的x和y，我已經提出了以下解決方案：

def forward_euler_vector(): 
    num_steps = 10000 
    h = 0.01 

    x = np.zeros([num_steps + 1, 2, 2]) # steps, variables, number of solutions 
    a = 1. 

    x[0, 0, 0] = 10. # initial conditions 1st solution 
    x[0, 1, 0] = 5. 

    x[0, 0, 1] = 0. # initial conditions 2nd solution 
    x[0, 1, 1] = 0.0000000001 

    def f(x): 
     return np.array([x[1], 
         -x[0] - a * x[1] * (x[0] ** 2 + x[1] ** 2 - 1)]) 

    for step in xrange(num_steps): 
     x[step + 1] = x[step] + h * f(x[step]) 

    return x 

的問題：forward_euler_vector（）的作品，但這是最好的方式來矢量化它？我問，因爲矢量版本的運行約20毫秒在我的筆記本電腦慢：

In [27]: %timeit forward_euler() 
1 loops, best of 3: 301 ms per loop 

In [65]: %timeit forward_euler_vector() 
1 loops, best of 3: 320 ms per loop 

Answer 1

@Ophion意見解釋非常清楚是怎麼回事。在f(x)內對array()的調用引入了一些開銷，這在表達式h * f(x[step])中消除了使用矩陣乘法的好處。

正如他所說，您可能有興趣查看scipy.integrate以獲得一組不錯的數值積分器。

爲了解決矢量化代碼的問題，您希望每次調用f時都不要重新創建數組。您想要初始化數組一次，並在每次調用時將其返回修改。這類似於C/C++中的static變量。

您可以在函數f(x)的定義時使用可變默認參數來實現此功能，該參數會被解釋一次，並且具有局部範圍。由於它必須是可變的，你將其封裝在一個單一的元素的列表：

def f(x,static_tmp=[empty((2,2))]): 
    static_tmp[0][0]=x[1] 
    static_tmp[0][1]=-x[0] - a * x[1] * (x[0] ** 2 + x[1] ** 2 - 1) 
    return static_tmp[0] 

有了這個修改你的代碼，陣列創建的開銷就消失了，我的機器上我收穫了小的改進：

%timeit forward_euler()  #258ms 
%timeit forward_euler_vector() #248ms 

這意味着優化與numpy的矩陣乘法的收益是相當小的，至少在手頭的問題。

您可能希望立即擺脫函數f，在for循環中執行其操作，從而消除通話開銷。然而，默認參數的這個技巧也可以用於更通用的時間積分器，其中您必須提供函數f。

編輯：如由Jaime指出，另一種方式去是治療static_tmp作爲函數f的屬性，並已宣佈後的功能，但稱這是之前來創建它：

def f(x): 
    f.static_tmp[0]=x[1] 
    f.static_tmp[1]=-x[0] - a * x[1] * (x[0] ** 2 + x[1] ** 2 - 1) 
    return f.static_tmp 
f.static_tmp=empty((2,2)) 

Answer 2

有總是瑣碎autojit解決方案：

def forward_euler(initial_x, initial_y, num_steps, h): 

    x = np.zeros([num_steps + 1, 2]) # steps, number of solutions 
    y = np.zeros([num_steps + 1, 2]) 
    a = 1. 

    x[0, 0] = initial_x[0] # initial condition 1st solution 
    y[0, 0] = initial_y[0] 

    x[0, 1] = initial_x[1] # initial condition 2nd solution 
    y[0, 1] = initial_y[1] 

    for step in xrange(int(num_steps)): 
     x[step + 1] = x[step] + h * y[step] 
     y[step + 1] = y[step] + h * (-x[step] - a * y[step] * (x[step] ** 2 + y[step] ** 2 - 1)) 

    return x, y 

時序：

from numba import autojit 
jit_forward_euler = autojit(forward_euler) 

%timeit forward_euler([10,0], [5,0.0000000001], 1E4, 0.01) 
1 loops, best of 3: 385 ms per loop 

%timeit jit_forward_euler([10,0], [5,0.0000000001], 1E4, 0.01) 
100 loops, best of 3: 3.51 ms per loop