Renewable energy use and the renewable energy sector’s development: public finance, environmental externalities and sustainability

1.   Introduction

The analysis of the dynamical behaviour of numerical methods for differential and integral equation epidemic models has been a topic of interest ^{[1,2,3,4,5,6,7,8]}, because it satisfies the need to have robust numerical methods that share the same properties of the analytical solution. While most of the existing papers are devoted to the simulation of specific models, here we consider the following general problem, which consists of a Volterra integro-differential system of $2M+1$ equations of the type

where $t\geq 0,$ $V_i(t) = \sum_{r = 1}^M\beta_{ir}\varphi_r(t)+\alpha_iP(t),$ and $i = 1, \ldots, M.$ Here $\varphi_{i0}(t),$ $A_i(t),$ $i = 1, \dots, M,$ $P_0(t)$ and $B(t)$ are given continuous functions and $\alpha_i, \beta_i, c_i, \beta_{ir}\geq 0,$ $i, r = 1, \dots, M,$ are given constants. At time $t = 0$ the value of $S_i(t),$ $i = 1, \ldots, M,$ is $S_i^0,$ given. The motivation for considering system (1.1) is that it represents a general framework that includes a variety of epidemic mathematical models that, taking into account the age of infection, involve memory terms. We report some examples where we give the biological definitions of the variables and parameters.

I. An age-of-infection epidemic model is described in [9,p.139] by the following system

$S(t)$ is the number of susceptibles and $\varphi(t)$ is the total infectivity at time $t,$ $\beta > 0$ is the rate of effective contacts and $A(\tau)$ is the mean infectivity of members of population with infection age $\tau.$ Moreover, $\varphi_0(t)$ represents the infectivity at time $t$ of people who were infected before the initial outbreak. The general model (1.1) reduces to (1.2), with $M = 1,$ $\alpha_1 = 0,$ $\beta_{11} = 1,$ $\beta_1 = \beta,$ $c_1 = 0,$ and $B(t) = P_0(t) = 0,$ for $t\geq0.$

II. A model with both symptomatic and asymptomatic infections is described in ^[10], by the following system

Model (1.3) traces the evolution of the epidemics by distinguishing the infectivity functions of symptomatic and asymptomatic people, being $\varphi^s(t)$ the former and $\varphi^a(t)$ the latter. Thus, an analogous distinction is made for the known functions $\varphi_0^s(t)$ and $\varphi_0^a(t),$ which describe the contribution to the total infectivity of people who got infected before the initial time. Here, $N > 0$ is the total size of the population, $a > 0$ is the average number of contacts made by a member per unit of time and $f(t)\in(0, 1)$ is the probability for an individual to become symptomatic after the infection. Furthermore, $A^s(\tau)$ and $A^a(\tau)$ represent the mean infectivity of symptomatic and asymptomatic individuals with infection age $\tau,$ respectively. System (1.3) corresponds to the general model (1.1) with $M = 2,$ $\alpha_i = 0,$ $\beta_{ij} = 1,$ $\beta_i = \frac{a}{N},$ $c_i = 0,$ for $i, j = 1, 2$ and $S_1(t) = S_2(t) = S(t),$ $\varphi_1(t) = \varphi^s(t),$ $\varphi_2(t) = \varphi^a(t),$ $A_1(t) = f(t)A^s(t),$ $A_2(t) = (1-f(t))A^a(t)$ and $B(t) = P_0(t) = 0,$ for $t\geq 0.$

III. An age of infection epidemic model in a multi-group heterogeneous populations is proposed in ^[11]. The system, when the state space is discrete and the population is divided into $M$ subgroups of sizes $N_1, \ldots, N_M,$ reads

$i = 1, \ldots, M.$ Here, in group $i,$ $S_i(t)$ represents the number of susceptible members and $\varphi_i(t)$ is the infectivity at time $t,$ $\varphi_{i0}(t)$ is the infectivity at time $t$ of members who were infected before time $0,$ and $A_i(\tau)$ is the mean infectivity of individuals of group $i$ with infection age $\tau.$ Furthermore, $a_i\geq 0$ is the contact rate for members of group $i$, and $p_{ij}$ is the fraction of contacts made by a member of group $i$ with a member of group $j.$ Referring to the notations in ^[11,12], the general model (1.1) reduces to (1.4), with $\alpha_i = 0,$ $\beta_{ij} = \frac{p_{ij}}{N_j},$ $\beta_i = a_i,$ $c_i = 0,$ for $i, j = 1, \dots, M$ and $B(t) = P_0(t) = 0,$ for $t\geq 0.$

IV. In ^[13] the virus shedding epidemic model is studied

