The general setting of this work is a dynamical system, on the set
, defined by the following system of autonomous ordinary differential equations:
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:
For both the continuous and the discrete systems, it is assumed that the functions and are as smooth as needed. Naturally, the consistency requirements,
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 second-order convergent, elementary stable, and preserve the property of having the set forward invariant with respect to the continuous system.
The construction of higher-order 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 , valid for a scalar differential equation, where the second-order 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]:
A numerical scheme (2) is called symmetric or time-reversible if
More generally, the scheme is said to be reversible for all whenever implies that .
We construct a NSFD scheme that is time-reversible and show that, in essence, this property makes our scheme of second-order 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 Kermack-McKendrik Susceptible-Infective (SI) epidemic model , 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  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 right-hand 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
We consider the numerical method,
which is a NSFD scheme , on which only the nonlocal approximation of the right-hand 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 , we will assume here that
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 fixed-point iterations, etc. However, in particular settings, as the ones considered in the sequel such conditions are relatively easy to formulate.
Assume that (7) holds. Then the
The NSFD scheme (6) is time-reversible;
be an eigenvalue of
with associated left eigenvector. Then multiplying both sides of (10) on the left by we obtain
It 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
If follows from (11) that
d) Let and let denote the solution of (1) with . Denote
Using the consistency conditions (3) one can show that
Then, for the local truncation error,
Then using (12) we obtain and second order accuracy follows from the standard theory for one-step 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
Let a scheme (2) be reversible for all . If
then the set is invariant under for every .
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 .
) is not easy to verify. It can be replaced by a stronger condition, requiring that the vectorbe 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
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 .
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
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
Considering that is compact, there exists such that
|are both (row or column) diagonally dominant matrices|
Then, these matrices are invertible and (19) can be written explicitly as
Further, we have
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.
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 .
Theorem 5 is a useful tool in determining the the invariance of under .
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,
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)]
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  is absorbed here in the parameters and .
where and are real numbers such that
It is easy so see, and it is also shown in , that for the populations , of vectors and hosts, respectively, we have
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.
with respective outward normal vectors and , we have
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:
Therefore, condition (21) holds for all and , where
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
It is easy to see from the structure of , and in (31) that
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.
The literature on high-order 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 time-reversible nonstandard finite difference (NSFD) scheme, with the classical denominator for the discrete derivative but a nonlocal discretization of the right-hand side of the differential equations. Our findings, directly applicable to mass action-type models in biology, are twofold. First, the NSFD scheme is second-order 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  as far as the stability for one-dimensional model is concerned.
-  R. Anguelov, T. Berge, M. Chapwanya, J.K. Djoko, P. Kama, J. M-S. 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), 818-854.
-  R. Anguelov, Y. Dumont, J.M-S. Lubuma and M. Shillor, Dynamically consistent nonstandard finite difference schemes for epidemiological models, Journal of Computational and Applied Mathematics, 255(2014), 161-182.
-  R. Anguelov, P. Kama and J.M.-S. Lubuma, On nonstandard finite difference models of reaction-diffusion equations, J. Computational and Applied Mathematics 175 (2005), 11–29.
R. Anguelov and J.M-S. Lubuma, Contributions to the
mathematics of the nonstandard finite difference method and applications,
Numerical Methods for Partial Differential Equations17 (2001), 518–543.
-  R. Anguelov, J.M-S. Lubuma, Nonstandard Finite Difference Method by Nonlocal Approximation, Mathematics and Computers in Simulation 61(3-6) (2003), 465–475.
-  R Anguelov, J M-S Lubuma, M. Shillor, Topological dynamic consistency of nonstandard finite difference schemes for dynamical systems, Journal of Difference Equations and Applications, 17 (12), (2011), 1769-1791
-  A Bermon, R J Plemmons, Nonnegative matrices in the mathematical sciences, SIAM, 1994.
-  D. T. Dimitrov and H.V. Kojouharov, Positive and elementary stable nonstandard numerical methods with applications to predator-prey models, Journal of Computational and Applied Mathematics 189 (2006), 98–108.
-  E. Hairer, C. Lubich, G. Wanner, Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations, Springer, Berlin, Heidelberg, New York, 2006.
-  H.W. Hethcote, The mathematics of infectious disease, SIAM Rev. 42 (2000), 599–653.
-  H.V. Kojouharov, S. Roy, M. Gupta, F. Alalhareth and J.M. Slezak, A second-order modified nonstandard theta method for autonomous differential equations, Applied Mathematics Letters 112, (2021), 106775.
-  J.M-S. Lubuma and A. Roux, An improved theta method for systems of ordinary differential equations, J. Difference Equations and Applications, 9 (2003), 1023–1035.
-  M Martcheva, Introduction to Mathematical Epidemiology, Springer, 2015.
-  R.E. Mickens, Nonstandard finite difference models of differential equations, World Scientific, Singapore, 1994.
-  R.E. Mickens (Ed.), Applications of Nonstandard Finite Difference Schemes, World Scientific, Singapore, 2000.
-  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.
-  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) 645-653.
-  A.M. Stuart and A.R. Humphries, Dynamical systems and numerical analysis, Cambridge University Press, New York, 1998.
-  W. Walter, Ordinary differential equations, New York, Springer, 1998.