Get the equilibrium value in coupled ode by python

I am dealing with a coupled ODE. I have already plotted the solutions out using odeint, but I want to get the value of equilibrium. The ode looks like this:

IN,EnN,c,l,XN,N_init,Bx_init=0.1,0.01,0.1,0.1,0.1,0.4,1
init = (N_init, Bx_init)
def model(init,t, IN, EnN, c, Qn, l, XN):
N_prime = IN - EnN - c*init[1]*init[0] + XN*l*Qn*init[1]
Bx_prime = c*init[1]*init[0]/Qn - l*init[1]
return [N_prime, Bx_prime]


where init, t, IN, EnN, c, Qn, l, XN are all constants. The plot I got looks like this: Bx and N. I tried using functions like solve_ivp, fsolve, root_scalar, however, I couldn't get the exact time t when Bx or N reaches the equilibrium i.e. get the t when dN/dt=0 and dBx/dt=0.

i followed Solving for a set of coupled ODEs to get correct variable values and tied to use solve_ivp

sol=solve_ivp(model, (0,40),init, method="LSODA", args=(0.1, 0.01, 0.1, 0.025, 0.1, 0.01))


However， I got the error

TypeError: 'float' object is not subscriptable


Could someone give me a hint? Thanks in advance.

• Please avoid cross-posts, or put at least links. Apr 26, 2022 at 18:53

The error you're seeing is because your model function has the first two arguments in the wrong order. solve_ivp takes the current time first, followed by y, then finally any extra arguments you want to pass.

As far as your question of finding the point where $$\frac{dN}{dt} = 0$$ and $$\frac{d B_x}{dt} = 0$$, mathematically the ODE asymptotically approaches this limit but never reaches it, so trying to check for this exactly is not going to give you any meaningful answer.

You can check for when $$\frac{dN}{dt}$$ and $$\frac{d B_x}{dt}$$ are smaller (by absolute value) than some defined limit, say $$10^{-2}$$ or $$10^{-3}$$ (this is up to you to figure out what a reasonable value is).

Here's an example script which does this:

from scipy.integrate import solve_ivp
from math import sqrt

def model(t, init, IN, EnN, c, Qn, l, XN):
N_prime = IN - EnN - c * init[1] * init[0] + XN * l * Qn * init[1]
Bx_prime = c * init[1] * init[0] / Qn - l * init[1]
return [N_prime, Bx_prime]

def detect_steady_state(t, init, IN, EnN, c, Qn, l, XN, tol):
diff = model(t, init, IN, EnN, c, Qn, l, XN)
l2_norm = sqrt(diff[0]**2+diff[1]**2)
return l2_norm - tol

args = (0.1,0.01,0.1,0.025,0.1,0.01)
tol=1e-3

model_cb = lambda t, y: model(t, y, *args)
detect_ss_cb = lambda t, y: detect_steady_state(t, y, *args, tol)
# stop solving once we've determined the "steady state" value and time
detect_ss_cb.terminal = True
# ensure we only stop when dy/dt falls below the tolerance, not starts below and rises above
detect_ss_cb.direction = -1

# the end time doesn't really matter as long as it's longer than the expected "steady state" time
sol = solve_ivp(model_cb, (0, 1e5), [0.4,1], method='LSODA', events=detect_ss_cb)
# time at which we reached "steady state"
print(sol.t_events[0])
# value of y when we reached "steady state"
print(sol.y_events[0][0])