$i = 1, 2,$ where $P(t)$ is the pathogen shed by infected individuals of each group at a rate $r_1$ and $r_2,$ respectively. Here, $\beta_i,$ $i = 1, 2,$ are the contact rates, $A_1(\tau)$ and $A_2(\tau)$ are the mean infectivity of individuals in group $1$ and $2$ at age of infection $\tau$ and $\Gamma(\tau)$ is the fraction of pathogen remaining $\tau$ time units after having been shed by an infectious individual (see also [14,p. 168] in case of homogeneous mixing). The general model (1.1) reduces to (1.5), with $M = 2,$ $\alpha_i = 1,$ $\beta_{ij} = 0,$ $c_i = r_i,$ $i, j = 1, 2$ and $B(t) = \Gamma(t),$ for $t\geq 0.$

Since numerical simulations are of fundamental importance in the process of understanding the dynamic and the asymptotic behaviour of the system, we focus on the construction of a numerical model that preserves the global properties of the continuous problem. In the literature, there have been several investigations looking at numerical approximations of integral and integro-differential equations whose solutions remain bounded or converge to zero, including ^{[15,16,17,18,19,20,21]}. Here we concentrate on the general system (1.1) where, under suitable assumptions on the known parameters, $S_i(t)$ tends, as $t\to +\infty,$ to a limit $S_i(\infty),$ which is the solution of a known nonlinear limiting equation. In case of the epidemic models listed above and in many other cases (see for example ^[22,23] and references therein), $S_i(\infty)$ is the final size, one of the most relevant parameters in the description of the epidemic. We focus on a numerical method for (1.1) and we are interested to prove that, for any positive value of the size $h$ of the discretization, the limit of the numerical approximation, $S_i^{\infty}(h),$ exists and turns out to be the solution of an analogous nonlinear limiting equation. Furthermore, we show that $S_i^{\infty}(h)$ tends to $S_i(\infty),$ as $h$ vanishes. The convergence as $h\to 0$ is an obvious property that a numerical method must satisfy while integrating over a limited range, but it is not at all guaranteed and, in general, difficult to prove for the asymptotic solution.

Then, an important feature of our investigation is that it directly shows, under suitable conditions on the integro-differential system, that the discrete model has properties which ensure its solutions to behave asymptotically as the continuous ones. The numerical investigation is performed on a general discrete system obtained approximating (1.1) by non-standard finite differences. Then, the stability analysis, consisting in a comparative study of the behaviour of the analytical and the numerical solutions, is carried out in significant cases (that is (1.2)–(1.5)) and can be extended to other models, that fall into the form (1.1).

In this viewpoint, we say that the numerical scheme is coherent with the continuous problem and system (1.1) can be considered as a test problem for proving such coherence.

Here we assume that the given functions in (1.1) belong to $C[0, +\infty)\cap L^1[0, +\infty)$, and that:

Assumptions A for $i = 1, \ldots, M:$

● $\alpha_i\geq 0, \; \beta_{ir}\geq 0, \; r = 1, \ldots, M,$

● $S_i^0 = S_i(0) > 0, \; \varphi_{i0}(t)\geq 0, \; A_i(t)\geq 0, \; t \geq 0,$

and

● $P_0(t)\geq 0, \; \; B(t)\geq 0, \; \; t\geq 0.$

Assumptions B there exist positive constants $\varphi_{0, max},$ $P_{0, max},$ $A_{max},$ and $\bar{B}$ such that:

● $\varphi_0(t)\leq \varphi_{0, max}, \; P_0(t)\leq P_{0, max},$ $t\geq 0,$

● $A_i(t)\leq A_{max},$ $i = 1, \ldots, M,$ $t\geq 0,$

● $h\sum_{n = 0}^{\infty}B(nh)\leq \bar{B},$ $h > 0.$

Assumptions C there exist positive constants $\bar{A},$ $\bar{\varphi}_0,$ $\bar{P}_0$ such that for $h > 0$:

● $h\sum_{n = 0}^{+\infty}A_i(nh)\leq \bar{A}, \; i = 1, \ldots, M,$

● $h\sum_{n = 0}^{+\infty}\varphi_{0i}(nh)\leq \bar{\varphi}_0, \; i = 1, \ldots, M,$

● $h\sum_{n = 0}^{+\infty}P_{0}(nh)\leq \bar{P}_0.$

Observe that the third of Assumptions B and the Assumptions C have the form $h\sum_{n = 0}^{+\infty}Q(nh)\leq \bar{Q},$ $h > 0,$ with $Q\in L^1[0, +\infty)$ and $\bar{Q} > 0,$ which is certainly accomplished when, for example, the function $Q^{\prime}(t)\in L^1[0, +\infty)$ (see ^[19]) or when (usually true in realistic situations) $Q(t)$ is definitely non-increasing.

A generalization to the results in ^[9,13,14] implies that:

● under the Assumptions A, $S_i(t)$ is positive and non-increasing, $\varphi_i(t)\geq 0,$ for $i = 1, \dots, M,$ and $P(t)\geq 0;$

