It sounds like you want y1(t) to remain constant (with the value 100) for the
equilibration phase. You can do this by ensuring that dy1(t)/dt = 0 during this
phase. There are (at least) two ways you can accomplish that. The first is to
modify the calculation of ydot[1]
in ode_rhs
as follows:
if t < 30000:
ydot[1] = 0.0
else:
ydot[1] = -p[7]*y[0]*y[1] + p[8]*y[8]
and use the intitial condition 100 for y[1]
.
Note that this introduces a discontinuity in the right-hand side of the system,
but the adaptive solver used by odeint
(the Fortran code LSODA) is usually robust enough to handle it.
Here's a self-contained example. I've made p
and t1
arguments to ode_rhs
.
t1
is the duration of the equilibration phase.
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()
A slight variation of the above method is to modify the system by adding a
parameter that is either 0 or 1. When the parameter is 0, the equlibration system is solved, and when the parameter is 1, the full system is solved. The code for ydot[1]
(in my smaller example) is then
ydot[1] = full * (-p[0]*y[0]*y[1] + p[1]*y[2])
where full
is the parameter.
To handle the equilibration phase, the system is solved once on 0 <= t < t1 with
full=0
. Then the final value of the equilibration solution is used as the
initial condition to the second solution, run with full=1
. The advantage of this method is that you are not forcing the solver to deal with the discontinuity.
Here's how it looks in code.
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):