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

Manh Tuan Hoang<EMAIL_ADDRESS><EMAIL_ADDRESS>Department of
Mathematics, FPT University, Hoa Lac Hi-Tech Park, Km29 Thang Long Blvd,
Hanoi, Viet Nam

###### Abstract

In this work, a novel second-order nonstandard finite difference (NSFD) method
that preserves simultaneously the positivity and local asymptotic stability of
one-dimensional autonomous dynamical systems is introduced and analyzed. This
method is based on novel non-local approximations for right-hand side
functions of differential equations in combination with nonstandard
denominator functions. The obtained results not only resolve the contradiction
between the dynamic consistency and high-order accuracy of NSFD methods but
also improve and extend some well-known results that have been published
recently in [Applied Mathematics Letters 112(2021) 106775], [AIP Conference
Proceedings 2302(2020) 110003] and [Applied Mathematics Letters 50(2015)
78-82]. Furthermore, as a simple but important application, we apply the
constructed NSFD method for solving the logistic, sine, cubic, and Monod
equations; consequently, the NSFD schemes constructed in the earlier work
[Journal of Computational and Applied Mathematics 110(1999) 181-185] are
improved significantly. Finally, we report some numerical experiments to
support and illustrate the theoretical assertions as well as advantages of the
constructed NSFD method.

###### keywords:

NSFD method , Second-order , Positivity , Elementary stable , Autonomous
dynamical system

###### MSC:

[2010] 37M05, 37M15, 65L05, 65P99

††journal: Elsevier

##  1 Introduction

The concept of nonstandard finite difference (NSFD) schemes was first
introduced by Mickens in 1980 to overcome drawbacks of standard finite
difference ones [2, 14, 15, 16, 17, 18, 20, 21]. Nowadays, NSFD schemes have
been widely used as one of the powerful and efficient approaches for solving
differential equation models arising in both theory and applications. The main
advantage of NSFD schemes over standard ones is that they can preserve
important qualitative properties of solutions of the corresponding
differential equations regardless of chosen step sizes. However, most of NSFD
schemes are only of first order accuracy. For this reason, there have been
many attempts to formulate higher-order NSFD schemes for ordinary differential
equations (ODEs) [4, 6, 7, 11, 12, 13].

In very recent works [6, 11], second-order NSFD methods preserving the local
asymptotic stability (LAS) of one-dimensional autonomous ODEs have been
introduced. These results partially resolved the contradiction between the
dynamic consistency and high-order accuracy of NSFD schemes; however, the
constructed methods did not preserve the positivity of solutions of
differential equations.

Motivated and inspired by this, in this work we construct a novel NSFD method
that not only preserves simultaneously the positivity and LAS of one-
dimensional autonomous dynamical systems but also is accurate of order $2$.
The method is based on novel non-local approximations for right-hand side
functions of differential equations in combination with nonstandard
denominator functions. The obtained results not only resolve the contradiction
between the dynamic consistency and high-order accuracy of NSFD schemes but
also improve and extend the well-known results constructed in [6, 11, 23].
Moreover, it should be emphasized that elementary stable NSFD schemes
introduced in some previous works [3, 8, 9] are only of first order accuracy
and not positive.

It is worthy noting that our constructed NSFD method is completely different
from NSFD schemes formulated in [23] and [14, 15, 16, 17, 18]. More clearly,
we use non-local approximations based on novel representations of right-hand
side functions and especially, denominator functions depend not only on the
step size but also on numerical approximations, meanwhile, the NSFD schemes
introduced in [14, 15, 16, 17, 18] and [23] used denominator functions only
depend on the step size. Moreover, the NSFD scheme formulated in [23] used a
non-local approximation based on the sign of the right-hand side functions at
each iterative step; however, our method does not require the determination of
the sign of the right-hand functions. In fact, the positive and elementary
stable NSFD method in [23] can be considered as a particular case of our
method, but it is only of first-order accuracy. This will be confirmed by a
numerical example performed in Subsection 4.1.

Importantly, the present idea used to construct the NSFD method can be
developed for formulating dynamically consistent NSFD methods having high-
order accuracy for general multi-dimensional autonomous and non-autonomous
dynamical systems, which are governed by ODEs, fractional-order ODEs and
partial differential equations. Furthermore, the constructed method is fully
explicit; hence, it is easy to be implemented and can be combined with
extrapolation techniques to improve the accuracy.

As a simple but important application, we apply the constructed NSFD method
for solving the logistic, sine, cubic, and Monod equations, which were
considered in [19], and general multi-dimensional dynamical systems. The
result is that we obtain positive and elementary stable NSFD schemes of
seconder-order accuracy for these equations. Meanwhile, the NSFD schemes
constructed in [19] are only accurate of order $1$. Therefore, the results
constructed in [19] are improved significantly.

This work is organized as follows. In Section 2, we provide some definitions
and preliminaries. The NSFD method is constructed and analyzed in Section 3.
Section 4 reports a set of numerical simulations and some applications of the
NSFD method. The last section presents some remarks and conclusions.

##  2 Definitions and preliminaries

Consider general autonomous ODEs of the form

$\dfrac{dy}{dt}=f(y),\quad t\geq 0,\quad y(0)=y_{0},$ (1)

where $y:\mathbb{R}\to\mathbb{R}$, $f\in C^{2}(\mathbb{R})$ and
$y_{0}\in\mathbb{R}$. Besides, we consider a general finite difference method
for Eq. (1)

$D_{h}(y_{n})=F_{h}(f;y_{n}),$ (2)

where $D_{h}(y_{n})\approx\dfrac{dy}{dt}|_{t=t_{n}}$, $F_{h}(f;y_{n})\approx
f(y(t_{n}))$, $y_{n}\approx y(t_{n})$, $t_{n}=nh$ and $h>0$ is the step size.

######  Definition 1 (Positivity-preserving numerical schemes [14, 15, 16, 17,
18, 23]).

The finite difference method (2) is called positive, if, for any value of the
step size $h$, and $y_{0}\in\mathbb{R}^{+}:=\\{y\in\mathbb{R}|\,y\geq 0\\}$
its solution remains positive, i.e., $y_{n}\in\mathbb{R}^{+}$ for all
$n\in\mathbb{N}$.

######  Definition 2 (Elementary stable numerical schemes [2, 8]).

The finite difference method (2) is said to be elementary stable if, for any
value of the step-size $h$, its only fixed points $y^{*}$ are the same as the
equilibria of the ODE (1) and the local stability properties of each $y^{*}$
are the same for both the differential equation and the difference method.

######  Definition 3 (Concept of NSFD schemes [2, 14, 15, 16, 17, 18]).

The one-step finite difference scheme (2) for solving Eq. (1) is called a
nonstandard finite difference (NSFD) method if at least one of the following
conditions is satisfied:

  * 1.

