使用 odeint 确定分叉点

Question

我正在尝试获取分叉图的数据，即触发器 [A] 的值（由 Aspace 的 linspace 指定）给出此集合中 [N] 的稳态（空斜线）值微分方程：

由于方程组太长，无法尝试通过重新排列获得零斜线，我尝试使用 odeint 在给定时间内积分，并获得所有 d[x]/d(t ) 到 0，并且 return 那些特定空斜线的 [A] 和 [N] 的值。从理论上讲，这段代码应该 return 预期的值，但是它似乎在没有进一步解释的情况下杀死了内核，使我无法追查问题。请找到下面的代码（不要害怕，其中大部分只是说明变量值）。

我的代码：

#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)

给出结果（缩短）

  0.00,   0.000099999991, [  9.99994911e-05   9.99994911e-05   9.99789864e-11   9.99999905e-05]
 10.00,   0.002883872464, [  7.25813895e+00   7.25813895e+00   5.26700470e-01   2.88387246e-03]
 20.00,   0.008780897525, [  1.21628435e+01   1.21628435e+01   1.47905181e+00   8.78089753e-03]
 30.00,   0.017500871488, [ 16.0785845   16.0785845    2.58469186   0.01750087]
 40.00,   0.030166487759, [ 19.40580443  19.40580443   3.76509943   0.03016649]
 50.00,   0.049390699930, [ 22.32996985  22.32996985   4.98527848   0.0493907 ]
 60.00,   0.081339098240, [ 24.95672023  24.95672023   6.22713342   0.0813391 ]
 70.00,   0.144113671629, [ 27.35699725  27.35699725   7.48255648   0.14411367]
100.00,  63.103453906658, [ 52.49789017  52.49789017  27.55477377  63.10345391]
110.00,  64.363708289042, [ 53.43024046  53.43024046  28.54219752  64.36370829]
120.00,  65.425943398185, [ 54.25586731  54.25586731  29.43110515  65.4259434 ]
130.00,  66.343605498371, [ 55.00078061  55.00078061  30.24480971  66.3436055 ]
140.00,  67.150809538919, [ 55.68201657  55.68201657  30.99866996  67.15080954]
150.00,  67.870747329663, [ 56.31143752  56.31143752  31.70343926  67.87074733]

使用 odeint 确定分叉点

Using odeint to determine bifurcation

python

scipy

differential-equations

odeint