Implementing scipy.integrate.RK45 for SIR model

Question

I looked at the documents for scipy.integrate.RK45, but I could not find good examples.

I'm trying to implement the ODE for an infectious disease, with pre-defined parameters beta, f, gamma and initial values for the groups of population (Susceptible, Exposed, Infectious, Recovered) S0, E0, I0, R0. (I will post these values if I need to, but I think these are irrelevant for my question.)

I know how to implement scipy.integrate.odeint:

def SIR_model(y, t, f, beta, gamma):
    S, E, I, R = y
    dS = -beta*S*I 
    dE = beta*S*I - f*E 
    dI = f*E - gamma*I 
    dR = gamma*I
    return [dS, dE, dI, dR]

sol = scipy.integrate.odeint(SIR_model, [S0,E0,I0,R0], t, args=(f, beta,gamma))

does the job. My goal right now is to use the RK45 method. I realized that this method needs to be implemented differently, but I got raised an error when I tried the following:

scipy.integrate.RK45(fun=SIR_model, t0 = t[0], y0 = [S0,E0,I0,R0],t_bound=t[-1])

# raised TypeError: SIR_model() missing 3 required positional arguments: 'f', 'beta', and 'gamma'

scipy.integrate.RK45(fun=SIR_model(y,t,f,beta,gamma), t0 = t[0], y0 = [S0,E0,I0,R0],t_bound=t[-1])

# raised TypeError: cannot unpack non-iterable NoneType object

y0 = [S0,E0,I0,R0]
scipy.integrate.RK45(fun=SIR_model(y0,t,f,beta,gamma), t0 = t[0], y0 = [S0,E0,I0,R0],t_bound=t[-1])

# raised TypeError: 'list' object is not callable

Any suggestions would be welcome.

Answer 1

Use solve_ivp with method="RK45", this is mostly similar to odeint in terms of automatization. Take care that the ODE function needs t as first argument.

RK45 itself is a stepper class, you would have to implement the main time loop yourself, which sometimes offers greater flexibility.

---

See Cannot import X problem. Stiff ODE solver for Oregonator Model for a longer example with the Radau stepper.