$D_{h}(y_{n})=\dfrac{y_{n+1}-y_{n}}{\varphi(h)}$, where
$\varphi(h)=h+\mathcal{O}(h^{2})$ is a non-negative function;

  * 2.

$F_{h}(f;y_{n})=g(y_{k},y_{k+1},h)$, where $g(y_{k},y_{k+1},h)$ is a non-local
approximation of the right-hand side of Eq. (1).

We recall that Eq. (1) is said to be positive if $y(t)\geq 0$ holds for all
$t\geq 0$ whenever $y_{0}\geq 0$. Thanks to [10, Lemma 1], we obtain the
following result on the positivity of solutions of Eq. (1).

######  Lemma 1.

The set $\mathbb{R}^{+}$ is a positively invariant set of Eq. (1) (i.e,
$y(t)\geq 0$ for all $t\geq 0$ whenever $y_{0}\geq 0$) if and only if
$f(0)\geq 0$.

##  3 Construction of the nonstandard finite difference method

###  3.1 Auxiliary results

We first introduce some notations

$\begin{split}&\mathcal{F}^{+}=\big{\\{}f:\mathbb{R}_{+}\to\mathbb{R}\big{|}\,f(y)\geq
0\quad\text{for all}\quad y\geq 0\big{\\}},\\\
&\mathcal{F}^{-}=\big{\\{}f:\mathbb{R}_{+}\to\mathbb{R}\big{|}\,f(y)\leq
0\quad\text{for all}\quad y\geq 0\big{\\}},\\\
&\mathcal{B}^{L}=\big{\\{}f:\mathbb{R}_{+}\to\mathbb{R}\big{|}\,f(y)\geq
L\quad\text{for all}\quad y\geq 0\big{\\}}\quad\text{for each}\quad M\geq
0.\end{split}$

It is easy to see that $\mathcal{B}^{L}$ is not empty for any $L\geq 0$, for
example, $f(y)=L$, $f(y)=y^{2}+L$ and $f(y)=Le^{y}$ belong to
$\mathcal{B}^{L}$. Moreover, the set $\mathcal{B}^{L}\cap
C^{k}(\mathbb{R}^{+})$ ($k\geq 1$) is also not empty.

######  Lemma 2.

Let $f:\mathbb{R}^{+}\to\mathbb{R}$ be a function satisfying:

  1. (i)

$f(y)\in C^{k}(\mathbb{R}^{+})$, $k\geq 1$.

  2. (ii)

$f$ has a finite number of zeros, i.e, the equation $f(y)=0$ has only a finite
number of solutions.

Then $f$ can be represented in the form

$f(y)=f^{+}(y)+f^{-}(y),\quad f^{+},f^{-}\in C^{k}(\mathbb{R}^{+}),\quad
f^{+}\in\mathcal{F}^{+},\quad f^{-}\in\mathcal{F}^{-}.$

###### Proof.

Assume that $0\leq y_{1}\leq y_{2}<\ldots\leq y_{m}$ with $m<\infty$ are zeros
of $f$. We set

$l=\min_{y\in[0,y_{m}]}f(y),\quad L=\max_{y\in[0,y_{m}]}f(y),\quad
M=max\\{|l|,\,\,|L|\\}.$

Since $f$ has only a finite number of zeros, $f(y)$ does not change sign when
$y>y_{m}$. We consider two following cases for the sign of $f(y)$ when
$y>y_{m}$.
Case 1: The case when $f(y)>0$ for all $y>y_{m}$. In this case, we take

$f^{+}(y)=f(y)+g(y),\qquad f^{-}(y)=-g(y),\quad\text{for some}\quad
g\in\mathcal{B}^{M}\cap C^{k}(\mathbb{R}^{+}).$

Obviously, $f^{-}(y)\leq 0$ for all $y\geq 0$. On the other hand, since $g\in
B^{M}$, $g(y)\geq M$ for all $y\geq 0$. Consequently, for $y\in[0,y_{m}]$

$f^{+}(y)=f(y)+g(y)\geq l+M\geq l+|l|\geq 0.$

Hence, $f^{+}(y)\geq 0$ for all $y\geq 0$.
Case 2: The case when $f(y)<0$ for all $y>y_{m}$. In this case, we choose

$f^{-}(y)=f(y)-g(y),\qquad f^{+}(y)=g(y),\quad\text{for some}\quad
g\in\mathcal{B}^{M}\cap C^{k}(\mathbb{R}^{+}).$

Obviously, $f^{+}(y)\geq 0$ for all $y\geq 0$. Furthermore

$f^{-}(y)=f(y)-f(y)\leq L-g(y)\leq L-M\leq L-|L|=0\quad\text{for all}\quad
y\in[0,y_{m}].$

Hence, $f^{-}(y)\geq 0$ for all $y\geq 0$.
The proof is completed by combining Cases 1 and 2. ∎

######  Remark 1.

Based on the proof of Lemma 2, we can write a procedure to determine
representations of the function $f$ as follows:
Step 1. Determine all zeros the function $f$, i.e., solve the equation
$f(y)=0$ (using numerical methods such as Newton’s method if needed). Since
$f$ has only only a finite number of zeros, this step is always done. Suppose
that $0\leq y_{1}\leq y_{2}\leq\ldots\leq y_{m}$ with $m<\infty$ are all the
zeros of $f$.
Step 2. Compute (by using numerical methods if needed)

$l=\min_{y\in[0,y_{m}]}f(y),\quad L=\max_{y\in[0,y_{m}]}f(y),\quad
M=max\\{|l|,\,\,|L|\\}.$

Step 3. Determine the sign (positive or negative) of the function $f$ on
$(y_{m},\infty)$. This step is always done because $f$ has only only a finite
number of zeros.
Step 4. If $f(y)>0$ for all $y>y_{m}$, we take

$f^{+}(y)=f(y)+g(y),\qquad f^{-}(y)=-g(y),\quad\text{for some}\quad
g\in\mathcal{B}^{M}\cap C^{k}(\mathbb{R}^{+}).$

If $f(y)<0$ for all $y>y_{m}$, we choose

$f^{-}(y)=f(y)-g(y),\qquad f^{+}(y)=g(y),\quad\text{for some}\quad
g\in\mathcal{B}^{M}\cap C^{k}(\mathbb{R}^{+}).$

Note that the set $\mathcal{B}^{M}\cap C^{k}(\mathbb{R}^{+})$ is always
nonempty for each $M>0$.

######  Theorem 1.

Let $f:\mathbb{R}_{+}\to\mathbb{R}$ be a function satisfying:

  1. (i)

$f(y)\in C^{2}(\mathbb{R}^{+})$.

  2. (ii)

$f$ has a finite number of zeros.

  3. (iii)

$f(0)\geq 0$.