● if in addition Assumptions B hold, then $S_i(t),$ $\varphi_i(t),$ $i = 1, \dots, M$ and $P(t)$ are bounded;

● $\lim_{t \to +\infty}\varphi_i(t) = 0$ and $\lim_{t \to + \infty}S_i(t) = S_i(\infty),$ where $S_i(\infty)$ satisfies the limiting algebraic system $R_i(x) = 0,$ $i = 1, \dots, M,$ where

with $x = (x_1, \ldots, x_M),$ if a solution to (1.6) exists.

The paper is organized as follows. In the next section we present a numerical method based on a nonlocal discretization to solve the general system (1.1) and we prove the basic non-negativity and boundedness properties of the numerical solution under the Assumptions A and B on the known functions. These results hold for any value of the size $h > 0$ of the discretization. In Section 3 we give some preliminary results that are used in Section 4, where the numerical asymptotic solution is shown to be the root of a nonlinear system of algebraic equations. Here we prove that this solution is unique and that it converges to its analytical counterpart as the stepsize $h\to 0.$ In Section 5 we present some cases of interest in epidemic models and apply the theory developed in the previous sections in order to study the coherence of the numerical method with the epidemic models. Furthermore, we show the results of some simulations performed. Finally, in Section 6 we conclude the paper with some remarks.

2.   Numerical method

Define a stepsize $h > 0,$ and an uniform mesh $\{t_n = nh, \; n = 0, 1, \ldots\}.$ Our numerical approach to solve the nonlinear system (1.1) is the following finite difference scheme

$n = 0, 1, \ldots,$ where $V_i^n = \sum_{r = 1}^M\beta_{ir}\varphi_r^n+\alpha_iP^n,$ and $i = 1, \ldots, M.$ Here $S_i^0 = S_i(0),$ $\varphi_{i}^0 = \varphi_{i0}(0),$ and $P^0 = P_{0}(0)$ are given. Furthermore, $S_i^{n}\approx S_i(t_n),$ $\varphi_i^{n}\approx \varphi_i(t_n),$ for $i = 1, \ldots, M,$ and $P^{n}\approx P(t_n),$ for $n = 0, 1, \ldots.$ We say that the method (2.1) is non-standard since the nonlinear terms in (1.1) are approximated in a nonlocal way. A pseudocode implementation of the numerical method (2.1) is reported with the Algorithm 1.

The following results, concerning non-negativity and boundedness, hold.

Theorem 2.1. Consider equation (2.1) with Assumptions A, then the solution sequences $\{S_i^n\}_{n\in\mathbb{N}_0},$ $\{\varphi_i^n\}_{n\in\mathbb{N}_0},$ $i = 1, \ldots, M$ and $\{P^n\}_{n\in\mathbb{N}_0},$ are non-negative, $\forall h > 0.$ Furthermore, the sequence $\{S_i^n\}_{n\in\mathbb{N}_0},$ $i = 1, \ldots, M,$ is positive and non-increasing.

Proof. We proceed by induction to prove that the statement $S_i^{n} > 0$, $\varphi_i^{n}\geq 0,$ $P^n\geq 0,$ holds for all $n\in \mathbb{N}_0,$ $h > 0$ and $1\leq i \leq M.$ The case $n = 0$ is true because the initial values are non-negative and $V_i^0 = \sum_{r = 1}^M\beta_{ir}\varphi_r^0+\alpha_iP^0\geq 0,$ $i = 1, \dots, M.$ Consider $n\geq1$ and assume that the properties are true for each $0\leq j \leq n-1.$ It follows that $V_i^j\geq 0,$ for $0\leq j \leq n-1,$ therefore it is

and then, from (2.1), also $\varphi_i^{n}\geq0,$ $i = 1, \dots, M$ and $P^{n}\geq 0.$ Furthermore, from (2.2), $S_i^n\leq S_i^{n-1},$ for all $n\geq 1$ and $i = 1, \dots, M,$ which completes the proof. □

From the first equation of (2.1) it is

and this leads to the following equation for the asymptotic numerical solution $S_i^{\infty}(h),$

where $S_i^{\infty}(h) = \lim_{n\to\infty}S_i^n,$ for any fixed stepsize $h > 0,$ which exists since the sequence $\{S_i^n\}_{n = 0}^{+\infty}$ is non-increasing, for each $1\leq i \leq M.$

Theorem 2.2. Consider equation (2.1) with Assumptions A and B, then the solution sequences $\{S_i^n\}_{n\in\mathbb{N}_0},$ $\{\varphi_i^n\}_{n\in\mathbb{N}_0},$ $i = 1, \ldots, M$ and $\{P^n\}_{n\in\mathbb{N}_0},$ are bounded, $\forall h > 0.$

