1 Introduction
The general setting of this work is a dynamical system, on the set
, defined by the following system of autonomous ordinary differential equations:
(1) 
The solutions, , of the continuous system are approximated at the discrete times, , , , by a sequence obtained recursively through a numerical method of the following form:
(2) 
For both the continuous and the discrete systems, it is assumed that the functions and are as smooth as needed. Naturally, the consistency requirements,
(3) 
are assumed to hold, and thus we suppose that is at least continuously differentiable on . Our aim is that the numerical method (2) be reliable in the three important directions below. The numerical method must be secondorder convergent, elementary stable, and preserve the property of having the set forward invariant with respect to the continuous system.
The construction of higherorder NSFD schemes that are dynamically consistent with the underlying features of the continuous differential equations models, particularly for those with transient dynamics, is a pending problem. Several authors have attempted to address this problem. These include the nonstandard versions of the classical method, with , developed in
[1, 2, 3, 12] where the positivity of the discrete solution is achieved by also using Mickens’ rule on a complex denominator function for the discrete derivative. In the same vein, we mention the recent work [11], valid for a scalar differential equation, where the secondorder accuracy and elementary stability are achieved by a modified nonstandard method, , in which the complex denominator function varies at each iteration.
In this paper, we deal with numerical methods characterized as follows [9, Section V.1]:
Definition 1
A numerical scheme (2) is called symmetric or timereversible if
(4) 
More generally, the scheme is said to be reversible for all whenever implies that .
We construct a NSFD scheme that is timereversible and show that, in essence, this property makes our scheme of secondorder accuracy and elementary stable (Section 2). The same property is used to establish a discrete analogue of the tangent condition for forward invariance of a convex set, which implies that the domain under consideration is positively invariant with respect to our NSFD scheme (Section 3). The case of mass action models of biological and chemical processes is considered in details with KermackMcKendrik SusceptibleInfective (SI) epidemic model [10], as an illustrative example (Section 4). Concluding remarks, with possible extensions of this work, are given in Section 5
2 Nonstandard finite difference method of order 2
In the construction of NSFD schemes, Mickens [14] made an important observation namely, that the discrete models of differential equations have a larger parameter space than the corresponding differential equations. Adding to this Mickens’ Rule 1, which says that the orders of the discrete derivatives must be exactly equal to the orders of the corresponding derivatives of the differential equations, it is convenient for nonlocal approximations to view the righthand side of the system (1) as a restriction, on the diagonal , of a certain function of two variables such that and hence the system takes the form
(5) 
We consider the numerical method,
(6) 
which is a NSFD scheme [4], on which only the nonlocal approximation of the righthand side of (5) is used, the rule on the complex denominator function of the derivatives being excluded. The NSFD scheme (6) is implicit. The existence and uniqueness of solution of the respective system of algebraic equation is an issue which requires attention on its own. To make this NSFD scheme a generalized dynamical system on [18], we will assume here that

(7) 
Finding general conditions for (7) to hold is challenging, as it is the case for general equations of the form solved by using for instance the implicit function theorem, sophisticated fixedpoint iterations, etc. However, in particular settings, as the ones considered in the sequel such conditions are relatively easy to formulate.
Proof. The assumption (7) implies that the NSFD scheme (6) can be written in the explicit form (2) where where the function satisfies the equation
(8) 
a) If in the equation (6) we interchange and and replace by , the equation remains the same. Therefore (4) holds, which means that the NSFD scheme is timereversible.
c) Let us use the explicit form (2) of (6). From (8), the Jacobian matrix (in the variable) of satisfies
(9)  
Let be an arbitrary fixed point of (6). At , using , equation (9) simplifies to
or. equivalently,
(10) 
Let
be an eigenvalue of
with associated left eigenvector
. Then multiplying both sides of (10) on the left by we obtainIt follows from (7) that the matrix is not singular. Then, due to the obvious impossibility for and to be both zero, neither of them is. Hence
Therefore, to every eigenvalue of with associated left eigenvector , there corresponds an eigenvalue of the matrix with the same left eigenvector . After some technical manipulation we obtain
(11)  
If follows from (11) that
which implies that a hyperbolic equilibrium point of (5) is asymptotically stable if, and only if, it is asymptotically stable as a fixedpoint of (6).
d) Let and let denote the solution of (1) with . Denote
Using the consistency conditions (3) one can show that
(12) 
Then, for the local truncation error,
of the NSFD scheme (2) & (8), Taylor expansions of about and of about yield
Then using (12) we obtain and second order accuracy follows from the standard theory for onestep numerical methods.
3 Domain preserving property of reversible schemes
The theory for continuous dynamical systems of the form (1) provides theorems proving that a set is forward invariant through conditions only on the boundary of the set, e.g. [19, Section 10]. In this section we derive similar theorems for reversible maps, that is appropriate assumptions on the boundary of a set imply that the set is forward invariant.
Let be a closed subset of . We denote by the interior of , by  the boundary of and assume that
(13) 
Theorem 4
Proof. Let . Assume that there exists such that . Since is continuous on , there exists such that . Equivalently, , which contradicts the condition (14). Therefore, . Considering that is an arbitrary element of , we have . Using the continuity of on and (13) we obtain .
Condition (14
) is not easy to verify. It can be replaced by a stronger condition, requiring that the vector
be outward directed or tangential to at . This condition is similar to the tangent condition in [19, Section 10.XV], which is the essential property for invariance of sets under flows. In the setting of the system (1) for the set the tangential condition states that(15) 
for every point and any outer normal vector at . Let us recall that a vector is called outer normal vector to at if the circle with center and radius has no intersection with . The straight line through perpendicular to is called a tangent. This general definition of normal vector and tangent does not require any smoothness of the boundary , which is rather convenient in applications. Let us note that the normal vector , similarly the associated tangent, need not exist and, if existing, need not be unique.
In the next theorem we make an additional assumption on to be convex. Convexity of implies that a normal vector and a tangent exist at every point . Further, a defining property of being convex is that is on one side of any tangent at any point . Specifically, on the side not containing .
Theorem 5
4 Schemes for mass action models
Modeling of biological and chemical processes by applying the mass action principle for representing the interaction of the involved species typically results in a system of the form (5) in , a convex and compact subset of , where the function is linear in both its arguments. Then can be represented as
(17) 
where and are matrix functions of and , respectively. The matrix and the vector are constant. In this particular setting the numerical method (6) can be written in the form
(18) 
or, equivalently,
(19) 
where is the identity matrix. Furthermore, and also expressing it reversibility, the equation (18) can be written in the form
(20) 
Considering that is compact, there exists such that
(21)  
are both (row or column) diagonally dominant matrices  
Then, these matrices are invertible and (19) can be written explicitly as
(22) 
where
(23) 
Further, we have
(24) 
Theorem 6
Proof. The scheme (18) is constructed in the general form (6). Therefore, it is sufficient to show that condition (7) holds. The solution is given in the explicit form (22). Therefore, it exists for every and it is unique. Further, it belongs to , due to the assumption that is invariant.
Remark 7
The existence of such (21) holds follows from the compactness of . In practice, it is determined from the condition of diagonal dominance using the boundedness of the variable .
Remark 8
Theorem 5 is a useful tool in determining the the invariance of under .
Remark 9
It very common in applications that and the nondiagonal entries of , , , are nonnegative, that is, these matrices are Metzler matrices. Then, assuming that (21) holds,
is an Mmatrix. Hence, its inverse is a nonnegative matrix, [7, Theorem 2.3]. Therefore, it follows from (22) that
(25) 
The property (25) is usually a first step towards proving the invariance of under , .
5 Application to a model in Mathematical Epidemiology
In this section we show how the theory and tools developed in Sections 2, 3 and 4 can be put together in deriving second order scheme with the properties a)–c) in Theorem 2. We note that properties a)–c) are of qualitative nature and do not follows from the standard numerical analysis theory involving stability, consistency and therefore convergence. While convergence as is expected, the computations are always done for some positive value of . Hence, the importance that the main characteristics of the dynamical system are preserved as given in properties a)–c) for . While the theoretical existence of is important, deriving constructively value of is critical for any specific application. We show how the value of can be computed for the considered model. It should be noted that, while the numerical approach proposed in this paper is exemplified on the model considered in this section, its realm of application is not limited to it.
We consider the model of vector borne diseases with temporal immunity as given in [13, Section 4.4, (4.39)–(4.40)]
(26)  
(27)  
(28)  
(29)  
(30) 
where and are the susceptible and infective vectors and , , are the susceptible, invective and recovered (with temporary immunity) hosts, e.g. humans. All parameters are positive. Bearing in mind the units , are the recruitments, ,  the death rates for vectors and hosts, respectively, and
are the probabilities of sufficient contact for transmission from host to vector and from vector to host, the average period of infectivity of host is
, while the average duration of immunity of host is . Let us note that an additional parameter in [13] is absorbed here in the parameters and .The set of equations (26)–(30) defines a dynamical system on
where and are real numbers such that
(32) 
It is easy so see, and it is also shown in [13], that for the populations , of vectors and hosts, respectively, we have
(33) 
Hence, the inequalities (32) are important. The limit properties (33) are used in [13]
to reduce the asymptotic analysis of (
26)–(30) to asymptotic analysis of a system to two equations. We consider here the full system (26)–(30), since the goal is to compute accurately the trajectories for all finite times and not only their limit at . Further, this illustrates that the applicability of the proposed numerical approach is not restricted in any way by the dimensionality of the system.The vector field , defined by right hand side of (26)–(30) satisfies the tangential condition in [19, Section 10.XV]. Specifically, on the planes
(34)  
(35) 
with respective outward normal vectors and , we have
(36)  
(37) 
We consider the numerical method (22) for the system (26)–(30) on . Our aim is apply Theorem 6. From the structure of the matrices , and it follows that the matrix is column diagonally dominant for every and . The column diagonal dominance of the matrix is equivalent to the following set of inequalities:
(38)  
Therefore, condition (21) holds for all and , where
(39) 
It remains to show the invariance of . Considering that , , are Metzler matrices, we have that (25) holds. Therefore, it is sufficient to apply Theorem 5 on the part of the boundary defined by planes (34)–(35). Using (24), we have
(40) 
It is easy to see from the structure of , and in (31) that
(41) 
Putting together (40) and (41) we obtain
Using that and (36), we have
for all on the plane (34) and . Similarly, we obtain
for all on the plane (35) and . Then it follows from Theorem 5 that is positively invariant under all maps , . Hence the numerical method (22) for the system (26)–(30) with given in (39) satisfies properties a)–d) in Theorem 2.
6 Conclusion
The literature on highorder numerical methods for dynamical systems defined by ordinary differential equations is rich. However, these classical schemes can exhibit spurious/ghost solutions or other elementary instability that do not correspond to the feature of the continuous equations. We overcome this difficulty by constructing a timereversible nonstandard finite difference (NSFD) scheme, with the classical denominator for the discrete derivative but a nonlocal discretization of the righthand side of the differential equations. Our findings, directly applicable to mass actiontype models in biology, are twofold. First, the NSFD scheme is secondorder convergent and elementary stable. Second, we formulate a discrete analogue of the tangent condition under which it is shown that the NSFD scheme preserves the propriety of the continuous model stating that its feasible region is forward invariant.
Our future research includes extending this study to a dynamical system with nonhyperbolic equilibrium points, an approach considered in [6] as far as the stability for onedimensional model is concerned.
References
 [1] R. Anguelov, T. Berge, M. Chapwanya, J.K. Djoko, P. Kama, J. MS. Lubuma and Y. Terefe, Nonstandard Finite difference method revisited and application to the Ebola Virus Disease dynamics transmission, Journal of Difference Equations and Applications 26, 2020(6), 818854.
 [2] R. Anguelov, Y. Dumont, J.MS. Lubuma and M. Shillor, Dynamically consistent nonstandard finite difference schemes for epidemiological models, Journal of Computational and Applied Mathematics, 255(2014), 161182.
 [3] R. Anguelov, P. Kama and J.M.S. Lubuma, On nonstandard finite difference models of reactiondiffusion equations, J. Computational and Applied Mathematics 175 (2005), 11–29.