Then $f$ can be represented in the form

$f(y)=f^{+}(y)+yf^{-}(y),\quad f^{+},f^{-}\in C^{1}(\mathbb{R}^{+}),\quad
f^{+}\in\mathcal{F}^{+},\quad f^{-}\in\mathcal{F}^{-}.$ (3)

###### Proof.

The proof includes two steps as follows.
Step 1. Consider the case when $f(0)=0$. We define

$g(y):=\begin{cases}&\dfrac{f(y)}{y},\qquad y>0,\\\ &f^{\prime}(0),\qquad
y=0.\end{cases}$

It is clear that $f(y)=yg(y)$ and $g\in C^{1}(\mathbb{R}^{+})$. From Lemma 2,
the function $g$ admits a representation

$g(y)=g^{+}(y)+g^{-}(y),\quad g^{+},g^{-}\in C^{1}(\mathbb{R}^{+}),\quad
g^{+}\in\mathcal{F}^{+},\quad g^{-}\in\mathcal{F}^{-}.$

Consequently,

$f(y)=f^{+}(y)+yf^{-}(y),\quad f^{+}(y):=yg^{+}(y),\quad f^{-}(y)=g^{-}(y).$

Step 2: If $f(0)>0$, we set $g(y)=f(y)-f(0)$. Then, $g(0)=0$. From Step 1 we
have that $g$ has a representation

$g(y)=g^{+}(y)+yg^{-}(y),\quad g^{+},g^{-}\in C^{1}(\mathbb{R}^{+}),\quad
g^{+}\in\mathcal{F}^{+},\quad g^{-}\in\mathcal{F}^{-}.$

Hence,

$f(y)=g(y)+f(0)=\big{[}f(0)+g^{+}(y)\big{]}+yg^{-}(y),$

where $f(0)+g^{+}(y)\in\mathcal{F}^{+}$ and $g^{-}(y)\in\mathcal{F}^{-}$. The
proof is completed by combining Step 1 and Step 2. ∎

######  Remark 2.

Similarly to Remark 1, the proof of Theorem 1 also provides a way to determine
representations of the function $f$. In many cases, it is easy to find
representations in the form (3) for a given function $f$, for instance

$\begin{split}&y-y^{2}=y+y(-y)=(y+y^{2})+2y(-y)=(y+ye^{y})+y(-y-e^{y}),\\\
&e^{y}-1=(e^{y}+y-1)+y(-1)=(e^{y}+y^{2}-1)+y(-y),\\\
&\sin(y)=\big{[}\sin(y)+y\big{]}+y(-1)=\big{[}\sin(y)+ye^{y}\big{]}+y(-e^{y}).\end{split}$

Hence, the representations in the form (3) of $f$ may not be unique. This
leads to an optimal problem for the representation of the function $f$.

######  Remark 3.

In [5], a representation theorem for real-valued functions was constructed
(see [5, Theorem 10]); however, this theorem did not affirm the
differentiability of negative and positive parts $f^{+}$ and $f^{-}$. It only
proved that $f^{+}$ and $f^{-}$ are continuous. Actually, as will be seen in
Lemma 4, the differentiability of $f^{+}$ and $f^{-}$ is very essential in the
construction of the elementary stable NSFD method.

###  3.2 Main results

Throughout this paper, we always assume that:

  1. (A1)

The right-hand side function $f$ of Eq. (1) satisfies $f(0)\geq 0$. This
condition implies that $y(t)\geq 0$ for all $t\geq 0$ whenever $y_{0}\geq 0$.

  2. (A2)

The equation (1) has a finite number of equilibria and each of which is
hyperbolic, i.e., $f^{\prime}(y_{*})\neq 0$ for any equilibrium point $y_{*}$
of (1). This condition means that the LAS of equilibria of (1) can be
investigated by the Lyapunov indirect method (see [1, 22]).

Our main objective is to construction an NSFD method preserving simultaneously
the positivity and LAS of Eq. (1). First, we utilize Theorem 1 to transform
the right-hand side function $f$ of Eq. (1) to the form:

$f(y)=f^{+}(y)+yf^{-}(y),\quad f^{+},f^{-}\in C^{1}(\mathbb{R}^{+}),\quad
f^{+}\in\mathcal{F}^{+},\quad f^{-}\in\mathcal{F}^{-}.$ (4)

Then we propose the following NSFD scheme

$\dfrac{y_{n+1}-y_{n}}{\varphi(h,y_{n})}=f^{+}(y_{n})+\alpha
y_{n}f^{-}(y_{n})+\beta y_{n+1}f^{-}(y_{n}),$ (5)

where $\varphi:\mathbb{R}^{+}\times\mathbb{R}^{+}\to\mathbb{R}^{+}$ satisfies

$\ \varphi(h,y)=h+\mathcal{O}(h^{2})\quad\mbox{as}\quad h\to 0,$ (6)

and $\alpha,\beta$ are real numbers satisfying

$\alpha+\beta=1,\quad\alpha\leq 0,\quad\beta\geq 0.$ (7)

For convenience, the variables in the function $\varphi(h,y)$ will be omitted
in some places; Furthermore, throughout the rest of this section, we always
assume that the conditions (6) and (7) are satisfied.

######  Remark 4.

  1. (i)

The denominator function $\varphi$ of the NSFD method (5) depends not only on
$h$ but also on $y_{n}$ to guarantee that it is accurate of order $2$ (see
Lemma 5). Furthermore, it is determined by the correct asymptotic behaviour
modelling.

  2. (ii)

Differently from most of NSFD schemes introduced in [14, 15, 16, 17, 18], the
NSFD scheme (5) uses a non-local approximation with the appearance of the
parameters $\alpha$ and $\beta$ as weights. These parameters not only ensure
the positivity (see Lemma 3) but also make the method more flexible. As will
be seen in Subsection 4.1, the NSFD method (5) can become an exact scheme with
suitable parameters $\alpha$ and $\beta$.

######  Lemma 3 (Positivity of the NSFD method).

Let $y_{0}\geq 0$ be any initial data for Eq. (1) and $\\{y_{n}\\}_{n\geq 1}$
be approximations generated by the NSFD method (5). Then we have $y_{n}\geq 0$
for all $n\geq 1$. In other words, the NSFD method (5) preserves the
positivity of Eq. (1) for all finite step sizes $h>0$.

###### Proof.

The theorem is proved by mathematical induction. We first transform the scheme
(5) to the explicit form

$y_{n+1}=\dfrac{y_{n}+\varphi f^{+}(y_{n})+\varphi\alpha
y_{n}f^{-}(y_{n})}{1-\varphi\beta f^{-}(y_{n})}$ (8)

Assume that $y_{n}\geq 0$. From

