# 使用odeint求解具有阶跃函数参数的ODE

##### 问题描述投票：0回答：2

``````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`的值。

python scipy ode differential-equations
##### 2个回答
1

``````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');
``````

``````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)
```现在的解决方案是：```

2

``````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
``````

0

``````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
``````