當前迭代在scipy odeint

Question

我正在使用Scipy的odeint（scipy.integrate.odeint）爲我解決一些ODE，所有工作都很好，很好。但是，現在我想在計算中包含另一個與時間相關的數據集，即t = [0, 1, 2, 3]我已將數據z = [0.1, 0.2, 0.25, 0.22]包含在計算中。我可以通過向量作爲參數，但是這給了我每個時間步的整個向量。有沒有一種有效的方法來獲取計算的當前步驟（迭代器）？這樣我可以獲得z[i]第i個時間步。請注意，z的長度爲t，並且兩者都可以包含數千個元素。

感謝

一個很簡單的例子：

import numpy as np 
from scipy.integrate import odeint 

def func(y, t, z): 
    # I'd like to get the i-th element 
    # of z, corresponding to t[i] 
    return y+z[i] 

result = odeint(func, [0], t, (z,)) 

Answer 1

此問題的變通解決方案是使用更通用scipy.integrate.ode功能。這個函數有幾個內置的集成方案，並且您可以更多地控制每次迭代過程中發生的情況。請參見下面的示例：

import numpy as np 
from scipy.integrate import ode 

def func(t, y, z): 
    return y+z 

t = np.linspace(0, 1.0, 100) 
dt = t[1]-t[0] 
z = np.random.rand(100) 
output = np.empty_like(t) 
r = ode(func).set_integrator("dop853") 
r.set_initial_value(0, 0).set_f_params(z[0]) 

for i in xrange(len(t)): 
    r.set_f_params(z[i]) 
    r.integrate(r.t+dt) 
    output[i] = r.y 

在每次迭代期間，的z的解算器的值被相應地更新。

Answer 2

對於求解器，您可以使用interpolation作爲時變輸入。在這種情況下：

import numpy as np 
from scipy.integrate import odeint 
from scipy.interpolate import interp1d 

t_arr = np.array([0,1,2,3]) 
z_arr = np.array([0.1, 0.2, 0.25, 0.22]) 
finterp = interp1d(t_arr,z_arr,fill_value='extrapolate') #create interpolation function 

def func(y,t,z): 
    print(t) 
    zt = finterp(t) # call interpolation at time t 
    return y+zt 

result = odeint(func, [0], t_arr, (z_arr,)) 

然而，解算器在odeint可以請求超出t_arr時刻（在你的情況下，在3.028），所以你必須超越T = 3 Z指定的值或允許與fill_value外推（如上所示）。選擇一種插值時要小心，以反映您期望的行爲。