Fourth Order Runge Kutta for SIR model

Question

I am implementing a fourth order Runge Kutta method to solve a system of three ODEs that describe the SIR model. Based on other simulations, I would expect a different solution (mine diverges very quickly). Here is my code for the runge kutta method:

def rk4(f, s0, x0, xf, n):
  x = np.linspace(x0, xf, n+1) #x grid
  s = np.array((n+1)*[s0]) #array of the state of the sytem for each x
  h = x[1] - x[0] #stepsize
  for i in range(n):       #Fourth Order Runge Kutta Method
    k0 = h * f(s[i]) 
    k1= h * f(s[i] + 0.5 * k0)
    k2 = h * f(s[i] + 0.5 * k1)
    k3 = h * f(s[i] + k2)
    s[i+1] = s[i] + (k0 + 2*(k1 + k2) + k3) / 6.
  return x, s

And here is the implementation for the SIR model:

def f(u):
  dS = -a*u[0]*u[1] #dY/dt = -aS(t)I(t)
  dI = a*u[0]*u[1] - b*u[1]  #dI/dt = aS(t)I(t)-bI(t)
  dR = b*u[1]  #dR/dt = bI(t)
  return np.array([dS, dI, dR])

for:

a = 0.2
b = 0.1
x, s = rk4(f, np.array([235.0, 14.0, 0.0]), 0., 109., 109000)

the solution I get is this although I expected to approach these data. What am I doing wrong?

Thank you in advance!

"the solution is different from what I expected. What am I doing wrong?": well, what did you expect? And what did you get? — 9769953, Jun 03 '22 at 22:38
Is there a reason you implemented your own RK4 method, instead of using a existing one? — 9769953, Jun 03 '22 at 22:38
Was it intentional that your derivative function `f` takes `x` as a parameter but does not use it? — ddejohn, Jun 03 '22 at 23:40
Can you check if the model you used is correct? It might be that you took the formulation for densities, but your calculation is with population numbers. The practical difference is that the infection term gets divided by the population number `N=sum(u)`, then `dS = -a*u[0]*u[1]/N`, the same in the second equation. — Lutz Lehmann, Jun 04 '22 at 12:21

score 0 · Accepted Answer · answered Jun 06 '22 at 11:45

With the given parameters a=1/5 and b=1/10 and the general context of infection models, it is reasonable to guess that the intention is that the model has one transmission per infected person every 5 days, and a resolution of every infection in 10 days on average.

The model to implement this via population numbers, whatever unit they may have, is

def f(u):
  S,I,R = u
  N = sum(u)
  dS = -a*S*I/N 
  dI = a*S*I/N - b*I  
  dR = b*I
  return np.array([dS, dI, dR])

Only if the variables are the population densities the sum N will always be one and can be omitted.

With this change the code as given produces the image

which looks very similar to the reference image

Thank you very much! That was the problem indeed. – George Ntino Jun 06 '22 at 15:34 — George Ntino, Jun 06 '22 at 15:34

Fourth Order Runge Kutta for SIR model

1 Answers1