我想要求解一个 ODE 系统,在前 30,000 秒内,我希望我的状态变量之一从相同的初始值开始。 30,000 秒之后,我想将该状态变量的初始值更改为不同的值,并在其余时间模拟系统。这是我的代码:
def ode_rhs(y, t):
ydot[0] = -p[7]*y[0]*y[1] + p[8]*y[8] + p[9]*y[8]
ydot[1] = -p[7]*y[0]*y[1] + p[8]*y[8]
ydot[2] = -p[10]*y[2]*y[3] + p[11]*y[9] + p[12]*y[9]
ydot[3] = -p[13]*y[3]*y[6] + p[14]*y[10] + p[15]*y[10] - p[10]*y[2]*y[3] + p[11]*y[9] + p[9]*y[8] - p[21]*y[3]
ydot[4] = -p[19]*y[4]*y[5] - p[16]*y[4]*y[5] + p[17]*y[11] - p[23]*y[4] + y[7]*p[20]
ydot[5] = -p[19]*y[4]*y[5] + p[15]*y[10] - p[16]*y[4]*y[5] + p[17]*y[11] + p[18]*y[11] + p[12]*y[9] - p[22]*y[5]
ydot[6] = -p[13]*y[3]*y[6] + p[14]*y[10] - p[22]*y[6] - p[25]*y[6] - p[23]*y[6]
ydot[7] = 0
ydot[8] = p[7]*y[0]*y[1] - p[8]*y[8] - p[9]*y[8]
ydot[9] = p[10]*y[2]*y[3] - p[11]*y[9] - p[12]*y[9] - p[21]*y[9]
ydot[10] = p[13]*y[3]*y[6] - p[14]*y[10] - p[15]*y[10] - p[22]*y[10] - p[21]*y[10] - p[23]*y[10]
ydot[11] = p[19]*y[4]*y[5] + p[16]*y[4]*y[5] - p[17]*y[11] - p[18]*y[11] - p[22]*y[11] - p[23]*y[11]
ydot[12] = p[22]*y[10] + p[22]*y[11] + p[22]*y[5] + p[22]*y[6] + p[21]*y[10] + p[21]*y[3] + p[21]*y[9] + p[24]*y[13] + p[25]*y[6] + p[23]*y[10] + p[23]*y[11] + p[23]*y[4] + p[23]*y[6]
ydot[13] = p[15]*y[10] + p[18]*y[11] - p[24]*y[13]
return ydot
pysb.bng.generate_equations(model)
alias_model_components()
p = np.array([k.value for k in model.parameters])
ydot = np.zeros(len(model.odes))
y0 = np.zeros(len(model.odes))
y0[0:7] = p[0:7]
t = np.linspace(0.0,1000000.0,100000)
r = odeint(ode_rhs,y0,t)
所以,换句话说,我想在前 30,000 秒内每次调用
odeint
时将 y0[1] 设置为相同的值(100)。在将信号输入系统之前,我实际上试图让系统平衡一段时间。我想过做像 if t < 30000: y0[1] = 100
这样的事情作为我的 ode_rhs()
函数的第一行,但我不太确定这是否有效。
听起来您希望 y1(t) 保持恒定(值为 100) 平衡阶段。您可以通过确保在此期间 dy1(t)/dt = 0 来做到这一点 阶段。有(至少)两种方法可以实现这一目标。第一个是 修改
ydot[1]
中ode_rhs
的计算如下:
if t < 30000:
ydot[1] = 0.0
else:
ydot[1] = -p[7]*y[0]*y[1] + p[8]*y[8]
并使用初始条件 100 作为
y[1]
。
请注意,这会在系统右侧引入不连续性, 但
odeint
(Fortran 代码 LSODA)使用的自适应求解器通常足够强大来处理它。
这是一个独立的示例。我对
p
提出了 t1
和 ode_rhs
论据。
t1
是平衡阶段的持续时间。
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
def ode_rhs(y, t, p, t1):
ydot[0] = -p[0]*y[0]*y[1] + p[1]*y[2] + p[2]*y[2]
if t < t1:
ydot[1] = 0.0
else:
ydot[1] = -p[0]*y[0]*y[1] + p[1]*y[2]
ydot[2] = p[0]*y[0]*y[1] - p[1]*y[2] - p[2]*y[2]
return ydot
ydot = np.zeros(3)
p = np.array([0.01, 0.25, 0.1])
y0 = [20.0, 100.0, 0.0]
t = np.linspace(0, 200, 2001)
t1 = 20.0
sol = odeint(ode_rhs, y0, t, args=(p, t1))
plt.figure(1)
plt.clf()
plt.subplot(3, 1, 1)
plt.plot(t, sol[:, 0])
plt.axvline(t1, color='r')
plt.grid(True)
plt.ylabel('y[0]')
plt.subplot(3, 1, 2)
plt.plot(t, sol[:, 1])
plt.axvline(t1, color='r')
plt.grid(True)
plt.ylabel('y[1]')
plt.ylim(0, 110)
plt.subplot(3, 1, 3)
plt.plot(t, sol[:, 2])
plt.axvline(t1, color='r')
plt.grid(True)
plt.ylabel('y[2]')
plt.xlabel('t')
plt.show()
上述方法的一个细微变化是通过添加一个来修改系统 参数为 0 或 1。当参数为 0 时,求解平衡系统;当参数为 1 时,求解整个系统。
ydot[1]
(在我的较小示例中)的代码是
ydot[1] = full * (-p[0]*y[0]*y[1] + p[1]*y[2])
其中
full
是参数。
为了处理平衡阶段,系统在 0 <= t < t1 with
full=0
上求解一次。然后将平衡解的最终值用作
第二个解的初始条件,使用 full=1
运行。这种方法的优点是您不会强迫求解器处理不连续性。
代码如下。
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
def ode_rhs(y, t, p, full):
ydot[0] = -p[0]*y[0]*y[1] + p[1]*y[2] + p[2]*y[2]
ydot[1] = full * (-p[0]*y[0]*y[1] + p[1]*y[2])
ydot[2] = p[0]*y[0]*y[1] - p[1]*y[2] - p[2]*y[2]
return ydot
ydot = np.zeros(3)
p = np.array([0.01, 0.25, 0.1])
y0 = [20.0, 100.0, 0.0]
t1 = 20.0 # Equilibration time
tf = 200.0 # Final time
# Solve the equilibration phase.
teq = np.linspace(0, t1, 100)
full = 0
soleq = odeint(ode_rhs, y0, teq, args=(p, full))
# Solve the full system, using the final point of the
# equilibration phase as the initial condition.
y0 = soleq[-1]
# Note: the system is autonomous, so we could just as well start
# at t0=0. But starting at t1 makes the plots (below) align without
# any additional shifting of the time arrays.
t = np.linspace(t1, tf, 2000)
full = 1
sol = odeint(ode_rhs, y0, t, args=(p, full))
plt.figure(2)
plt.clf()
plt.subplot(3, 1, 1)
plt.plot(teq, soleq[:, 0], t, sol[:, 0])
plt.axvline(t1, color='r')
plt.grid(True)
plt.ylabel('y[0]')
plt.subplot(3, 1, 2)
plt.plot(teq, soleq[:, 1], t, sol[:, 1])
plt.axvline(t1, color='r')
plt.grid(True)
plt.ylabel('y[1]')
plt.ylim(0, 110)
plt.subplot(3, 1, 3)
plt.plot(teq, soleq[:, 2], t, sol[:, 2])
plt.axvline(t1, color='r')
plt.grid(True)
plt.ylabel('y[2]')
plt.xlabel('t')
plt.show()
这是它生成的图(第一个示例的图是 几乎完全相同):
这个答案对我不起作用,因为我需要定期改变我的初始条件。因此,我想提出替代解决方案,即在函数本身内部改变微分方程的条件:
我们查看t的值并调整它:
if int(t) > 20:
full = 1
else:
full = 0
这里是一个函数内部:
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
def ode_rhs(y, t, p, full):
if int(t) > 20:
full = 1
else:
full = 0
ydot[0] = -p[0]*y[0]*y[1] + p[1]*y[2] + p[2]*y[2]
ydot[1] = full * (-p[0]*y[0]*y[1] + p[1]*y[2])
ydot[2] = p[0]*y[0]*y[1] - p[1]*y[2] - p[2]*y[2]
return ydot
ydot = np.zeros(3)
# intial conditions
p = np.array([0.01, 0.25, 0.1])
y0 = [20.0, 100.0, 0.0]
t = np.linspace(0, 200, 100)
full = 0
# solve equation
solution = odeint(ode_rhs, y0, t, args=(p, full))
plt.figure()
plt.clf()
plt.subplot(3, 1, 1)
plt.plot(t, solution[:, 0])
plt.axvline(20, color='r') # vertical line
plt.grid(True)
plt.ylabel('y[0]')
plt.subplot(3, 1, 2)
plt.plot(t, solution[:, 1])
plt.axvline(20, color='r') # vertical line
plt.grid(True)
plt.ylabel('y[1]')
plt.ylim(0, 110)
plt.subplot(3, 1, 3)
plt.plot(t, solution[:, 2])
plt.axvline(20, color='r') # x=20 vertical line
plt.grid(True)
plt.ylabel('y[2]')
plt.xlabel('t')
plt.show()
这允许调用函数一次来求解方程。
更重要的是,您现在可以定期调整方程内的参数。例如,假设您有 t = [0:200] 并且您想要每 20 步更改 full 的值,您可以这样做:
if int(t/20) % 2 == 0:
full = 1
else:
full = 0