Positive Numerical Approximation of Integro-Differential Epidemic Model

Abstract

In this paper, we study a dynamically consistent numerical method for the approximation of a nonlinear integro-differential equation modeling an epidemic with age of infection. The discrete scheme is based on direct quadrature methods with Gregory convolution weights and preserves, with no restrictive conditions on the step-length of integration h, some of the essential properties of the continuous system. In particular, the numerical solution is positive and bounded and, in cases of interest in applications, it is monotone. We prove an order of convergence theorem and show by numerical experiments that the discrete final size tends to its continuous equivalent as h tends to zero.

1. Introduction

Epidemic models based on nonlinear integral and integro-differential equations [1,2,3,4,5,6] allow to include the contribution of time since becoming infected to the total infectivity. This property makes these (age-of-infection) models suitable to describe diseases such as smallpox [7], cholera [8], SARS [9], COVID-19 [10,11,12] and AIDS [13,14]. Numerical simulations support the comprehension of the phenomenon and provide a tool for its quantitative description. Therefore, attention should be paid to set up a numerical framework that allows to supply real-time and reliable answers. We consider the integro-differential equation representing the Kermack and McKendrick age-of-infection epidemic model:

$$S^{\prime }\left(t\right)=\beta S\left(t\right){\int }_0^{\infty }{S}^{\prime }(t-s)A\left(s\right)ds,$$

(1)

for which we refer to [5,15] and references therein. Here, $S\left(t\right)$ is the number of susceptibles at time t, $\beta$ represents the rate of effective contacts and $A\left(s\right)\ge 0$ is the mean infectivity of members of the population with infection age s ( $A\in L^1\left({\mathbb{R}}^{+}\right)$). Equation (1) represents a special nonlinear Volterra integral equation and we refer to [16,17] for the theory about these problems. From the numerical point of view, there is a growing interest in the development and analysis of methods for the approximation of the solution to integral and integro-differential equations of Volterra type (see, for example, [16,18] and references cited therein). Recently, great attention is paid to the construction of methods that reflect and preserve the essential properties of the problem and that efficiently integrate the memory term (see, for example, [19,20]). In particular, for epidemic models, non-standard discretization methods allow the design of numerical models that qualitatively and quantitatively behave like the reference problem [21,22,23,24,25]. In [26], a non-standard discretization is used to integrate problem (1). Such a scheme is explicit, it preserves positivity of solutions, boundedness and the final size relation for the epidemic, and therefore it assumes the role of a discrete model. On the other hand, the method is only linearly convergent, and this can be a disadvantage when you have to integrate over long times, balancing accuracy and computational cost. Our aim here is to develop a higher order method. For numerical solutions obtained by direct discretizations of (1), results that guarantee the unconditional dynamic consistency properties of non-standard schemes are not straightforward. These considerations motivated us to use a different approach. The idea is to exploit the exponential form of the evolution operator in (1) and integrate it by direct quadrature methods of any order. The derived discrete equation has the advantage of automatically ensuring the positivity of the solution and, at the same time, allowing an asymptotic analysis that leads to results consistent with the continuous model. This paper is organized as follows: In Section 2, we report the main facts about the age-of-infection Model (1) as developed in [15,27]. Then, in Section 3, we formulate the numerical method, we study the convergence and prove the positivity and boundedness of the numerical solution, for any value of the stepsize. Furthermore, we analyze the asymptotic dynamics of the numerical model by deriving the numerical final size relation. In Section 4, we obtain additional properties of the numerical solution in cases of infectivity functions of interest in applications. Finally, numerical experiments are reported in Section 5, to show the theoretical results obtained. Some remarks, in Section 6, conclude the paper.

2. The Continuous Model

We refer to [15,27] for a complete description of Model (1). Here, we mention the main features that will be useful in the rest of the paper. In the following, we assume that the population has a constant size N and that

$${\int }_t^{\infty }{S}^{\prime }(t-s)A\left(s\right)ds=(S_0-N)A\left(t\right),$$

(2)

where $S_0=S\left(0\right)\le N$ (see, for example, [28]). The basic reproduction number is defined as