[4]
R. Anguelov and J.MS. Lubuma, Contributions to the
mathematics of the nonstandard finite difference method and applications,
Numerical Methods for Partial Differential Equations
17 (2001), 518–543.  [5] R. Anguelov, J.MS. Lubuma, Nonstandard Finite Difference Method by Nonlocal Approximation, Mathematics and Computers in Simulation 61(36) (2003), 465–475.
 [6] R Anguelov, J MS Lubuma, M. Shillor, Topological dynamic consistency of nonstandard finite difference schemes for dynamical systems, Journal of Difference Equations and Applications, 17 (12), (2011), 17691791
 [7] A Bermon, R J Plemmons, Nonnegative matrices in the mathematical sciences, SIAM, 1994.
 [8] D. T. Dimitrov and H.V. Kojouharov, Positive and elementary stable nonstandard numerical methods with applications to predatorprey models, Journal of Computational and Applied Mathematics 189 (2006), 98–108.
 [9] E. Hairer, C. Lubich, G. Wanner, Geometric Numerical Integration: StructurePreserving Algorithms for Ordinary Differential Equations, Springer, Berlin, Heidelberg, New York, 2006.
 [10] H.W. Hethcote, The mathematics of infectious disease, SIAM Rev. 42 (2000), 599–653.
 [11] H.V. Kojouharov, S. Roy, M. Gupta, F. Alalhareth and J.M. Slezak, A secondorder modified nonstandard theta method for autonomous differential equations, Applied Mathematics Letters 112, (2021), 106775.
 [12] J.MS. Lubuma and A. Roux, An improved theta method for systems of ordinary differential equations, J. Difference Equations and Applications, 9 (2003), 1023–1035.
 [13] M Martcheva, Introduction to Mathematical Epidemiology, Springer, 2015.
 [14] R.E. Mickens, Nonstandard finite difference models of differential equations, World Scientific, Singapore, 1994.
 [15] R.E. Mickens (Ed.), Applications of Nonstandard Finite Difference Schemes, World Scientific, Singapore, 2000.
 [16] R.E. Mickens, Nonstandard finite difference methods, In: R.E. Mickens (Ed), Advances in the applications of nonstandard finite difference schemes, World Scientific, Singapore, 2005, pp. 1–9.
 [17] R.E. Mickens, Dynamic consistency: a fundamental principle for constructing nonstandard finite difference schemes for differential equations, Journal of Difference Equations and Applications 11 (2005) 645653.
 [18] A.M. Stuart and A.R. Humphries, Dynamical systems and numerical analysis, Cambridge University Press, New York, 1998.
 [19] W. Walter, Ordinary differential equations, New York, Springer, 1998.