Solving ODE set with a step-function parameter using odeint

Question

I'm trying to solve the following system numerically with odeint in Python::

My problem is that k1 depends on the value of p. When p<=0.01 k1=0 and otherwise k1=5.

I'm not sure how to write the system with the above constraint. I can start with:

def sys(x,t,flip):
S,T,p = x

dx = [1-T-k1*T,
     0.1*(1-S)-k1*S,
     -0.2*(1-T-k1*T)+0.1*(1-S)-k1*S]
return dx

Where flip can be the 0.01 value, but I'm not sure how to implement the if statement for this case, because I need the integration process to check the value of p after each iteration.

Answer 1

Here is a solution using solve_ivp:

import numpy as np
import matplotlib.pylab as plt

from scipy.integrate import solve_ivp

def ode(t, TS):
    T, S = TS

    p = -0.2*T + S
    k1 = 0 if p <= 0.01 else 5

    dTdt = 1 - T - k1*T
    dSdt = 0.1*(1 - S) - k1*S

    return dTdt, dSdt

# Solve
TS_start = (0.7, 0.4)
t_span = [0, 3]
sol = solve_ivp(ode, t_span, TS_start, method='RK45',
                rtol=1e-5)
print(sol.message)
print('nbr eval', sol.nfev)

# Graph
plt.plot(sol.t, sol.y[0, :], label='T')
plt.plot(sol.t, sol.y[1, :], label='S'); plt.legend();
plt.xlabel('time');

It seems to be a difficult problem to solve numerically. Other tested solvers don't converge, and the use of events is not helpful as they are more and more frequent near the steady state solution

Edit: indeed changing the flipping value to -0.01 instead of +0.01 leads to a limit cycle solution, for which the events method of solve_ivp is working:

import numpy as np
import matplotlib.pylab as plt

from scipy.integrate import solve_ivp

def ode(t, TS):
    T, S = TS

    p = -0.2*T + S
    k1 = 0 if sign_change(t, TS) <= 0 else 5

    dTdt = 1 - T - k1*T
    dSdt = 0.1*(1 - S) - k1*S

    return dTdt, dSdt

def sign_change(t, TS):
    T, S = TS
    p = -0.2*T + S
    flip = -0.01
    return p - flip

# Solve
TS_start = [0.3, 0.1]
t_span = [0, 5]
sol = solve_ivp(ode, t_span, TS_start,
                method='RK45', events=sign_change)
print(sol.message)
print('nbr eval', sol.nfev)

# Graph
plt.plot(sol.t, sol.y[0, :], '-x', label='T')
plt.plot(sol.t, sol.y[1, :], label='S'); plt.legend();
plt.xlabel('time');

for t in sol.t_events[0]:
    plt.axvline(x=t, color='gray', linewidth=1)

now the solution is: