A non-standard ﬁnite-di ﬀ erence-method for a non-autonomous epidemiological model: analysis, parameter identiﬁcation and applications

: In this work, we propose a new non-standard ﬁnite-di ﬀ erence-method for the numerical solution of the time-continuous non-autonomous susceptible-infected-recovered model. For our time-discrete numerical solution algorithm, we prove preservation of non-negativity and show that the unique time-discrete solution converges linearly towards the time-continuous unique solution. In addition to that, we introduce a parameter identiﬁcation algorithm for the susceptible-infected-recovered model. Finally, we provide two numerical examples to stress our theoretical ﬁndings


Introductory motivation and literature overview
Differential equations, especially those possessing non-negative solutions, play an important role in many areas of sciences such as biology [1][2][3][4][5], chemistry [6][7][8], epidemiology [9][10][11] or population dynamics [12][13][14] to name a few of these fields. Since we focus on examining a time-discrete, mathematical model of epidemiology in our work, we limit ourselves to a concise review of literature in this area.
In this article, we want to consider the time-continuous non-autonomous susceptible-infected-recovered (SIR) model S (0) = S 0 , I (0) = I 0 , in epidemiology thoroughly investigated by Wacker and Schlüter in 2020 [15]. We provide further assumptions regarding this time-continuous model in Section 2. In 1927, one of the most famous papers in mathematical epidemiology, written by Kermack and McKendrick, was published [16]. Previously, Ross and Hudson proposed important foundations for epidemiological models [17,18]. Additional well-written books and review articles on mathematical epidemiology are authored by Hethcote [19] and Brauer and Castillo-Chávez [20]. Since the establishment of the classical SIR model, many modifications have been suggested and, due to the on-going COVID-19 pandemic, many articles have been published with respect to this topic and we can only summarize a brief selection of these works.
To this day, there are works which provide closed solution formulas for simple models in epidemiology [21][22][23]. However, more articles consider numerical approximations due to appearing nonlinearities. Especially during the COVID-19 pandemic, different approaches have been published [15,24]. In 2022, Chen and co-authors suggested a susceptible-exposed-infected-unreported-recovered (SEIUR) model for the spread of COVID-19 with piecewise constant dynamical parameters [25] as one possible extension of classical SIR models. Taking age-and sex-structure of a population into account, Wacker and Schlüter proposed an extension of classical SIR models by multiple groups which can be considered as, for example, age groups or locally related groups [26]. Using fractional differential operators, Xu and co-authors analyzed the spread of COVID-19 in 2020 with respect to the United States by means of a fractional-order susceptible-exposed-infected-recovered (SEIR) model [27]. Regarding parameter identification problems in epidemiological models, Marinov and Marinova designed adaptive SIR models to solve these inverse problems [28,29]. In context of inverse problems, Comunian, Gaburro and Giudici concluded that high-quality data are needed to calibrate constant parameters of epidemiological models [30]. Other interesting lines of research concern the consideration of epidemics on networks [31][32][33][34][35] or application of deep-learning approaches to, for example, forecast time-series data of epidemics by these methods [36][37][38].

Article's motivation and contributions
One important aspect in epidemiological models is preservation of positivity. It plays a huge role in different scientific problems such as wave equations [39], simulations of chemical reactions [7], population dynamics [12][13][14] or other biological systems [3]. However, many classical time discretization schemes suffer from conserving non-negativity for arbitraty time step sizes [40]. Hence, designing non-negativity preserving numerical solution schemes for such problems can be regarded as an art as stated by Mickens [41]. For these reasons, Wacker and Schlüter examined a numerical solution scheme for (1.1) based on implicit Eulerian time discretization in [15] whose algorithm is presented in Algorithm 1. Define h j+1 := t j+1 − t j as current time step size 11: Define A j+1 := 1 + β j+1 · h j+1 · α j+1 · h j+1 12: 14: Compute R j+1 := R j + β j+1 · h j+1 · I j+1 16: end for 17: Outputs 18: Calculated vectors S, I, R ∈ R M Such non-standard finite-difference-methods for time discretization are also important in other application areas such as magnetohydrodynamics [42]. Here, we use f t j := f j as an abbreviation for time-continuous functions at time points t j and f h t j := f h, j as an abbreviation for time-discrete functions at time time points t j on possibly non-equidistant meshes t 1 < t 2 < . . . < t M−1 < t M . Although Algorithm 1 preserves all desired properties of the time-continuous model (1.1) as shown in [15], there are still some disadvantages: 1) From its formulation (1.1), there is to solve a quadratic equation where the a-priori knowledge of non-negative solutions from the time-continuous model is used to rule out one of both possible solutions for Algorithm 1 to obtain unique solvability; 2) In each time step, we must recalculate the auxiliary quantities A j+1 and B j+1 ; 3) By both aforementioned arguments, we conclude that Algorithm 1 does not intrinsically calculate biologically correct solutions and, from an algorithm point of view, recalculation of these auxiliary quantities is undesirable.
We want to propose a different discretization by non-standard finite-difference-methods based on non-local approximations. To explain the idea of non-standard finite-difference-methods based on non-local approximations, we consider a time-continuous dynamical system Standard time-discretizations normally approximate the right-hand-side function by function values of y (t) at one previous or the current time point. Non-standard finite-difference-methods based on nonlocal approximations, on the contrary, approximate f (t, y (t)) by function values from more then solely one previous or the current time point. For further details, we refer to the book [41] by Mickens. For the aforementioned three reasons, we suggest a non-standard finite-difference-method (NSFDM) given by in Section 3 as our first main contribution, after summarizing analytical properties of the timecontinuous problem formulation (1.1) in Section 2. In Section 3, we additionally prove unconditional preservation of non-negativity and linear convergence of (1.3) towards the time-continuous solution of (1.1). Additionally, (1.3) transfers to a uniquely solvable algorithm which preserves non-negativity intrinsically. Based on (1.3), we design a numerical algorithm to solve the inverse problem of finding the time-varying parameter functions in Section 4 as our next main contribution. Finally, we strengthen our theoretical findings in Section 5 by providing some numerical examples, especially one numerical example where linear convergence is stressed, in contrast to [15], where linear convergence was not shown numerically.

Time-continuous problem formulation
Following [15], we take the following assumptions into consideration with regard to our timecontinuous model (1.1): (1). The total population size N is fixed for all time t ≥ 0, i.e., it holds for all t ≥ 0 and we assume I 0 > 0 to really get an evolution in time of the infected subgroup of the population; (2). We divide the total population into three homogeneous subgroups, namely susceptible people (S), infected people (I) and recovered people (R). The equation is valid for all t ≥ 0; (3). The time-varying transmission rate α : [0, ∞) −→ [0, ∞) is at least continuously differentiable once and there exist positive constants α min > 0 and α max > 0 such that holds for all t ≥ 0; (4). The time-varying recovery rate β : [0, ∞) −→ [0, ∞) is at least continuously differentiable once and there exist positive constants β min > 0 and β max > 0 such that holds for all t ≥ 0.
We want to note that S and R depend on the initial conditions.
while the same supremum norm on t p , t p+1 reads

Conservation of total population size N
As our first result, our time-discrete problem formulation (1.3) conserves total population size N for all time points t j ∈ [0, T ] with j ∈ {1, 2, . . . , M}.
Proof. This statement can be proven by induction. For j = 1, it is clear by our initial conditions that holds. Let us assume that is valid for an arbitrary j ∈ {1, 2, . . . , M − 1}. Adding all three equations of (1.3), we obtain and this implies S h, j+1 + I h, j+1 + R h, j+1 = S h, j + I h, j + R h, j which proves our assertion (3.1).

Boundedness and non-negativity
As our next result, we show that approximate solutions of our time-continuous model (1.1) from our discretization model (1.3) remain non-negative and bounded.

Unique solvability
It directly follows that our numerical solution scheme (1.3) is uniquely solvable for all j ∈ {1, 2, . . . , M}. Proof. (3.3)-(3.5) prove that our numerical solution scheme (1.3) is well-defined. Unique solvability is a direct consequence of these equations as well. Hence, our assertion is proven.

