# A generalized nonstandard finite difference method for a class of autonomous
dynamical systems and its applications

Manh Tuan Hoang  Department of Mathematics, FPT University, Hoa Lac Hi-Tech
Park, Km29 Thang Long Blvd, Hanoi, Viet Nam<EMAIL_ADDRESS><EMAIL_ADDRESS>This paper is dedicated in honor of Ronald E.
Mickens’ 80th birthday.

###### Abstract.

In this work, a class of continuous-time autonomous dynamical systems
describing many important phenomena and processes arising in real-world
applications is considered. We apply the nonstandard finite difference (NSFD)
methodology proposed by Mickens to design a generalized NSFD method for the
dynamical system models under consideration. This method is constructed based
on a novel non-local approximation for the right-side functions of the
dynamical systems. It is proved by rigorous mathematical analyses that the
NSFD method is dynamically consistent with respect to positivity, asymptotic
stability and three classes of conservation laws, including direct
conservation, generalized conservation and sub-conservation laws. Furthermore,
the NSFD method is easy to be implemented and can be applied to solve a broad
range of mathematical models arising in real-life. Finally, a set of numerical
experiments is performed to illustrate the theoretical findings and to show
advantages of the proposed NSFD method.

###### Key words and phrases:

Mickens’ methodology, nonstandard finite difference method, dynamic
consistency, positivity, stability, conservation laws

###### 2020 Mathematics Subject Classification:

Primary 34A45, 65L05

##  1\. Introduction

Many important problems in real-world situations can be mathematically modeled
by autonomous dynamical systems of the general form:

(1.1) $\dfrac{dy(t)}{dt}=f(y(t)),\quad t\geq 0,\quad\quad
y(0)=y_{0}\in\mathbb{R}^{n},$

where $y$ is a vector function of $t$ with $n$ components; $f$ is a function
of $y$ and it is assumed to satisfy suitable conditions that guarantee
solutions of the model (1.1) exist and are unique (see, for instance, [4, 14,
32, 49, 54, 57, 84, 85]). The solutions of the model (1.1) often possess
several important properties, typically positivity, asymptotic stability,
boundedness, periodicity and physical properties.

In this work, we consider the model (1.1) under the following two conditions:

  * •

(C1): There exists a real number $\alpha$ such that

(1.2) $f(y)+\alpha y\geq 0\quad\mbox{for all}\quad y\geq 0.$

Here the notation $\geq$ is meant entry-wise for vectors.

  * •

(C2): The set of equilibria is finite and all the equilibria are hyperbolic.

The condition (C1) implies that the set
$\mathbb{R}^{n}_{+}=\big{\\{}y\in\mathbb{R}^{n}|y\geq 0\big{\\}}$ is a
positively invariant set of the model (1.1), that is, $y(t)\geq 0$ whenever
$y_{0}\geq 0$ (see [45, 84]), meanwhile, the condition (C2) means that the
local asymptotic stability (LAS) of all the equilibria can be determined by
the Lyapunov’s indirect method (see [49, 52, 53, 85]). Also, in many cases the
model (1.1) can satisfy three classes of conservation laws, which include
direct conservation, generalized conservation and sub-conservation laws [65].
It is easy to find a large number of mathematical models in real-world
applications with the above-mentioned properties. Below, we list some notable
models that satisfy the conditions (C1) and (C2) and/or the conservation laws

  * •

mathematical models in biology, ecology and epidemiology [4, 14, 32, 54];

  * •

models for interacting species [14];

  * •

chemostat systems [84];

  * •

epidemic models for infectious diseases [39, 58];

  * •

metapopulation models [5, 48];

  * •

single-species population models [18];

  * •

predator-prey models with linear prey growth and Beddington-DeAngelis
functional response [27, 31];

  * •

a simple vaccination model with multiple endemic states [51];

  * •

a vaccination model with non-linear incidence [38];

  * •

a two-stage epidemic model with generalized non-linear incidence [77];

  * •

a generalization of the Kermack-McKendrick epidemic model [15].

It should be emphasized that dynamically consistent NSFD schemes for most of
the above models have not been studied.

Our main objective is to construct a numerical method preserving the
positivity, LAS and conservation laws of the model (1.1) for all the values of
the step size. To achieve this goal, we will apply the Mickens’ methodology
[59, 60, 61, 62, 63] to propose a dynamically consistent nonstandard finite
difference (NSFD) method for the model (1.1). We recall that the concept of
NSFD schemes was first introduced by Mickens to overcome usual numerical
instabilities caused by standard finite difference (SFD) schemes as well as
drawbacks of SFD schemes [59, 60, 61, 62, 63, 64]. One of the main and
outstanding advantages of NSFD schemes is that they have the ability to
preserve essential mathematical properties of the corresponding differential
equation models regardless of the values of the step size; furthermore, they
are easy to be implemented. Therefore, NSFD schemes are suitable and effective
to simulate dynamics of differential equation models over very long time
periods. During the past several decades, NSFD schemes have been strongly
developed by mathematicians and engineers and have become one of the most
powerful methods for solving differential equations. Nowadays, NSFD schemes
have been widely used for ordinary differential equations, partial
differential equations, delay differential equations, fractional differential
equations and integro-differential equations (see, for instance, [1, 2, 6, 7,
8, 9, 12, 16, 19, 20, 28, 29, 30, 33, 35, 36, 41, 50, 59, 60, 61, 62, 63, 64,
65, 66, 67, 68, 69, 70, 71, 71, 72, 73, 74, 75, 78, 79, 80, 81, 82, 83, 86,
87, 88, 89] and references therein). Here, we refer the readers to [62, 78,
79] for notable reviews of NSFD schemes.

It is well-known that the positivity and LAS are essential properties of most
systems in biology, ecology and epidemiology (see, for example, [4, 14, 32,
49, 54, 84, 85]). For this reason, the construction of NSFD schemes preserving
these properties is very important and has attracted the attention of many
researchers. Below, we mention some notable results in this topic.
In an early work [6], Anguelov and Lubuma introduced an NSFD method preserving
monotonic properties of solutions of differential equations and the LAS of
hyperbolic equilibrium points. In 2005 and 2007, Dimitrov and Kojouharov in
[28, 29] analyzed stability-preserving NSFD schemes based on the
$\theta$-methods and the second-order Runge-Kutta methods for general
autonomous dynamical systems. Note that the NSFD schemes constructed in [6,
28, 29] are not positive. In [30, 87], dynamically consistent NSFD schemes for
general productive-destructive systems were formulated. In 2015, Wood and
Kojouharov in [88] designed an NSFD method preserving the positivity of
solutions and the local behavior of dynamical systems near equilibria. In
[19], Cresson and Pierret studied the construction of NSFD schemes for a
general class of two-dimensional dynamical systems including many models in
population dynamics. In 2020, Dang and Hoang constructed and analyzed high-
order nonstandard Runge-Kutta methods for a class of autonomous dynamical
systems possessing the positivity and LAS. Very recently, Hoang in [40]
proposed a second-order NSFD method that simultaneously preserves the
positivity and LAS of one-dimensional autonomous dynamical systems. Besides,
many results focusing on NSFD schemes preserving the positivity and LAS of
specific models in biology, ecology, biology and other fields are also
remarkable (see, for example, [1, 2, 8, 80, 81, 89] and references therein).

Motivated and inspired by the above reasons, we will apply the NSFD
methodology proposed by Mickens to construct a dynamically consistent NSFD
scheme for the model (1.1). This scheme is constructed based on a novel non-
local approximation for the right-side functions. It is proved by rigorous
mathematical analyses that the NSFD scheme preserves the positivity,
asymptotic stability and conservation laws for all the values of the step
size. Moreover, the constructed NSFD scheme is easy to be implemented and can
be used to solve a broad range of mathematical problems arising in real-life.
This claim is also confirmed by a set of illustrative numerical examples.

The plan of this work is as follows:
Section 2 provides some basic concepts and preliminaries. The NSFD scheme is
introduced and analyzed in Section 3. A set of numerical experiments is
investigated and reported in Section 4. Finally, conclusions, remarks and some
open problems are given in Section 5.

##  2\. Definitions and Preliminaries

In this section, we recall from literature some essential concepts and
preliminaries related to qualitative theory of continuous and discrete
dynamical systems, conservation laws and NSFD schemes.

###  2.1. Stability of equilibria of continuous-time dynamical systems

Consider the continuous-time dynamical system (1.1). We recall that a point
$y^{*}\in\mathbb{R}^{n}$ is called an equilibrium point of System (1.1) if
$f(y^{*})=0$ (see [49, 85]). Without loss of generality, we can always assume
that $y^{*}=0$.

######  Definition 2.1 ([49, 85]).

The equilibrium point $y^{*}=0$ of (1.1) is said to be:

  * •

stable, if, for each $\epsilon>0$, there is $\delta=\delta(\epsilon)$ such
that

$\|y(0)\|<\delta\quad\mbox{implies that}\quad\|y(t)\|<\epsilon,\quad\forall
t\geq 0;$

  * •

unstable if it is not stable;

  * •

(locally) asymptotically stable if it is stable and $\delta$ can be chosen
such that

$\|y(0)\|<\delta\quad\mbox{implies that}\quad\lim_{t\to\infty}y(t)=0.$

The following theorem is known as the Lyapunov’s indirect method [49, 52, 53,
85].

######  Theorem 2.2.