Proof. From Theorem 2.1, $\{S_i^n\}_{n\in \mathbb{N}_0}$ is a non-increasing sequence, thus it is bounded from above by $S_i^0,$ $i = 1, \dots, M.$ Consider $\varphi_i^n,$ from the first equation and the second equation of (2.1), it is

Finally, the inequality

directly follows from the assumptions and the third equation in (2.1). □

As already pointed out in the previous section, the time-continuous solution to (1.1) is non-negative and bounded and the proposed numerical solution algorithm preserves boundedness and nonnegativity unconditionally with respect to time step size.

In order to study the convergence of the numerical method, we assume that the known functions are continuously differentiable on $[0, T]$, with $T < +\infty,$ and we investigate the behaviour of the local truncation error (see ^[24] for the definition). Here, we consider the case of $M = 1,$ since the generalization to $M > 1$ is straightforward. In this case, the local truncation error of the discretization in (2.1) reads, for $T = \bar{n}h$ and $n = 0, \dots, \bar{n},$

By the mean value theorem

with $\theta_j\in(0, 1),$ $j = 0, \dots, n-1.$ Moreover, due to the convergence properties of the rectangular quadrature rule (see, for instance, ^[25]), the bound

holds for each $0\leq j \leq n-1,$ with $C > 0,$ independent of $h.$ It then follows from (2.4) that

for each $n = 0, \dots, \bar{n}.$ Hence

being $\tilde{C}$ a positive constant not depending on $h.$ Denote by

the global error of the discretization (2.1). The following theorem, that we prove by standard techniques, provides sufficient conditions for the convergence of the numerical method.

Theorem 2.3. Assume that the given functions $A_i(t),$ $i = 1, \ldots, M$ and $B(t),$ describing problem (1.1), are continuously differentiable on an interval $[0, T]$ and that $\{S_i^n\}_{n \in \mathbb{N}_0},$ $\{\varphi_i^n\}_{n \in \mathbb{N}_0},$ $\{P^n\}_{n \in \mathbb{N}_0}$ are the approximations to the solution of (1.1), defined by (2.1). Let $h = T/\bar{n},$ with $\bar{n}$ positive integer, then

Furthermore, the order of convergence is $1.$

Proof. We prove the result for $M = 1,$ since the generalization to $M > 1$ is straightforward. In this case, the time-continuous system (1.1) reads

and subtracting (2.1) from it leads to the global error

for $n = 0, \dots, \bar{n}.$ Furthermore, for each $j = 0, \dots, n-1,$

hence, from (2.6),

with $K$ and $\tilde{K}$ positive constants depending on the parameters of the problem but not on $h$. Therefore, for a sufficiently small $h,$

Finally, the Gronwall discrete inequality yields (see, for instance [24,p.101]),

for $n = 0, \dots, \bar{n}.$ Then, from (2.5), the result follows. □

3.   Preliminary results

In this section we prove some preliminary results needed to describe the asymptotic behaviour of the solution to system (2.1).

Lemma 3.1. Consider system (2.1), for $i = 1, \ldots, M,$ it is

Proof. Summing from $0$ to $\infty$ in the last equation of (2.1), interchanging the order of summation and adding to both members $P(0) = P^0,$ it is

The same can be done for the second equation of (2.1), where taking into account that, from the first of (2.1) it is $h\beta_r\sum_{j = 0}^{+\infty}S_r^{j+1}V_r^j = S_r^0-S_r^{\infty}(h),$ we have

for any $r = 1, \ldots, M.$ Combining the previous expressions with

we get the result. □

Theorem 3.1. Consider equation (2.1) with Assumptions A, B and C, then there exists $0 < \bar V < +\infty,$ such that $h\sum_{n = 0}^{+\infty}V_r^n < \bar V.$

Proof. From Lemma 3.1, for $i = 1, \ldots, M,$ it is

□

4.   Asymptotic behavior of the numerical solution

The results of Section 2 ensure that the approximation of the solution to (1.1), obtained by (2.1), unconditionally retains the basic properties of the model. Here, we investigate the asymptotic properties and prove that the limit of the numerical solution behaves exactly as its continuous counterpart.

First of all, we observe that from Assumptions C and (3.1), for each positive $h,$

Define, for $x = (x_1, \dots, x_M),$ $h > 0$ and $i = 1, \ldots, M,$

Here, $U_i(h)$ is given by

We observe that, due to (4.3), the relation (2.3) is equivalent to

Hence, by substituting in it the expression of $\sum_{n = 0}^{+\infty}V_i^n$ from Lemma 3.1, it is clear that for any fixed $h > 0,$ the asymptotic numerical solution $S_i^{\infty}(h),$ $i = 1, \ldots, M,$ is a root of the nonlinear system of equations

if a solution to (4.4) exists.

