Efficient relaxation scheme for the SIR and related compartmental models

It is well known that the SIR model can be solved from a transcendental differential equation. Here, we revisit how to get such an equation to complete our analysis of the proposed scheme below. Let $\mu=\beta/\gamma$ be the reciprocal relative removal rate. By the second and third equations of system (I.1), we have

$$\frac{dS}{dR}=-\mu S,$$

which leads to $\ln\left(S\right)=-\mu R+c$. Therefore, we arrive at

$$S\left(t\right)=e^{c-\mu R\left(t\right)}.$$

(II.1)

To find $c$, we set $t=0$ in (II.1) and use (A3). Indeed, by $n=S\left(0\right)=e^{c-\mu R\left(0\right)}=e^{c}$, we get $c=\ln\left(n\right)$ and thus, deduce that

$$S\left(t\right)=ne^{-\mu R\left(t\right)}.$$

(II.2)

Combining this and (A1), we derive the following nonlinear differential equation for $R\left(t\right)$:

$$R^{\prime}\left(t\right)=\gamma\left(N-ne^{-\mu R\left(t\right)}-R\left(t\right)\right).$$

(II.3)

Observe that if one can approximate $R\left(t\right)$ well in (II.3), then $S\left(t\right)$ and $I\left(t\right)$ will be well approximated via (II.2) and (I.1), respectively. Define $g\left(r\right)=\gamma ne^{-\mu r}+\gamma r$ for $r\in\mathbb{R}$. We can rewrite (II.3) as

$$R^{\prime}\left(t\right)=\gamma N-g\left(R\left(t\right)\right).$$

Let $\left\{R_{k}\right\}_{k=0}^{\infty}$ be a time-dependent sequence satisfying, for $k=1,2,3,\ldots$,

$$R_{k}^{\prime}\left(t\right)+MR_{k}\left(t\right)=\gamma N-g\left(R_{k-1}\left(t\right)\right)+MR_{k-1}\left(t\right).$$

(II.5)

The sequence $\left\{R_{k}\right\}_{k=0}^{\infty}$ aims to approximate $R\left(t\right)$ in (II.3) in the sense that $R_{k}$ will be close to $R$ as $k\to\infty$ uniformly in time. The accompanying initial condition for equation (II.5) is $R_{k}\left(0\right)=0$ for any $k\geq 0$. Since our approximate model performs as an iterative scheme, we need a starting point, $R_{0}\left(t\right)$. Here, we choose $R_{0}\left(t\right)=0$ based on the initial condition of $R\left(t\right)$ (cf. (A3)), which is the best information given to the sought $R\left(t\right)$.

Also, in (II.5), we introduce a $k$-independent constant $M\geq\gamma>0$ for the so-called relaxation process. This relaxation term plays a very important role. It allows us to prove the non-negativity of the relaxation scheme, while many numerical approaches, including the regular linearization method, do not have or cannot prove this feature. Herewith, the regular linearization method we meant is the scheme $\left\{R_{k}\right\}_{k=0}^{\infty}$ in (II.5) with either $M=0$ or only the exponential term in $g\left(r\right)$ being linearized.

Formulated below is the theorem showing that the scheme $\left\{R_{k}\right\}_{k=0}^{\infty}$ preserves the non-negativity of the removals over the relaxation process.

In the following, we formulate the strong convergence result for the scheme $\left\{R_{k}\right\}_{k=0}^{\infty}$. For ease of presentation, proof of this result is deliberately placed in the Appendix. It is worth mentioning that proof of the strong convergence of the scheme relies so much on the strict estimation of $g^{\prime}$. Such an estimation can merely be obtained by the aid of the non-negativity of the scheme.

Our theoretical finding below shows that when $n\mu<1$, the scheme $\left\{R_{k}\left(t\right)\right\}_{k=0}^{\infty}$ converges faster than the case $n\mu\geq 1$. The proof of the following corollary is also found in the Appendix.

As readily expected, for every step $k$, we obtain a non-homogeneous differential equation that can be effectively approximated in the discrete framework. Mimicking the proof of Theorem II.3 and using Theorem 2, we have