$$R_0=\beta N{\int }_0^{\infty }A\left(s\right)ds.$$

(3)

We want to point out that, for the solution $S\left(t\right)$ of (1), it holds:

• $S:\mathbb{R}\to {\mathbb{R}}^{+}$ is such that $t_1\le t_2\Rightarrow S\left(t_1\right)\ge S\left(t_2\right)$;
• ${lim}_{t\to +\infty }S\left(t\right)=S_{\infty }\in (0,S_0]$;
• $S_{\infty }$ is the unique solution of the final size relation.

  $$log\left(\frac{S_0}{S_{\infty }}\right)=R_0\left(1-\frac{S_{\infty }}{N}\right).$$

  (4)

Our starting point is the following equation

$$S\left(t\right)=S_{0}exp\left(-\beta {\int }_0^{t}A(t-s)\left(N-S\left(s\right)\right)ds\right),$$

(5)

which is equivalent to (1) with (2). This equivalence is obtained by using the identity

$${\int }_0^{x}A\left(s\right)S^{\prime }(x-s)ds={\displaystyle \frac{\partial }{\partial x}}\left({\int }_0^{x}A\left(s\right)S(x-s)ds\right)-A\left(x\right)S_0,$$

into Equation (1) with (2). It comes out that

$$\frac{S^{\prime }\left(x\right)}{S\left(x\right)}=\beta \left(-\left((N-S_0)A\left(x\right)+A\left(x\right)S_0\right)+\frac{\partial }{\partial x}\left({\int }_0^{x}A\left(s\right)S(x-s)ds\right)\right).$$

Then, integration in $[0,t],$ and straightforward calculations lead to the form (5).

3. The Numerical Model

Let $h>0$ and ${\left\{t_n\right\}}_{n\ge 0}$ be a uniform mesh such that $t_n=nh$. The direct discretization of (5) by direct quadrature ($DQ$) methods with Gregory convolution weights (see, for example, [18]) reads

$$S_n=S_{0}exp\left(-h\beta \left(\sum _{j=0}^{n_0-1}{w}_{nj}A\left(t_{n-j}\right)(N-S_j)+\sum _{j=n_0}^n{\omega }_{n-j}A\left(t_{n-j}\right)(N-S_j)\right)\right),$$

(6)

where $S_n\approx S\left(t_n\right)$, $n\ge n_0\ge 1,$ and $w_{nj}$ and ${\omega }_j$ are the non-negative quadrature weights. The starting values $S_0,S_1,\dots ,S_{n_0-1}$ are given. When $n_0=1,$${\omega }_0=0$ and ${\omega }_j=1,$ the numerical method above corresponds to the discrete-time Kermack–McKendrick model introduced by Diekmann in [29]. This motivates us to deepen the analysis of the dynamic behavior of (6) since, as it will be clear later in this section, from a numerical point of view it possesses good convergence properties and, in many cases, is easy to implement.

3.1. Basic Properties

At each step and for any fixed $h>0,$ Equation (6) implicitly defines $S_n$ as the solution of the nonlinear equation ${\mathcal{F}}_n(x;h)=0,$ where

$$\begin{array}{cc}\hfill {\mathcal{F}}_n(x;h)& =x-S_{0}exp\left(-h\beta {\omega }_{0}A\left(0\right)(N-x)\right)\hfill \\ \hfill & ·exp\left(-h\beta \left(\sum _{j=0}^{n_0-1}{w}_{nj}A\left(t_{n-j}\right)(N-S_j)+\sum _{j=n_0}^{n-1}{\omega }_{n-j}A\left(t_{n-j}\right)(N-S_j)\right)\right).\hfill \end{array}$$

(7)

The nonlinear function (7) is concave and ${\mathcal{F}}_n(0;h)<0.$ If we assume that $S_0,\dots ,S_{n-1}\in [0,N],$ then ${\mathcal{F}}_n(N;h)>0.$ Thus, $\mathcal{F}$ has exactly one root in $[0,N].$ Furthermore, this root is less than $S_0$. The application of this result leads to state the following theorem.

Theorem 1.

