2

私が扱っている問題 (示されている例は非常に単純化されています) は一般的な問題のようですが、まだ解決策を見つけていません。v1、v2、v3 の 3 つの異なる反応があり、次のように定義されています。

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 を使用して、最初の 2 つの反応はそれぞれ A と C および B と C を生成しますが、3 番目の反応は A と B を生成物 P に変換します (k1、k2、k3 は定数で、ここでは 1 に設定されています)。

3 番目の反応は、C が Cthr (ここでは Cthr = 25) と呼ばれる特定のしきい値を超えた場合にのみ発生すると想定されています。それ以外の場合は v3 は 0 です。製品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

方程式はわかりませんが、 がC0 にならない理由はわかります。

Cは のローカル変数であるため、値をthresholdmodel()変更Cしても のステータスは変わりませんodeint。によって返されるodeint常に統合されます。だからどんどん増えていく。dCthresholdmodel()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-19T07:56:44.677 に答える