In order to show that $S_i^{\infty}(h)$ tends to $S_i(\infty),$ as $h\to 0,$ we need the following result.

Lemma 4.1. Consider equation (2.1) with Assumptions A, B and C, then, for $U_i(h)$ defined in (4.3), it is

Proof. From its definition and Lemma 3.1, $V_i^n,$ $i = 1, \dots, M,$ $n\geq 0,$ is bounded by a constant independent of $h.$ This implies that $\lim_{h\to 0}hV_i^n = 0,$ uniformly with respect to $h.$ In this situation, we are allowed to proceed as in the proof of [19,Theorem 4.3] to get (4.5). □

Now we want to prove that system (4.4) has a unique solution for $h > 0,$ so we consider a general nonlinear algebraic system,

for which the following result holds.

Theorem 4.1. Denote by $D = \prod_{i = 1}^MD_i,$ with $D_i = (a_i, b_i],$ $i = 1, \dots, M.$ Let $\Gamma_i:D\to \mathbb R,$ $i = 1, \dots, M,$ be a twice continuously differentiable function and assume that, for each $i = 1, \dots, M$ and any fixed $(x_1, \dots, x_{i-1}, x_{i+1}, \dots, x_M)\in \prod\limits_{j = 1, j\neq i}^M D_j,$

a) $\Gamma_i(x_1, \dots, x_{i-1}, \zeta, x_{i+1}, \dots, x_M)$ admits at least one zero in $D_i;$

b) $\lim_{\zeta \to a_i^+}\Gamma_i(\dots, x_{i-1}, \zeta, x_{i+1}, \dots) > 0$ and $\Gamma_i(\dots, x_{i-1}, b_i, x_{i+1}, \dots) < 0;$

c) $\partial_{x_j} \Gamma_i(x) > 0,$ for all $x\in D,$ $i, j = 1, \dots, M$ and $i\neq j;$

d) $\partial_{x_jx_k}^2 \Gamma_i(x)\geq0,$ for all $x\in D$ and $j, k = 1, \dots, M.$

Then, system (4.6) has a unique solution in $D.$

Proof. We proceed by induction on $M.$ Choose $M = 2.$ Let $x_1\in D_1,$ then from a), b) and d), $\Gamma_2(x_1, \zeta)$ has a unique zero $\xi_2 = \xi_2(x_1)\in D_2,$ with $\partial_{x_2}\Gamma_2(x_1, \xi_2) < 0.$ Therefore, taking in account the arbitrariness of $x_1$ in $D_1,$ from the implicit function theorem, $\forall x \in D_1,$ there exists a unique $z(x)\in D_2$ such that $u_2(x): = \Gamma_2(x, z(x))\equiv0.$ Since $u_2^\prime(x)\equiv0$ and $u_2^{\prime\prime}(x)\equiv0,$ using assumption c) and the fact that a), b) and d) imply $\partial_{x_2}\Gamma_2(x, z(x)) < 0,$ we conclude that $z^\prime(x) > 0$ and $z^{\prime\prime}(x)\geq0.$ Now, we exploit the function $z(x)$ to build a solution to system (4.6) with $M = 2.$\\ For an arbitrary $\alpha_0\in D_1,$ let $\alpha_1\in D_1$ be the unique root of the function $\Gamma_1(\zeta, z(\alpha_0)).$ If $\alpha_1 = \alpha_0,$ then $(\alpha_0, z(\alpha_0))$ is a solution to (4.6) with $M = 2.$ Otherwise we can suppose, with no loss of generality, that $\alpha_0 < \alpha_1.$ It follows that $z(\alpha_0) < z(\alpha_1)$ and, from c), $0 = \Gamma_1(\alpha_1, z(\alpha_0)) < \Gamma_1(\alpha_1, z(\alpha_1)).$ Thus $\Gamma_1(\zeta, z(\alpha_1))$ has a unique zero $\alpha_2\in D_1,$ for which $\alpha_0 < \alpha_1 < \alpha_2.$ Similar arguments lead to an increasing sequence $\{\alpha_{n}\}_{n\in \mathbb{N}}\subset D_1,$ such that

Hence, $(\alpha, z(\alpha))\in D$ is a solution to system (4.6), with $M = 2.$

To prove its uniqueness, we consider another solution $(\beta, \gamma)\in D$ and define the function $u_1: x\in D_1 \to \Gamma_1(x, z(x))$. Since $\Gamma_2(\beta, \gamma) = 0,$ it follows that $\gamma = z(\beta)$ and $u_1(\beta) = 0 = u_1(\alpha).$ Since $u_1$ is convex and, from b), it admits a unique zero $\beta = \alpha,$ then the solution to (4.6) for $M = 2$ is unique.