$$R_{k}^{\prime}\left(t\right)+MR_{k}\left(t\right)\leq\gamma N-\gamma n+\left(M-\gamma\right)R_{k-1}\left(t\right),$$

(II.7)

leading to

	$\displaystyle R_{k}\left(t\right)$	$\displaystyle\leq e^{-Mt}\int_{0}^{t}e^{Ms}\left(\gamma a+\left(M-\gamma\right)R_{k-1}\left(s\right)\right)ds$
		$\displaystyle\leq t\gamma a+\left(M-\gamma\right)\int_{0}^{t}R_{k-1}\left(s\right)ds.$

By induction and by the choice $R_{0}\left(t\right)=0$, we can show that

$$\left|R_{k}\left(t\right)\right|\leq\gamma a\sum_{i=1}^{k}\left(M-\gamma\right)^{i-1}\frac{t^{i}}{i!}.$$

Therefore, if we choose $\gamma\leq M\leq\gamma+\frac{1}{T}$, then $\left|R_{k}\left(t\right)\right|\leq T\gamma a\left(e-1\right)$. Note that this bound is independent of $k$. Thus, by Theorem 1, we get $R_{k}\in C^{1}$ for any $k$ with, cf. (II.7), Theorem 2 and the choice $M\leq\gamma+\frac{1}{T}$,

$$\left|R_{k}^{\prime}\left(t\right)\right|\leq\gamma a+\left(M-\gamma\right)T\gamma a\left(e-1\right)\leq\gamma ae.$$

Furthermore, by differentiating both sides of (II.5) with respect to time, we can demonstrate that $R_{k}\in C^{2}$ for any $k$. Indeed,

	$\displaystyle R_{k}^{\prime\prime}\left(t\right)$	$\displaystyle+MR_{k}^{\prime}\left(t\right)$
		$\displaystyle=\left(M-\gamma+\gamma n\mu e^{-\mu R_{k-1}\left(t\right)}\right)R_{k-1}^{\prime}\left(t\right)$
		$\displaystyle\leq\left(M-\gamma+\gamma n\mu\right)\left|R_{k-1}^{\prime}\left(t\right)\right|.$

This yields that $\left|R_{k}^{\prime\prime}\left(t\right)\right|\leq 2M-\gamma+\gamma n\mu$, which is a $k$-independent upper bound. This ensures the Euler method’s global error by leveraging its existing convergence theory. For completeness, we present below the discrete solution to our proposed scheme.

Consider the time increment $\Delta t=T/P$ for $P\geq 2$ being a fixed integer. Then, we set the mesh-point in time by $t_{p}=p\Delta t$ for $0\leq p\leq P$. We seek $R_{k}^{p}\approx R_{k}\left(t_{p}\right)$ as a discrete solution to equation (II.5). By the standard Euler method, $R_{k}^{p}$ is determined by the following equation:

	$\displaystyle R_{k}^{p}$	$\displaystyle+\Delta tMR_{k}^{p}$
		$\displaystyle=R_{k}^{p-1}+\Delta t\left(\gamma N-g\left(R_{k-1}^{p}\right)+MR_{k-1}^{p}\right).$		(II.8)

By this way, the global error of the Euler method is attained in the sense that for every $k$, there exists a constant $\tilde{C}>0$ such that

$$\max_{0\leq p\leq P}\left|R_{k}^{p}-R_{k}\left(t_{p}\right)\right|\leq\tilde{C}\Delta t.$$

We accentuate that by the above analysis of $R_{k}$, i.e. $R_{k}\in C^{2}$ for any $k$, the constant $\tilde{C}$ is independent of $k$. Thus, by Theorem 3, we can estimate the distance between the discrete (approximate) solution $R_{k}^{p}$ and the true solution $R$ at each mesh-point,

$$\max_{0\leq p\leq P}\left|R_{k}^{p}-R\left(t_{p}\right)\right|\leq\tilde{C}\Delta t+\sqrt{\frac{C^{k}}{k!}}\max_{0\leq t\leq T}\left|R\left(t\right)\right|.$$

(II.9)

The SIRD model extends the SIR model by distinguishing between recovered and deceased individuals. In this framework, the removals in the SIR model no longer encompass the number of infected individuals who have passed away. To account for mortality, a mortality rate $\sigma>0$ is introduced, representing the rate at which infected individuals succumb to death. Consequently, the death rate per unit of time is calculated as the product of the mortality rate and the number of infected individuals. Additionally, as the number of deceased individuals is excluded from the removals, the rate of change of infections over time is adjusted to reflect the loss caused by mortality. Mathematically, the SIRD model reads as

$$\begin{cases}I^{\prime}\left(t\right)=\beta S\left(t\right)I\left(t\right)-\left(\gamma+\sigma\right)I\left(t\right),\\ S^{\prime}\left(t\right)=-\beta S\left(t\right)I\left(t\right),\\ R^{\prime}\left(t\right)=\gamma I\left(t\right),\\ D^{\prime}\left(t\right)=\sigma I\left(t\right),\end{cases}$$

(III.1)

where $D\left(t\right)$ stands for the number of deceased people (after infection) at time $t$. In (III.1), we make use of the following assumptions.

(B1) The total population size is always conserved with $N>0$, meaning that $S\left(t\right)+I\left(t\right)+R\left(t\right)+D\left(t\right)=N$ for all $t$.

(B2) We know the infection rate $\beta>0$ from the infection process, the removal rate $\gamma>0$ from the removal process, and the death rate $\sigma>0$ from the mortality process.

(B3) The initial conditions are $S\left(0\right)=n>1$, $I\left(0\right)=a\geq 1$, $R\left(0\right)=0$ and $D\left(0\right)=0$.

Now, we detail the transcendental equation for $R\left(t\right)$ and the application of the relaxation scheme. From the second and third equations of system (III.1), we deduce that

$$S\left(t\right)=ne^{-\mu R\left(t\right)},$$

(III.3)

where we have recalled the reciprocal relative removal constant $\mu=\beta/\gamma$. Using the same way, the third and last equations of system (III.1) give

$$D\left(t\right)=\frac{\sigma}{\gamma}R\left(t\right)$$

(III.4)

by virtue of $R\left(0\right)=D\left(0\right)=0$ (cf. (B3)). Then, plugging (III.3), (III.4) and (B1) into the third equation of (III.1), we obtain the following differential equation for $R\left(t\right)$:

$$R^{\prime}\left(t\right)=\gamma\left[N-ne^{-\mu R\left(t\right)}-\left(1+\frac{\sigma}{\gamma}\right)R\left(t\right)\right].$$

(III.5)

Henceforth, our relaxation scheme in this case becomes

$$R_{k}^{\prime}\left(t\right)+\overline{M}R_{k}\left(t\right)=\gamma N-\overline{g}\left(R_{k-1}\left(t\right)\right)+\overline{M}R_{k-1}\left(t\right),$$

(III.6)

where $\overline{g}\left(r\right)=\gamma ne^{-\mu r}+\left(\gamma+\sigma\right)r$. Similar to the SIR model, here we rely on (B3) to choose $R_{k}\left(0\right)=0$ for any $k\geq 0$ as the initial condition and $R_{0}\left(t\right)=0$ as the starting point.

By choosing $\overline{M}\geq\gamma+\sigma$, our sequence $\left\{R_{k}\right\}_{k=0}^{\infty}$ (defined in (III.6)) is non-negativity-preserving and globally strongly convergent to $R$ of the transcendental equation (III.5). These findings are analogous to our central Theorems 2 and 3, and therefore, we omit their formulations. Besides, Theorem 1 is applied to (III.5), guaranteeing the global existence and uniqueness of the $C^{1}$ solutions to the SIRD model (III.1).

The SIR model, along with its variants SIRD and SIRX, assumes a constant population size. These models, known as epidemiological SIR-type models without vital dynamics, are limited in their representation of population changes; see (A1), (B1) and (C1). The SIR model with vital dynamics addresses this limitation by incorporating birth and death rates to account for population size fluctuations.

In the present work, we explore that the transcendental system governing the SIR model with background mortality takes the form of (II.5). With $\sigma$ being the death rate, the population experiences changes over time. Here, individuals from all compartments can exit through deaths, allowing for a more realistic representation of population dynamics. Mathematically, the SIR model with background mortality can be expressed as follows:

$$\begin{cases}I^{\prime}\left(t\right)=\beta S\left(t\right)I\left(t\right)-\gamma I\left(t\right)-\sigma I\left(t\right),\\ S^{\prime}\left(t\right)=-\beta S\left(t\right)I\left(t\right)-\sigma S\left(t\right),\\ R^{\prime}\left(t\right)=\gamma I\left(t\right)-\sigma R\left(t\right).\end{cases}$$

(III.7)

In this perspective, we make use of the following assumptions.

(C1) The total population size is dependent of $t$, i.e. $N=N\left(t\right)>0$. It can be computed that $N\left(t\right)=S\left(t\right)+I\left(t\right)+R\left(t\right)=e^{-\sigma t}N_{0}$ for some fixed $N_{0}>0$.

(C2) We know the infection rate $\beta>0$ from the infection process, the removal rate $\gamma>0$ from the removal process, and the death rate $\sigma>0$ from the mortality process.

(C3) The initial conditions are $S\left(0\right)=n>1$, $I\left(0\right)=a\geq 1$ and $R\left(0\right)=0$. This implies that $N_{0}=n+a$.

Similar to the classical SIR model, we seek the transcendental equation for $R$ prior to the application of our proposed numerical scheme. When doing so, we define $\overline{R}\left(t\right)=e^{\sigma t}R\left(t\right)$ and $\overline{S}\left(t\right)=e^{\sigma t}S\left(t\right)$. By the second and third equations of system (III.7), we find that

	$\displaystyle\overline{S}^{\prime}\left(t\right)$	$\displaystyle=-\beta\overline{S}\left(t\right)I\left(t\right),$		(III.8)
	$\displaystyle\overline{R}^{\prime}\left(t\right)$	$\displaystyle=\gamma e^{\sigma t}I\left(t\right).$		(III.9)

Therefore, we deduce that

$\displaystyle\frac{d\overline{S}}{d\overline{R}}$

$\displaystyle=-\mu e^{-\sigma t}\overline{S},$

(III.10)

or equivalently, $\ln\left(\overline{S}\right)=-\mu e^{-\sigma t}\overline{R}+\tilde{c}\left(t\right)$. Herewith, we have recalled the reciprocal relative removal rate $\mu=\beta/\gamma$. Henceforth, we have

$$\overline{S}\left(t\right)=e^{\tilde{c}\left(t\right)-\mu e^{-\sigma t}\overline{R}\left(t\right)}.$$

(III.11)

Since $\overline{S}\left(0\right)=n$ and $\overline{R}\left(0\right)=R\left(0\right)=0$ by (D3), we find that $\tilde{c}\left(0\right)=\ln\left(n\right)$. Moreover, by taking the derivative in time of (III.11), we arrive at

$$\overline{S}^{\prime}\left(t\right)=e^{\tilde{c}\left(t\right)-\mu e^{-\sigma t}\overline{R}\left(t\right)}\left[\tilde{c}^{\prime}\left(t\right)-\mu e^{-\sigma t}\left(-\sigma+\overline{R}^{\prime}\left(t\right)\right)\right].$$

(III.12)

Dividing both sides of (III.12) by $\overline{R}^{\prime}\left(t\right)$, we find that

	$\displaystyle\frac{\overline{S}^{\prime}\left(t\right)}{\overline{R}^{\prime}\left(t\right)}$
	$\displaystyle=\frac{e^{\tilde{c}\left(t\right)-\mu e^{-\sigma t}\overline{R}\left(t\right)}\left[\tilde{c}^{\prime}\left(t\right)-\mu e^{-\sigma t}\left(-\sigma+\overline{R}^{\prime}\left(t\right)\right)\right]}{\overline{R}^{\prime}\left(t\right)}.$

Then combining this with (III.10) and (III.11), we derive the following differential equation for $\tilde{c}\left(t\right)$:

	$\displaystyle e^{\tilde{c}\left(t\right)-\mu e^{-\sigma t}\overline{R}\left(t\right)}\left[\tilde{c}^{\prime}\left(t\right)-\mu e^{-\sigma t}\left(-\sigma+\overline{R}^{\prime}\left(t\right)\right)\right]$
	$\displaystyle=-\mu e^{-\sigma t}e^{\tilde{c}\left(t\right)-\mu e^{-\sigma t}\overline{R}\left(t\right)}\overline{R}^{\prime}\left(t\right),$

or equivalently, $\tilde{c}^{\prime}\left(t\right)=-\mu\sigma e^{-\sigma t}$. Thus, we obtain $\tilde{c}\left(t\right)=\ln\left(n\right)+\mu\left(e^{-\sigma t}-1\right)$ and

$$\overline{S}\left(t\right)=ne^{\mu\left(e^{-\sigma t}-1\right)}e^{-\mu e^{-\sigma t}\overline{R}\left(t\right)}.$$

Together with the back-substitution $e^{-\sigma t}\overline{R}\left(t\right)=R\left(t\right)$, we thereby get $\overline{S}\left(t\right)=ne^{\mu\left(e^{-\sigma t}-1\right)}e^{-\mu R\left(t\right)}$. Now, we note that by (C1) and (C3), $e^{\sigma t}I\left(t\right)=N_{0}-\overline{S}\left(t\right)-\overline{R}\left(t\right)$ holds true for any $t$. Plugging this into (III.9) and using the fact that $\overline{S}\left(t\right)=ne^{\mu\left(e^{-\sigma t}-1\right)}e^{-\mu R\left(t\right)}=ne^{\mu\left(e^{-\sigma t}-1\right)}e^{-\mu e^{-\sigma t}\overline{R}\left(t\right)}$, we derive the transcendental equation for $\overline{R}$ as follows:

$$\overline{R}^{\prime}\left(t\right)=\gamma\left[N_{0}-ne^{\mu\left(e^{-\sigma t}-1\right)}e^{-\mu e^{-\sigma t}\overline{R}\left(t\right)}-\overline{R}\left(t\right)\right].$$

(III.13)

By setting $\hat{g}\left(r\right)=\gamma ne^{\mu\left(e^{-\sigma t}-1\right)}e^{-\mu e^{-\sigma t}r}+\gamma r$, our relaxation scheme for the SIR model with background mortality is structured by

	$\displaystyle\overline{R}_{k}^{\prime}\left(t\right)$	$\displaystyle+\hat{M}\overline{R}_{k}\left(t\right)$
		$\displaystyle=\gamma N_{0}-\hat{g}\left(\overline{R}_{k-1}\left(t\right)\right)+\hat{M}\overline{R}_{k-1}\left(t\right).$		(III.14)

Similar to the above-mentioned SIR-based models, we choose $\overline{R}_{k}\left(0\right)=0$ for any $k\geq 0$ as the initial condition and $\overline{R}_{0}\left(t\right)=0$ as the starting point, based on the fact that $\overline{R}\left(0\right)=R\left(0\right)=0$. Also, here we take $\hat{M}\geq\gamma$ to ensure the non-negativity preservation and global strong convergence of the sequence $\left\{\overline{R}_{k}\right\}_{k=0}^{\infty}$ (defined in (III.14)) to the sought $\overline{R}$ of the transcendental equation (III.13). As another analog of Theorems 2 and 3, we omit details of the formulations of the theoretical results for the sequence $\left\{\overline{R}_{k}\right\}_{k=0}^{\infty}$. It is also worth mentioning that Theorem 1 remains true in this case, providing the global existence and uniqueness of $C^{1}$ solutions to the SIR model with mortality (III.7).