$f^{+}(y_{n})\geq 0,\quad f^{-}(y_{n})\leq 0,\quad\alpha\leq 0,\quad\beta\geq
0,\quad\varphi\geq 0,$

we deduce that $y_{n+1}\geq 0$. The proof is complete. ∎

To determine the Jacobian of the scheme (5) at its equilibria conveniently, we
rewrite the scheme (8) in the form

$y_{n+1}=y_{n}+\dfrac{\varphi f(y_{n})}{1-\varphi\beta f^{-}(y_{n})}.$ (9)

Assume that the set of equilibria of Eq. (1) is
$\Omega=\\{y_{1}^{*},y_{2}^{*},\ldots,y_{m}^{*}\\}$. We denote by $\Omega^{s}$
is the set of locally asymptotically stable equilibria of Eq. (1), that is

$\Omega^{s}=\\{y^{*}\in\Omega|f^{\prime}(y^{*})<0\\}.$

######  Lemma 4 (Stability of the NSFD method).

Assume that the function $\varphi(h,y)$ satisfies

$\varphi(h,y^{*})<\dfrac{2}{2\beta
f^{-}(y^{*})-f^{\prime}(y^{*})},\qquad\forall y^{*}\in\Omega^{s}.$ (10)

Then the NSFD method (5) is elementary stable.

###### Proof.

It is easy to verify that equilibria of (1) and (9) are identical. Let $y^{*}$
be any equilibrium point of (1). The Jacobian of (9) at $y^{*}$ is given by

$J(y^{*})=1+\dfrac{\varphi(h,y^{*})f^{\prime}(y^{*})}{1-\varphi(h,y^{*})\beta
f^{-}(y^{*})}$

We consider two following cases:
(i) If $f^{\prime}(y^{*})>0$, then $y^{*}$ is an unstable equilibrium point of
(1). We need to prove that $y^{*}$ is also an unstable equilibrium point of
(9). Indeed, since

$J(y^{*})=1+\dfrac{\varphi(h,y^{*})f^{\prime}(y^{*})}{1-\varphi(h,y^{*})\beta
f^{-}(y^{*})}>1,$

$y^{*}$ is an unstable equilibrium point of (9).
(ii) If $f^{\prime}(y^{*})<0$, then $y^{*}$ is a locally asymptotically stable
equilirium point of (1). We will prove that $y^{*}$ is a locally
asymptotically stable equilirium point of (9). Indeed, we have

$J(y^{*})=1+\dfrac{\varphi(h,y^{*})f^{\prime}(y^{*})}{1-\varphi(h,y^{*})\beta
f^{-}(y^{*})}<1.$

On the other hand, $J(y^{*})>-1$ if and only if

$2>\varphi(h,y^{*})\big{[}2\beta f^{-}(y^{*})-f^{\prime}(y^{*})\big{]}.$

This inequality is satisfied if (10) holds. Hence, if (10) then
$|J(y^{*})|<1$, or equivalently, $y^{*}$ is a stable equilibrium point of (9).
The proof is completed. ∎

######  Lemma 5 (The second-order NSFD method).

Let $\varphi(h,y)$ be a function satisfying

$\dfrac{\partial^{2}\varphi}{\partial h^{2}}(0,y)=f^{\prime}(y)-2\beta
f^{-}(y),\quad\forall y\geq 0.$ (11)

Then the NSFD method (5) is accurate of order $2$.

###### Proof.

Using the Taylor expansion for $y(t)$ we obtain

$y(t_{n+1})=y(t_{n}+h)=y(t_{n})+hy^{\prime}(t_{n})+\dfrac{h^{2}}{2}y^{\prime\prime}(t_{n})+\mathcal{O}(h^{3}).$
(12)

Let us denote

$G(h,y)=y+\dfrac{\varphi(h,y)f(y)}{1-\varphi(h,y)\beta f^{-}(y)}.$

It is easy to verify that

$G(0,y)=y,\quad\dfrac{\partial G}{\partial
h}(0,y)=f(y),\quad\dfrac{\partial^{2}G}{\partial
h^{2}}(0,y)=\dfrac{\partial^{2}\varphi}{\partial h^{2}}(0,y)+2\beta f^{-}(y).$

Consequently,

$\begin{split}y(t_{n})+\dfrac{\varphi f(y(t_{n}))}{1-\varphi\beta
f^{-}(y(t_{n}))}&=G(0,y(t_{n}))+h\dfrac{\partial G}{\partial
h}(0,y(t_{n}))+\dfrac{h^{2}}{2}\dfrac{\partial^{2}G}{\partial
h^{2}}(0,y(t_{n}))+\mathcal{O}(h^{3})\\\
&=y(t_{n})+hf(y(t_{n}))+\dfrac{h^{2}}{2}\bigg{[}\dfrac{\partial^{2}\varphi}{\partial
h^{2}}(0,y(t_{n}))+2\beta
f^{-}(y(t_{n}))\bigg{]}+\mathcal{O}(h^{3}).\end{split}$ (13)

From (12) and (13) we deduce that if (11) holds then

$y(t_{n+1})-\bigg{[}y(t_{n})+\dfrac{\varphi f(y(t_{n}))}{1-\varphi\beta
f^{-}(y(t_{n}))}\bigg{]}=\mathcal{O}(h^{3}),$

which means that the NSFD scheme is accurate of order 2. ∎

From Lemmas 3, 4 and 5, we obtain

######  Theorem 2.

The NSFD method (5) is positive, elementary stable and accurate of order $2$
if

$\begin{split}(H1)\qquad&\varphi(h,y)=h+\mathcal{O}(h^{2})\quad\mbox{as}\quad
h\to 0,\quad\varphi(h,y)\geq 0\quad\text{for
all}\quad(h,y)\in\mathbb{R}^{+}\times\mathbb{R}^{+},\\\
(H2)\qquad&\varphi(h,y^{*})<\dfrac{2}{2\beta
f^{-}(y^{*})-f^{\prime}(y^{*})},\qquad\forall y^{*}\in\Omega^{s},\\\
(H3)\qquad&\dfrac{\partial^{2}\varphi}{\partial
h^{2}}(0,y)=f^{\prime}(y)-2\beta f^{-}(y),\quad\forall y\geq 0.\\\
(H4)\qquad&\alpha+\beta=1,\quad\alpha\leq 0,\quad\beta\geq 0.\end{split}$ (14)

Next, we introduce a function $\varphi(h)$ satisfying the conditions
$(H1)-(H3)$. Let us denote

$\lambda(y):=-f^{\prime}(y)+2\beta f^{-}(y),$
$\varphi(h,y)=\begin{cases}\dfrac{1-e^{-h\lambda(y)}}{\lambda(y)},\quad&\lambda(y)\neq
0,\\\ h,\qquad&\lambda(y)=0.\end{cases}$ (15)

######  Proposition 1.

