Python:使用 odeint 实现阈值模型

Question

我正在处理的问题(所示的示例是高度简化的)似乎是一个常见问题，但我还没有找到解决方案。我有三种不同的 react ，v1、v2 和 v3，定义如下:

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，前两个 react 分别产生 A 和 C 以及 B 和 C，而第三个 react 将 A 和 B 转化为产物 P(k1、k2、k3 是常数，这里将它们设置为 1 )。

只有当 C 超过称为 Cthr 的特定阈值(此处:Cthr = 25)时，第三个 react 才会发生，否则 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，我不知道为什么。
我想要得到的是:C产生直到阈值，下降到0，再次产生，下降到0等等，看起来像锯齿波。
P 的时间进程应该看起来像楼梯(仅在超过 Cthr 时产生)。

有谁知道我必须做什么才能将 C 设置回 0 并收到预期的输出？非常感谢!

Answer 1

我不明白你的方程，但我可以告诉为什么 C 不会下降到 0。

由于 C 是 thresholdmodel() 中的局部变量，因此将 C 更改为任何值都不会更改 odeint 中的状态。 odeint 将始终集成 thresholdmodel() 返回的 dC。所以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)