Let ${\left\{S_n\right\}}_{n\ge 0}$ be the solution to the discrete Equation (6). Assume that $0<S_j\le N$ for $j=0,\cdots ,n_0-1,$ then for each $h>0$:

• The sequence ${\left\{S_n\right\}}_{n\ge 0}$ is positive;
• The sequence ${\left\{S_n\right\}}_{n\ge 0}$ is bounded from above by $S_0$.

3.2. Convergence

From now on, we assume that the quadrature weights in (6) satisfy (see, for example, [18] (p. 79, Theorem 2.6.10))

$$\underset{n\ge 0}{sup}\underset{0\le j<n_0}{max}{w}_{nj}\le W<+\infty and\underset{n\ge 0}{sup}{\omega }_n\le \mathrm{\Omega }<+\infty .$$

(8)

Consider $T>0$ and $t\in [0,T],$ and let $h=T/M,$ with $M>0.$ Define, for $n=0,\dots ,M,$

$$\begin{array}{cc}\hfill \delta (h;t_n)& =exp\left(-\beta {\int }_0^{t_n}A(t_n-s)(N-S\left(s\right))ds\right)\hfill \\ \hfill & -exp\left(-h\beta \left(\sum _{j=0}^{n_0-1}{w}_{nj}A\left(t_{n-j}\right)(N-S\left(t_j\right))+\sum _{j=n_0}^n{\omega }_{n-j}A\left(t_{n-j}\right)(N-S\left(t_j\right))\right)\right).\hfill \end{array}$$

(9)

$\delta (h;t_n)$ is the local error for (6). It is easy to prove that $p=n_0+1$ is the largest integer so that

$$\underset{0\le n\le M}{max}\left|\delta (h;t_n)\right|\le c{h}^p.$$

(10)

The following convergence result holds:

Theorem 2.

Consider the approximate solution $S_n,$$n=0,\dots ,M,$ of (5) obtained by (6), with $h=T/M,$ and $M>0.$ Let the given function $A\left(t\right)$ be sufficiently smooth and assume that (10) holds and:

• the starting errors ${\eta }_j\left(h\right)=S\left(t_j\right)-S_j,$ for $0\le j<n_0-1,$ satisfy

  $$|{\eta }_j\left(h\right)|=\mathcal{O}\left(h^{p-1}\right);$$

  (11)
• the weights $w_{nj}$ and ${\omega }_n,$$j=0,\dots ,n_0-1,$$n=0,\dots ,M,$ satisfy (8).

Then the method (6) is convergent of order $p.$

Proof. 

For $n=n_0,\dots ,M,$

$$\begin{array}{cc}\hfill \frac{S\left(t_n\right)-S_n}{S_0}=& \delta (h;t_n)\hfill \\ & +exp\left(-h\beta \left(\sum _{j=0}^{n_0-1}{w}_{nj}A\left(t_{n-j}\right)(N-S\left(t_j\right))+\sum _{j=n_0}^n{\omega }_{n-j}A\left(t_{n-j}\right)(N-S\left(t_j\right))\right)\right)\hfill \\ & -exp\left(-h\beta \left(\sum _{j=0}^{n_0-1}{w}_{nj}A\left(t_{n-j}\right)(N-S_j)+\sum _{j=n_0}^n{\omega }_{n-j}A\left(t_{n-j}\right)(N-S_j)\right)\right),\hfill \end{array}$$

(12)

where $\delta (h;t_n)$ is the local truncation error defined in (9). Let $e(h;t_n)=S\left(t_n\right)-S_n,$ be the global error. Using standard manipulations and choosing $h<\frac{1}{\beta {\omega }_{0}A\left(0\right)S_0},$ from (12) we have

$$\left|e(h;t_n)\right|\le \frac{S_0\left|\delta (h;t_n)\right|}{1-h\beta {\omega }_{0}A\left(0\right)S_0}+h\frac{S_0\beta \underset{0\le t\le T}{max}A\left(t\right)·max\{W,\mathrm{\Omega }\}}{1-h\beta {\omega }_{0}A\left(0\right)S_0}\sum _{j=0}^{n-1}\left|e_j(h;t_n)\right|,n=n_0,\dots ,M.$$

