我正在尝试使用python中的odeint数字地解决以下系统::
我的问题是k1取决于p的值。当p<=0.01
k1=0
,否则为k1=5
。
我不确定如何使用上述约束编写系统。我可以开始:
def sys(x,t,flip):
S,T,p = x
dx = [1-T-k1*T,
0.1*(1-S)-k1*S,
-0.2*(1-T-k1*T)+0.1*(1-S)-k1*S]
return dx
flip
可以是0.01
的值,但是我不确定这种情况下如何实现if语句,因为我需要集成过程来在每次迭代后检查p
的值。
这里是使用solve_ivp的解决方案:
import numpy as np
import matplotlib.pylab as plt
from scipy.integrate import solve_ivp
def ode(t, TS):
T, S = TS
p = -0.2*T + S
k1 = 0 if p <= 0.01 else 5
dTdt = 1 - T - k1*T
dSdt = 0.1*(1 - S) - k1*S
return dTdt, dSdt
# Solve
TS_start = (0.7, 0.4)
t_span = [0, 3]
sol = solve_ivp(ode, t_span, TS_start, method='RK45',
rtol=1e-5)
print(sol.message)
print('nbr eval', sol.nfev)
# Graph
plt.plot(sol.t, sol.y[0, :], label='T')
plt.plot(sol.t, sol.y[1, :], label='S'); plt.legend();
plt.xlabel('time');
这似乎很难解决。其他经过测试的求解器无法收敛,并且事件的使用无济于事,因为它们在稳态解附近越来越频繁[]
编辑:的确将翻转值更改为-0.01而不是+0.01会导致极限循环解,对此events
的solve_ivp
方法起作用:
import numpy as np import matplotlib.pylab as plt from scipy.integrate import solve_ivp def ode(t, TS): T, S = TS p = -0.2*T + S k1 = 0 if sign_change(t, TS) <= 0 else 5 dTdt = 1 - T - k1*T dSdt = 0.1*(1 - S) - k1*S return dTdt, dSdt def sign_change(t, TS): T, S = TS p = -0.2*T + S flip = -0.01 return p - flip # Solve TS_start = [0.3, 0.1] t_span = [0, 5] sol = solve_ivp(ode, t_span, TS_start, method='RK45', events=sign_change) print(sol.message) print('nbr eval', sol.nfev) # Graph plt.plot(sol.t, sol.y[0, :], '-x', label='T') plt.plot(sol.t, sol.y[1, :], label='S'); plt.legend(); plt.xlabel('time'); for t in sol.t_events[0]: plt.axvline(x=t, color='gray', linewidth=1)
现在的解决方案是:
要在积分后重建k1
和p
的值,请将此计算放入附加函数中>>
def kp(T,S):
p = -0.2*T + S
k1lo, k1hi = 0, 5
flip = -0.01
k1 = 0.5*(k1lo+k1hi + (k1hi-k1lo)*np.tanh((p-flip)/1e-8))
return k1,p
要在积分后重建k1
和p
的值,请将此计算放入附加函数中>>
def kp(T,S):
p = -0.2*T + S
k1lo, k1hi = 0, 5
flip = -0.01
k1 = 0.5*(k1lo+k1hi + (k1hi-k1lo)*np.tanh((p-flip)/1e-8))
return k1,p