Let $y^{*}=0$ be an equilibrium point for the system (1.1), where
$f:D\to\mathbb{R}^{n}$ is continuously differentiable and $D$ is a
neighborhood of the origin. Let

$A=\dfrac{\partial f}{\partial y}(y)\bigg{|}_{y=0}.$

Then,

  * •

$y^{*}$ is LAS if $Re\lambda_{i}<0$ for all eigenvalues of $A$;

  * •

$y^{*}$ is unstable if $Re\lambda_{i}>0$ for one or more the eigenvalues of
$A$.

######  Remark 2.3.

An equilibrium point $y^{*}$ of (1.1) is said to be hyperbolic if none of the
eigenvalues of the matrix $A$ lies on the imaginary axis and non-hyperbolic
otherwise (see [85]). Hence, Theorem 2.2 is only applicable for hyperbolic
equilibrium points.

######  Theorem 2.4 ([45, 84]).

Suppose that $f$ in (1.2) has the property that solutions of initial value
problems $y(0)=y_{0}\geq 0$ are unique and, for all $i$, $f_{i}(y)\geq 0$
whenever $y\geq 0$ satisfies $y_{i}=0$. Then $y(t)\geq 0$ for all $t\geq 0$
for which it is defined, provided $y(0)\geq 0$.

######  Remark 2.5.

We infer from Theorem 2.4 that the set $\mathbb{R}^{n}_{+}$ is a positively
invariant set of the model (1.1) if the condition (1.2) holds.

###  2.2. Conservation laws

We consider three classes of conservation laws for the model (1.1), that are
direct conservation, generalized conservation and sub-conservation laws [65].

For the system (1.1), let us denote by
$P_{T}(t)=y_{1}(t)+y_{2}(t)+\ldots+y_{n}(t)$. We say that the system (1.1)
satisfies

  * •

a direct conservation law (DCL) if

$\dfrac{dP_{T}}{dt}=0,$

which implies that $P_{T}$ is constant;

  * •

a generalized conservation law (GCL) if

$\dfrac{dP_{T}}{dt}=a_{1}-b_{1}P_{T},\quad a_{1},b_{1}>0,$

which implies that $P_{T}(t)$ monotonically converges to $a_{1}/b_{1}$ as
$t\to\infty$;

  * •

a sub-conservation law (SCL) if a sub-set of the populations, $P_{S}(t)$,
where

$P_{S}(t)=P_{1S}(t)+P_{2S}(t)+\ldots+P_{mS}(t),\quad 2\leq m\leq n,$

has the property

$\dfrac{dP_{T}}{dt}=0,$

or

$\dfrac{dP_{T}}{dt}=a_{1}-b_{1}P_{T}.$

In [65], Mickens and Washington introduced NSFD discretizations of some
mathematical models satisfying the above conservation laws.

###  2.3. Stability of equilibria of discrete-time dynamical systems

Consider a general discrete-time dynamical system governed by first-order
difference equations of the form

(2.1) $y_{k+1}=g(y_{k}),\quad y_{0}\in\mathbb{R}^{n},$

where $g:D\to\mathbb{R}$ and $D\subset\mathbf{R}^{n}$ is the domain of the
function $g$. The vector $y_{0}$ is called the initial value or initial data.
Note that the uniqueness of the sequence $\\{y_{k}\\}_{k=0}^{\infty}$ is
automatic if it exists.

A point $y^{*}\in\mathbf{R}^{n}$ is called a fixed point or an equilibrium
point of the system (2.1) if $g(y^{*})=y^{*}$ (see [34, 85]).

######  Definition 2.6 ([34, 85]).

Let $y^{*}$ be an equilibrium point of the system (2.1). Then, $y^{*}$ is said
to be:

  * •

stable if given $\epsilon>0$ there exists $\delta>0$ such that

$\|y_{-}y^{*}\|<\delta\quad\mbox{implies}\quad\|y_{k}-y^{*}\|<\epsilon,$

for all $k>0$.

  * •

unstable if it is not stable;

  * •

attracting if there exists $\eta>0$ such that

$\|y_{0}-y^{*}\|<\eta\quad\mbox{implies}\quad\lim_{k\to\infty}y_{k}=y^{*};$

  * •

(locally) asymptotically stable if it is stable and attracting.

######  Theorem 2.7 ([34, 85]).

Let $f\in\mathcal{C}^{2}(\mathbb{R}^{n},\mathbb{R}^{n})$. Then a fixed point
$y^{*}$ of the system (2.1) is locally asymptotically stable if the
eigenvalues of $df(y^{*})$ lie strictly inside the unit circle. If any of the
eigenvalues lie outside the unit circle the fixed point is unstable.

######  Remark 2.8.

A fixed point $y^{*}$ is said to be hyperbolic if none of the eigenvalues of
the matrix $df(y^{*})$ lies on the unit circle and non-hyperbolic if otherwise
(see [85]). Hence, Theorem 2.7 only applicable for hyperbolic equilibrium
points.

###  2.4. Nonstandard finite difference schemes

Consider a finite difference scheme that approximates solutions of (1.1) in
the form

(2.2) $D_{\Delta t}(y_{k})=F_{\Delta t}(f;y_{k}),$

where $D_{\Delta t}(y_{k})\approx dy/dt$, $F_{\Delta t}(f;y_{k})\approx f(y)$
and $t_{k}=k\Delta t$, $\Delta t$ is the step size.

According to the Mickens’ methodology [59, 60, 61, 62, 63, 64], an NSFD scheme
for the equation (1.1) is a discrete model constructed based on a set of six
rules. In particular, NSFD schemes for first-order differential equations can
be defined as follows [6, 7, 29].

######  Definition 2.9.

The finite difference scheme (2.2) is called an NSFD scheme if at least one of
the following conditions is satisfied:

  * •

$D_{\Delta t}(y_{k})=\dfrac{y_{k+1}-y_{k}}{\phi(\Delta t)}$, where
$\phi(\Delta t)=\Delta t+\mathcal{O}(\Delta t^{2})$ is a non-negative function
and is called a nonstandard denominator function;

  * •

$F_{\Delta t}(f;y_{k})=g(y_{k},y_{k+1},\Delta t)$, where
$g(y_{k},y_{k+1},\Delta t)$ is a non-local approximation of the right-hand
side of the system (1.1).

######  Definition 2.10 ([7, 8]).

Assume that the solutions of the equation (1.1) satisfy some property
$\mathcal{P}$. The numerical scheme (2.2) is called (qualitatively) stable
with respect to property $\mathcal{P}$ (or $\mathcal{P}$-stable), if for every
value of $\Delta t>0$ the set of solutions of (2.2) satisfies property
$\mathcal{P}$.

######  Definition 2.11 ([6, 7, 64]).

Consider the differential equation $dy/dt=f(y)$. Let a finite difference
scheme for the equation be $y_{k+1}=F(y_{k};\Delta t)$. Let the differential
equation and/or its solutions have property $\mathcal{P}$. The discrete model
equation is dynamically consistent with the differential equation if it and/or
its solutions also have property $\mathcal{P}$.

######  Definition 2.12.

(see [88]) The NSFD method (2.2) is called positive or positivity-preserving,
if, for any value of the step size $\Delta t$, and
$y_{0}\in\mathbb{R}_{+}^{n}$ its solution remains positive, i.e.,
$y_{k}\in\mathbb{R}_{+}^{n}$ for all $k\in\mathbb{N}$.

##  3\. Construction of the NSFD scheme

In this section, we construct a generalized NSFD scheme for the model (1.1)
and investigate its qualitative dynamical properties.

###  3.1. Derivation of the NSFD scheme and its positivity

First, suppose that the function $f$ satisfies the condition (C1) with a real
number $\alpha$. We use the following non-local approximation for the function
$f$

(3.1) $f(y(t_{k}))=f(y(t_{k}))+\big{(}my(t_{k})-my(t_{k})\big{)}\approx
f(y_{k})+my_{k}-my_{k+1},$

where $m$ is a real number. By combining (3.1) with the nonstandard
approximation for the first-order derivative $dy/dt$, we obtain the following
NSFD model for the system (1.1):

(3.2) $\dfrac{y_{k+1}-y_{k}}{\phi(\Delta t)}=f(y_{k})+my_{k}-my_{k+1},$

where $\phi(\Delta t)=\Delta t+\mathcal{O}(\Delta t^{2})$ as $\Delta t\to 0$.
The explicit form of the scheme (3.2) is given by

(3.3) $y_{k+1}=\dfrac{y_{k}+\phi f(y_{k})+\phi my_{k}}{1+m\phi},$

or equivalently,

(3.4) $y_{k+1}=y_{k}+\dfrac{\phi}{1+m\phi}f(y_{k}).$

######  Remark 3.1.

Although the rule of nonlocal approximation of the NSFD methodology states
that nonlinear terms in the right-hand side $f(y)$ should be approximated in a
nonlocal way, in many cases linear terms with negative coefficients, namely,
$-cy$ ($c>0$) can be also nonlocally approximated by $-cy_{k+1}$ to guarantee
the positivity of the corresponding NSFD schemes (see, for instance, [1, 8,
12, 22, 43, 44, 55, 56, 65]). Hence, the approximation (3.1) can be considered
as a nonlocal approximation, which follows the NSFD methodology.

######  Lemma 3.2.

The NSFD sheme (3.2) preserves the positivity of the model (1.1) for all the
values of the step size $\Delta t$ whenever

(3.5) $m\geq m_{P}:=\max\\{\alpha,\,0\\}.$

###### Proof.