The Gronwall discrete inequality (see, for example, [30] (p. 101)) yields

$$\begin{array}{cc}\hfill \left|e(h;t_n)\right|\le & \frac{1}{1-h\beta {\omega }_{0}A\left(0\right)S_0}·exp\left(\frac{\beta S_{0}T\underset{0\le t\le t_n}{max}A\left(t\right)·max\{W,\mathrm{\Omega }\}}{1-h\beta {\omega }_{0}A\left(0\right)S_0}\right)\hfill \\ \hfill & ·\left(S_0\underset{n_0\le n\le M}{max}\left|\delta (h;t_n)\right|+h{S}_0\beta \underset{0\le t\le T}{max}A\left(t\right)·max\{W,\mathrm{\Omega }\}\sum _{j=0}^{n_0-1}\left|{\eta }_j\left(h\right)\right|\right).\hfill \end{array}$$

Since the local error $\delta (h;t_n)$ and the starting errors ${\eta }_j\left(h\right)$ satisfy (10) and (11), respectively, it follows that

$$\underset{n=n_0,\dots ,M}{max}\left|e(h;t_n)\right|\le C{h}^p,$$

with C positive constant not depending on $h.$    □

3.3. The Numerical Final Size

In this section, we assume that

$$\forall h>0,\exists S_{\infty }\left(h\right)\in (0,S_0]:S_{\infty }\left(h\right)=\underset{n\to +\infty }{lim}{S}_n\left(h\right).$$

(13)

Moreover, we also assume that $h{\sum }_{n=0}^{\infty }{\omega }_{n}A\left(t_n\right)<+\infty$ (true if $A^{\prime }\left(t\right)\in L^1\left({\mathbb{R}}^{+}\right)$, see [31] (Lem.1), or if $A\left(t\right)$ is ultimately monotonic, see [32] (p. 208)). We define the numerical basic reproduction number

$$R_0\left(h\right)=h\beta N\sum _{n=0}^{+\infty }{\omega }_{n}A\left(t_n\right),$$

as the direct discretization (see [18]) of (3). In the Theorem below, we obtain a relation between $S_{\infty }\left(h\right)$ and $R_0\left(h\right)$ which represents the discrete version of (4).

Theorem 3.

Let $A\left(t\right)$ in (5) be sufficiently smooth and let ${\left\{S_n\right\}}_{n\ge 0}$ be the solution to the discrete Equation (6). If the starting values $S_j\le S_0$, for $j=0,\dots ,n_0-1$, are positive and (8) and (13) hold, then for each $h>0$

$$log\left(\frac{S_0}{S_{\infty }\left(h\right)}\right)=R_0\left(h\right)\left(1-\frac{S_{\infty }\left(h\right)}{N}\right).$$

(14)

Proof. 

We rewrite (6) as

$$log\left(\frac{S_0}{S_n}\right)=h\beta \left(\sum _{j=0}^{n_0-1}{w}_{nj}A\left(t_{n-j}\right)(N-S_j)+\sum _{j=n_0}^n{\omega }_{n-j}A\left(t_{n-j}\right)(N-S_j)\right),n=n_0,\dots ,$$

(15)

then

$$\begin{array}{cc}\hfill log\left(\frac{S_{n-1}}{S_n}\right)=& \beta h\left(\sum _{j=n_0+1}^n{\omega }_{n-j}A\left(t_{n-j}\right)(S_{j-1}-S_j)+{\omega }_{n-n_0}A\left(t_{n-n_0}\right)(N-S_{n_0})\right)\hfill \\ & +\beta h\sum _{j=0}^{n_0-1}\left(w_{nj}A\left(t_{n-j}\right)-w_{n-1j}A\left(t_{n-j-1}\right)\right)(N-S_j).\hfill \end{array}$$

Summing for n ranging from $n_0+1$ to $+\infty ,$ we obtain

