2

我正在处理的问题（所示示例高度简化）似乎是一个常见问题，但我还没有找到解决方案。我有三种不同的反应，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，前两个反应分别产生 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，我不知道为什么。
我想得到的是：C一直产生到阈值，下降到0，再次产生，下降到0等等，看起来像一个锯齿波。
P 的时间过程应该看起来像楼梯（仅在超过 Cthr 时产生）。

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

4

1 回答 1

2

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

因为C是局部变量thresholdmodel()，所以C对任何值的更改都不会改变 中的状态odeint。odeint将始终集成dC返回的thresholdmodel(). 所以C会不断增加。

编辑：

C 会不断增加，所以每次 C > threshold 都需要改变阈值，如：

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-19T07:56:44.677 回答