We use mathematical induction to prove the positivity of the NSFD scheme
(3.2). Indeed, assume that $y_{k}\geq 0$. Then,

$f(y_{k})+my_{k}\geq f(y_{k})+\alpha y_{k}\geq 0.$

Combining this estimate with (3.3), we obtain

$y_{k+1}=\dfrac{y_{k}+\phi(f(y_{k})+my_{k})}{1+m\phi}\geq 0.$

This is the desired conclusion. The proof is complete. ∎

######  Remark 3.3.

The discretization (3.1) guarantees that the NSFD scheme (3.2) is consistent.
Similarly to convergence analysis presented in [19], we can conclude that the
NSFD scheme (3.2) is convergent of order $1$.

######  Remark 3.4.

It is easy to see from (3.4) that the NSFD scheme (3.2) can be rewritten in
the form

$\dfrac{y_{k+1}-y_{k}}{\Phi(\Delta t)}=f(y_{k}),$

where

$\Phi(\Delta t)=\dfrac{\phi(\Delta t)}{1+m\phi(\Delta t)}=\phi(\Delta
t)+\mathcal{O}(\phi^{2}(\Delta t))=\Delta t+\mathcal{O}(\Delta t^{2}).$

Hence, the NSFD scheme (3.2) can be equivalently derived by using
renormalization of the denominator. Here, it is important to note that
conditions for the NSFD scheme (3.2) to be dynamically consistent with the
model (1.1) will be set to $m$.

###  3.2. Equilibria and their asymptotic stability

Let us denote by $\mathcal{E}_{C}$ and $\mathcal{E}_{D}$ the sets of
equilibrium points of the continuous-time dynamical system (1.1) and the
discrete-time model (3.2), respectively. It should be emphasized that
$\mathcal{E}_{C}\subset\mathcal{E}_{D}$ for most standard finite difference
schemes such as the Taylor’s and Runge-Kutta schemes since they often generate
spurious equilibrium points, which depend on the step size [59, 60, 61, 62,
63, 64]. However, we infer from the system (3.4) that:

######  Lemma 3.5.

$\mathcal{E}_{C}=\mathcal{E}_{D}$ for all the finite values of the step size
$\Delta t$. In other words, the NSFD scheme (3.2) preserves the set of
equilibrium points of the model (1.1).

Let us denote by $\mathcal{E}_{C}^{S}$ the set of all asymptotically stable
equilibrium points of the system (1.1). To analyze the stability of the NSFD
scheme (3.2), we introduce the following notation

(3.6)
$\begin{split}&\mathcal{F}=\bigcup_{y^{*}\in\mathcal{E}_{C}^{S}}\sigma(J_{C}(y^{*})),\\\
&m_{S}=\max_{\lambda\in\mathcal{F}}\bigg{\\{}-\dfrac{|\lambda|^{2}}{2Re\lambda}\bigg{\\}},\end{split}$

where $J_{C}(y^{*})$ is the Jacobian matrix of the system (1.1) evaluating at
$y^{*}$ and $\sigma(J_{C}(y^{*}))$ is the set of eigenvalues of
$J_{C}(y^{*})$. Note that $Re\lambda<0$ for all $\lambda\in\mathcal{F}$ (see
Theorem 2.2).

######  Theorem 3.6.

Let $m$ be a real number with the property that

(3.7) $m\geq m_{S},$

where $m_{S}$ is given in (3.6). Then, the NSFD scheme (3.2) preserves the LAS
of the system (1.1) for all the values of the step size $\Delta t$.

###### Proof.

Let $y^{*}$ be any equilibrium point of the continuous system (1.1). Let us
denote by ${J}_{C}(y^{*})$ and ${J}_{D}(y^{*})$ the Jacobian matrices of the
systems (1.1) and (3.2) evaluating at $y^{*}$, respectively. Then, it follows
from the equation (3.4) that

(3.8) ${J}_{D}=I+\dfrac{\phi}{1+m\phi}{J}_{C},$

where $I$ is the identity matrix. The relation (3.8) implies that

$\begin{split}\det(xI-{J}_{C})&=\det\bigg{[}\bigg{(}x+\dfrac{1+m\phi}{\phi}\bigg{)}I-\frac{1+m\phi}{\phi}J_{D}\bigg{]}\\\
&=\bigg{(}\dfrac{1+m\phi}{\phi}\bigg{)}^{n}\det\bigg{[}\bigg{(}1+\dfrac{\phi
x}{1+m\phi}\bigg{)}I-{J}_{D}\bigg{]}.\end{split}$

Hence, $\lambda$ is an eigenvalue of ${J}_{C}$ if and only if
$\mu=1+\dfrac{\phi}{1+m\phi}\lambda$ is an eigenvalue of ${J}_{D}$. This
implies that

(3.9)
$\begin{split}|\mu_{i}|^{2}&=\bigg{|}1+\dfrac{\phi}{1+m\phi}\lambda_{i}\bigg{|}^{2}=\bigg{(}1+\dfrac{\phi}{1+m\phi}Re\lambda_{i}\bigg{)}^{2}+\bigg{(}\dfrac{\phi}{1+m\phi}Im\lambda_{i}\bigg{)}^{2}\\\
&=1+\dfrac{2\phi}{1+m\phi}Re\lambda_{i}+\bigg{(}\dfrac{\phi}{1+m\phi}\bigg{)}^{2}|\lambda_{i}|^{2}.\end{split}$

We now consider the following two cases of the stability of $y^{*}$.
Case 1. $y^{*}$ is an unstable equilibrium point of the system (1.1). In this
case, it follows from Theorem 2.2 that there exists
$\lambda_{i}\in\sigma({J}_{C}(y^{*}))$ such that $Re\lambda_{i}>0$. Then, we
infer from (3.9) that the corresponding eigenvalue $\mu_{i}$ of
${J}_{D}(y^{*})$ satisfies

$|\mu_{i}|^{2}=1+\dfrac{2\phi}{1+m\phi}Re\lambda_{i}+\bigg{(}\dfrac{\phi}{1+m\phi}\bigg{)}^{2}|\lambda_{i}|^{2}>1,$

which implies that $y^{*}$ is an unstable equilibrium point of the system
(3.2) (see Theorem 2.7).
Case 2. $y^{*}$ is an asymptotically stable equilibrium point of the system
(1.1). In this case, $Re\lambda<0$ for all $\lambda\in\sigma({J}_{C}(y^{*}))$.
By using (3.9), we obtain $|\mu|^{2}<1$ if and only if

$\big{(}2mRe\lambda+|\lambda|^{2}\big{)}\phi<-2Re\lambda.$

This condition is always satisfied if $m\geq-{|\lambda|^{2}}/{(2Re\lambda)}$.
Hence, we obtain from the condition $m\geq m_{s}$ and Theorem 2.7 that $y^{*}$
is an asymptotically stable equilibrium point of the NSFD scheme (3.2).

Finally, the proof is completed by combining Cases $1$ and $2$. ∎

######  Remark 3.7.

In many cases, it is not easy to obtain conditions for NSFD schemes to be
dynamically consistent with respect to the asymptotic stability of dynamical
systems (see, for instance, [19, 22, 24, 26, 30, 43, 44, 87]). However, we can
obtain the stability threshold $m_{S}$ of the NSFD scheme (3.2) by simple
calculations or a simple numerical algorithm.

###  3.3. Conservation laws

Assume that the system (1.1) satisfies conservation laws. We show that the
NSFD scheme (3.2) has the ability to preserve this property.

######  Lemma 3.8.

The NSFD scheme (3.2) preserves the DCL for all finite step sizes.

###### Proof.

Assume that the system (1.1) satisfies a DCL, that is the total population
satisfies

$\dfrac{dP_{T}}{dt}=\dfrac{dy_{1}}{dt}+\dfrac{dy_{2}}{dt}+\ldots+\dfrac{dy_{n}}{dt}=0,$

which implies that $\sum_{i=1}^{n}y_{i}(t)$ is constant, and hence,
$\sum_{i=1}^{n}f_{i}(y(t))=0$ for all $t\geq 0$. By using (3.4) we obtain

$P_{T,k+1}:=\sum_{i=1}^{n}y_{i,k+1}=\sum_{i=1}^{n}y_{i,k}+\dfrac{\phi}{1+m\phi}\sum_{i=1}^{n}f_{i}(y_{k})=\sum_{i=1}^{n}y_{i,k+1}:=P_{T,k}.$

This implies that the DCL is preserved automatically by the NSFD scheme (3.2).
The proof is complete. ∎

Next, assume that a GCL exists for the system (1.1), that is

(3.10) $\dfrac{dP_{T}}{dt}=a_{1}-b_{1}P_{T},\quad a_{1},b_{1}>0.$

Note that (3.10) has an exact finite difference scheme of the form (see [65])

(3.11) $\dfrac{P_{T,k+1}-P_{T,k}}{\phi^{*}(\Delta
t)}=a_{1}-b_{1}P_{T,k},\quad\phi^{*}(\Delta t)=\dfrac{1-e^{-b_{1}\Delta
t}}{b_{1}}.$

Applying the NSFD scheme (3.2) we obtain

(3.12) $P_{T,k+1}=P_{T,k}+\dfrac{\phi}{1+m\phi}(a_{1}-b_{1}P_{T,k}).$

######  Lemma 3.9.

The following assertions hold for the difference equation (3.12):

  1. (1)

the sequence $\\{P_{T,k}\\}_{k}$ defined in (3.12) monotonically converges to
${a_{1}}/{b_{1}}$ as $k\to\infty$, provided that