Monotonicity properties and long-time behavior
We show that monotonicity properties and long-time behavior of the time-continuous problem formulation (1.1) translates to our time discretization (1.3). is non-negative and bounded above by the total population size N, the inequality (3.6) implies that this sequence is also monotonically decreasing which proves our first claim. (2). Since R h, j M j=1 is non-negative and bounded above by the total population size N, the inequality (3.5) yields that this sequence is monotonically increasing which shows our second assertion.
(3). We conclude (3.5) and this shows that our third claim is also true.
Remark 2. Both Theorems 1 and 5 imply that our time-continuous problem formulation (1.1) and our time-discrete problem formulation (1.3) converge towards the disease-free equilibrium state in the long run.

Error analysis
We start this subsection by listing all our assumptions for our main result regarding linear convergence of the time-discrete solution towards the time-continuous one: there are non-negative constants α min , α max , β min , β max such that Proof. We shortly explain our proof's strategy because it is relatively technical. If we first assume that our time-discrete solution coincides with the time-continuous solution at a certain time point t p ∈ [0, T ] for an arbitrary p ∈ {1, 2, . . . , M − 1}, then we investigate the local error propagation to the next time point t p+1 . Secondly, we consider error propagation in time and finally, we examine accumulation of these errors. We denote time-discrete solutions, for example, by S h,p and time-continuous solutions by S t p at the same time point t p .
1) For investigation of local errors, we assume that t p , S h,p = t p , S t p ; t p , I h,p = t p , I t p and t p , R h,p = t p , R t p hold for an arbitrary p ∈ {1, 2, . . . , M − 1} on the time-interval t p , t p+1 . Since we only take one time step into account, we denote corresponding time-discrete solutions by S h,p+1 , I h,p+1 and R h,p+1 . 1.1) We have by (3.6) from Remark 1 that holds. This implies (3.10) By the mean value Theorem, there exists ξ S ,1 ∈ t p , t p+1 such that we obtain This yields (3.12) For II S ,1 , we notice that and hence, holds due to boundedness of time-continuous solutions. Using again the mean value theorem for III S ,1 , we see that there is ξ α,1 ∈ t p , t p+1 such that we obtain (3.14) Hence, we conclude for III S ,1 . Plugging (3.15) into (3.13) yields (3.16) Thus, we finally obtain S t p+1 − S h,p+1 (3.10) (3.17) 1.2) By (3.7) from Remark 1, we get .

(3.23)
By boundedness and the mean value theorem, we receive the following estimates for V I I,1 , V I,1 and IV I,1 . By boundedness, the triangle inequality and the mean value theorem, we can deduce III I,1 = α p+1 N · I t p+1 · S t p+1 − I t p · S t p+1 and consequently, this implies (3.26) Hence, we conclude that ≤ I I,1 + II I,1 holds. By defining C loc,I : + h p+1 · β max · I t p+1 − I h,p+1 . (3.31) By using similar argument as in previous steps with the help of the mean value theorem and (3.29), we can deduce (3.32) 1.4) By setting C loc := max C loc,S , C loc,I , C loc,R , we conclude (3.34) 2) Normally, the points t p , S p , t p , I p and t p , R p do not exactly lie on the graph of the time-continuous solution. We then must investigate how procedural errors such as S h,p − S t p , holds where I S ,2 is given and estimated by N · S h,p · I h,p · I t p − S t p · I t p · I h,p =:III S ,2 . (3.36) We have and III S ,2 = S h,p · I h,p · I t p − S t p · I t p · I h,p (3.39) Using (3.35), we obtain (3.40)

2.2) We see that
holds and get (3.42)

2.4) For our summary, we define
which proves linear convergence of the time-discrete solution towards the time-continuous solution.

Algorithmic summary
Our algorithm for our proposed non-standard finite-difference-method for (1.1) is portrayed in Algorithm 2 below.
from (3.6) 12: from (3.7) 13: Our idea now is to reduce (1.3) to for all p ∈ {1, . . . , M − 1}. Hence, we obtain the sought time-varying transmission rate and the sought time-varying recovery rate by This finishes our claim's proof.

Algorithmic summary
Our algorithmic summary for our proposed parameter identification problem (4.2) reads:

Numerical examples
In this section, we want to strengthen our theoretical findings by two numerical examples. We apply GNU Octave for our implementation [43].

