如何從文本文件中讀取微分方程組來解決scipy.odeint的系統問題？

Question

我有一個龐大（> 2000方程）的ODE系統，我想用python scipy的odeint解決。

我有三個問題，我想解決（也許我會問3個不同的問題？）。 爲了簡單起見，我會在這裏用玩具模型來解釋它們，但請記住我的系統很大。 假設我有ODE的以下系統：

dS/dt = -beta*S 
dI/dt = beta*S - gamma*I 
dR/dt = gamma*I 

使用β= C p我

，其中C，P和γ是我想傳遞給odeint參數。

odeint期待一個文件是這樣的：

def myODEs(y, t, params): 
    c,p, gamma = params 
    beta = c*p 
    S = y[0] 
    I = y[1] 
    R = y[2] 
    dydt = [-beta*S*I, 
      beta*S*I - gamma*I, 
      - gamma*I] 
    return dydt 

那則可以傳遞給odeint這樣的：

myoutput = odeint(myODEs, [1000, 1, 0], np.linspace(0, 100, 50), args = ([c,p,gamma],)) 

我產生了數學的文本文件，說myOdes.txt，其中文件的每一行對應於我的ODE系統的RHS，所以看起來像這樣

#myODEs.txt 

-beta*S*I 
beta*S*I - gamma*I 
- gamma*I 

我的文本文件看起來類似於odeint的期望，但我還沒有完成。 我有三個主要問題：

1. 我怎樣才能把我的文本文件，以便odeint明白，這是我的系統的RHS？
2. 如何以一種智能的方式，即系統地定義我的變量？由於它們有2000多個，我不能手動定義它們。理想情況下，我會在一個單獨的文件中定義它們，並讀取它們。
3. 如何將參數（也有很多）作爲文本文件傳遞？

我讀this question接近於我的問題1和2，並試圖複製它（我直接把值參數，這樣，我沒有擔心我的3點以上）：

systemOfEquations = [] 
    with open("myODEs.txt", "r") as fp : 
     for line in fp : 
      systemOfEquations.append(line) 

    def dX_dt(X, t): 
     vals = dict(S=X[0], I=X[1], R=X[2], t=t) 
     return [eq for eq in systemOfEquations] 

    out = odeint(dX_dt, [1000,1,0], np.linspace(0, 1, 5)) 

，但我得到了錯誤：

odepack.error: Result from function call is not a proper array of   floats. 
    ValueError: could not convert string to float: -((12*0.01/1000)*I*S), 

編輯：我修改了代碼：

systemOfEquations = [] 
    with open("SIREquationsMathematica2.txt", "r") as fp : 
     for line in fp : 
       pattern = regex.compile(r'.+?\s+=\s+(.+?)$') 
       expressionString = regex.search(pattern, line) 
       systemOfEquations.append(sympy.sympify(expressionString)) 


    def dX_dt(X, t): 
     vals = dict(S=X[0], I=X[1], R=X[2], t=t) 
     return [eq for eq in systemOfEquations] 

    out = odeint(dX_dt, [1000,1,0], np.linspace(0, 100, 50),) 

和這個工程（我不完全得到for循環的前兩行正在做什麼）。但是，我想要更自動地定義變量的過程，我仍然不知道如何使用此解決方案並將參數傳遞到文本文件中。沿着同樣的路線，我怎樣才能在dX_dt函數中定義參數（取決於變量）？

在此先感謝！

Answer 1

這不是一個完整的答案，而是一些觀察/問題，但它們太長而無法評論。

dX_dt被odeint多次調用，其中1D陣列y和元組t。您通過args參數提供t。 y由odeint生成並隨每個步驟而變化。 dX_dt應該簡化，以便運行速度更快。

通常像[eq for eq in systemOfEquations]這樣的表達式可以簡化爲systemOfEquations。 [eq for eq...]沒有做任何有意義的事情。但可能有些需要它的systemOfEquations。

我建議你打印出systemOfEquations（對於這個小小的3線情況），無論是爲了你和我們的利益。您正在使用sympy將字符串從文件翻譯成等式。我們需要看看它產生了什麼。

請注意myODEs是一個函數，而不是文件。它可以從模塊導入，當然這是一個文件。

指向vals = dict(S=X[0], I=X[1], R=X[2], t=t)的意思是生成sympy表達式可以使用的字典。更直接的（我想更快）dX_dt功能看起來像：

def myODEs(y, t, params): 
    c,p, gamma = params 
    beta = c*p 
    dydt = [-beta*y[0]*y[1], 
      beta*y[0]*y[1] - gamma*y[1], 
      - gamma*y[1]] 
    return dydt 

我猜想，運行sympy產生表達式dX_dt會比一個「硬編碼」一個像這樣的慢了許多。

我要添加sympy標籤，因爲正如所寫的那樣，這是將文本文件翻譯爲odeint可以使用的功能的關鍵。

我傾向於將t參數中的方程變化性，而不是sympy表達式列表。

也就是說取代：

dydt = [-beta*y[0]*y[1], 
      beta*y[0]*y[1] - gamma*y[1], 
      - gamma*y[1]] 

的東西，如

arg12=np.array([-beta, beta, 0]) 
    arg1 = np.array([0, -gamma, -gamma]) 
    arg0 = np.array([0,0,0]) 
    dydt = arg12*y[0]*y[1] + arg1*y[1] + arg0*y[0] 

一旦這是正確的，那麼argxx定義可以是移動之外dX_dt，並通過args通過。現在dX_dt只是一個簡單而快速的計算。

這整個sympy方法可能正常工作，但恐怕在實踐中它會很慢。但有更多sympy經驗的人可能會有其他見解。