(3.13) $m\geq m_{GCL}:=b_{1};$

  2. (2)

The scheme (3.12) is an exact finite difference scheme for (3.10) if

(3.14) $\phi(\Delta t)=\phi_{E}(\Delta t):=\dfrac{\phi^{*}(\Delta
t)}{1-m\phi^{*}(\Delta t)},$

where $\phi^{*}$ is given in (3.11).

###### Proof.

Proof of Part (1). It follows from (3.12) that

$\begin{split}P_{T,k+1}&=P_{T,k}+\dfrac{\phi}{1+m\phi}(a_{1}-b_{1}P_{T,k})=\dfrac{a_{1}\phi}{1+m\phi}+\bigg{(}1-\dfrac{b_{1}\phi}{1+m\phi}\bigg{)}P_{T,k}\\\
&=\dfrac{a_{1}\phi}{1+m\phi}+\bigg{(}1-\dfrac{b_{1}\phi}{1+m\phi}\bigg{)}\Bigg{[}\dfrac{a_{1}\phi}{1+m\phi}+\bigg{(}1-\dfrac{b_{1}\phi}{1+m\phi}\bigg{)}P_{T,k-1}\Bigg{]}\\\
&=\dfrac{a_{1}\phi}{1+m\phi}+\dfrac{a_{1}\phi}{1+m\phi}\bigg{(}1-\dfrac{b_{1}\phi}{1+m\phi}\bigg{)}+\bigg{(}1-\dfrac{b_{1}\phi}{1+m\phi}\bigg{)}^{2}P_{T,k-1}.\end{split}$

By continuing the above process, we have

$P_{T,k+1}=\dfrac{a_{1}\phi}{1+m\phi}{\sum_{i=0}^{k}}\bigg{(}1-\dfrac{b_{1}\phi}{1+m\phi}\bigg{)}^{i}+\bigg{(}1-\dfrac{b_{1}\phi}{1+m\phi}\bigg{)}^{k}P_{T,0}.$

Since $m\geq b_{1}$, $1-\dfrac{b_{1}\phi}{1+m\phi}\in(-1,1)$. This implies
that

(3.15)
$\lim_{t\to\infty}P_{T,k}=\dfrac{a_{1}\phi}{1+m\phi}\bigg{(}\dfrac{b_{1}\phi}{1+m\phi}\bigg{)}^{-1}=\dfrac{a_{1}}{b_{1}}.$

On the other hand, we infer from $m\geq b_{1}$ that

$\dfrac{\partial P_{T,k+1}}{\partial
P_{T,k}}=1+\dfrac{\phi}{1+m\phi}(-b_{1})\geq 0,$

which implies the monotonocity of $\\{P_{T,k}\\}_{k}$. The proof is completed.
Proof of Part (2). The exact scheme (3.11) can be written in the form

(3.16) $P_{T,k+1}=P_{T,k}+\phi^{*}(a_{1}-b_{1}P_{T,k}).$

We infer from (3.16) and (3.12) that (3.12) is an exact scheme for (3.10) if

$\phi^{*}=\dfrac{\phi}{1+m\phi}.$

This condition is equivalent to (3.14). The proof is complete. ∎

From Lemma 3.9 we have:

######  Theorem 3.10.

The NSFD scheme (3.2) satisfies the GCL if one of the following conditions
holds:

  1. (1)

$m\geq m_{GCL}:=b_{1}$.

  2. (2)

$\phi(\Delta t)=\phi_{E}(\Delta t)$, where $\phi_{E}$ is given in (3.14).

######  Remark 3.11.

The function $\phi_{E}$ given in (3.14) satisfies $\phi(\Delta t)=\Delta
t+\mathcal{O}(\Delta t^{2})$ and is defined for all $\Delta t>0$ if $b_{1}\geq
m$. When $b_{1}<m$, it is defined only if $\Delta
t\neq-\dfrac{\ln(m-b_{1})}{b_{1}}$.

######  Remark 3.12.

The sub-conservation law will be satisfied if the DCL and GCL are satisfied.

###  3.4. Dynamics consistency

For the sake of convenience, we summarize the results constructed in
Subsections 3.1-3.4 in the following theorem.

######  Theorem 3.13 (Dynamic consistency).

The NSFD scheme (3.2) is

  * •

dynamically consistent with the positivity of the model (1.1) if the condition
(3.5) is satisfied;

  * •

dynamically consistent with the LAS of the model (1.1) if the condition (3.7)
is satisfied;

  * •

dynamically consistent with the DCL of the model (1.1) for all the values of
$m$;

  * •

dynamically consistent the GCL the model (1.1) if

(3.17) $m\geq m_{GCL}:=b_{1}$

or

(3.18) $\phi(\Delta t)=\phi_{E}(\Delta t),$

where $\phi_{E}$ is given in (3.14).

######  Remark 3.14.

  * •

The positivity, LAS and conservation laws are independent of each other. So,
they may not occur simultaneously. Therefore, the threshold parameter $m$
should be selected depending on properties that the model (1.1) can have.

  * •

The determination of $m_{P}$ in (3.5) and $m_{GCL}$ in (3.17) is trivial while
it is also easy to write a numerical algorithm to compute $m_{S}$.

  * •

Since the conditions for the NSFD scheme (3.2) to be dynamically consistent
with respect to the positiivity, LAS and conservation laws of the model (1.1)
are set to the parameter $m$, we can choose the denominator function
$\phi(\Delta t)$ to be $\Delta t$. This shows the role of the parameter $m$
for the dynamic consistency of the NSFD scheme (3.2).

###  3.5. A note on the nonstandard Euler scheme

As emphasized in Remark 3.4, the NSFD scheme (3.2) can be considered as a
nonstandard Euler method. For the sake of convenience, we now consider a
special case of the parameter $m$, namely, $m=0$. In this case, the NSFD
scheme (3.2) becomes

(3.19) $y_{k+1}=y_{k}+\phi(\Delta t)f(y_{k}).$

This can be considered as a ”weakly” NSFD scheme [28].

Because of the absence of parameter $m$, conditions for the dynamic
consistency of the scheme (3.19) will be set to the denominator function
$\phi$. Below, we will give conditions for the scheme (3.19) to be dynamically
consistent with the model (1.1).

We infer from (3.19) that $y_{k+1}\geq 0$ for all $\phi>0$ if $\alpha\leq 0$.
In this case, the nonstandard Euler scheme is unconditionally positive.
Otherwise, if $\alpha>0$, then it follows from (1.2) that

$y_{k+1}=y_{k}+\phi f(y_{k})\geq y_{k}+(-\alpha y_{k})=(1-\phi\alpha)y_{k}.$

Consequently, $y_{k+1}\geq 0$ if $\phi<(\alpha)^{-1}$. So, we have the
following result on the positivity of the nonstandard Euler scheme.

######  Lemma 3.15.

The nonstandard Euler scheme (3.19) preserves the positivity of the model
(1.1) for all the values of the step size $\Delta t$ whenever the denominator
function satisfies the following condition for all $\Delta t>0$

(3.20) $\phi(\Delta
t)<\phi_{P}:=\begin{cases}&\infty,\,\,\,\quad\qquad\alpha\leq 0,\\\
&(\alpha)^{-1},\qquad\alpha>0.\end{cases}$

The following result is a direct consequence of the results presented in [28,
29].

######  Corollary 3.16.

The nonstandard Euler scheme (3.19) preserves the LAS of the model (1.1) for
all the values of the step size $\Delta t$ if the denominator function
satisfies the following condition for all $\Delta t>0$

(3.21) $\phi(\Delta
t)<\phi_{S}:=\min_{\lambda\in\mathcal{F}}\bigg{\\{}-\dfrac{2Re\lambda}{|\lambda|^{2}}\bigg{\\}},$

where the set $\mathcal{F}$ is defined in (3.6).

It is easily seen that the nonstandard Euler scheme always preserves the DCL.
On the other hand, when applying the nonstandard Euler scheme for systems
satisfying the GLC, we obtain

(3.22) $\dfrac{P_{T,k+1}-P_{T,k}}{\varphi(\Delta t)}=a_{1}-b_{1}P_{T,k}.$

Hence, the GLC is preserved if

(3.23) $\phi<\phi_{GCL}=(b_{1})^{-1}.$

Furthermore, the scheme (3.22) becomes an exact scheme for (3.10) if

(3.24) $\phi(\Delta t)=\dfrac{1-e^{-b_{1}\Delta t}}{b_{1}}.$

Hence, we obtain the following assertion.

######  Corollary 3.17.

The nonstandard Euler scheme (3.19) is always dynamically consistent with
respect to the DCL and is dynamically consistent with respect to the GCL if
the condition (3.23) or the condition (3.24) holds.

######  Remark 3.18.

  * •

By choosing suitable denominator functions, the nonstandard Euler scheme
(3.19) can preserve the positivity, LAS and conservation laws as the NSFD
scheme (3.2). Although the schemes (3.2) and (3.19) are all convergent of
order $1$, in the next section we will provide a numerical example in which
errors of Scheme (3.2) are better than errors of Scheme (3.19).

  * •

The condition (3.24) can conflict with the conditions (3.20) and (3.21).

##  4\. Numerical experiments

In this section, we perform some numerical examples to illustrate the
theoretical findings and to show advantages of the proposed NSFD scheme.

######  Example 4.1 (Numerical simulation of single-species population
models).