Now assume that the result holds for any $M-1$ dimensional system satisfying a)–d). Consider $M > 2,$ proceeding as in the previous case, for each $(x_1, \dots, x_{M-1})\in \prod_{j = 1}^{M-1}D_j,$ there exists a unique function $z(x_1, \dots, x_{M-1})$ such that

Define, for $i = 1, \dots, M-1,$ the functions

Now we want to prove that assumptions a)-d) are true for the $M-1$ dimensional system

which is equivalent to (4.6). Regarding the assumption a), we need to prove that, for each fixed $(x_1, \dots, x_{i-1}, x_{i+1}, \dots, x_{M-1})\in \prod\limits_{j = 1, j\neq i}^{M-1}D_j$, the function

has a zero in $D_i,$ $i = 1, \dots, M-1.$ For the theorem assumptions, given $\alpha_0\in D_i,$ there exists a unique $\alpha_1\in D_i,$ such that

If $\alpha_0 = \alpha_1,$ we get the result. Otherwise, proceeding as in the $M = 2$ case, we construct a monotone sequence $\{\alpha_n\}_{n\in \mathbb{N}}$ such that $\alpha = \lim_{n\to +\infty}\alpha_n\in D_i$ and $u_i(x_1, \dots, x_{i-1}, \alpha, x_{i-1}, \dots, x_{M-1}) = 0$ holds. Furthermore, b), c) and d) immediately follow from the hypotheses of the theorem, from (4.7) and from (4.8). Therefore, induction hypotheses assure that there exists a unique $\lambda\in \prod_{i = 1}^{M-1}D_i,$ root of the system (4.9). It follows that $(\lambda, z(\lambda))$ is a solution to the original system (4.6).

Finally, if $(\eta, \theta)$ is another solution of (4.6), with $\eta\in \prod_{i = 1}^{M-1}D_i$ and $\theta \in D_M,$ then $\theta = z(\eta).$ Thus $\eta$ solves (4.9) and the uniqueness of its solution, that we have just proved, yields $\eta = \lambda,$ which completes the proof. □

Consider the nonlinear system (4.4) for $h > 0,$ with $R_i(x, h)$ defined in (4.2). The assumptions of Theorem 4.1 apply to this system with $\Gamma_i(x) = R_i(x, h)$ and $D_i = (0, S_i^0],$ $i = 1, \dots, M.$ As a matter of fact we analyze, for $h > 0,$ $i, j = 1, \dots, M$ and $i\neq j,$ the twice continuously differentiable function $R_i(x_1, \dots, x_{i-1}, \zeta, x_{i+1}, \dots, x_M, h),$ with $0 < \zeta \leq S_i^0$ and $0 < x_j\leq S_j^0$ fixed. Since

there exists $0 < S_i^\infty(h)\leq S_i^0,$ such that $R_i(x_1, \dots, x_{i-1}, S^\infty_i(h), x_{i+1}, \dots, x_M, h) = 0$. Furthermore, for each $i, j = 1, \dots, M,$ $i\neq j,$ it is

Thus, according to Theorem 4.1, system (4.4) has, for each $h > 0,$ a unique solution $(S^\infty_1(h), \ldots, S^\infty_M(h)),$ with $0 < S^\infty_i(h)\leq S_i^0,$ $i = 1, \dots, M,$ which is the asymptotic numerical solution.

Regarding the time-continuous model, we consider the nonlinear system

with $R_i(x)$ defined in (1.6). It can be easily seen that the assumptions of Theorem 4.1 are accomplished and then the asymptotic analytical solution $S_i(\infty),$ $i = 1, \dots, M,$ to (1.1) is the unique solution of (4.10), with $0 < S_i(\infty)\leq S_i^0,$ $i = 1, \dots, M.$

In order to investigate the relation between the asymptotic properties of the numerical solution as $h\to 0$ and of the time-continuous solution, we assume that

for each $h > 0$ and $i = 1, \dots, M.$ The conditions in (4.11) are true, for instance, if the involved functions are ultimately non-increasing, or if their derivatives belong to $L^1[0, +\infty)$ (see ^[19]).

As the stepsize $h$ vanishes, we expect that the asymptotic numerical solution converges to the continuous one. In fact, we apply the following result.

Theorem 4.2. Consider a bounded subset $D$ of $\mathbb R^M,$ and a function

satisfying:

a) $\Phi(w, h),$ continuous for each $(w, h)\in D\times [0, +\infty);$

b) equation $\Phi(w, h) = 0,$ has at least one solution $\bar w(h)\in D,$ $\forall h\in[0, +\infty);$

c) equation $\Phi(w, h) = 0$ has a unique solution for $h = 0,$ namely $\bar w = \bar w(0)\in D.$

Then $\lim_{h\to 0}\bar w(h) = \bar w.$

Proof. Straightforward extension of [26,Theorem 3.5] to the case of functions defined on open sets. □