The function $\varphi(h,y)$ given by (15) satisfies the conditions $(H1)-(H3)$
in Theorem 2.

###### Proof.

It is easy to verify that $\varphi(h,y)$ satisfies $(H1)$ and $(H3)$. Since
$\lambda(y^{*})>0$ for all $y^{*}\in\Omega^{s}$, we have

$\varphi(h,y^{*})<\dfrac{1}{\lambda(y^{*})}=\dfrac{1}{-f^{\prime}(y^{*})+2\beta
f^{-}(y^{*})}<\dfrac{2}{-f^{\prime}(y^{*})+2\beta f^{-}(y^{*})},$

which implies that $(H2)$ is also satisfied. The proof is completed. ∎

######  Remark 5.

  1. (i)

Lemma 2 in [6] provides an important result for selecting the denominator
functions satisfying the conditions $(H1)-(H3)$ in Theorem 2.

  2. (ii)

The convergence of the NSFD scheme (5) can be established similarly to the
results constructed in [5, Section 6].

##  4 Some applications of the NSFD method

In this section, we present some numerical examples and applications of the
NSFD method (5). In all the following examples, the denominator functions
$\varphi$ for the NSFD method (5) are determined by using Proposition 1.

###  4.1 The logistic differential equation

We consider the well-known logistic equation (see [1])

$\dfrac{dy}{dt}=2\bigg{(}y-\dfrac{y^{2}}{2}\bigg{)}:=f(y),\quad 0\leq t\leq
1,\quad y(0)=0.5.$ (16)

The exact solution of this equation is given by

$y(t)=\dfrac{1}{0.5+1.5e^{-2t}}.$

The equation has two equilibria $y_{0}^{*}=0$ (unstable) and $y_{1}^{*}=2$
(stable).

From the NSFD method (5), we obtain the following second-order NSFD schemes
for the Eq. (16), which correspond to three different representations of the
right-hand side function $f(y)=2y-y^{2}$ given in Table 1:

Table 1: Second-order NSFD schemes for E.q (16). Scheme | $f^{+}$ | $f^{-}$ | $\beta$ | Denominator function | Remark
---|---|---|---|---|---
sNSFD1 | $2y$ | $-y$ | $1.25$ | $\varphi(h,y)=\dfrac{1-e^{(2+0.5y)h}}{2+0.5y}$ | Second-order scheme
sNSFD2 | $2y+y^{2}$ | $-2y$ | $1.25$ | $\quad\varphi(h,y)=\dfrac{1-e^{-h(2+3y)}}{2+3y}$ | Second-order scheme
sNSFD3 | $2y$ | $-y$ | $1$ | $\varphi(h,y)=\dfrac{1-e^{-2h}}{2}$ | Exact NSFD scheme

We now observe absolute errors and convergence rates of the schemes sNSFD1 and
sNSFD2 at the end of the integration interval, namely, $T=1$. Numerical
results are presented in Table 2. In this table,

$error=|y(t_{N})-y_{N}|,\quad
rate=:log_{(h_{1}/h_{2})}\bigg{(}\dfrac{error(h_{1})}{error(h_{2})}\bigg{)}.$

Table 2: Absolute errors and rates of the NSFD schemes for Eq. (16). Step h | sNSFD1 error | rate | sNSFD2 error | rate | Wood & Kojourharov’s scheme (17) | rate |
---|---|---|---|---|---|---|---
$10^{-1}$ | $0.0014$ |  | $0.0127$ |  | $0.0470$ |  |
$10^{-2}$ | $1.4678e-005$ | $1.9795$ | $1.3823e-004$ | $1.9632$ | $0.0045$ | $1.0189$ |
$10^{-3}$ | $1.4749e-007$ | $1.9979$ | $1.3910e-006$ | $1.9973$ | $4.4841e-004$ | $1.0015$ |
$10^{-4}$ | $1.4756e-009$ | $1.9998$ | $1.3918e-008$ | $1.9998$ | $4.4820e-005$ | $1.0002$ |
$10^{-5}$ | $1.4742e-011$ | $2.0004$ | $1.3917e-010$ | $2.0000$ | $4.4818e-006$ | $1.0000$ |

It is clear that the schemes sNSFD1 and sNSFD2 are accurate of order 2.
However, the scheme sNSFD1 provides better errors than the scheme sNSFD2. This
means that the representations of the function $f$ affect the errors of the
NSFD schemes. This observation leads to an optimal problem for the
representation of the function $f$.

It is important to note that in the case $\beta=1$, the NSFD method (5) is
accurate of order $2$ for Eq. (16) if

$\dfrac{\partial^{2}\varphi}{\partial h^{2}}(0,y)=f^{\prime}(y)-2\beta
f^{-}(y)\equiv 2.$

In this case, the function $\varphi(h)=(1-e^{-2h})/{2}$ satisfies this
condition. For this function, we do not need to update the value of the
denominator function in each iterative step; consequently, computational
volume and computation time are reduced. This show the role of the parameter
$\beta$. Actually, in the case when $\varphi(h)=(1-e^{-2h})/{2}$ we obtain an
exact scheme for the logistic equation [14, 15, 16, 17, 18].

We now use compare the second NSFD schemes with positive and elementary stable
NSFD (PSENSFD) method constructed by Wood and Kojouharov in [23]. Note that
the Wood and Kojourharov’s scheme is given by

$y_{n+1}=\begin{cases}&y_{n}+\phi(h)\big{(}2y_{n}-y_{n}^{2}\big{)},\quad\text{if}\quad
2y_{n}-y_{n}^{2}\geq 0,\\\
&\dfrac{y_{n}^{2}}{y_{n}-\phi(h)\big{(}2y_{n}-y_{n}^{2}\big{)}},\quad\text{if}\quad
2y_{n}-y_{n}^{2}<0,\end{cases}$ (17)

where $\varphi(h)$ satisfies $\phi(h)<1$ for all $h>0$. Therefore, we take
$\phi(h)=1-e^{-h}$. The errors and convergence rates of the scheme (17) are
also presented in Table 2. It is clear that it is only accurate of order $1$.
However, it requires the computation of the sign of the right-hand side
function at each iterative step.

To end this subsection, we consider dynamics of numerical solutions generated
by the standard Euler scheme, the standard second-order Runge-Kutta (RK2)
scheme and the NSFD method (5) (scheme sNSFD1) with the step size $h=1.25$.
The numerical approximations are depicted in Figures 1 and 2. It is easy to
observe that the standard Euler scheme generated the numerical approximation
oscillates around the equilibrium position, the approximation generated by the
RK2 scheme converges to a spurious equilibrium point; however, the NSFD scheme
preserves the dynamics of the logistic equation. This result is completely
consistent with recognized results on NSFD schemes for differential equations
[2, 11, 14, 15, 16, 16, 17, 18, 20, 21, 23]. On the other hand, it is easy to
observe that the NSFD method (5) is more accurate than the PESNSFD method.