Consider single-species population models, which describe the population size
changes in the absence of disease and of maturation delay [18]

(4.1) $\dfrac{dN}{dt}=B(N)N-dN,$

where $N(t)$ stands for the population size at the time $t$, $d>0$ is the
death rate constant, and $B(N)N$ is a birth rate function with $B(N)$ has the
following properties for all $N\in(0,\infty)$
(A1) $B(N)>0$;
(A2) $B(N)$ is continuously differentiable with $B^{\prime}(N)<0$;
(A3) $B(0^{+})>d>B(\infty)$.
The following functions are examples of birth functions $B(N)$ that are found
in biological systems and satisfy (A1)-(A3) (see [18])
(B1) $B_{1}(N)=be^{-aN}$, with $a>0$ and $b>d$;
(B2) $B_{2}(N)=\dfrac{p}{q+N^{n}}$, with $p,q,n>0$ and $p>qd$;
(B3) $B_{3}(N)=\dfrac{A}{N}+c$, with $A>0$, $d>c>0$.
It is clear that the model (4.1) satisfies the condition (C1) with $\alpha=d$.
On the other hand, it possesses two equilibrium points, that are $N^{0}=0$ and
$N^{*}=B^{-1}(d)$. Moreover, $N^{0}$ is always unstable while $N^{*}$ is
always asymptotically stable. Therefore, the condition (C2) also satisfies. If
$B(N)=B_{3}(N)$, then the model (4.1) becomes linear and satisfies the GCL.

We consider the model (4.1) with

$B(N)=B_{1}(N)=be^{-aN}N-dN,\quad a=d=1,\quad b=12$

as a test problem for numerical schemes. In this case, the model has a unique
equilibrium point

$N^{*}=-{\ln(d/b)}/{a}\approx 2.4849,$

which is locally asymptotically stable. Note that the conservation laws are
absent in this example.

For the NSFD scheme (3.2), we compute the parameter $m$ as

$m\geq\max\big{\\{}m_{P},\,\,m_{S}\big{\\}}=\max\\{1,\,\,1.2425\\}=1.2425.$

For the nonstandard Euler scheme, the condition for the denominator function
can be computed as

$\phi(\Delta t)<\min\\{\phi_{P},\,\,\phi_{S}\\}=\min\\{1,\,\,0.8049\\}=0.8049$

Numerical solutions generated by the NSFD scheme (3.2) with $m=2.5$ and
$\Delta t\in\\{2.0,\,0.1\\}$ and by the nonstandard Euler scheme with
$\phi(\Delta t)=({1-e^{-2.5\Delta t}})/{2.5}$ and $\Delta t\in\\{1,\,0.1\\}$
are depicted in Figure 1. In this figure, each blue curve represents a
solution corresponding to a specific initial value and the green arrows
represents the evolution of the continuous model. We observe from Figure 1
that the numerical solutions generated by the NSFD model (3.2) are positive
and stable. Furthermore, the NSFD model replicates the dynamics of the
continuous model and its behaviour does not depend on the chosen step sizes.
Similarly, the nonstandard Euler scheme (3.19) also preserves the positivity
and stability of the continuous model regardless of the chosen step sizes.

We now use the standard Euler scheme and a well-known two-stage Runge-Kutta
(RK2) scheme, namely, the explicit trapezoidal scheme (see [13]) for the
purpose of comparing them with the NSFD schemes. Approximations obtained from
these standard schemes are sketched in Figure 2. It is clear that the standard
Euler scheme generates an unstable numerical solution, which oscillates around
the equilibrium position. Hence, the dynamics of the continuous model cannot
be preserved. The similar unusual phenomena occurs when using the standard RK2
scheme with step sizes $\Delta t=1$ and $\Delta t=1.2$. More clearly, the RK2
scheme using $\Delta t=1.2$ generates the approximation oscillating around the
equilibrium point and when $\Delta t=1.0$, the corresponding numerical
solution converges to a spurious equilibrium point. This result is completely
consistent with analysis performed by Mickens in [59, 60, 61, 62, 63].

Before ending this example, we compute errors obtained from the NSFD scheme
(3.2), the nonstandard Euler scheme (3.19) and the standard Euler scheme. Note
that theses schemes are all convergent of order $1$ since $\phi(\Delta
t)=\Delta t+\mathcal{O}(\Delta t^{2})$. For this purpose, we consider $N(0)=2$
and use the approximation generated by the classical four-stage Runge-Kutta
scheme (see [13]) with $\Delta t=10^{-5}$ as a benchmark (reference) solution
since it is impossible to find the exact solution. Errors provided by the
three numerical schemes are shown in Tables 1 and 2. In these tables, $error$
and $error_{T}$ are computed by

$error=\max_{1\leq k\leq N}|N(t_{k})-N_{k}|,\quad
error_{T}=\sum_{i=1}^{N}|N(t_{k})-N_{k}|.$

From Tables 1 and 2, we see that the errors provided by the NSFD scheme (3.2)
are better. The errors of the nonstandard Euler scheme are better than the
ones provided by the standard Euler scheme. Hence, advantages of the
nonstandard schemes in term of errors are shown.

(a)

(b)

(c)

(d)

Figure 1. Numerical approximations of the model (4.1) generated by the NSFD
scheme (3.2) with $m=2.5$ and $\Delta t\in\\{2.0,\,0.1\\}$ and by the
nonstandard Euler scheme (3.19) with $\phi(\Delta t)=(1-e^{-2.5\Delta t})/2.5$
and $\Delta t\in\\{1,0.1\\}$.

(a)

(b)

Figure 2. Numerical approximation of the model (4.1) generated by the standard Euler and RK2 schemes. Table 1. The values of the $error$ of the NSFD scheme (3.2) with $m_{s}=1.5$, the nonstandard Euler scheme with $\phi(h)=(1-e^{-1.5\Delta t})/1.5$ and the standard Euler scheme $\Delta t$ | NSFD $error$ | N. Euler $error$ | S. Euler $error$
---|---|---|---
$1$ | $0.053229803837052$ | $0.200391028658851$ | $1.895818024330166$
$10^{-1}$ | $0.002547712779198$ | $0.010718671314132$ | $0.025728835277831$
$10^{-2}$ | $3.474294297345359e-004$ | $9.794823732640623e-004$ | $0.002320869273385$
$10^{-3}$ | $3.568404980169859e-005$ | $9.709160784732163e-005$ | $2.298145180823497e-004$
$10^{-4}$ | $3.577852512925972e-006$ | $9.700655808853043e-006$ | $2.295902635651714e-005$
$10^{-5}$ | $3.578797715952931e-007$ | $9.699803174001431e-007$ | $2.295678681374369e-006$
Table 2. The values of the $error_{T}$ of the NSFD scheme (3.2) with $m_{s}=1.5$, the nonstandard Euler scheme with $\phi(h)=(1-e^{-1.5\Delta t})/1.5$ and the standard Euler scheme $\Delta t$ | NSFD $error_{T}$ | N. Euler $error_{T}$ | S. Euler $error_{T}$
---|---|---|---
$1$ | $0.056862533841534$ | $0.257861518971072$ | $15.530365282265436$
$10^{-1}$ | $0.030063654051754$ | $0.109862865600443$ | $0.256967636704650$
$10^{-2}$ | $0.039734990275286$ | $0.103390564043856$ | $0.247233541153000$
$10^{-3}$ | $0.040703612013705$ | $0.102744510888389$ | $0.246264374931792$
$10^{-4}$ | $0.040800485618428$ | $0.102679929285270$ | $0.246167506109657$
$10^{-5}$ | $0.040810166884441$ | $0.102673437145564$ | $0.246157825305468$

######  Example 4.2 (Numerical simulation of a predator-prey system).

Consider predator-prey model with a Beddington-DeAngelis functional response
and linear intrinsic growth rate of the prey population, which is derived from
a general predator-prey model proposed by DeAngelis et al in [27]. The model
under consideration is given by

(4.2) $\begin{split}&\dfrac{dx}{dt}=x-\dfrac{Axy}{1+x+y},\\\
&\dfrac{dy}{dt}=\dfrac{Exy}{1+x+y}-Dy,\end{split}$

where $x$ and $y$ are the prey and predator population sizes, respectively,
and all the parameters are assumed to be positive. We refer the readers to
[31] for details of this model.

Clearly, the model (4.2) satisfies the condition (C1) with
$\alpha=\max\\{A-1,D\\}$. On the other hand, stability analysis in [31]
implies that (C2) holds if $A\neq E$. Numerical simulation of the model (4.2)
by NSFD schemes was also examined in [29].

Clearly, the model (4.2) satisfies the condition (C1) with
$\alpha=\max\\{A-1,D\\}$. On the other hand, stability analysis in [31]
implies that (C2) holds if $A\neq E$. Numerical simulation of the model (4.2)
by NSFD schemes was also examined in [29].

We now consider the predator-prey model (4.2) with the following parameters

$A=6,\quad D=5,\quad E=7.$

In this case, the model has a unique positive equilibrium point
$(x^{*},y^{*})=(6,\,\,1.4)$ and the Jacobian matrix evaluating at
$(x^{*},y^{*})$ is

$J(x^{*},y^{*})=\begin{pmatrix}0.7143&-3.5714\\\ &\\\
0.3333&-0.8333\end{pmatrix}.$

Consequently, its eigenvalues are

$\lambda_{1,2}=-0.0595\pm 0.7692i,$

which implies that $(x^{*},y^{*})$ is locally asymptotically stable.