The vector function $R(x, h) = (R_1(x, h), \ldots, R_M(x, h)),$ with $R_i$ defined in (4.2), satisfies all the assumptions of Theorem 4.2, in $D = (0, S_1^0]\times\ldots\times(0, S_M^0].$ Now, if (4.11) holds, the system corresponding to $h = 0$ is given by (4.10). Therefore, this result establishes a connection between the asymptotic numerical solution $S^{\infty}_i(h),$ $i = 1, \ldots, M$ and the solution $S_i(\infty),$ $i = 1, \ldots, M$ of the nonlinear system (4.10), thus emphasizing that the limit of $S^{\infty}_i(h),$ $i = 1, \ldots, M,$ for $h\to 0$ exists and

5.   Case studies and numerical experiments

As we have emphasized in Section 1, system (1.1) represents a general theoretical setting which includes a variety of epidemic models in the literature. In this context, the nonlinear system (4.10) is the final size relation for the epidemic and system (4.4) represents the discrete final size relation corresponding to the numerical solution. From this point of view, Sections 3 and 4 analyze how the qualitative properties of the model are preserved when the system is integrated by the numerical scheme (2.1). We focus here on four cases of interest, which act as test cases in our analysis. Due to Theorems 2.1 and 2.2, the numerical solution remains non-negative and bounded for any value of the stepsize $h > 0,$ this guarantees that the epidemic is correctly simulated by method (2.1). The asymptotic analysis of the previous section applies to show that also the long time behaviour is preserved since the numerical final size converges, as $h\to 0,$ to the final size of the epidemic. This is also clear in the figures, which represent the results of the numerical experiments for each test model considered. For our experiments we choose illustrative test equations and we use the non-standard method (2.1).

● The age-of-infection epidemic model described in [9,p. 139] by the system (1.2). Here, in case $\varphi_0(t) = (N-S^0)A(t),$ being $N$ the constant size of the population, the final size of the epidemic $S(\infty)$ is the unique root of the nonlinear equation (see for example ^[9])

where $R_0 = \beta N\int_0^{+\infty}A(t)dt$ is the basic reproduction number. Given the identifications of the parameters in I, the nonlinear system (4.10) corresponds to system (5.1) (see ^[27,28] for a general discussion on the meaning and on the analytical and numerical computation of $R_0$). For this model the numerical method (2.1) reduces to the non-standard numerical scheme proposed in ^[19] and the nonlinear equation (4.4) represents the relation for the numerical final size of the epidemic $S^\infty(h)$ as obtained in that paper. In ^[19] it is proved that $\lim_{h\to 0}S^{\infty}(h) = S(\infty),$ provided that this limit exists. In Section 4, Theorem 4.2, we have proved that, in fact, this limit exists since, under the above mentioned assumptions on $A(t),$ (4.11) holds. So, Theorem 4.2 completes the analysis made in ^[19] and we can assert that the numerical method is asymptotically coherent with (1.2). Figure 1 shows the result when integrating problem (1.2) for $t\in[0, 100]$, when the infectivity function has a low regularity (see ^[29] for reference on this problem)

We have set the parameters $\beta = 10^{-4},$ $\tau_e = 12,$ $\tau_a = 14,$ $\tau_b = 16,$ $\tau_h = 19,$ and we have used $S^0 = 49950,$ $N = 50000,$ and stepsize $h = 0.1.$ We see that the numerical solution behaves coherently with the theoretical findings. The dot at the left end point of the integration interval is the value for $S(\infty)$ obtained by solving the nonlinear equation (5.1) through the Matlab routine fzero. This value is $S(\infty) = 148.83.$ The value $S^{\infty}(h)$ obtained by running the method (2.1) in the interval $[0, 1000],$ with $h = 10^{-3},$ is $149.97$ and the accuracy in the approximation improves linearly as the stepsize $h\to 0.$ Furthermore, in compliance with (4.1), the endpoint approximation of $\varphi$ is numerically zero. Other numerical tests are reported in ^[19].

● The model (1.3) with both symptomatic and asymptomatic infections in ^[10]. In this case the final size of the epidemic $S(\infty)$ is the root of the nonlinear equation

where $\mathcal{F} = a\int_{0}^{+\infty}f(t)A^s(t)\ dt+a\int_{0}^{+\infty}(1-f(t))A^a(t)\ dt,$ is considered to be finite. System (5.3) is equivalent to (4.10) with the parameters specified in II. Hence, the findings of Section 4 provide a uniqueness result for $S(\infty)$ as the only root in $(0, S^0]$ of (5.3). Furthermore, the convergence relation (4.12) confirms the numerical method (2.1) as a reliable tool to predict the asymptotic behaviour of the epidemic and the total number of symptomatic and asymptomatic patients. These considerations extend the investigation of ^[10]. As an example, we integrate problem (1.3) with

