Simulation of SIR model for COVID-19

In the present page we numerically solve the initial value proble of the SIR model for COVID-19 of the form \begin{alignat*}{4} \frac{dS}{dt} & = -\beta IS, & \quad \frac{dI}{dt} & = \beta IS - \gamma I, & \quad \frac{dR}{dt} & = \gamma I, \quad t\in[0,\infty), \\ S(0) & = S_0, & \quad I(0) & = I_0, & \quad R(0) & = R_0, & \quad \end{alignat*} where \(S(t)\), \(I(t)\) and \(R(t)\) are real-value unknown functions of \(t \in \mathbb{R}\), \(S_0\), \(I_0\) and \(R_0\) are given initial data, and \(\beta\) and \(\gamma\) are positive constants. In mathematical epidemiology, \(t\) is time, \(S(t)\) is the number of susceptible people, \(I(t)\) is the number of people infected, \(R(t)\)s the number of people who have recovered and developed immunity to the infection, \(\beta\) the infection rate, and \(\gamma\) is the recovery rate. This systems of ordinary differential equations is well-known as the SIR model for infection disease.

We briefly mention the basic properties of the SIR model. We are concerned with the simulation of SIR model for COVID-19 below. Set \(T(t)=I(t)+R(t)\) which is the total cases at the time \(t\). A unit of time here is supposed to be day, and set the initial date, say, today as \(t=0\). Visit the government website dealing with the statistics about COVID-19, and get the total cases \(T(0)\), the daily new cases \(T(0)-T\), the total number of discharged individuals \(R(0)\) and the daily new numbers of discharged individuals \(R(0)-R(-1)\). Set \(S(0):=N-T(0)\) and \(I(0):=T(0)-R(0)\). See, e.g., South Korea, Hong Kong, Okinawa Osaka Tokyo. Note that a daily increase of a quantity can be seen as a differentiation of the quantity on the day. So we set the infection and recovery rates by \[ \beta := \frac{T(0)-T(-1)}{(T(0)-R(0))(N-T(0))}, \quad \gamma := \frac{R(0)-R(-1)}{T(0)-R(0)}. \] since \[ \frac{dT}{dt} = \beta(T-R)(N-T), \quad \frac{dR}{dt} = \gamma(T-R). \] Then we have \[ \rho_e(0) = \frac{T(0)-T(-1)}{R(0)-R(-1)}. \] By using MATLAB we draw a solution to the initial value problem for the SIR model with data including area, population \(N\), staring date, ending date, \(T(0)\), \(T(0)-T(-1)\), \(R(0)\) and \(R(0)-R(-1)\). See sir1.m for the detail.

The behavior of solutions to the SIR model depends sensitively on \(\rho_e(0)\), and the above simulation follows the daily change of \(\rho_e(0)\). So we take the average of some \(d\) days for \(\beta\) and \(\gamma\) instead of the values determined only by the data of the initial date. We replace them by \[ \beta := \frac{T(0)-T(-d)}{d\bigl(T(0)-R(0)\bigr)(N-T(0))}, \quad \gamma := \frac{R(0)-R(-d)}{d\bigl(T(0)-R(0)\bigr)}. \] for some \(d=1,2,3,\dotsc\). Then the effective reproduction number of the initial date becomes \[ \rho_e(0) = \frac{T(0)-T(-d)}{R(0)-R(-d)}. \] By using Julia Language we draw a solution to the initial value problem for the SIR model with data including area, population \(N\), staring date \(d0\), ending date \(d1\), the date of the begining of taking average \(d2=d0-d\), \(T(0)\), \(T(-d)\), \(R(0)\) and \(R(-d)\). If we fix the area, \(N\), \(d_1\) (or \(d1-d0\)), \(d2\), \(T(-d)\) and \(R(-d)\), we can repeat the simulation only by inputting \(d0\), \(T(0)\) and \(R(0)\). See sir2.m and sir3.m for the detail.