By simple calculations, we determine the condition for the parameter $m$ of
the NSFD scheme (3.2) as

$m\geq\max\\{m_{P},m_{S}\\}=5.$

Therefore, we can take $m=5.1$.

Numerical approximations provided by the Euler and RK2 schemes and the NSFD
scheme (3.2) using $(x(0),y(0))=(4.5,1)$ are depicted in Figure 3. We see that
the phase plane generated by the Euler scheme is unstable and spirals out from
the equilibrium position; the RK2 scheme provides an approximation which
oscillates around the equilibrium point; meanwhile, the approximations
obtained from the NSFD scheme (3.2) correctly preserves the dynamics of the
predator-prey model.

Note that the condition for the denominator function $\phi$ of the nonstandard
Euler scheme (3.19) is

$\phi<\min\\{\phi_{P},\phi_{S}\\}=\dfrac{1}{5}.$

Hence, we can take $\phi(\Delta t)=(1-e^{-\tau\Delta t})/\tau$, where
$\tau>5$. Then, the dynamic consistency of the nonstandard Euler scheme is
guaranteed.

(a)

(b)

(c)

(d)

Figure 3. Numerical approximations of the model (4.2) generated by the NSFD
scheme (3.2) and the standard Euler and RK2 schemes.

######  Example 4.3 (Numerical simulation of a vaccination model).

Consider an infectious disease with a preventive (prophylactic) vaccine and an
effective therapeutic treatment proposed by Gumel and Moghadas in [38]. This
model monitors the temporal dynamics of susceptible individuals ($S$),
vaccinated individuals ($V$) and infected individuals $(I)$ and is described
by

(4.3) $\begin{split}&\dfrac{dS}{dt}=\Pi-\dfrac{c\beta_{1}I}{1+I}S-\xi
S+\alpha^{*}I-\mu S,\\\ &\dfrac{dV}{dt}=\xi S-\dfrac{c\beta_{2}I}{1+I}V-\mu
V,\\\
&\dfrac{dI}{dt}=\dfrac{c\beta_{1}I}{1+I}S+\dfrac{c\beta_{2}I}{1+I}V-\alpha^{*}I-\mu
I.\end{split}$

Here, all the parameters are positive. Mathematical formulation and
qualitative study for the model (4.3) were performed in [38].

It is easy to check that the model (4.3) satisfies the condition (C1) with
$\alpha=\max\\{c\beta_{1}+\xi+\mu,c\beta_{2}+\mu,\alpha^{*}+\mu\\}$.
Meanwhile, the condition (C2) holds if the basic reproduction number (BRN) of
the model satisfies

$\mathcal{R}_{0}=\dfrac{c\beta_{1}\Pi}{(\mu+\alpha^{*})(\mu+\xi)}+\dfrac{c\beta_{2}\Pi\xi}{\mu(\mu+\alpha^{*})(\mu+\xi)}\neq
1.$

Also, the GCL is satisfied because

$\dfrac{d(S+V+I)}{dt}=\Pi-\mu(S+V+I).$

In this example, we consider the vaccination model (4.3) with the following
parameters

$\begin{split}&\Pi=700,\quad\beta_{1}=10^{-4},\quad\beta_{2}=10^{-6},\quad\mu=0.03\qquad(\mbox{see
\cite[cite]{[\@@bibref{}{Gumel1}{}{}]}}),\\\ &\alpha^{*}=0.95,\qquad
c=8,\qquad\xi=0.95,\qquad\qquad\qquad\quad\,\,(\mbox{assumed}).\end{split}$

In this case, the BRN is $\mathcal{R}_{0}=0.7677<1$ and the model has a unique
disease-free (uninfected) equilibrium point
$E^{0}=(S^{0},V^{0},I^{0})=(714.2857,\,\,2.2619\times 10^{4},\,\,0)$, which is
locally asymptotically stable (see [38]).

We now compute the parameter $m$ for the NSFD scheme (3.2). First, it is easy
to see that

$m_{P}=\max\\{c\beta_{1}+\xi+\mu,c\beta_{2}+\mu,\alpha+\mu\\}=0.98.$

Since the Jacobian matrix of the model evaluating at $E^{*}$ has three
eigenvalues, that are

$\lambda_{1}=-0.98,\quad\lambda_{2}=-0.03,\quad\lambda_{3}=-0.23,$

we obtain

$m_{S}=0.49.$

Note that the model satisfying a GCL, namely,

$\dfrac{d(S+V+I)}{dt}=\Pi-\mu(S+V+I).$

Consequently, the condition for the preservation of the GCL is

$m\geq m_{GCL}={\mu}=0.03.$

Therefore, the condition for the parameter $m$ is

$m\geq\max\\{m_{P},m_{S},m_{GCL}\\}=0.98.$

In a similar way, we can determine the condition for the nonstandard Euler
scheme as

$\phi(\Delta t)=\dfrac{1-e^{-\tau\Delta t}}{\tau},\quad\tau>(0.98)^{-1}\approx
1.02.$

Numerical solutions obtained from the standard Euler, RK2 and NSFD scheme are
depicted in Figure 4, respectively. In these figures, each blue curve
represents a phase space corresponding to a specific initial value, the green
arrows show the evolution of the continuous model and the red circle indicates
the position of the stable equilibrium point.
We observe from Figures 4(a) and 4(b) that the standard Euler scheme generates
a negative and unstable numerical approximation, which oscillates around the
equilibrium point with an increasing amplitude; the solution provided by the
standard RK2 scheme is not positive and converges to a spurious equilibrium
point; conversely, the NSFD scheme correctly preserves the dynamics of the
continuous model regardless of the chosen step sizes (see Figures 4(c) and
4(d)).

Note that the vaccination model satisfies a GCL. For this reason, we now
consider approximations of the total population generated by the Euler, RK2
and NSFD scheme (3.2). From numerical approximations in Figure 5, we see that
the NSFD schemes preserves the GCL but this property is not preserved by the
standard schemes.

(a)

(b)

(c)

(d)

Figure 4. Numerical approximations of the model (4.3) generated by the
standard Euler and RK2 schemes and by the NSFD scheme (3.2).

(a)

(b)

(c)

Figure 5. Approximations of the total population generated by the standard
Euler and RK2 schemes and the NSFD scheme (3.2) using
$(S(0),\,V(0),\,I(0))=(5000,\,20000,\,1000)$.

##  5\. Conclusions, discussions and some open problems

In this work, we have applied the NSFD methodology proposed by Mickens to
construct a generalized NSFD method for a class of autonomous dynamical
systems that satisfy positivity, stability and conservation laws. This NSFD
method is constructed based on a novel non-local approximation for the right-
side functions of the dynamical systems. By rigorous mathematical analyses, we
have shown that the NSFD method is dynamically consistent with respect to the
positivity, asymptotic stability and three classes of conservation laws,
including direct conservation, generalized conservation and sub-conservation
laws. Additionally, a set of numerical experiments has been also performed to
illustrate the theoretical findings and to show advantages of the proposed
NSFD method over some well-known standard numerical schemes. The numerical
results are not only consistent with the theoretical assertions but also
indicate that the NSFD method is easy to be implemented and can be used to
solve a broad range of mathematical models arising in real-life.

It worth noting that the threshold parameters $(m_{S},m_{P},m_{GCL})$ for the
NSFD method (3.2) and $(\phi_{P},\phi_{S},\phi_{GCL})$ for the nonstandard
Euler method (3.19) can be computed easily. More clearly, these threshold
parameters can be obtained by simple calculations or simple numerical
algorithms.

Similarly to the NSFD methods formulated in [7, 28, 29], the NSFD method (3.2)
also provides a stability-preserving numerical scheme for autonomous dynamical
systems in general contexts. However, the positivity of the NSFD schemes in
[7, 28, 29] is not guaranteed automatically while the NSFD method (3.2) has
the ability to preserve not only the positivity but also the conservation laws
of the system (1.1).

In [88], Wood and Kojouharov introduced a nonstandard numerical method
preserving the positivity and LAS of autonomous dynamical systems. This method
is computationally efficient, easy to be implemented and can be used to solve
a broad range of problems in science and engineering. However, it requires the
determination of the sign of the right-hand side functions at each iteration,
and hence, the computational time is increased. For the NSFD method (3.2), we
only need to compute the parameter $m$ exactly once and then use it for all
iterations. Furthermore, the computation of $m$ is easy and can be done by a
simple algorithm.

Before ending this section, we mention some open problems and future
directions that research in NSFD methods can take.

  1. (1)

The NSFD method (3.2) requires the determination of the number $\alpha$ in the
condition (1.2). Although $\alpha$ can be determined easily for most of the
models in literature, it is necessary to construct an algorithm for computing
$\alpha$ in general cases. Even though many important mathematical models in
real-word applications satisfy the condition (1.2), this condition leads to a
big restriction of the application of the proposed NSFD method. Hence, it is
important to study extensions of the NSFD method (3.2) for mathematical models
having the positivity and asymptotic stability but not satisfying Condition
(1.2).

  2. (2)

NSFD methods preserving the asymptotic stability of non-hyperbolic equilibrium
points: Most stability-preserving NSFD methods including the NSFD method (3.2)
require equilibrium points of dynamical systems to be hyperbolic. The
construction of generalized NSFD methods preserving the LAS of dynamical
systems even when they have non-hyperbolic equilibrium points has been studied
in some previous work (see, for instance, [11]) but is still an important open
problem.

  3. (3)