where

is the gamma probability density function. The number of symptomatic and asymptomatic individuals at time $t,$ $I^s(t)$ and $I^a(t)$, respectively, satisfy the system

Because of the absence of demographic turnover, the number of recovered people, at time $t,$ is $R(t) = N-(S(t)+I^a(t)+I^s(t)).$ Figure 2 shows the approximation of $S(t)$ by (1.1), as well as the approximations of $I^s(t)$, $I^a(t)$ and $R(t)$ computed, for $n\geq 0,$ by the same nonlocal technique employed in (2.1).

● The multi-group heterogeneous populations model (1.4) proposed in ^[11]. The nonlinear system for the final sizes $S_i(\infty),$ of the epidemic in group $i$, $i = 1, \ldots, M,$ is

System (4.10), with the parameter specifications in III, corresponds to the final size system (5.6). In Figure 3 we report the numerical solution to equation (1.4), when the population of size $N = 1000$ is divided in $2$ subgroups of sizes $N_1 = 0.1N$ and $N_2 = 0.9N,$ respectively. We consider $\varphi_{i0}(t) = (N_i-S_i^0)A_i(t),$ $i = 1, 2$ and we choose $S_1^0 = 99,$ $S_2^0 = 899,$ $a_1 = 5,$ $a_2 = 10,$ $p_{11} = 0.4,$ $p_{12} = 0.6,$ $p_{21} = 0.5,$ $p_{22} = 0.5,$ and

In order to check the asymptotic properties of the numerical solution, we have solved problem (1.4) by (2.1) on a large interval many times, each time halving the value of the stepsize $h.$ Then we have used the numerical solution at the end point of the integration interval to evaluate the expression (5.6). This procedure gives a measure for the errors in approximating $S_i(\infty),$ which are listed in Table 1, where it is clear that the numerical final size converges linearly to the corresponding continuous one.

● The virus shedding epidemic model (1.5) studied in ^[13]. Here we assume that $\varphi_{i0} = \int_{-\infty}^tA_i(t-s)S(s)P(s)ds,$ $i = 1, 2$ and $P_0(t) = \int_{-\infty}^t\Gamma(t-s)(r_1\varphi_1(s)+r_2\varphi_2(s))ds,$ are known functions respectively of the form $\varphi_{i0}(t) = (N_i-S_{i}^0)A_i(t),$ and $P_0(t) = P_0\Gamma(t).$ In this case the final size relation for $S_{i}(\infty),$ $i = 1, 2,$ is the following

$i = 1, 2.$ For each group, $N_i$ ($i = 1, 2$) are the sizes and

are the basic reproduction numbers. Moreover, under the conditions in IV, system (4.10) corresponds to the nonlinear system (5.8), which has a unique solution because of Theorem 4.1. The numerical results are shown in Figure 4 for

$N_1 = 200,$ $N_2 = 300,$ $S_{1}^0 = 199,$ $S_{2}^0 = 298,$ $P_0 = 2,$ $\beta_1 = 0.015,$ $\beta_2 = 0.03,$ $r_1 = 0.1,$ $r_2 = 1.$ Again a convergence of order one of the numerical final size $S_i^{\infty}(h),$ $i = 1, 2,$ to the one given by the model (1.5), can be experimentally observed.

6.   Concluding remarks

Problem (1.1) and the corresponding discrete system (2.1) represent a general framework to study and compare the qualitative behaviour of the analytical solution of some epidemic models, currently of interest in the scientific community, and its numerical approximation. The numerical solution is obtained by a non-standard discretization of the nonlinear terms in the system, and agrees with the analytical solution in many important qualitative aspects. Both the behaviour at finite time and the asymptotic properties of the solution are preserved for any value of the stepsize $h > 0.$ Furthermore, it is proved that the asymptotic numerical solution is the unique root of a nonlinear system, and that it converges to the analytical one as $h\to 0.$ The system of equations (1.1) includes, in its general form, some of the well known age-of-infection epidemic models in the literature. In general, the study carried out in this paper can be regarded as a stability numerical investigation on a class of epidemic problems that act as test equations. The role of the test equations here can be considered from two different points of view. First of all, since the proposed method preserves the qualitative peculiarities of the continuous solution, it can be considered reliable also for its quantitative evaluation and so the numerical solution can answer crucial questions such as the peak of the epidemic represented by (1.1). From another point of view, the aforementioned characteristics of the method make it reasonable to expect that it will behave well even on more complex problems, not included in (1.1), such as models with death for disease or with time-varying coefficients ^[30], where the theory needs simulation to describe the dynamics of the epidemic. The issue of time-varying coefficients represents an interesting point to investigate in a future work within the framework of the general model (1.1).

Acknowledgments

This research was supported by GNCS-INDAM.

Conflict of interest

The authors declare there is no conflict of interest.