听起来您希望 y1(t) 保持恒定(值为 100)
平衡阶段。您可以通过确保在此期间 dy1(t)/dt = 0 来做到这一点
阶段。有(至少)两种方法可以实现这一目标。第一个是
修改计算ydot[1]
in ode_rhs
如下:
if t < 30000:
ydot[1] = 0.0
else:
ydot[1] = -p[7]*y[0]*y[1] + p[8]*y[8]
并使用初始条件 100y[1]
.
请注意,这会在系统右侧引入不连续性,
但使用的自适应求解器odeint
(Fortran 代码 LSODA)通常足够强大来处理它。
这是一个独立的示例。我做了p
and 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])
where full
是参数。
为了处理平衡阶段,系统在 0 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()
And here's the plot that it generates (the plot from the first example is
almost exactly the same):