Global asymptotic stability (GAS) of NSFD methods: It is well-known that
analyzing the GAS of NSFD schemes for mathematical models having the GAS
property is very important but not a simple task in general. This problem has
been mentioned and considered in many previous works (see, for instance [3, 8,
22, 23, 24, 25, 26, 42, 81, 82]). We observe from the numerical examples in
Section 4 that the NSFD method can preserve not only the LAS but also the GAS
of the continuous models under consideration. Therefore, it is suitable to
predict that there exists a positive number $m_{GAS}$ that plays as a GAS
threshold of the NSFD method (3.2), i.e., the NSFD method (3.2) preserves the
GAS of the system (1.1) whenever $m\geq m_{GAS}$. This prediction will be
useful in studying the construction of NSFD methods preserving the GAS of
dynamical systems.

  4. (4)

High-order NSFD methods: High-order NSFD methods for differential equations
have attracted the attention of many researchers with various approaches (see,
for example [17, 21, 37, 40, 50, 55, 56]); hence, the improvement accuracy of
the NSFD method (3.2) is essential with many practice applications. A simple
approach for improving accuracy of the NSFD method (3.2) is to combine it with
Richardson’s extrapolation method or variable step size strategies. This is
easy to do because (3.2) is simple and fully explicit.

  5. (5)

New NSFD methods for differential equation models arising in real-world
applications: Although NSFD methods have been strongly developed and achieved
many important results, fields of theory and practice always pose new and
complex problems; hence, constructing effective NSFD methods for these
problems is very essential. Recently, NSFD methods for fractional differential
equations have been widely studied (see, for instance [12, 20]). It was proved
that NSFD methods for fractional-order systems are efficient and powerful;
however, except for the positivity, other properties of fractional-order NSFD
methods have not been widely discussed. It is easy to extend the NSFD method
(3.2) for systems in the context of fractional derivative operators. Then, the
asymptotic stability in particular and qualitative study in general for the
fractional-order NSFD schemes will be posed. It is appropriate to predict that
the parameter $m$ under suitable conditions will guarantee the dynamic
consistency of the proposed NSFD schemes.

## Acknowledgments

We would like to thank the editor and anonymous referees for useful and
valuable comments that led to a great improvement of the paper.

