Python：使用odeint實現閾值模型

我正在處理的問題（顯示的示例高度簡化）似乎是一個常見問題，但我還沒有找到解決方案。我有三種不同的反應，v1，v2和v3，定義如下：Python：使用odeint實現閾值模型

v1：R < - > A + C; v1 = k1 *（R-A * C/5000）

v2：R <→B + C; v2 = k2 *（R-B * C/5000）

v3：A + B→P; v3 = k3 * A * B

使用資源R，前兩個反應分別產生A和C以及B和C，而第三個反應將A和B轉換爲產物P（k1，k2，k3）常量，在這裏它們被設置爲1）。

只有當C超過了一個叫做Cthr的閾值（這裏：Cthr = 25）時，第三個反應纔會發生，否則v3是0.所以這個想法是C積累了，然後一旦達到一定的濃度，導致生產的產品P.

我實施瞭如下：

def thresholdmodel (yn,tvec,allpara,R): 

    (A, B, C, P) = yn 

    k1, k2, k3 = allpara['kv']  

    Cthr = allpara['Cthresh'] 

    if C <= Cthr: 
     v3 = 0 
    else:  
     v3 = k3*A*B 
     C = 0 #does not(!) affect the ouput, why? 

    v1 = k1*(R - A*C/5000.) 
    v2 = k2*(R - B*C/5000.) 

    dA = v1 - v3 
    dB = v2 - v3 
    dC = v1 + v2 
    dP = v3 

    return (dA, dB, dC, dP)

模擬輸出看起來像這樣：http://i50.tinypic.com/apdvkj.png

因此很明顯，它的工作原理，直到達到該門檻第一次（v3是0，P沒有產生），但是之後C沒有被設置爲0，我不知道爲什麼。
我想要得到的是：直到閾值降到0，再次產生C，降到0等等，看起來像鋸齒波。
P的時間過程應該看起來像樓梯（只有在超過Cthr時纔會產生）。

有誰知道我必須做什麼才能將C設置回0並獲得預期的輸出？非常感謝！

來源

2012-07-19 Cleb

我不明白你的公式，但我可以告訴爲什麼C不下降到0

因爲C是一個局部變量thresholdmodel()，改變C爲任意值不會改變odeint狀態。 odeint將始終合併dC由thresholdmodel()返回。所以C會不斷增加。

編輯：

C的溫度將持續增加，所以你需要改變閾值每一次C>閾值，如：

lastC = 0 

def thresholdmodel (yn,tvec,allpara,R): 
    global lastC 
    (A, B, C, P) = yn 

    k1, k2, k3 = allpara['kv']  

    Cthr = allpara['Cthresh'] 

    if C <= lastC + Cthr: 
     v3 = 0 
    else:  
     v3 = k3*A*B 
     lastC = C 

    v1 = k1*(R - A*C/5000.) 
    v2 = k2*(R - B*C/5000.) 

    dA = v1 - v3 
    dB = v2 - v3 
    dC = v1 + v2 
    dP = v3 

    return (dA, dB, dC, dP)

來源

2012-07-19 07:56:44 HYRY

好，非常感謝您的快速反應！那麼，這只是一個玩具的例子來說明我的問題。你不明白什麼？只有簡單的質量作用動力學應用，沒有什麼特別的。所以，謝謝澄清問題！你知道一個可行的方法來解決它嗎？ – Cleb 2012-07-19 08:19:14

我不知道<->和 - >的意思。如果你的代碼中C可以下降到0，那麼v3將幾乎爲0，並且當C == 0時只有k3 * A * B。這是你的意圖嗎？ – HYRY 2012-07-19 10:50:24

再次感謝你的回覆！ '<->'表示反應是可逆的（雙向），' - >'表示它是不可逆的（只有一個方向）。好吧，對於小於閾值（Cthr）的所有Cs，v3幾乎爲零。一旦達到該閾值，v3變爲k3 * A * B並且可以產生P.之後，C再次設置爲零，因此v3變爲零，C再次建立，直到下一次達到閾值，依此類推。但是因爲C只是一個局部變量，所以在總體輸出中將其設置爲零似乎很困難，正如您所指出的那樣。任何建議來解決這個問題？ – Cleb 2012-07-20 02:48:09

Python：使用odeint實現閾值模型

回答

相關問題