Example 1: Theoretical results
The initial conditions for our first example read S (0) = 5000, I (0) = 5000, R (0) = 0, N = 10, 000, α = 0.5, β = 0.1 and our simulation interval is given by [0, T ] with T = 70. We compare three numerical algorithms, namely our proposed non-standard finite-difference-method (NSFDM), the explicit Eulerian scheme (EE) and an explicit Runge-Kutta scheme with two stages (RK2) and use equidistant time step sizes h for our computations.      As it can be clearly seen in Figures 1 and 2, only our novel non-standard finite-difference-method is unconditionally non-negativity-preserving and stable. Both explicit standard time-discretizations suffer from time-step restrictions to obtain non-negativity preservation. In Figure 3, we observe that all three numerical algorithms possess same properties for smaller time step sizes h. On the positive side for all algorithms, they conserve total population size N as depicted in Figure 4.
Finally, underlining our convergence result from Theorem 6, we compare results of our nonstandard finite-difference-method (NSFDM) for different time step sizes h with a fine scale Runge-Kutta scheme of order 4 (RK4) at T = 70 where a fine time step size of 0.0001 is applied for (RK4). The errors  Table 1 that the theoretical convergence order is achieved by our non-standard finite-difference-method. We consider real-world data from Spain [44]. N = 47, 370, 000 is the total population size of Spain. The cumulative cases of infected people are modified such hold from the given data. We take α (t) = 0.52 · e −0.03·t and β (t) = 0.045. For our simulation, we additionally assume h = 0.75 as an equidistant time step size.   Our estimated transmission rate α (t) and constant recovery rate β (t) are user-chosen and adapted to the data in Figures 6 and 7, while the modified data for infected and recovered people are depicted in Figure 5. We solely consider the beginning of the pandemic because, later, vaccines have been available and reinfections have occurred after virus modifications. In Figure 8, we see results of our numerical simulations with h = 0.75. Results are relatively sensitive with respect to time step sizes and estimates rates, e.g., time-varying transmission rates. Hence, we can only conclude that good data acquisition and different predictions in different settings are crucial for future pandemics to make predictions for courses of such future epidemics.

Concluding remarks
In this article, we first summarized some analytical properties of a time-continuous epidemiological model (1.1) as given in [15]. Afterwards, we designed a new non-standard finite-difference-method (1.3) as a time-discrete version of the continuous model. In Section 3, we contributed our main results on our proposed non-standard finite-difference-method (1.3). Namely, we investigated unconditionally non-negativity preservation, conservation of total population size, convergence towards a disease-free equilibrium point and linear convergence by letting the time size h → 0. Conclusively for our theoretical findings, we developped a simple parameter identification algorithm in Section (4). In Section 5, we gave two numerical examples to underline our theoretical findings. Comparing Algorithms 1 and 2, we conclude that the latter one is definitely easier to implement and it can serve as a starting point for constructing higher order time discretization schemes for (1.1). Regarding practical applications, our proposed time-discretization by a non-standard finite-difference-method preserves all desired properties from the time-continuous model which are expected from our real-word intuition. Additionally, the simple parameter identification algorithm produce fast graphs from given data where trends in time-varying transmission rates such as contact restrictions might be seen from.
There might be different directions to extend our work: • It seems straightforward to integrate further compartments into our model such that we can additionally investigate influences such as, for example, vaccines or the simple extension to the SEIR model. For example, we readily can construct as one variant of a non-standard finite-difference-method for the SEIR model and it follows directly that this scheme conserves total population size N and can be written in an explicit fashion as our proposed non-standard finite-difference-method for our considered SIR model by , R h, j+1 = R h, j + h j+1 · β j+1 · I h, j+1 (6.2) to see that it also preserves non-negativity assuming positive initial conditions. One disadvantage is that there are different possibilities to construct such non-standard finite-difference-methods by non-local approximations and a suitable choice depends on the time-continuous properties which should be transferred to the time-discrete case. However, as seen by this example, it is flexible to construct different possibilities of algorithms and they are easy to implement; • It might be of interest to construct higher order non-standard finite-difference-methods from our first-order non-standard finite-difference-method in order to get higher accuracy although, especially, our second example demonstrated that predictions of such models need to be considered with care. One straightforward way to improve accuracy, and hence, the convergence order of numerical solution schemes is application of Richardson's extrapolation [45]; • One can use our discretization model for the discretization of optimal control problems for minimizing certain epidemic aspect such as eradication time [46].