$$\begin{array}{cc}\hfill & log\left(\frac{S_{n_0}}{S_{\infty }\left(h\right)}\right)=h(N-S_{n_0})\beta \sum _{n=1}^{+\infty }{\omega }_{n}A\left(t_n\right)+h\beta \sum _{n=n_0+1}^{+\infty }\sum _{j=n_0+1}^n{\omega }_{n-j}A\left(t_{n-j}\right)(S_{j-1}-S_j)\hfill \\ \hfill & +h\beta \sum _{j=0}^{n_0-1}\left(\sum _{n=n_0+1}^{+\infty }(w_{nj}A\left(t_{n-j}\right)-w_{n-1,j}A\left(t_{n-1-j}\right))\right)(N-S_j)\hfill \\ \hfill & =h\beta \left((N-S_{\infty }\left(h\right))\sum _{n=0}^{+\infty }{\omega }_{n}A\left(t_n\right)-(N-S_{n_0}){\omega }_{0}A\left(t_0\right)-\sum _{j=0}^{n_0-1}{w}_{n_{0}j}A\left(t_{n_0-j}\right)(N-S_j)\right).\hfill \end{array}$$

(16)

Consider Equation (15) with $n=n_0,$

$$log\left(\frac{S_0}{S_{n_0}}\right)=h\beta \left(\sum _{j=0}^{n_0-1}{w}_{n_{0}j}A\left(t_{n_0-j}\right)(N-S_j)+{\omega }_{0}A\left(0\right)(N-S_{n_0})\right).$$

Then, taking into account (16), it is

$$log\left(\frac{S_0}{S_{\infty }\left(h\right)}\right)=\left(N-S_{\infty }\left(h\right)\right)h\beta \sum _{n=0}^{+\infty }{\omega }_{n}A\left(t_n\right).$$

   □

Thus (14) can be interpreted as the discrete counterpart of (4) and, just like (4) (see [28]), it possesses a unique solution in $\left(0,S_0\right)$.

4. Trapezoidal Discretization in Some Realistic Cases

In this section, we consider particular choices for $A\left(t\right)$ in (5) which appear in applications (see, for example, [7], [27] (sec. 2.1), [11,33]). Indeed, we focus on cases where the infectivity function $A\left(t\right)$ is identically zero on an initial time interval,

$$A\left(t\right)=0,for t\le \tau ,with \tau \ge 0.$$

(17)

The condition $A\left(0\right)=0$ makes the numerical scheme (6) explicit, simplifying its analysis and reducing the computational cost at each step. We choose to approximate this type of model using the trapezoidal DQ method. In this case, the numerical solution not only agrees with the basic properties described in Theorem 1, without restrictions on the stepsize, but it is also monotone, and therefore retains the asymptotic dynamic of the analytical solution. The resulting explicit scheme has the form,

$$S_n=S_0·exp\left(-\frac{h}{2}\beta \left(A\left(t_n\right)(N-S_0)+2\sum _{j=1}^{n-1}A\left(t_{n-j}\right)(N-S_j)\right)\right),$$

(18)

for $n>0$, with $n_0=1,$ for which we state the following theorem.

Theorem 4.

Let ${\left\{S_n\right\}}_{n\ge 0}$ be the solution to the discrete Equation (18). Let $S_0>0$ then, for each $h>0,$ the sequence ${\left\{S_n\right\}}_{n\ge 0}$ is non-increasing.

Proof. 

Theorem 1 assures non-negativity and boundedness of the sequence. Here, we prove that it is also non-increasing by induction, starting from

$$S_1=S_{0}exp\left(-\frac{h}{2}\beta A\left(t_1\right)(N-S_0)\right)\le S_0.$$

Consider $n>1$ and assume that the assertion is true for each $j\le n-1.$ Since we have

$$\begin{array}{cc}\hfill log\left(\frac{S_{n-1}}{S_n}\right)=& \frac{h}{2}\beta \left(A\left(t_{n-1}\right)(N+S_0-2S_1)\right)\hfill \\ \hfill +& \frac{h}{2}\beta \left(2\sum _{j=2}^{n-1}A\left(t_{n-j}\right)(S_{j-1}-S_j)+A\left(t_n\right)(N-S_0)\right)\ge 0,\hfill \end{array}$$