Figure 1: The approximations generated by the Euler scheme, the RK2 scheme and
the NSFD scheme with $h=1.25$. Figure 2: The approximations generated by the
new NSFD method and PENSFD scheme (scheme (17)) with $h=1.25$.

###  4.2 The cubic differential equation

Consider the cubic differential equation

$\dfrac{dy}{dt}=y(1-y^{2}),\quad y(0)\geq 0.$ (18)

In $\mathbb{R}^{+}$, the cubic equation (18) has two equilibrium point
$y^{*}=0$ (unstable) and $y^{*}=1$ (stable).

Now, Theorem 2 will be applied to obtain a positive and elementary stable NSFD
scheme of second order for Eq. (18). For this purpose, we choose

$f^{+}(y)=y,\quad
f^{-}(y)=y^{2},\quad\beta=\dfrac{3}{2},\quad\alpha=-\dfrac{1}{2}.$

Then we obtain

$\dfrac{y_{n+1}-y_{n}}{\varphi}=y_{n}+\dfrac{1}{2}(y_{n})^{3}-\dfrac{3}{2}(y_{n})^{2}y_{n+1},$

or equivalently,

$y_{n+1}=\dfrac{(2+\varphi)y_{n}+\varphi(y_{n})^{3}}{2+3\varphi(y_{n})^{2}}.$
(19)

By Theorem 2, the scheme (19) is positive, elementary stable and accurate of
order $2$ if

$\dfrac{\partial^{2}\varphi}{\partial h^{2}}(0,y)=f^{\prime}(y)-2\beta
f^{-}(y)=1.$

For this condition, the function $1-e^{-h}$ is most suitable; consequently, we
obtain

$y_{n+1}=\dfrac{(2+\varphi)y_{n}+\varphi(y_{n})^{3}}{2+3\varphi(y_{n})^{2}},\quad\varphi(h)=1-e^{-h}.$
(20)

In [19], Mickens proposed the following NSFD scheme that incorporates the
maximum symmetry in modeling the nonlinear term:

$y_{n+1}=\Bigg{[}\dfrac{(2+\phi)+\phi(y_{n})^{2}}{(2-\phi)+3\phi(y_{n})^{2}}\Bigg{]}y_{n},\quad\phi=\dfrac{1-e^{-2h}}{2}.$
(21)

Clearly, the scheme (20) is different from (21). However, it is easy to show
that the NSFD scheme (21) is a particular case of the NSFD method (5). Indeed,
for the NSFD scheme (5) we choose

$f^{+}(y)=\dfrac{1}{2}y+(y)^{3},\quad
f^{-}(y)=\dfrac{1}{2}-\dfrac{3}{2}y^{2},\quad\alpha=0,\quad\beta=1.$

Note that in this form $f^{-}\notin\mathcal{F}^{-}$. Then the NSFD method (5)
for Eq. (18) becomes

$\dfrac{y_{n+1}-y_{n}}{\varphi(h)}=\dfrac{1}{2}y_{n}+\dfrac{1}{2}(y_{n})^{3}+\dfrac{1}{2}y_{n+1}-\dfrac{3}{2}(y_{n})^{2}y_{n+1},$

or equivalently,

$y_{n+1}=\Bigg{[}\dfrac{(2+\varphi)+\varphi(y_{n})^{2}}{(2-\varphi)+3\varphi(y_{n})^{2}}\Bigg{]}y_{n}.$
(22)

It is clear that the scheme (22) will coincident with (21) if
$\varphi=({1-e^{-2h})}/{2}$. Since
$f^{-}(y)=(1/2)-(3/2)y^{2}\notin\mathcal{F}^{-}$, the scheme (22) is not
positive for any function $\varphi(h)>0$. A simple condition for the scheme to
be positive is $\varphi(h)<1/2$ for all $h>0$. Fortunately, the function
$\varphi=({1-e^{-2h})}/{2}$ satisfies this condition; therefore, the scheme
(21) is positive. Importantly, the function $\varphi=({1-e^{-2h})}/{2}$
satisfies Theorem 2; consequently, the scheme (21) is also of second order
accuracy.

The results presented in this subsection and Subsection 4.1 can be generalized
as follows.

######  Proposition 2.

The ODE

$\dfrac{dy}{dt}=ay-by^{m},\quad a,b>0,\quad m\geq 2,\quad y(0)\geq 0$ (23)

has a positive and elementary stable NSFD scheme of second order accuracy
defined by

$\dfrac{y_{n+1}-y_{n}}{\varphi(h)}=ay_{n}-b\bigg{(}1-\dfrac{m}{2}\bigg{)}(y_{n})^{m}-b\dfrac{m}{2}{(y_{n})}^{m-1}y_{n+1},\quad\varphi(h)=\dfrac{1-e^{-ah}}{a}.$
(24)

###  4.3 The modified Monod equation

In [19], Mickens considered the modified Monod equation

$\dfrac{dy}{dt}=\dfrac{(\lambda-1)y-(\lambda+1)y^{2}}{1+y},\quad\lambda>1,$
(25)

and proposed the following NSFD scheme

$\dfrac{y_{n+1}-y_{n}}{\phi}=\Bigg{[}\dfrac{(\lambda-1)-(1+\lambda)y_{n+1}}{1+y_{n}}\Bigg{]}y_{n},\quad\phi=\dfrac{1-e^{-Rh}}{R},\quad
R=\lambda-1.$ (26)

Clearly, the scheme (26) is a particular case of the NSFD method (6) with

$f^{+}(y)=\dfrac{(\lambda-1)y}{1+y},\,\,f^{-}(y)=\dfrac{(\lambda+1)y}{1+y},\,\,\beta=1,\,\,\varphi(h,y)=\dfrac{1-e^{-Rh}}{R},\,\,R=\lambda-1.$

Although the scheme (26) is positive and elementary stable, it is only of
first order accuracy. However, thanks to Theorem 2 it is easy to determine
conditions for this scheme to be accurate of order $2$. More clearly, the
scheme (26) is of second-order accuracy if

$\dfrac{\partial^{2}\varphi}{\partial
h^{2}}(0,y)=-\dfrac{(\lambda+3)+4(\lambda+1)y+3(\lambda+1)y^{2}}{(1+y)^{2}}.$

A denominator function satisfying this above condition is

$\varphi(h,y)=\dfrac{1-e^{-R(y)h}}{R(y)},\quad
R(y)=\dfrac{(\lambda+1)+4(\lambda+1)y+3(\lambda+1)y^{2}}{(1+y)^{2}}.$

This function ensures that the scheme (26) is of second order accuracy.

###  4.4 The sine equation

