用odeint确定分岔

Question

我试图得到分岔图的数据，即触发器[A]的值（由Aspace的linspace指定），在这组微分方程中给出了[N]的稳态（nullcline）值：

由于方程组太长，无法尝试通过重新排列来获得空斜，因此我尝试使用odeint在给定时间内进行积分，并将所有d[x]/d（t）的值设为0，并为这些特定的空斜返回[A]和[N]的值。理论上，这段代码应该返回预期的值，但是它似乎在没有进一步解释的情况下杀死了内核，阻止了我追踪这个问题。请找到下面的代码（不要害怕，大部分代码只是声明变量值）。在

我的代码：

#basic conditions as specified in the material
OS = 0
O = 0
S= 0
N = 0
B = 0

n1 = 1*(10**-4)
a1 = 1
a2 = 0.01
a3 = 0.2
n2 = 1*(10**-7)
b1 = 0.0011
b2 = 0.001
b3 = 0.0007
n3 = 1*(10**-4)
c1 = 1
c2 = 0.01
c3 = 0.2
n4 = 1*(10**-7)
d1 = 0.0011
d2 = 0.001
d3 = 0.0007
n5 = 1*(10**-4)
e1 = 0.005
e2 = 0.1
n6 = 1*(10**-7)
f1 =0.001
f2 = 9.95*(10**-4)
f3 = 0.01
k1c = 0.05
k2c = 0.001
k3c = 5
g1 = 1
g2 = 1
g3 = 1
def differential_eq(y, t, A, B):
    O, S, OS, N = y
    dydt = np.empty(4)
    dydt[0] = ((n1 + a1*A + a2*OS + a3*OS*N)/(1 + n2 + b1*A + b2*OS + b3*OS*N)) - g1*O - k1c*O*S + k2c*OS
    dydt[1] = (n3 + c1*A + c2*OS + c3*OS*N)/(1 + n4 + d1*A + d2*OS + d3*OS*N) - g2*S - k1c*O*S + k2c*OS
    dydt[2] = k1c*O*S -k2c*OS - k3c*OS
    dydt[3] = ((n5 + e1*OS + e2*OS*N)/(1 + n6 +f1*OS + f2*OS*N + f3*B)) - g3*N
    return dydt

timepool =np.linspace(0,100,10)

def simulate(A):
A = A
y = (0,0,0,0)
B =0
solution = sp.odeint (differential_eq(y, timepool, A, B), (A, B), timepool)
    solution = sp.odeint (differential_eq, initial, timepool)
    if dydt[0] == 0 and dydt[1] ==0 and dydt[2] ==0 and dydt[3] ==0:
        print (A,N, solution)
        return (A, N)

Aspace = np.linspace(0,150,150)
for A in Aspace:
    simulate(A)

现在，这似乎行不通，而且我看不出原因的任何迹象，因为它只会杀死内核：如果有人知道为什么会这样，请告诉我。在

Answer 1

结合注释中的注释和其他内容，您的simulate函数应该如下所示

def simulate(A):
    y = (0,0,0,0)
    B =0
    # parameter for the accuracy level
    eps = 1e-9
    # pass the correct arguments: first the function reference, then
    # initial point and sample times, the additional arguments and
    # the absolute and relative tolerances (they are combined as atol+rtol*(norm(y))
    solution = sp.odeint (differential_eq,y, timepool, args = (A, B), atol=1e-4*eps, rtol=1e-3*eps)
    y, t = solution[-1], timepool[-1]
    # compute the slopes at the last sample point by calling the ODE function
    dydt = differential_eq(y,t,A,B)
    # variables not explicitly declared global are local to the ODE function, 
    # they do not influence the global variables. 
    O, S, OS, N = y
    # never compare floating point numbers with "==" in numerics
    # always account for the floating point noise accumulated in the computation
    if max(abs(dydt)) < eps :
        #here print only the last sample point 
        print "%6.2f, %16.12f, %s"%(A, N, y)

    return (A, N)

得出结果（缩短）

^{pr2}$