it is $S_{n-1}\ge S_n$.    □

Theorem 4 assures that the solution ${\left\{S_n\left(h\right)\right\}}_{n\ge 0}$ to the discrete Equation (18) is non-increasing for each positive value of h and then it admits a finite limit $S_{\infty }\left(h\right)$, as n goes to $+\infty$. Hence, condition (13) holds for each $h>0,$ and the solution of Equation (14) is indeed the numerical final size $S_{\infty }\left(h\right).$

For the sake of comparison, we consider the reformulation of problem (1) given in [15]

$$\begin{array}{cc}\hfill S\left(t\right)& =S_0-\beta {\int }_0^{t}S\left(s\right)\phi \left(s\right)ds,\hfill \\ \hfill \phi \left(t\right)& =(N-S_0)A\left(t\right)+\beta {\int }_0^{t}S\left(s\right)\phi \left(s\right)A(t-s)ds,\hfill \end{array}$$

(19)

and compare the behavior of the numerical solutions obtained by (18) and the following trapezoidal discretization of (19)

$$\begin{array}{cc}\hfill S_n& =S_0-h\beta \sum _{j=1}^{n-1}{S}_j{\phi }_j-\frac{h}{2}\beta S_n{\phi }_n,\hfill \\ \hfill {\phi }_n& =(N-S_0)A\left(t_n\right)+h\beta \sum _{j=1}^{n-1}A\left(t_{n-j}\right)S_j{\phi }_j,\hfill \end{array}$$

(20)

for $n>0.$ Both methods are convergent of order $2.$ However, in (20), positiveness and monotonicity of the solution can be proved only under the following constraint on the discretization stepsize, $h<2/(\beta N·supA(t\left)\right)$, leading to a severe limitation in case of large values of the population size $N.$

5. Numerical Experiments

In this section, we report some numerical examples in order to show experimentally the theoretical results proved in Section 3 and Section 4. For our experiments, we choose illustrative test equations and we use method (6). Our first example consists of problem (5) with

$$A\left(t\right)=\frac{1.5}{\sigma \sqrt{2\pi }}{e}^{-\frac{{(t-\mu )}^2}{2{\sigma }^2}},\sigma =0.6,\mu =0.4,N={10}^2,S_0=90,\beta ={10}^{-3},$$

(21)

for $t\in [0,1].$ Here, we have integrated problem (21) by DQ method with trapezoidal, first and second Gregory weights (the orders of convergence are 2, 3 and 4, respectively). Since $A\left(t\right)$ does not meet condition (17), the methods are implicit and require at each step the resolution of a nonlinear equation for which a fixed-point iteration process is used.

Table 1. Error values and experimental order of convergence for examples (5)–(21).

NSFD in [26]	h	$\parallel {\mathbf{E}}_{\mathbf{n}}{\parallel }_2$	Experimental Order
	${10}^{-1}$	$2.478\times {10}^{-2}$	-
	${10}^{-2}$	$2.322\times {10}^{-3}$	$1.028$
	${10}^{-3}$	$2.308\times {10}^{-4}$	$1.003$
	${10}^{-4}$	$2.307\times {10}^{-5}$	$1.000$
Quadrature rule in (6)	h	$\parallel {\mathbf{E}}_{\mathbf{n}}{\parallel }_{\mathbf{2}}$	Experimental order
Trapezoidal Rule	${10}^{-1}$	$3.341\times {10}^{-3}$	-
	${10}^{-2}$	$3.326\times {10}^{-5}$	$2.002$
	${10}^{-3}$	$3.326\times {10}^{-7}$	$2.000$
	${10}^{-4}$	$3.326\times {10}^{-9}$	$2.000$
First Gregory Rule	${10}^{-1}$	$4.377\times {10}^{-4}$	-
	${10}^{-2}$	$4.082\times {10}^{-7}$	$3.030$
	${10}^{-3}$	$4.052\times {10}^{-10}$	$3.003$
	${10}^{-4}$	$3.941\times {10}^{-13}$	$3.012$