In [19], Mickens considered the since equation

$\dfrac{dy}{dt}=\sin(\pi y)$ (27)

and formulated the following NSFD scheme

$\dfrac{y_{n+1}-y_{n}}{\phi}=sin(\pi y_{n}),\quad\phi=\dfrac{1-e^{-\pi
h}}{\pi}.$ (28)

It is easy to verify that the scheme (28) is of first order accuracy. We now
use Theorem 2 to obtain a second-order NSFD scheme for the sine equation. For
this purpose, we rewrite Eq. (27) in the form

$\dfrac{dy}{dt}=\sin(\pi y)=f^{+}(y)+yf^{-}(y)=\big{[}sin(\pi y)+\pi
y\big{]}+y(-\pi).$

Note that $y+sin(y)\geq 0$ for all $y\geq 0$. Applying the NSFD method (5)
with $\beta=1$ we obtain

$y_{n+1}=\dfrac{y_{n}+\varphi\big{[}sin(\pi y_{n})+\pi
y_{n}\big{]}}{1+\varphi\pi}.$ (29)

The scheme is elementary stable and positive, and it is accurate of order $2$
if

$\dfrac{\partial^{2}\varphi}{\partial h^{2}}(0,y)=\pi\cos(\pi y)-2\pi.$

Hence, the function

$\varphi(h,y)=\dfrac{1-e^{-\lambda(y)h}}{\lambda(y)},\quad\lambda(y)=\pi\cos(\pi
y)-2\pi$

guarantees that the NSFD scheme (29) is positive, elementary stable and
accurate of order $2$.

###  4.5 An extension of the NSFD method for systems of ODEs

In this subsection, we generalize the NSFD method (5) for systems of ODEs. For
this purpose, we first consider the following general class of two-dimensional
differential equations including several models in population dynamics

$\begin{split}&\dfrac{dx}{dt}=x\big{(}f_{+}(x,y)-f_{-}(x,y)\big{)}:=f(x,y),\quad
x(0)=x_{0}\geq 0,\\\
&\dfrac{dy}{dt}=y\big{(}g_{+}(x,y)-g_{-}(x,y)\big{)}:=g(x,y),\quad
y(0)=y_{0}\geq 0,\end{split}$ (30)

where $f_{+},f_{-}$ and $g_{+},g_{-}$ are positive for all
$(x,y)\in\mathbb{R}_{+}\times\mathbb{R}_{+}$ and of class $\mathcal{C}^{1}$.
In [5], Cresson and Pierret constructed a nonstandard finite difference method
preserving dynamical properties including positivity and stability of the
model (30). It was proved that the constructed NSFD method is convergent and
of order one (see [5, Theorem 5.2]).

We now utilize the approach is this work to construct a second-order NSFD
method for the model (30). For this reason, we propose the following NSFD
method

$\begin{split}&\dfrac{x_{n+1}-x_{n}}{\varphi_{1}(x_{n},y_{n},h)}=x_{n}f_{+}(x_{n},y_{n})-\alpha_{1}x_{n}f_{-}(x_{n},y_{n})-\beta_{1}x_{n+1}f_{-}(x_{n},y_{n}),\\\
&\dfrac{y_{n+1}-y_{n}}{\varphi_{2}(x_{n},y_{n},h)}=y_{n}g_{+}(x_{n},y_{n})-\alpha_{2}y_{n}g_{-}(x_{n},y_{n})-\beta_{2}y_{n+1}g_{-}(x_{n},y_{n}),\end{split}$
(31)

where $\varphi_{1}$ and $\varphi_{2}$ are denominator functions determined by
the correct asymptotic behaviour modelling and

$\begin{split}&\alpha_{1}+\beta_{1}=1,\quad\alpha_{1}\leq 0,\quad\beta_{1}\geq
0,\\\ &\alpha_{2}+\beta_{2}=1,\quad\alpha_{2}\leq 0,\quad\beta_{2}\geq
0.\end{split}$ (32)

Note that when $\varphi_{1}$ and $\varphi_{2}$ do not depend on $x_{n}$ and
$y_{n}$, we will obtain the NSFD method constructed in [5].

The method (31) can be rewritten in the explicit form

$\begin{split}&x_{n+1}=\dfrac{x_{n}+\varphi_{1}x_{n}f_{+}(x_{n},y_{n})-\varphi_{1}\alpha_{1}x_{n}f_{-}(x_{n},y_{n})}{1+\varphi_{1}\beta_{1}f_{-}(x_{n},y_{n})},\\\
&y_{n+1}=\dfrac{y_{n}+\varphi_{2}y_{n}g_{+}(x_{n},y_{n})-\varphi_{2}\alpha_{2}y_{n}g_{-}(x_{n},y_{n})}{1+\varphi_{2}\beta_{2}g_{-}(x_{n},y_{n})}.\end{split}$
(33)

Therefore, the NSFD method (33) is positive for all $h>0$ if the condition
(32) holds. Assume that (32) is satisfied. Similarly to Theorem 5, we obtain:

######  Theorem 3.

The NSFD method (33) is positive and accurate of order $2$ if

$\begin{split}&\dfrac{\partial^{2}\varphi_{1}}{\partial
h^{2}}(x,y,0)=\dfrac{\partial f}{\partial x}(x,y)f(x,y)+\dfrac{\partial
f}{\partial y}(x,y)g(x,y)-2\beta_{1}f_{-}(x,y),\\\
&\dfrac{\partial^{2}\varphi_{2}}{\partial h^{2}}(x,y,0)=\dfrac{\partial
g}{\partial x}(x,y)f(x,y)+\dfrac{\partial g}{\partial
y}(x,y)g(x,y)-2\beta_{2}g_{-}(x,y).\end{split}$ (34)

Using the approach in [5], thanks to the Lyapunov indirect method [1, 22], we
can determine a stability thresholds $h^{*}=(h_{1}^{*},h_{2}^{*})$ for the
NSFD method (33), i.e., the method (33) preserves the local asymptotic
stability of the model (30) if

$0<\varphi_{1}<h_{1}^{*},\quad 0<\varphi_{2}<h_{2}^{*}.$ (35)

Hence, we only need to choose denominator functions satisfying the conditions
(34) and (35). Such denominator functions can be determined based on
Proposition 1.

Generally, we now consider the following general class of dynamical systems

$\dfrac{dy_{i}}{dt}=f_{i}(y_{1},y_{2},\ldots,y_{n}),\quad
y_{i}(0)=y_{i}^{0}\geq 0,\quad i=1,2,\ldots,n,$ (36)

where the right-hand function $f$ satisfies appropriate conditions, which
guarantee that solutions of (36) exist, are unique and positive for $t>0$.

Suppose that the functions $f_{i}$ ($i=1,2,\ldots,n$) admit a representation
of the form

