At a fixed input the production terms are constant, so both linear ordinary differential equations can also be solved exactly:Each displacement from its equilibrium concentration decays with time constant one. The Jacobian matrix is , with two eigenvalues equal to , so the fixed point is globally asymptotically stable for the nonnegative state space.
Here is a numerical solver using the classical fourth-order Runge-Kutta method; it uses only the repository's existing NumPy and matplotlib dependencies. Continue the definitions from part (b):
def integrate(rhs, initial, times):
out = np.empty((len(times), len(initial)))
out[0] = initial
for j, dt in enumerate(np.diff(times)):
state = out[j]
k1 = rhs(state)
k2 = rhs(state + dt*k1/2)
k3 = rhs(state + dt*k2/2)
k4 = rhs(state + dt*k3)
out[j+1] = state + dt*(k1 + 2*k2 + 2*k3 + k4)/6
return out
def direct_rhs(state, X):
Y, Z = state
return np.array([2*hill(X, 1, 2) - Y,
3*hill(X, 2, 4) - Z])
times = np.linspace(0, 8, 801)
plt.figure(figsize=(7, 4))
for initial, style in [([0, 0], "-"), ([3, 4], "--")]:
trajectory = integrate(lambda s: direct_rhs(s, 2), initial, times)
exact = np.array([1.6, 1.5]) + (np.array(initial) - [1.6, 1.5]) \
* np.exp(-times[:, None])
assert np.max(np.abs(trajectory - exact)) < 1e-8
plt.plot(times, trajectory[:, 0], style, label=f"Y, initial {initial}")
plt.plot(times, trajectory[:, 1], style, label=f"Z, initial {initial}")
plt.xlabel("Time")
plt.ylabel("Concentration")
plt.legend()
plt.tight_layout()
plt.show()For , . The initial state approaches those values from below, while approaches from above. Over eight time units the initial displacement is reduced by , explaining why both trajectories nearly coincide with the equilibria from part (b). The two-panel figure above displays these simulations.
Articles by others on the same topic
There are currently no matching articles.