Second Gregory Rule	${10}^{-1}$	$1.976\times {10}^{-4}$	-
	${10}^{-2}$	$6.785\times {10}^{-9}$	$4.464$
	${10}^{-3}$	$6.857\times {10}^{-13}$	$3.995$

and Figure 1 show the numerical error behavior for different values of the stepsize $h.$ The numerical errors

$$E_n=S_n-{\overline{S}}_n,$$

are computed against the reference solution ${\overline{S}}_n,$ obtained by second Gregory rule with $h={10}^{-6}$. For the sake of comparison, we also report the errors obtained by integrating (21) with the first order non-standard finite difference (NSFD) method used in [26]. It is clear that the experimental order of convergence of DQ methods is coherent with the theoretical one predicted in Theorem 2.

Figure 2 represents the work precision diagram of the number of function evaluations with respect to the accuracy of the numerical solutions by the four methods considered here. It shows better performances in higher order methods. Indeed, even though we are comparing an explicit method with implicit ones, we see that, at the same computational cost, the accuracy of DQ methods is higher than NSFD. In addition, for given accuracies, even low ones, the NSFD scheme requires much more effort.

In order to show the long time behavior of the numerical solution, we consider problem (5) with:

$$A\left(t\right)=2t{e}^{-t},T=80,N={10}^5,S_0=9\times {10}^4,\beta ={10}^{-5}.$$

(22)

In Figure 3, the behavior of the numerical solution, computed by (6) and second Gregory weights, with $h={10}^{-3}$, is reported. In this case, the theoretical value for the final state, obtained by solving the nonlinear Equation (4) by the MATLAB routine fzero (see [34]), is $S_{\infty }=1717$.

We introduce the truncated numerical basic reproduction number ${\tilde{R}}_0\left(h\right)$ and the truncated numerical final size ${\tilde{S}}_{\infty }(h;T)$ as follows

$$M=T/h,{\tilde{R}}_0\left(h\right)=\beta Nh\sum _{n=0}^M{\omega }_{n}A\left(t_n\right),{\tilde{S}}_{\infty }\left(h\right)=S_M,$$

(23)

and define the residual on the numerical final size relation (14) as

$$r(h;T)=\left|log\left(\frac{S_0}{{\tilde{S}}_{\infty }\left(h\right)}\right)-{\tilde{R}}_0\left(h\right)\left(1-\frac{{\tilde{S}}_{\infty }\left(h\right)}{N}\right)\right|.$$

In this case, since (17) holds, Theorem 4 guarantees the existence of ${lim}_{n\to +\infty }{S}_n$ when the DQ method has trapezoidal weights. What is more, from our tests we observe that, for fixed $h>0,$ the value of ${\tilde{S}}_{\infty }\left(h\right)$ converges to a finite limit as T grows for any choice of the quadrature weights. Here, we have used $T=80$ for computing ${\tilde{S}}_{\infty }\left(h\right),$ and the effectiveness of this choice is confirmed by the values of $r(h;T)$ in

Table 2. Long time solution for problems (5)–(22): numerical truncated $R_0$ and final size.

Quad. Rule in (6)	h	Rel. Err. on ${\tilde{\mathbf{R}}}_0(\mathbf{h};\mathbf{T})$	Rel. Err. on ${\tilde{{S}}}_{\mathit{\infty }}(\mathbf{h};\mathbf{T})$	$\mathbf{r}(\mathbf{h};\mathbf{T})$
Trap. Rule	${10}^{-1}$	$8.330\times {10}^{-4}$	$2.104\times {10}^{-3}$	$3.975\times {10}^{-14}$
	${10}^{-2}$	$8.333\times {10}^{-6}$	$2.103\times {10}^{-5}$	$3.863\times {10}^{-14}$
	${10}^{-3}$	$8.333\times {10}^{-8}$	$2.103\times {10}^{-7}$	$4.174\times {10}^{-14}$
I Greg. Rule	${10}^{-1}$	$7.889\times {10}^{-5}$	$1.991\times {10}^{-4}$	$3.819\times {10}^{-14}$
	${10}^{-2}$	$8.288\times {10}^{-8}$	$2.091\times {10}^{-7}$	$3.797\times {10}^{-14}$
	${10}^{-3}$	$8.329\times {10}^{-11}$	$2.102\times {10}^{-10}$	$3.886\times {10}^{-14}$
