DeepAI

# Second-order nonstandard finite difference schemes for a class of models in bioscience

We consider a dynamical system, defined by a system of autonomous differential equations, on Ω⊂ℝ^n. By using Mickens' rule on the nonlocal approximation of nonlinear terms, we construct an implicit Nonstandard Finite Difference (NSFD) scheme that, under an existence and uniqueness condition, is an explicit time reversible scheme. Apart from being elementary stable, we show that the NSFD scheme is of second-order and domain-preserving, thereby solving a pending problem on the construction of higher-order nonstandard schemes without spurious solutions, and extending the tangent condition to discrete dynamical systems. It is shown that the new scheme applies directly for mass action-based models of biological and chemical processes.

• 3 publications
• 1 publication
10/29/2022

### Structure preserving schemes for Fokker-Planck equations of irreversible processes

In this paper, we introduce second order and fourth order space discreti...
05/31/2021

### A novel second-order nonstandard finite difference method for solving one-dimensional autonomous dynamical systems

In this work, a novel second-order nonstandard finite difference (NSFD) ...
04/22/2021

### Explicit, time-reversible and symplectic integrator for Hamiltonians in isotropic uniformly curved geometries

The kinetic term of the N-body Hamiltonian system defined on the surface...
01/20/2020

### TPFA Finite Volume Approximation of Wasserstein Gradient Flows

Numerous infinite dimensional dynamical systems arising in different fie...
12/28/2020

### A High-Order Harmonic Balance Method for Systems With Distinct States

A pure frequency domain method for the computation of periodic solutions...
07/12/2021

### Extending nonstandard finite difference schemes rules to systems of nonlinear ODEs with constant coefficients

In this paper, we present a reformulation of Mickens' rules for nonstand...
04/14/2021

### Existence, Uniqueness and Numerical Modeling of Wine Fermentation Based on Integro-Differential Equations

Predictive modeling is the key factor for saving time and resources with...