## References

  * [1] E. M. Adamu, K. C. Patidar, and A. Ramanantoanina, An unconditionally stable nonstandard finite difference method to solve a mathematical model describing Visceral Leishmaniasis, Math. Comput. Simul. 187 (2021) 171-190.
  * [2] O. Adekanye, and T. Washington, Nonstandard finite difference scheme for a Tacoma Narrows Bridge model, Appl. Math. Model. 62 (2018) 223-236.
  * [3] H. Al-Kahby, F. Dannan and S. Elaydi, Non-standard Discretization Methods for Some Biological Models, Applications of Nonstandard Finite Difference Schemes, pp. 155-180 (2000).
  * [4] L. J. S. Allen, Introduction to Mathematical Biology, Pearson Education, Inc, Prentice Hall, 2007.
  * [5] P. Amarasekare and H. Possingham, Patch dynamics and metapopulation theory: the case of successional species, J. Theor. Biol. 209 (2001)333-344.
  * [6] R. Anguelov and J. M. -S. Lubuma, Nonstandard finite difference method by nonlocal approximation, Math. Comput. Simul. 61 (2003) 465-475.
  * [7] R. Anguelov and J. M.-S. Lubuma, Contributions to the mathematics of the nonstandard finite difference method and applications, Numer. Methods Partial Differ. Equ. 17 (2001) 518-543.
  * [8] R. Anguelov, Y. Dumont, J. M. -S. Lubuma, and M. Shillor, Dynamically consistent nonstandard finite difference schemes for epidemiological models, J. Comput. Appl. Math. 255 (2014) 161-182.
  * [9] R. Anguelov, J. M.-S. Lubuma, and S. K. Mahudu, Qualitatively stable finite difference schemes for advection-reaction equations, J. Comput. Appl. Math. 158 (2003) 19-30.
  * [10] R. Anguelov, P. Kama, and J. M. -S. Lubuma, On non-standard finite difference models of reaction-diffusion equations, J. Comput. Appl. Math. 175 (2005) 11-29.
  * [11] R. Anguelov, J. M. -S. Lubuma, and M. Shillor, Topological dynamic consistency of non-standard finite difference schemes for dynamical systems, J. Differ. Equ. Appl. 17 (2011) 1769-1791.
  * [12] A. J. Arenas, G. Gonzalez-Parra, and B. M. Chen-Charpentier, Construction of nonstandard finite difference schemes for the SI and SIR epidemic models of fractional order, Math. Comput. Simul. 121(2016) 48-63.
  * [13] U. Ascher and L. Petzold, Computer methods for ordinary differential equations and differential-algebraic equations, Society for Industrial and Applied Mathematics, Philadelphia, 1998.
  * [14] F. Brauer and C. Castillo-Chavez, Mathematical Models in Population Biology and Epidemiology, Second Edition, Springer New York, NY, 2012.
  * [15] V. Capasso and G. Serio, A generalization of the Kermack-McKendrick deterministic epidemic model, Math. Biosci. 42 (1978) 43-61.
  * [16] M. Chapwanya, J. M.-S. Lubuma, and R. E. Mickens, Positivity-preserving nonstandard finite difference schemes for cross-diffusion equations in biosciences, Comput. Math. Appl. 68 (2014) 1071-1082.
  * [17] B. M. Chen-Charpentier, D. T. Dimitrov, and H. V. Kojouharov, Combined nonstandard numerical methods for ODEs with polynomial right-hand sides, Math. Comput. Simul. 73 (2006) 105-113.
  * [18] K. Cooke, P. van den Driessche, and X. Zou, Interaction of maturation delay and nonlinear birth in population and epidemic models, J. Math. Biol. 39 (1999) 332-352.
  * [19] J. Cresson and F. Pierret, Non standard finite difference scheme preserving dynamical properties, J. Comput. Appl. Math. 303 (2016) 15-30.
  * [20] J. Cresson and A. Szafrańska, Discrete and continuous fractional persistence problems – the positivity property and applications, Commun. Nonlinear Sci. Numer. Simul. 44 (2017) 424-448.
  * [21] Q. A. Dang and M. T. Hoang, Positive and elementary stable explicit nonstandard Runge-Kutta methods for a class of autonomous dynamical systems, Int. J. Comput. Math. 97 (2020) 2036-2054.
  * [22] Q. A. Dang and M. T. Hoang, Nonstandard finite difference schemes for a general predator-prey system, J. Comput. Sci. 36 (2019) 101015.
  * [23] Q. A. Dang and M. T. Hoang, Complete Global Stability of a Metapopulation Model and Its Dynamically Consistent Discrete Models, Qual. Theory Dyn. Syst. 18 (2019) 461-475.
  * [24] Q. A. Dang, M. T. Hoang, Lyapunov direct method for investigating stability of nonstandard finite difference schemes for metapopulation models, J. Differ. Equ. Appl. 24 (2018) 15-47.
  * [25] Q. A. Dang and M. T. Hoang, Positivity and global stability preserving NSFD schemes for a mixing propagation model of computer viruses, J. Comput. Appl. Math. 374 (2020) 112753.
  * [26] Q. A. Dang and M. T. Hoang, Dynamically consistent discrete metapopulation model, J. Differ. Equ. Appl. 22 (2016) 1325-1349.
  * [27] D. L. DeAngelis, R. A. Goldstein, and R. V. O’Neill, A Model for Tropic Interaction, Ecology 56 (1975) 881-892.
  * [28] D. T. Dimitrov and H. V. Kojouharov, Nonstandard finite-difference schemes for general two-dimensional autonomous dynamical systems, Appl. Math. Lett. 18 (2005) 769-774.
  * [29] D. T. Dimitrov and H. V. Kojouharov, Stability-Preserving Finite-Difference Methods for General Multi-Dimensional Autonomous Dynamical Systems, Int. J. Numer. Anal. Model. 4 (2007) 280-290.
  * [30] D. T. Dimitrov and H. V. Kojouharov, Dynamically consistent numerical methods for general productive-destructive systems, J. Differ. Equ. Appl. 17 (2011) 1721-1736.
  * [31] D. T. Dimitrov and H. V. Kojouharov, Complete mathematical analysis of predator-prey models with linear prey growth and Beddington-DeAngelis functional response, Comput. Appl. Math. 162 (2005) 523-538.
  * [32] L. Edelstein-Keshet, Mathematical Models in Biology, Society for Industrial and Applied Mathematics, Philadelphia, 2005.
  * [33] M. Ehrhardt and R. E.Mickens, A nonstandard finite difference scheme for convection-diffusion equations having constant coefficients, Appl. Math. Comput. 219 (2013) 6591-6604.
  * [34] S. Elaydi, An Introduction to Difference Equations, Springer New York, 2005.
  * [35] H. Fatoorehchi and M. Ehrhardt, Numerical and semi-numerical solutions of a modified Thévenin model for calculating terminal voltage of battery cells, J Energy Storage 145(2022) 103746.
  * [36] S. M. Garba, A. B. Gumel, A. S. Hassan, and J. M. -S. Lubuma, Switching from exact scheme to nonstandard finite difference scheme for linear delay differential equation, Appl. Math. Comput. 258 (2015) 388-403.
  * [37] G. González-Parra, A. J. Arenas, and B. M. Chen-Charpentier, Combination of nonstandard schemes and Richardson’s extrapolation to improve the numerical solution of population models, Math Comput Model 52 (2010) 1030-1036.
  * [38] A. B. Gumel and S. M. Moghadas, A qualitative study of a vaccination model with non-linear incidence, Appl. Math. Comput. 143 (2003) 409-419.
  * [39] W. H. Hethcote, M. Zhien, and L. Shengbing, Effects of quarantine in six endemic models for infectious diseases, Math. Biosci. 180 (2002) 141-160.
  * [40] M. T. Hoang, A novel second-order nonstandard finite difference method for solving one-dimensional autonomous dynamical systems, Commun. Nonlinear Sci. Numer. Simul. 114 (2022) 106654.
  * [41] M. T. Hoang, Positivity and boundedness preserving nonstandard finite difference schemes for solving Volterra’s population growth model, Math. Comput. Simul. 199 (2022) 359-373.
  * [42] M. T. Hoang, Dynamically consistent nonstandard finite difference schemes for a virus-patch dynamic model, J Appl Math Comput 68 (2022) 3397-3423.
  * [43] M. T. Hoang, Reliable approximations for a hepatitis B virus model by nonstandard numerical schemes, Math. Comput. Simul. 193(2022) 32-56.
  * [44] M. T. Hoang, Dynamical analysis of a generalized hepatitis B epidemic model and its dynamically consistent discrete model, Math. Comput. Simul. 205(2023) 291-314
  * [45] Z. Horváth, On the positivity step size threshold of Runge-Kutta methods, Appl. Numer. Math. 53 (2005) 341-356.
  * [46] W. O. Kermack and A. G. McKendrick, Contributions to the mathematical theory of epidemics - I, Proc. R. Soc. Lond. 115 (1927) 700-721.
  * [47] W. O. Kermack and A. G. McKendrick, Contributions to the mathematical theory of epidemics - II, Proc. R. Soc. Lond. 138 (1932) 55-83.
  * [48] J. E. Keymer, P.A. Marquet, J. X. Velasco-Hernandez, and S. A. Levin, Extinction thresholds and metapopulation persistence in dynamic landscapes, Am. Nat. 156 (2000) 478-494.
  * [49] H. K. Khalil,  Nonlinear systems, Third Edition, Prentice Hall, 2002.
  * [50] H. V. Kojouharov, S. Roy, M. Gupta, F. Alalhareth, and J. M.Slezak, A second-order modified nonstandard theta method for one-dimensional autonomous differential equations, Appl. Math. Lett. 112 (2021) 106775.
  * [51] C. M.Kribs-Zaleta and J. X. Velasco-Hernández, A simple vaccination model with multiple endemic states, Math. Biosci. 164 (2000) 183-201.
  * [52] J. P. Lasalle, The Stability of Dynamical Systems, SIAM, Philadelphia, 1976.
  * [53] A. M. Lyapunov, General Problem of the Stability Of Motion, 1st Edition, CRC Press, 1992.
  * [54] M. Martcheva, An Introduction to Mathematical Epidemiology, Springer New York, NY, 2015.
  * [55] J. Martín-Vaquero, A. Martín del Rey, A. H. Encinas, J. D. Hernández Guillén, A. Queiruga-Dios, and G. Rodríguez Sánchez, Higher-order nonstandard finite difference schemes for a MSEIR model for a malware propagation, J. Comput. Appl. Math. 317 (2017) 146-156.
  * [56] J. Martín-Vaquero, A. Queiruga-Dios, A.Martín del Rey, A. H. Encinas, J. D. Hernández Guillén, and G. Rodríguez Sánchez, Variable step length algorithms with high-order extrapolated non-standard finite difference schemes for a SEIR model, J. Comput. Appl. Math. 330 (2018) 848-854.
  * [57] R. Mattheij and J. Molenaar, Ordinary Differential Equations in Theory and Practice, Society for Industrial and Applied Mathematics, Philadelphia, 2002.
  * [58] J. Mena-Lorca and H. W. Hethcote, Dynamic models of infectious diseases as regulators of population sizes, J. Math. Biol. 30 (1992) 693-716.
  * [59] R. E. Mickens, Nonstandard Finite Difference Models of Differential Equations, World Scientific, 1993.
  * [60] R. E. Mickens, Applications of Nonstandard Finite Difference Schemes, World Scientific, 2000.
  * [61] R. E. Mickens, Advances in the Applications of Nonstandard Finite Difference Schemes, World Scientific, 2005.
  * [62] R. E. Mickens, Nonstandard Finite Difference Schemes for Differential Equations, J. Differ. Equ. Appl. 8 (2002) 823-847.
  * [63] R. E. Mickens, Nonstandard Finite Difference Schemes: Methodology and Applications, World Scientific, 2020.
  * [64] R. E. Mickens, Dynamic consistency: a fundamental principle for constructing nonstandard finite difference schemes for differential equations, J. Differ. Equ. Appl. 11 (2005) 645-653.
  * [65] R. E. Mickens and T. M. Washington, NSFD discretizations of interacting population models satisfying conservation laws, Comput. Math. with Appl. 66 (2013) 2307-2316.
  * [66] R. E. Mickens, A nonstandard finite difference scheme for a PDE modeling combustion with nonlinear advection and diffusion, Math. Comput. Simul. 69 (2005) 439-446.
  * [67] R. E. Mickens, A nonstandard finite difference scheme for the diffusionless Burgers equation with logistic reaction, Math. Comput. Simul. 62 (2003) 117-124.
  * [68] R. E. Mickens, A nonstandard finite difference scheme for a Fisher PDE having nonlinear diffusion, Comput. Math. with Appl. 45 (2003) 429-436.
  * [69] R. E. Mickens, A nonstandard finite difference scheme for a nonlinear PDE having diffusive shock wave solutions, Math. Comput. Simul. 55 (2001) 549-555.
  * [70] R. E. Mickens, A nonstandard finite-difference scheme for the Lotka-Volterra system, Appl. Numer. Math. 45 (2003) 309-314
  * [71] R. E. Mickens and I. H. Herron, Approximate rational solutions to the Thomas-Fermi equation based on dynamic consistency, Appl. Math. Lett. 116 (2021) 106994.
  * [72] R. E. Mickens, K. Oyedeji, and S. Rucker, Exact finite difference scheme for second-order, linear ODEs having constant coefficients, J. Sound Vib. 287 (2005) 1052-1056.
  * [73] R. E. Mickens and I. Ramadhani, Finite-difference schemes having the correct linear stability properties for all finite step-sizes III, Comput. Math. with Appl. 27 (1994) 77-84.
  * [74] R.E. Mickens and A. B. Gumel, Construction and analysis of a non-standard finite difference scheme for the Burgers-Fisher equation, J. Sound Vib.  257 (2002) 791-797.
  * [75] R. E. Mickens, Numerical study of a non-standard finite difference scheme for the Val der Pol equation, J. Sound Vib. 250 (2002) 955-963.
  * [76] R. E. Mickens, NSFD scheme for acoustic propagation with the linearized Euler equations, Math. Comput. Simul. 127 (2016) 189-193.
  * [77] S. M. Moghadas and A. B. Gumel, Global stability of a two-stage epidemic model with generalized non-linear incidence, Math. Comput. Simul. 60 (2002) 107-118.
  * [78] K. C. Patidar, On the use of nonstandard finite difference methods, J. Differ. Equ. Appl. 11 (2005) 735-758.
  * [79] K. C. Patidar, Nonstandard finite difference methods: recent trends and further developments, J. Differ. Equ. Appl. 22 (2016) 817-849.
  * [80] Lih-Ing W. Roeger, Dynamically consistent discrete Lotka-Volterra competition models derived from nonstandard finite-difference schemes, Discrete Contin. Dyn. Syst. - B. 9 (2008) 415-429.
  * [81] Lih-Ing W. Roeger and G. L. Jr, Dynamically consistent discrete Lotka-Volterra competition systems, J. Differ. Equ. Appl. 19 (2013) 191-200
  * [82] Lih-Ing W. Roeger, Dynamically consistent discrete-time SI and SIS epidemic models, Conference Publications, 2013(special): 653-662. doi: 10.3934/proc.2013.2013.653.
  * [83] Lih-Ing W. Roeger, Exact finite-difference schemes for two-dimensional linear systems with constant coefficients, J. Comput. Appl. Math. 219 (2008) 102-109.
  * [84] H. L. Smith and P. Waltman, The Theory of the Chemostat: Dynamics of Microbial Competition, Cambridge University Press, 2009.
  * [85] A. Stuart and A. R. Humphries, Dynamical Systems and Numerical Analysis, Cambridge University Press, 1998.
  * [86] A. K. Verma, M. K. Rawani, and C. Cattani, A numerical scheme for a class of generalized Burgers’ equation based on Haar wavelet nonstandard finite difference method, Appl. Numer. Math. 168 (2021) 41-54.
  * [87] D. T. Wood, D. T. Dimitrov, and H. V. Kojouharov, A nonstandard finite difference method for $n$-dimensional productive-destructive systems, J. Differ. Equ. Appl. 21 (2015) 240-254.
  * [88] D. T. Wood and H. V. Kojouharov, A class of nonstandard numerical methods for autonomous dynamical systems, Appl. Math. Lett. 50 (2015) 78-82.
  * [89] D. T. Wood, H. V. Kojouharov, and D. T. Dimitrov, Universal approaches to approximate biological systems with nonstandard finite difference methods, Math. Comput. Simul. 133 (2017) 337-350.