II Greg. Rule	${10}^{-1}$	$7.130\times {10}^{-6}$	$1.799\times {10}^{-5}$	$3.819\times {10}^{-14}$
	${10}^{-2}$	$7.834\times {10}^{-10}$	$1.977\times {10}^{-9}$	$3.819\times {10}^{-14}$
	${10}^{-3}$	$7.905\times {10}^{-14}$	$2.581\times {10}^{-13}$	$3.864\times {10}^{-14}$

. In particular, Table 2 shows, for different values of h, the relative errors for the approximation of $R_0$ and $S_{\infty }$ by (23), as well as the corresponding numerical final size residuals.

Numerical results represent a validation of Theorem 3 but also give evidence of unproven theoretical properties. Indeed, it is clear from Table 2 that, for $h\to 0$, the numerical final size converges to the continuous one, with the same order $p\ge 1$ of the quadrature rule employed in (6).

Here, we consider, referring to [33], a kernel function suitable for modeling influenza

$$\begin{array}{cc}\hfill & t_1=0.4\mathrm{days},\hfill \\ \hfill & t_2=2.8\mathrm{days},\hfill \\ \hfill & t_3=3.8\mathrm{days},\hfill \\ \hfill & t_4=7.4\mathrm{days},\hfill \end{array}A\left(t\right)=\left\{\begin{array}{cc}\frac{1}{T_I}\frac{t-t_1}{t_2-t_1}\hfill & if{t}_1\le t\le t_2,\hfill \\ \frac{1}{T_I}\hfill & if{t}_2\le t\le t_3,\hfill \\ \frac{1}{T_I}\frac{t_4-t}{t_4-t_3}\hfill & if{t}_3\le t\le t_4,\hfill \\ 0\hfill & elsewhere,\hfill \end{array}\right.$$

(24)

based on the assumptions that no member of the population is infectious before $t_1$ days or after $t_4$ days post-exposure and maximum infectivity occurs between $t_2$ and $t_3$ days after exposure. The latent period $T_L$ and the infectious period $T_I$ are

$$T_L=\frac{t_1+t_2}{2}=1.6\mathrm{days},T_I=\frac{t_4+t_3-t_2-t_1}{2}=4.0\mathrm{days}.$$

We integrate problems (5)–(24) with

$$T=10,N={10}^5,S_0=9\times {10}^4,\beta ={10}^{-2},$$

by methods (18) and (20), respectively. The results, reported in Figure 4, have been obtained with two different values of the stepsize. We remind that, according to Theorem 4, the numerical solution computed by the trapezoidal scheme (18) is positive and non-increasing, regardless of the stepsize h. Conversely, we observe in Figure 4 that the numerical solution computed by means of the direct trapezoidal scheme (20), with $h=0.5,$ is increasing or even negative in some points.

6. Conclusions

We have introduced a numerical method specially designed for age-of-infection epidemic models represented by integro-differential equations of the type (1), reformulated as an integral equation in exponential form. The proposed method discretizes the integral term by DQ methods with Gregory convolution weights, from which it inherits the order of convergence. We show that the numerical solution is positive and bounded, which are crucial properties taking into account the biophysical meaning of variables. In some cases of realistic infectivity functions (for example, incorporating a non-zero finite incubation period and a finite infectious period), we prove that the proposed scheme, based on the trapezoidal method, also retains the monotonicity of the solution. In general situations, we provide numerical evidence of this behavior. A comprehensive analysis allows us to state that the numerical solution is dynamically consistent with the continuous one. We introduce a discretization of the basic reproduction number and prove that the numerical final size satisfies a nonlinear equation which is the discrete counterpart of the continuous final size relation. Numerical simulations show convergence, with respect to h, of the numerical final size to the continuous one, at a rate that reflects the order of the employed quadrature rule. A numerical comparison clearly shows that the numerical approach described here outperforms the first order dynamic-preserving non-standard explicit scheme introduced in [26]. This fact is advantageous in long-time integrations.