$f_{i}(y)=f_{i}^{+}(y)+y_{i}f_{i}^{-}(y),$ (37)

where $f^{+},f^{-}$ are of $\mathcal{C}^{1}$ class, $f^{+}$ is positive and
$f^{-}$ is negative for all
$y=(y_{1},y_{2},\ldots,y_{n})\in\mathbb{R}^{n}_{+}$. Then, we propose the
following NSFD method

$\dfrac{y_{i}^{n+1}-y_{i}^{n}}{\varphi_{i}(y^{n},h)}=f_{i}^{+}(y^{n})+\alpha_{i}y^{n}f_{i}^{-}(y^{n})+\beta_{i}y_{i}^{n+1}f^{-}(y^{n}),$
(38)

with

$\alpha_{i}+\beta_{i}=1,\quad\alpha_{i}\leq 0,\quad\beta_{i}\geq 0.$

Using the approach used in Section 3, we can determine conditions for the
function $\varphi_{i}$ such that the NSFD scheme (38) is positive, elementary
stable and of second-order accuracy.

In general, it is not an easy task to verify whether the function $f$ admits a
representation of the form (37). However, in some cases, it is easy to find a
such representation of the function $f$, for example, for the classical Lotka-
Volterra system [1], we can write

$\begin{split}&\dfrac{dx}{dt}=ax-bxy=ax+bx(-y),\quad a,b>0,\\\
&\dfrac{dy}{dt}=-cy+exy=exy+y(-c),\quad c,e>0\end{split}$ (39)

or for the classical SIRS epidemic model [1], we have

$\begin{split}&\dfrac{dS}{dt}=-\dfrac{\beta}{N}SI+\mu R=\mu
R+S\bigg{(}-\dfrac{\beta}{N}I\bigg{)},\quad\beta,\mu>0,\\\
&\dfrac{dI}{dt}=\dfrac{\beta}{N}SI-\gamma
I=\dfrac{\beta}{N}SI+I(-\gamma),\quad\gamma>0,\\\ &\dfrac{dR}{dt}=\gamma I-\mu
R=\gamma I+R(-\mu).\end{split}$ (40)

##  5 Conclusions and remarks

In this work, we have constructed the NSFD method that is not only positive
and elementary stable but also accurate of order $2$. The theoretical
assertions are supported by a set of numerical simulations. The obtained
results not only resolve the contradiction between the dynamic consistency and
high-order accuracy of NSFD schemes but also improve the well-known results
constructed in [6, 11, 19, 23]. Furthermore, the constructed method is easy to
be implemented and can be used to solve a larger class of differential
equations arising in both theory and applications.

Our future research directions include the construction of high-order NSFD
methods that are not only elementary stable and positive but also preserve
other important qualitative features of solutions of differential equations,
such as global asymptotic stability and periodicity. Also, the construction of
dynamically consistent NSFD methods having high-order accuracy for general
multi-dimensional autonomous and non-autonomous dynamical systems governed by
ODEs and fractional-order ODEs will be studied.

## References

## References

  * [1] L. J. S. Allen, An Introduction to Mathematical Biology, Prentice Hall, 2007.
  * [2] R. Anguelov, 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.
  * [3] R. Anguelov, J. M. -S. Lubuma, Nonstandard finite difference method by nonlocal approximation, Mathematics and Computers in Simulation 61(2003) 465-475.
  * [4] B. M. Chen-Charpentier, D. T. Dimitrov, H. V. Kojouharov, Combined nonstandard numerical methods for ODEs with polynomial right-hand sides, Mathematics and Computers in Simulation 73(2006) 105-113.
  * [5] J. Cresson, A. Szafranska, Discrete and continuous fractional persistence problems – the positivity property and applications, Communications in Nonlinear Science and Numerical Simulation 44(2017) 424-448.
  * [6] M. Gupta, J. M. Slezak, F. Alalhareth, S. Roy, H. V. Kojouharov, Second-order Nonstandard Explicit Euler Method, AIP Conference Proceedings 2302(2020) 110003.
  * [7] Quang A Dang, Manh Tuan Hoang, Positive and elementary stable explicit nonstandard Runge-Kutta methods for a class of autonomous dynamical systems, International Journal of Computer Mathematics 97 (2020) 2036-2054.
  * [8] D. T. Dimitrov, H. V. Kojouharov, Nonstandard finite-difference schemes for general two-dimensional autonomous dynamical systems, Applied Mathematics Letters 18(2005) 769-774.
  * [9] D. T. Dimitrov, H. V. Kojouharov, Stability-preserving Finite Difference Methods for General Multi-dimensional Autonomous Dynamical Systems, International Journal of Numerical Analysis and Modeling 4(2007) 280-290.
  * [10] Z. Horváth, Positivity of Runge-Kutta and diagonally split Runge-Kutta methods, Applied Numerical Mathematics 28(1998) 309-326.
  * [11] H. V. Kojouharov, S. Roy, M. Gupta, F. Alalhareth, J. M. Slezak, A second-order modified nonstandard theta method for one-dimensional autonomous differential equations, Applied Mathematics Letters 112(2021) 106775.
  * [12] 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, Journal of Computational and Applied Mathematics 317(2017) 146-156.
  * [13] 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, Variable step length algorithms with high-order extrapolated non-standard finite difference schemes for a SEIR model, Journal of Computational and Applied Mathematics 330(2018) 848-854.
  * [14] R. E. Mickens, Nonstandard Finite Difference Models of Differential Equations, World Scientific, 1993.
  * [15] R. E. Mickens, Applications of Nonstandard Finite Difference Schemes, World Scientific, 2000.
  * [16] R. E. Mickens, Advances in the Applications of Nonstandard Finite Difference Schemes, World Scientific, 2005.
  * [17] R. E. Mickens, Nonstandard Finite Difference Schemes for Differential Equations, Journal of Difference Equations and Applications 8(2002) 823-847.
  * [18] R. E. Mickens, Nonstandard Finite Difference Schemes: Methodology and Applications, World Scientific, 2020.
  * [19] R. E. Mickens, Discretizations of nonlinear differential equations using explicit nonstandard methods, Journal of Computational and Applied Mathematics 110(1999) 181-185.
  * [20] K. C. Patidar, On the use of nonstandard finite difference methods, Journal of Difference Equations and Applications 11(2005) 735-758.
  * [21] K. C. Patidar, Nonstandard finite difference methods: recent trends and further developments, Journal of Difference Equations and Applications 22(2016) 817-849.
  * [22] A. Stuart, A. R. Humphries, Dynamical systems and numerical analysis, Cambridge University Press, 1998.
  * [23] D. Wood, H. V. Kojouharov, A class of nonstandard numerical methods for autonomous dynamical systems, Applied Mathematics Letters 50(2015) 78-82.