## 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:

 dxdt≡˙x=f(x). (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:

 xk+1=F(h,xk)≡F(h)(xk). (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,

 F(0)(x)=x and dF(0)dh(x)=f(x), (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 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 [11], 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]:

###### Definition 1

A numerical scheme (2) is called symmetric or time-reversible if

 xk=F(−h,xk+1). (4)

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 [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 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

 ˙x=φ(x,x),  x∈Ω⊆Rn. (5)

We consider the numerical method,

 xk+1−xkh=12(φ(xk+1,xk)+φ(xk,xk+1)), (6)

which is a NSFD scheme [4], 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 [18], we will assume here that

 (i) given xk∈Ω, equation (6) has a unique solution xk+1 in Ω; (ii) different values of xk result in different solutions for xk+1.
(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 fixed-point iterations, etc. However, in particular settings, as the ones considered in the sequel such conditions are relatively easy to formulate.

###### Theorem 2

Assume that (7) holds. Then the

• The NSFD scheme (6) is time-reversible;

• The fixed points of the NSFD scheme (6) are precisely the equilibrium points of the system (5);

• The NSFD scheme (6) replicates the local stability of all hyperbolic equilbrium points of (5);

• The NSFD scheme (6) is a second-order scheme for the equation (5).

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

 F(h,x)−x=h2[φ(F(h,x),x)+φ(x,F(h,x))]. (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 time-reversible.

b) It is clear that is a fixed point of (6) if, and only if, , i.e. is an equilibrium point of (5).

c) Let us use the explicit form (2) of (6). From (8), the Jacobian matrix (in the variable) of satisfies

 ∂F(h,x)∂x−I = h2(∂φ∂y(F(h,x),x)∂F(h,x)∂x+∂φ∂z(F(h,x),x) (9) +∂φ∂y(x,F(h,x))+∂φ∂z(x,F(h,x))∂F(h,x)∂x)

Let be an arbitrary fixed point of (6). At , using , equation (9) simplifies to

or. equivalently,

 (I−h2∂f(¯¯¯x)∂x)∂F(h,¯¯¯x)∂x=I+h2∂f(¯¯¯x)∂x. (10)

Let

be an eigenvalue of

with associated left eigenvector

. Then multiplying both sides of (10) on the left by we obtain

 (1−h2λ)v∂F(h,¯¯¯x)∂x=(1+h2λ)v

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

 v∂F(h,¯¯¯x)∂x=1+h2λ1−h2λv.

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

 |μ|2 = (1+h2R(λ))2+h24I2(λ)(1−h2R(λ))2+h24I2(λ) (11) = 1+h24|λ|2+hR(λ)1+h24|λ|2−hR(λ)

If follows from (11) that

 R(λ)<0  ⟺  |μ|<1,

which implies that a hyperbolic equilibrium point of (5) is asymptotically stable if, and only if, it is asymptotically stable as a fixed-point of (6).

d) Let and let denote the solution of (1) with . Denote

 ξ=12(F(h)(x)+x).

Using the consistency conditions (3) one can show that

 ξ=12(x+F(0,x)+dFdh(0,x)+O(h2))=x+h2f(x)+O(h2)=u(h2)+O(h2) (12)

Then, for the local truncation error,

 E(h)=u(h)−F(h,x)=u(h)−x−h2φ(F(h)(x),x)−h2φ(x,F(h)(x)),

of the NSFD scheme (2) & (8), Taylor expansions of about and of about yield

 E(h) = h˙u(h2)−h2(φ(ξ,ξ)+∂φ(ξ,ξ)∂y(F(h)(x)−ξ)+∂φ(ξ,ξ)∂z(x−ξ)) −h2(φ(ξ,ξ)+∂φ(ξ,ξ)∂y(x−ξ)+∂φ(ξ,ξ)∂z(F(h)(x)−ξ))+O(h3) = h(f(u(h2))−f(ξ))+O(h3).

Then using (12) we obtain and second order accuracy follows from the standard theory for one-step numerical methods.

###### Remark 3

Since the scheme (6) is symmetric, the statement d) of the Theorem can be derived from the general theory of symmetric schemes [9, Theorem II.3.2].

## 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

 D=closure(∘D)=∘D∪∂D. (13)
###### Theorem 4

Let a scheme (2) be reversible for all . If

 F(−h,∂D)∩∘D=∅, h∈(0,¯¯¯h], (14)

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 .

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

 n(x)⋅f(x)≤0, (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

Let a scheme (2) be reversible for all and let be a closed and convex subset of satisfying (13). If for every , and an outer normal vector we have

 n(x)⋅(F(−h,x)−x)≥0, (16)

then the set is invariant under for every .

Proof. Let and . The inequality (16) implies that the point is on the same side as with respect to the tangent to at and associated with . Considering that is convex, this implies that . Then the conclusion follows from Theorem 4.

## 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

 φ(x,y)=P(x)y+A(x+y)+b=Q(y)x+A(x+y)+b (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

 xk+1−xkh=12(P(xk)xk+1+Q(xk)xk+1)+A(xk+1+xk)+b (18)

or, equivalently,

 (I−h2(P(xk)+Q(xk))−hA)xk+1=(I+hA)xk+hb, (19)

where is the identity matrix. Furthermore, and also expressing it reversibility, the equation (18) can be written in the form

 (I+h2(P(xk+1)+Q(xk+1))+hA)xk=(I−hA)xk+1−hb. (20)

Considering that is compact, there exists such that

 I−h2(P(x)+Q(x))−hA  and  I+h2(P(x)+Q(x))+hA (21) are both (row or column) diagonally dominant matrices for all x∈(0,¯¯¯h] and x∈Ω.

Then, these matrices are invertible and (19) can be written explicitly as

 xk+1=F(h,xk) (22)

where

 F(h,x)=(I−h2(P(x)+Q(x))−hA)−1((I+hA)x+hb). (23)

Further, we have

 F−1(h,x)=F(−h,x)=(I+h2(P(x)+Q(x))+hA)−1((I−hA)x−hb). (24)
###### Theorem 6

Let there exist such that (21) holds and is an invariant set of the map for all . Then the scheme (18) satisfies the properties a)–d) in Theorem 2.

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,

 I−h2(P(x)+Q(x))−hA

is an M-matrix. Hence, its inverse is a nonnegative matrix, [7, Theorem 2.3]. Therefore, it follows from (22) that

 Ω⊆Rn+ ⟹ F(h,Ω)⊆Rn+, h∈(0,¯h]. (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)]

 dSvdt = Λv−pSvI−μvSv, (26) dIvdt = pSvI−μvIv, (27) dSdt = Λ−qSIv−μS+γR, (28) dIdt = qSIv−(μ+α)I, (29) dRdt = αI−(μ+γ)R, (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 .

Let . Then the system (26)–(30) can be stated in the form (5) with given by (17), where

 P(x)=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝−pI0000pI000000−qIv0000qIv0000000⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,  Q(x)=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝000−pSv0000pSv00−qS0000qS00000000⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠, (31) A=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝−μv00000−μv00000−μ0γ000−(μ+α)0000α−(μ+γ)⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,  b=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝Λv0Λ00⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠.

The set of equations (26)–(30) defines a dynamical system on

 Ω={x∈R5+:Sv+Iv≤Mv,S+I+R≤M},

where and are real numbers such that

 Mv≥Λvμv,  M≥Λμ. (32)

It is easy so see, and it is also shown in [13], that for the populations , of vectors and hosts, respectively, we have

 limt→+∞Nv(t)=Λvμv,  limt→+∞N=Λμ. (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

 Sv+Iv=Mv, (34) S+I+R=M, (35)

with respective outward normal vectors and , we have

 u⋅f(x) = Λv−μv(Sv+Iv)=Λv−μMv≤0, (36) w⋅f(x) = Λ−μ(S+I+R)=Λ−μM≤0. (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:

 h(pI+μv) ≤ 1, (38) h(qS+μv ≤ 1, h(qIv+μ) ≤ 1, h(pSv+μ+α) ≤ 1, h(μ+2γ) ≤ 1.

Therefore, condition (21) holds for all and , where

 ¯¯¯h=min{1pM+μv,1qM+μv,1qMv+μ,1pMv+μ+α,1μ+2γ}. (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

 (I+h2(P(x)+Q(x))+hA)(F(−h,x)−x) =(I−hA)x−hb−(I+h2(P(x)+Q(x))+hA)x =h2(P(x)x+2Ax)+h2(Q(x)x+2Ax =−h2φ(x,x)−h2φ(x,x) = −hf(x) (40)

It is easy to see from the structure of , and in (31) that

 uTP=wTP=uTQ=wTQ=(0,0,0,0,0)T,  uTA=−μvu,  wTA=−μw. (41)

Putting together (40) and (41) we obtain

 (1−hμv)u⋅(F(−h,x)−x) = uT(I+h2(P(x)+Q(x))+hA)(F(−h,x)−x) = −hu⋅f(x)≥0.

Using that and (36), we have

 u⋅(F(−h,x)−x)=−h1−hμvu⋅f(x)≥0

for all on the plane (34) and . Similarly, we obtain

 w⋅(F(−h,x)−x)=−h1−hμvu⋅f(x)≥0

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 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 [6] as far as the stability for one-dimensional model is concerned.

## References

• [1] 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.
• [2] 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.
• [3] 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.
• [4] R. Anguelov and J.M-S. 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.M-S. Lubuma, Nonstandard Finite Difference Method by Nonlocal Approximation, Mathematics and Computers in Simulation 61(3-6) (2003), 465–475.
• [6] 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
• [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 predator-prey models, Journal of Computational and Applied Mathematics 189 (2006), 98–108.
• [9] E. Hairer, C. Lubich, G. Wanner, Geometric Numerical Integration: Structure-Preserving 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 second-order modified nonstandard theta method for autonomous differential equations, Applied Mathematics Letters 112, (2021), 106775.
• [12] 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.
• [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) 645-653.
• [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.