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实现阈值模型

回答

相关问题