Nonstandard finite difference schemes for a general predator-prey system

Quang A Dang; Manh Tuan Hoang

arXiv:1701.05663·math.NA·January 23, 2017·J. Comput. Sci.

Nonstandard finite difference schemes for a general predator-prey system

Quang A Dang, Manh Tuan Hoang

PDF

TL;DR

This paper develops a nonstandard finite difference scheme for predator-prey models, ensuring dynamic consistency and stability, with numerical simulations validating the theoretical findings.

Contribution

It introduces a novel NSFD scheme for predator-prey systems that preserves stability and dynamics across any step size, including non-hyperbolic equilibria.

Findings

01

NSFD scheme maintains dynamic consistency with continuous models

02

Global stability of non-hyperbolic equilibrium proved

03

Numerical simulations confirm theoretical stability results

Abstract

In this paper we transform a continuous-time predator-prey system with general functional response and recruitment for both species into a discrete-time model by nonstandard finite difference scheme (NSFD). The NSFD model shows complete dynamic consistency with its continuous counterpart for any step size. Especially, the global stability of a non-hyperbolic equilibrium point in a particular case of parameters is proved by the Lyapunov stability theorem. The performed numerical simulations confirmed the validity of the obtained theoretical results.

Tables1

Table 1. Table 1 : The sufficient conditions for dynamical consistency

$(m_{1}, m_{2})$	Set of equilibria	Positivity	Stability
$m_{1} \geq r (0)$ and $m_{2} \geq s (0)$	*	(4)	(14)
$m_{1} < r (0)$ and $m_{2} > s (0) + c K ϕ (K)$	*	(4)	(15)
$m_{1} > r (0) - M ϕ (0)$ and $m_{2} < s (0)$	*	(4)	(16)
$m_{1} < r (0) - M ϕ (0)$ and $m_{2} < s (0)$	*	(4)	(17)
$m_{1} < r (0)$ and $s (0) < m_{2} < s (0) + c K ϕ (K)$	*	(4)	(17)

Equations98

\frac{d x}{d t} = b x (x - 1) - a g (x) x y, \frac{d y}{d t} = g (x) x y - d y,

\frac{d x}{d t} = b x (x - 1) - a g (x) x y, \frac{d y}{d t} = g (x) x y - d y,

\frac{d x}{d t} = p (x) - a f (x, y) y, \frac{d y}{d t} = f (x, y) y - μ (y),

\frac{d x}{d t} = p (x) - a f (x, y) y, \frac{d y}{d t} = f (x, y) y - μ (y),

\frac{d x}{d t} = x (x - 1) - \frac{α x y}{1 + β x + μ y}, \frac{d y}{d t} = \frac{E x y}{1 + β x + μ y} - D y,

\frac{d x}{d t} = x (x - 1) - \frac{α x y}{1 + β x + μ y}, \frac{d y}{d t} = \frac{E x y}{1 + β x + μ y} - D y,

\begin{split}\dot{x}(t)&=x(t)f(x(t),y(t))=x(t)\big{[}r(x(t))-y(t)\phi(x(t))-m_{1}\big{]},\\ \dot{y}(t)&=y(t)g(x(t),y(t))=y(t)\big{[}s(y(t))+cx(t)\phi(x(t))-m_{2}\big{]},\end{split}

\begin{split}\dot{x}(t)&=x(t)f(x(t),y(t))=x(t)\big{[}r(x(t))-y(t)\phi(x(t))-m_{1}\big{]},\\ \dot{y}(t)&=y(t)g(x(t),y(t))=y(t)\big{[}s(y(t))+cx(t)\phi(x(t))-m_{2}\big{]},\end{split}

\forall x \geq 0, r (x) > 0, r^{'} (x) < 0, [x r (x)]^{'} \geq 0, and x \to \infty lim r (x) = 0, \forall y \geq 0, s (y) > 0, s^{'} (y) < 0, [y s (y)]^{'} \geq 0, and y \to \infty lim s (y) = 0, \forall x \geq 0, ϕ (x) > 0, ϕ^{'} (x) \leq 0, and [x ϕ (x)]^{'} \geq 0.

\forall x \geq 0, r (x) > 0, r^{'} (x) < 0, [x r (x)]^{'} \geq 0, and x \to \infty lim r (x) = 0, \forall y \geq 0, s (y) > 0, s^{'} (y) < 0, [y s (y)]^{'} \geq 0, and y \to \infty lim s (y) = 0, \forall x \geq 0, ϕ (x) > 0, ϕ^{'} (x) \leq 0, and [x ϕ (x)]^{'} \geq 0.

cx^{*}\phi(x^{*})+s\Big{(}\dfrac{r(x^{*})-m_{1}}{\phi(x^{*})}\Big{)}-m_{2}=0,

cx^{*}\phi(x^{*})+s\Big{(}\dfrac{r(x^{*})-m_{1}}{\phi(x^{*})}\Big{)}-m_{2}=0,

y^{*} = \frac{r ( x ^{*} ) - m _{1}}{ϕ ( x ^{*} )},

y^{*} = \frac{r ( x ^{*} ) - m _{1}}{ϕ ( x ^{*} )},

\frac{x _{k + 1} - x _{k}}{φ ( h )} \frac{y _{k + 1} - y _{k}}{φ ( h )} α_{j} + α_{j + 1} = α_{1} x_{k} r (x_{k}) + α_{2} x_{k + 1} r (x_{k}) - α_{3} x_{k} y_{k} ϕ (x_{k}) - α_{4} x_{k + 1} y_{k} ϕ (x_{k}) - α_{5} m_{1} x_{k} - α_{6} m_{1} x_{k + 1}, = β_{1} y_{k} s (y_{k}) + β_{2} y_{k + 1} s (y_{k}) + c β_{3} x_{k} y_{k} ϕ (x_{k}) + c β_{4} x_{k} y_{k + 1} ϕ (x_{k}) - β_{5} m_{2} y_{k} - β_{6} m_{2} y_{k + 1}, = β_{j} + β_{j + 1} = 1, j = 1, 3, 5; φ = h + O (h^{2}), h \to 0.

\frac{x _{k + 1} - x _{k}}{φ ( h )} \frac{y _{k + 1} - y _{k}}{φ ( h )} α_{j} + α_{j + 1} = α_{1} x_{k} r (x_{k}) + α_{2} x_{k + 1} r (x_{k}) - α_{3} x_{k} y_{k} ϕ (x_{k}) - α_{4} x_{k + 1} y_{k} ϕ (x_{k}) - α_{5} m_{1} x_{k} - α_{6} m_{1} x_{k + 1}, = β_{1} y_{k} s (y_{k}) + β_{2} y_{k + 1} s (y_{k}) + c β_{3} x_{k} y_{k} ϕ (x_{k}) + c β_{4} x_{k} y_{k + 1} ϕ (x_{k}) - β_{5} m_{2} y_{k} - β_{6} m_{2} y_{k + 1}, = β_{j} + β_{j + 1} = 1, j = 1, 3, 5; φ = h + O (h^{2}), h \to 0.

α_{1} \geq 0, α_{2} \leq 0, α_{3} \leq 0, α_{4} \geq 0, α_{5} \leq 0, α_{6} \geq 0, β_{1} \geq 0, β_{2} \leq 0, β_{3} \geq 0, β_{4} \leq 0, β_{5} \leq 0, β_{6} \geq 0.

α_{1} \geq 0, α_{2} \leq 0, α_{3} \leq 0, α_{4} \geq 0, α_{5} \leq 0, α_{6} \geq 0, β_{1} \geq 0, β_{2} \leq 0, β_{3} \geq 0, β_{4} \leq 0, β_{5} \leq 0, β_{6} \geq 0.

x_{k + 1} y_{k + 1} = F (x_{k}, y_{k}) := \frac{x _{k} + φ α _{1} x _{k} r ( x _{k} ) - φ α _{3} x _{k} y _{k} ϕ ( x _{k} ) - φ α _{5} m _{1} x _{k}}{1 - φ α _{2} r ( x _{k} ) + φ α _{4} y _{k} ϕ ( x _{k} ) + φ α _{6} m _{1}}, = G (x_{k}, y_{k}) := \frac{y _{k} + φ β _{1} y _{k} s ( y _{k} ) + φ β _{3} c x _{k} y _{k} ϕ ( x _{k} ) - φ β _{5} m _{2} y _{k}}{1 - φ β _{2} s ( y _{k} ) - φ β _{4} c x _{k} ϕ ( x _{k} ) + φ β _{6} m _{2}} .

x_{k + 1} y_{k + 1} = F (x_{k}, y_{k}) := \frac{x _{k} + φ α _{1} x _{k} r ( x _{k} ) - φ α _{3} x _{k} y _{k} ϕ ( x _{k} ) - φ α _{5} m _{1} x _{k}}{1 - φ α _{2} r ( x _{k} ) + φ α _{4} y _{k} ϕ ( x _{k} ) + φ α _{6} m _{1}}, = G (x_{k}, y_{k}) := \frac{y _{k} + φ β _{1} y _{k} s ( y _{k} ) + φ β _{3} c x _{k} y _{k} ϕ ( x _{k} ) - φ β _{5} m _{2} y _{k}}{1 - φ β _{2} s ( y _{k} ) - φ β _{4} c x _{k} ϕ ( x _{k} ) + φ β _{6} m _{2}} .

x^{*} = F (x^{*}, y^{*}), y^{*} = G (x^{*}, y^{*}),

x^{*} = F (x^{*}, y^{*}), y^{*} = G (x^{*}, y^{*}),

x^{*}\big{[}r(x^{*})-y^{*}\phi(x^{*})-m_{1}\big{]}=0,\qquad y^{*}\big{[}s(y^{*})+cx^{*}\phi(x^{*})-m_{2}\big{]}=0.

x^{*}\big{[}r(x^{*})-y^{*}\phi(x^{*})-m_{1}\big{]}=0,\qquad y^{*}\big{[}s(y^{*})+cx^{*}\phi(x^{*})-m_{2}\big{]}=0.

x_{k + 1} = f (x_{k}, y_{k}), y_{k + 1} = g (x_{k}, y_{k}),

x_{k + 1} = f (x_{k}, y_{k}), y_{k + 1} = g (x_{k}, y_{k}),

(i) det (J) < 1, (ii) 1 - T r (J) + d e t (J) > 0, (iii) 1 + T r (J) + det (J) > 0.

(i) det (J) < 1, (ii) 1 - T r (J) + d e t (J) > 0, (iii) 1 + T r (J) + det (J) > 0.

(i) det (J) > 1, (ii) 1 - T r (J) + d e t (J) < 0, (iii) 1 + T r (J) + det (J) < 0.

(i) det (J) > 1, (ii) 1 - T r (J) + d e t (J) < 0, (iii) 1 + T r (J) + det (J) < 0.

J (P_{0}^{*}) = \frac{1 + φ α _{1} r ( 0 ) - φ α _{5} m _{1}}{1 - φ α _{2} r ( 0 ) + φ m _{1} α _{6}} 0 0 \frac{1 + φ β _{1} s ( 0 ) - φ β _{5} m _{2}}{1 - φ β _{2} s ( 0 ) + φ β _{6} m _{2}} .

J (P_{0}^{*}) = \frac{1 + φ α _{1} r ( 0 ) - φ α _{5} m _{1}}{1 - φ α _{2} r ( 0 ) + φ m _{1} α _{6}} 0 0 \frac{1 + φ β _{1} s ( 0 ) - φ β _{5} m _{2}}{1 - φ β _{2} s ( 0 ) + φ β _{6} m _{2}} .

λ_{1} = \frac{1 + φ α _{1} r ( 0 ) - φ α _{5} m _{1}}{1 - φ α _{2} r ( 0 ) + φ m _{1} α _{6}}, λ_{2} = \frac{1 + φ β _{1} s ( 0 ) - φ β _{5} m _{2}}{1 - φ β _{2} s ( 0 ) + φ β _{6} m _{2}} .

λ_{1} = \frac{1 + φ α _{1} r ( 0 ) - φ α _{5} m _{1}}{1 - φ α _{2} r ( 0 ) + φ m _{1} α _{6}}, λ_{2} = \frac{1 + φ β _{1} s ( 0 ) - φ β _{5} m _{2}}{1 - φ β _{2} s ( 0 ) + φ β _{6} m _{2}} .

∣ λ_{1} ∣ - 1 = φ \frac{r ( 0 ) - m _{1}}{1 - φ α _{2} r ( 0 ) + φ m _{1} α _{6}}, ∣ λ_{2} ∣ - 1 = φ \frac{s ( 0 ) - m _{2}}{1 - φ β _{2} s ( 0 ) + φ β _{6} m _{2}},

∣ λ_{1} ∣ - 1 = φ \frac{r ( 0 ) - m _{1}}{1 - φ α _{2} r ( 0 ) + φ m _{1} α _{6}}, ∣ λ_{2} ∣ - 1 = φ \frac{s ( 0 ) - m _{2}}{1 - φ β _{2} s ( 0 ) + φ β _{6} m _{2}},

V (x_{k}, y_{k}) := α x_{k} y_{k} + β x_{k}^{2} + γ x_{k} + δ y_{k}, (x_{k}, y_{k}) \in R_{+}^{2},

V (x_{k}, y_{k}) := α x_{k} y_{k} + β x_{k}^{2} + γ x_{k} + δ y_{k}, (x_{k}, y_{k}) \in R_{+}^{2},

Δ V (x_{k}, y_{k}) = V (x_{k + 1}, y_{k + 1}) - V (x_{k}, y_{k}) < 0 \forall (x_{k}, y_{k}) \in R_{+}^{2} \ {(0, 0)} .

Δ V (x_{k}, y_{k}) = V (x_{k + 1}, y_{k + 1}) - V (x_{k}, y_{k}) < 0 \forall (x_{k}, y_{k}) \in R_{+}^{2} \ {(0, 0)} .

Δ V (x_{k}, y_{k}) = α (x_{k + 1} y_{k + 1} - x_{k} y_{k}) + β (x_{k + 1}^{2} - x_{k}^{2}) + γ (x_{k + 1} - x_{k}) + δ (y_{k + 1} - y_{k}) .

Δ V (x_{k}, y_{k}) = α (x_{k + 1} y_{k + 1} - x_{k} y_{k}) + β (x_{k + 1}^{2} - x_{k}^{2}) + γ (x_{k + 1} - x_{k}) + δ (y_{k + 1} - y_{k}) .

x_{k + 1} y_{k + 1} = x_{k} + φ \frac{x _{k} [ r ( x _{k} ) - m _{1} ] - x _{k} y _{k} ϕ ( x _{k} )}{1 - φ α _{2} r ( x _{k} ) + φ α _{4} y _{k} ϕ ( x _{k} ) + φ α _{6} m _{1}}, = y_{k} + φ \frac{y _{k} [ s ( y _{k} ) - m _{2} ] + c x _{k} y _{k} ϕ ( x _{k} )}{1 - φ β _{2} s ( y _{k} ) - φ β _{4} c x _{k} ϕ ( x _{k} ) + φ β _{6} m _{2}} .

x_{k + 1} y_{k + 1} = x_{k} + φ \frac{x _{k} [ r ( x _{k} ) - m _{1} ] - x _{k} y _{k} ϕ ( x _{k} )}{1 - φ α _{2} r ( x _{k} ) + φ α _{4} y _{k} ϕ ( x _{k} ) + φ α _{6} m _{1}}, = y_{k} + φ \frac{y _{k} [ s ( y _{k} ) - m _{2} ] + c x _{k} y _{k} ϕ ( x _{k} )}{1 - φ β _{2} s ( y _{k} ) - φ β _{4} c x _{k} ϕ ( x _{k} ) + φ β _{6} m _{2}} .

x_{k + 1} y_{k + 1} - x_{k} y_{k} = ξ_{y_{k}} (x_{k + 1} - x_{k}) + ξ_{x_{k}} (y_{k + 1} - y_{k}) \leq y_{k} (x_{k + 1} - x_{k}) + x_{k} (y_{k + 1} - y_{k}),

x_{k + 1} y_{k + 1} - x_{k} y_{k} = ξ_{y_{k}} (x_{k + 1} - x_{k}) + ξ_{x_{k}} (y_{k + 1} - y_{k}) \leq y_{k} (x_{k + 1} - x_{k}) + x_{k} (y_{k + 1} - y_{k}),

x_{k + 1}^{2} - x_{k}^{2} \leq x_{k} x_{k + 1} - x_{k}^{2} = x_{k} (x_{k + 1} - x_{k}) .

x_{k + 1}^{2} - x_{k}^{2} \leq x_{k} x_{k + 1} - x_{k}^{2} = x_{k} (x_{k + 1} - x_{k}) .

Δ V (x_{k}, y_{k}) \leq (α y_{k} + β x_{k} + γ) (x_{k + 1} - x_{k}) + (α x_{k} + δ) (y_{k + 1} - y_{k}) = (α y_{k} + β x_{k} + γ) φ \frac{x _{k} [ r ( x _{k} ) - m _{1} ] - x _{k} y _{k} ϕ ( x _{k} )}{1 - φ α _{2} r ( x _{k} ) + φ α _{4} y _{k} ϕ ( x _{k} ) + φ α _{6} m _{1}} + (α x_{k} + δ) φ \frac{y _{k} [ s ( y _{k} ) - m _{2} ] + c x _{k} y _{k} ϕ ( x _{k} )}{1 - φ β _{2} s ( y _{k} ) - φ β _{4} c x _{k} ϕ ( x _{k} ) + φ β _{6} m _{2}} \leq φ x_{k} y_{k} ϕ (x_{k}) Q (x_{k}, y_{k}),

Δ V (x_{k}, y_{k}) \leq (α y_{k} + β x_{k} + γ) (x_{k + 1} - x_{k}) + (α x_{k} + δ) (y_{k + 1} - y_{k}) = (α y_{k} + β x_{k} + γ) φ \frac{x _{k} [ r ( x _{k} ) - m _{1} ] - x _{k} y _{k} ϕ ( x _{k} )}{1 - φ α _{2} r ( x _{k} ) + φ α _{4} y _{k} ϕ ( x _{k} ) + φ α _{6} m _{1}} + (α x_{k} + δ) φ \frac{y _{k} [ s ( y _{k} ) - m _{2} ] + c x _{k} y _{k} ϕ ( x _{k} )}{1 - φ β _{2} s ( y _{k} ) - φ β _{4} c x _{k} ϕ ( x _{k} ) + φ β _{6} m _{2}} \leq φ x_{k} y_{k} ϕ (x_{k}) Q (x_{k}, y_{k}),

Q (x_{k}, y_{k}) := \frac{c ( α x _{k} + δ )}{1 - φ β _{2} s ( y _{k} ) - φ β _{4} c x _{k} ϕ ( x _{k} ) + φ β _{6} m _{2}} - \frac{α y _{k} + β x _{k} + γ}{1 - φ α _{2} r ( x _{k} ) + φ α _{4} y _{k} ϕ ( x _{k} ) + φ α _{6} m _{1}} .

Q (x_{k}, y_{k}) := \frac{c ( α x _{k} + δ )}{1 - φ β _{2} s ( y _{k} ) - φ β _{4} c x _{k} ϕ ( x _{k} ) + φ β _{6} m _{2}} - \frac{α y _{k} + β x _{k} + γ}{1 - φ α _{2} r ( x _{k} ) + φ α _{4} y _{k} ϕ ( x _{k} ) + φ α _{6} m _{1}} .

Q (x_{k}, y_{k}) [1 - φ α_{2} r (x_{k}) + φ α_{4} y_{k} ϕ (x_{k}) + φ α_{6} m_{1}] [1 - φ β_{2} s (y_{k}) - φ β_{4} c x_{k} ϕ (x_{k}) + φ β_{6} m_{2}] = τ_{1} x_{k} + τ_{2} φ x_{k} + τ_{3} φ c α x_{k} y_{k} ϕ (x_{k}) + τ_{4} + φ τ_{5} + τ_{6} φ y_{k} - Q_{1} (x_{k}, y_{k}),

Q (x_{k}, y_{k}) [1 - φ α_{2} r (x_{k}) + φ α_{4} y_{k} ϕ (x_{k}) + φ α_{6} m_{1}] [1 - φ β_{2} s (y_{k}) - φ β_{4} c x_{k} ϕ (x_{k}) + φ β_{6} m_{2}] = τ_{1} x_{k} + τ_{2} φ x_{k} + τ_{3} φ c α x_{k} y_{k} ϕ (x_{k}) + τ_{4} + φ τ_{5} + τ_{6} φ y_{k} - Q_{1} (x_{k}, y_{k}),

τ_{1} := c α - β, τ_{2} := - α_{2} c α r (x_{k}) + α_{6} c α m_{1} - β β_{6} m_{2}, τ_{3} := α_{4} + β_{4}, τ_{4} := cδ - γ, τ_{5} := - α_{2} cδ r (x_{k}) + cδ α_{6} m_{1} - β_{6} γ m_{2}, τ_{6} := α_{4} cδ ϕ (x_{k}) - α β_{6} m_{2}, Q_{1} (x_{k}, y_{k}) := β x_{k} [- φ β_{2} s (y_{k}) - φ β_{4} c x_{k} ϕ (x_{k})] + α y_{k} [1 - φ β_{2} s (y_{k})] + γ [- φ β_{2} s (y_{k}) - φ β_{4} c x_{k} ϕ (x_{k})] .

τ_{1} := c α - β, τ_{2} := - α_{2} c α r (x_{k}) + α_{6} c α m_{1} - β β_{6} m_{2}, τ_{3} := α_{4} + β_{4}, τ_{4} := cδ - γ, τ_{5} := - α_{2} cδ r (x_{k}) + cδ α_{6} m_{1} - β_{6} γ m_{2}, τ_{6} := α_{4} cδ ϕ (x_{k}) - α β_{6} m_{2}, Q_{1} (x_{k}, y_{k}) := β x_{k} [- φ β_{2} s (y_{k}) - φ β_{4} c x_{k} ϕ (x_{k})] + α y_{k} [1 - φ β_{2} s (y_{k})] + γ [- φ β_{2} s (y_{k}) - φ β_{4} c x_{k} ϕ (x_{k})] .

\begin{split}&\max\Big{\{}c,\dfrac{-\alpha_{2}cr(0)+\alpha_{6}cm_{1}}{\beta_{6}m_{2}}\Big{\}}<\dfrac{\beta}{\alpha},\quad\max\Big{\{}c,\dfrac{-\alpha_{2}cr(0)+c\alpha_{6}m_{1}}{\beta_{6}m_{2}}\Big{\}}<\dfrac{\gamma}{\delta},\\ &\dfrac{c\alpha_{4}\phi(0)}{\beta_{6}m_{2}}<\dfrac{\alpha}{\delta},\qquad\alpha_{4}+\beta_{4}<0.\end{split}

\begin{split}&\max\Big{\{}c,\dfrac{-\alpha_{2}cr(0)+\alpha_{6}cm_{1}}{\beta_{6}m_{2}}\Big{\}}<\dfrac{\beta}{\alpha},\quad\max\Big{\{}c,\dfrac{-\alpha_{2}cr(0)+c\alpha_{6}m_{1}}{\beta_{6}m_{2}}\Big{\}}<\dfrac{\gamma}{\delta},\\ &\dfrac{c\alpha_{4}\phi(0)}{\beta_{6}m_{2}}<\dfrac{\alpha}{\delta},\qquad\alpha_{4}+\beta_{4}<0.\end{split}

α_{4} + β_{4} < 0,

α_{4} + β_{4} < 0,

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Nonstandard finite difference schemes

for a general predator-prey system

Quang A Dang

[email protected]

Manh Tuan Hoang

[email protected]

Center for Informatics and Computing, Vietnam Academy of Science and Technology (VAST),

18 Hoang Quoc Viet, Cau Giay, Hanoi, Vietnam

Institute of Information Technology, Vietnam Academy of Science and Technology (VAST),

18 Hoang Quoc Viet, Cau Giay, Hanoi, Vietnam

Abstract

In this paper we transform a continuous-time predator-prey system with general functional response and recruitment for both species into a discrete-time model by nonstandard finite difference scheme (NSFD). The NSFD model shows complete dynamic consistency with its continuous counterpart for any step size. Especially, the global stability of a non-hyperbolic equilibrium point in a particular case of parameters is proved by the Lyapunov stability theorem. The performed numerical simulations confirmed the validity of the obtained theoretical results.

keywords:

Predator-Prey system; Nonstandard finite-difference scheme; Dynamically consistent; Lyapunov stability theorem; Global stability.

††journal: Elsevier

1 Introduction

Predator-prey systems are among the most interested topics in mathematical biology and ecology. Their dynamics continues to attract attention from both applied mathematicians and ecologists because of its universal existence and importance [1, 5, 6, 18]. The majority of the results in this topic is mainly concentrated on the study of the qualitative aspects of the continuous systems described by systems of differential equations. It is worth to mention some recent typical works such as [16, 17, 20, 29, 33, 34]…In that time the conversion of the continuous systems to discrete systems preserving the properties of the original continuous systems is of great importance. This problem has attracted attention from many researchers and a lot of results are reached for some dynamical systems. Nevertheless, some continuous models are studied fully from theoretical point of view but their discrete counterparts still are not investigated.

Biological systems including predator-prey systems often are described by ordinary or partial differential equations. There are many ways for converting continuous models to discrete counterparts. The most popular way for this purpose is to use standard difference methods such as Euler, Runge-Kutta schemes. However, in many nonlinear problems the standard difference schemes reveal a serious drawback which is called ”numerical instabilities” [23, 24, 26]. Under this concept we have in mind the phenomena when the discrete models, for example, the difference schemes, do not preserve properties of the corresponding differential equations. In [23, 24, 25, 26] Mickens showed many examples and analysed the numerical instabilities when using standard difference schemes. In general, standard difference schemes preserve the properties of the differential equations only in the case if the discretization parameter $h$ is sufficiently small. Therefore, when studying dynamical models in large time intervals the selection of small time steps will requires very large computational effort, so these discrete models are inefficient. Besides, for some special dynamical problems standard difference schemes cannot preserve the properties of the problems for any step sizes.

In order to overcome the numerical instabilities phenomena in 1989 Mickens [22] introduced the concept Nonstandard Finite Difference (NSDF) schemes and after that has developed NSDF methods in many works, such as [23, 24, 25, 26]. According to Mickens, NSDF schemes are those constructed following a set of five basic rules. The NSDF schemes preserve main properties of the differential counterparts, such as positivity, monotonicity, periodicity, stability and some other invariants including energy and geometrical shapes. It should be emphasized that NSFD schemes can preserve all properties of the continuous models for any discretization parameters. The discrete models with these properties are called dynamically consistent [2, 3, 8, 10, 11, 12, 13, 14, 26, 27, 30].

Up to date NSFD schemes became a power and efficient tool for simulating dynamical systems, especially in converting continuous models to dynamically consistent discrete counterparts [2, 3, 7, 9, 10, 11, 12, 13, 25, 28, 31, 32, 35]. The majority of these models are met in physics, mechanics, chemistry, biology, epidemiology, finance, …with complicated dynamical behaviour. For predator-prey systems, some NSFD schemes, which are dynamically consistent with them, are constructed, such as [4, 9, 12, 14]. Below we mention some models.

In 2006 Dimitrov and Kojouharov [12] used NSFD schemes for constructing discrete model for the general Rosenzweig-MacArthur predator-prey model with a logistic intrinsic growth of the prey population. The model has the form

[TABLE]

where $x$ and $y$ represent the prey and predator population sizes, respectively, $b>0$ represents the intrinsic growth rate of the prey, $a>0$ stands for the capturing rate and $d>0$ is the predator death rate. Under some assumptions for the function $g(x)$ the obtained discrete model preserves the positivity of solutions and locally asymptotical stability of the set of equilibrium points. The finite difference method is called positive and elementary stable nonstandard (PESN) method.

Next, in 2008 Dimitrov and Kojouharov [14] constructed a PESN scheme for predator-prey models with general functional response of the form

[TABLE]

where functions $p(x)$ and $\mu(y)$ describe the intrinsic growth rate of the prey and the mortality rate of the predator, respectively, the function $f(x,y)$ is called ”functional response” and represents the per capita predator ”feeding rate” per unit time. Recently, in 2016, Bairagi and Biswas [4] used PESN scheme for converting a predator-prey model with Beddington-DeAngelis Functional Response to a dynamically consistent discrete model. The model under consideration is of the form

[TABLE]

where $\alpha,\beta,\mu,E,D$ are positive constants. This model is a particular case of the model considered by Dimitrov and Kojouharov in [14]. An another interesting predator-prey model, which should be mentioned is the harvesting Leslie-Gower predator-prey model. In the recent paper [9] the authors also constructed NSFD scheme preserving the positivity of solutions and asymptotical stability for this model.

It should be emphasized that all the equilibrium points in the mentioned above predator-prey models are hyperbolic, therefore, for establishing the stability, it suffices to consider eigenvalues of the Jacobian matrix of the linearized system around equilibrium points.

Aiming at the conversion of continuous predator-prey models studied fully from the theoretical point of view to dynamically consistent discrete models, in the present paper we consider a mathematical model for a predator-prey system with general functional response and recruitment, which also includes capture on both species [21]. The results of qualitative aspects of the model are given in [21]. One important difference of this model in comparison with the above models is that the model [21] has one non-hyperbolic equilibrium point.

The paper is organized as follows. In Section 2 we recall from [21] the predator-prey system under consideration with the theoretical results of the existence of equilibrium points and their stability properties. In Section 3 we propose NSFD schemes and study their positivity and existence of equilibrium points. Next, Section 4 is devoted to the stability analysis of the equilibrium points. Some numerical simulations for demonstrating the validity of the theoretical results obtained in the previous section are given in Section 5. Finally, Section 6 is Conclusion.

2 Mathematical model

We consider a mathematical model for a predator-prey system with general functional response and recruitment for both species [21]. The model is described by the system of nonlinear differential equations

[TABLE]

where:

$x(t)$ and $y(t)$ are prey population and predator population, respectively;

2.

$r(x)$ and $s(y)$ are the per capita recruitment rates of prey and predators, respectively;

3.

$xy\phi(x)$ is the predator response, and $x\phi(x)$ is the number of prey consumed per predator in a unit of time;

4.

$c$ is a constant named conversion efficiency of prey into predators, generally $0<c<1$ , and $cxy\phi(x)$ is the predator numerical response;

5.

$m_{1}>0$ and $m_{2}>$ are the total mortality rates of prey and predators, respectively.

From the biological significance we have

[TABLE]

It is easy to check that the region $\Omega=\mathbb{R}_{+}^{2}$ is a positive invariant set for the system (1). Other qualitative properties of the model (1) including the existence of equilibrium points, their stability, are studied in detail in [21].

In order to easily track the paper, below we mention these properties.

Theorem 1

*(Existence of equilibrium points) [21, Proposition 1]

System (1) has four distinct kinds of possible equilibrium points in $\Omega$ :*

(i)

A trivial equilibrium point $P_{0}^{*}=(x_{0}^{*},y_{0}^{*})=(0,0)$ , for all the values of the parameter. 2. (ii)

An equilibrium point of the form $P_{1}^{*}=(x_{1}^{*},y_{1}^{*})=(K,0)$ , with $r(K)=m_{1}$ , if and only if $m_{1}<r(0)$ . 3. (iii)

An equilibrium point of the form $P_{2}^{*}=(x_{2}^{*},y_{2}^{*})=(0,M)$ , with $s(M)=m_{2}$ , if and only if $m_{2}<s(0)$ . 4. (iv)

An equilibrium point of the form $P_{3}^{*}=(x_{3}^{*},y_{3}^{*})=(x^{*},y^{*})$ , where $x^{*}$ satisfies the equation

[TABLE]

and $y^{*}$ is given, as a function of $x^{*}$ , by

[TABLE]

if and only if $(m_{1},m_{2})$ verifies $m_{1}<r(0)-M\phi(0)$ and $m_{2}<s(0)$ or $m_{1}<r(0)$ and $s(0)<m_{2}<s(0)+cK\phi(K)$ .

Theorem 2

(Stability analysis)

(i)

If $m_{1}>r(0)$ and $m_{2}>s(0)$ , then the extinction equilibrium point $P_{0}^{*}=(0,0)$ is locally asymptotically stable, and unstable otherwise. 2. (ii)

If $m_{1}\geq r(0)$ and $m_{2}\geq s(0)$ then the extinction equilibrium point $P_{0}^{*}$ is globally asymptotically stable. 3. (iii)

If $m_{1}<r(0)$ and $m_{2}>s(0)+cK\phi(K)$ , then the equilibrium point of the form $P_{1}^{*}=(K,0)$ is locally asymptotically stable, and unstable otherwise.The equilibrium point $(K,0)$ shall be called the equilibrium point of extinction of the predator species. 4. (iv)

If $m_{1}>r(0)-M\phi(0)$ and $m_{2}<s(0)$ , then the equilibrium point of the form $P_{2}^{*}=(0,M)$ is locally asymptotically stable, and unstable otherwise. The equilibrium point $(0,M)$ shall be called the equilibrium point of extinction of the prey species. 5. (v)

If an equilibrium point of the form $P_{3}^{*}=(x^{*},y^{*})$ belongs to $\Omega$ , then it is locally asymptotically stable. This equilibrium point $(x^{*},y^{*})$ shall be named the ecological stability equilibrium.

In general, the property of stability of the set of equilibria of differential equations plays an essential role in the study of asymptotical behavior of the solutions of differential equations. The construction of difference schemes, which preserve the stability of the equilibrium points, is important in numerical simulation of differential equations. The difference schemes with this stability property is called elementary stable schemes [2, 10, 11, 35]. There are many works concerning the elementary stable schemes. The typical results are for general dynamical systems [10, 11, 19] and for other specific systems [14, 30, 31, 32, 35] …. One popular approach to the elementary stability is the investigation of Jacobian matrices of the discrete models at the equilibria, namely, determination of conditions ensuring that all eigenvalues of Jacobian matrices have moduli less or equal 1. This is the necessary and sufficient condition for hyperbolic equilibrium points to be locally stable [1, 18]. The mentioned above approach has the following weaknesses and limitations:

It is applicable when all the equilibrium points are hyperbolic. To our best knowledge, at present no results on NSFD schemes preserving the stability of non-hyperbolic equilibrium points are available. 2. 2.

The consideration of Jacobian only guarantees the local stability meanwhile many models have the global stability.

Return to the model (1). A difficulty is that when $m_{1}=r(0)$ or $m_{2}=s(0)$ the equilibrium point $P_{0}^{*}=(0,0)$ becomes non-hyperbolic. Therefore, it is impossible to study its stability via the eigenvalues of Jacobian $J(P_{0}^{*})$ . Consequently, we cannot directly apply the results concerning the elementary stability [10, 11, 19] of the system. For the continuous model, the Lyapunov stability theorem should be used for investigating the global stability [21, Theorem 1].

The models considered in [4, 9, 12, 14] have equilibrium points, which are all hyperbolic. Meanwhile the system (1) has a non-hyperbolic equilibrium point in a particular case of the parameters. Moreover, it is globally asymptotical stable. Therefore, for the corresponding discrete system we shall use the Lyapunov stability theorem for proving the stability of this equilibrium point. This is the important contribution of our paper. Besides, we propose a more general model with many iterative parameters in the discretization of the right-hand sides. The combination of appropriate selection of denominator with these parameters will give sufficient conditions for dynamical consistency of the discrete model with the continuous system.

3 Construction of NSFD scheme

In this section we construct NSFD scheme for the system (1) preserving all dynamic properties of the original continuous model for any discretization parameter or step size $h>0$ . Recall that, according to Mickens, a finite difference scheme is called nonstandard if at least one of the following conditions is satisfied [23, 24, 25, 26]:

A nonlocal approximation is used.

2.

The discretization of the derivative is not traditional and uses a function $0<\varphi(h)=h+\mathcal{O}(h^{2})$ .

We propose a general NSFD scheme for the model (1) in the form

[TABLE]

Concerning the set of equilibrium points and positivity of the scheme (3) there hold the following assertions.

Proposition 1

The region $\Omega=\big{\{}(x,y)\in\mathbb{R}^{2}\big{|}x\geq 0,\quad y\geq 0\big{\}}$ is a positive invariant set for the scheme (3) if

[TABLE]

Proof

It is easy to reduce the scheme (3) to the explicit form

[TABLE]

Since the parameters $\alpha_{j}$ and $\beta_{j}$ ( $j=\overline{1,6}$ ) satisfy (4) from the above formulas we obtain the fact to be proved.

Proposition 2

The scheme (3) preserves the set of equilibrium points of system (3).

Proof

The equilibrium $(x^{*},y^{*})$ of the scheme (3) is the solution of the equation

[TABLE]

where $F(x,y)$ and $G(x,y)$ are defined by (5). Since the parameters $\alpha_{i},\beta_{i}$ satisfy (3), if the scheme (3) is defined then it is equivalent to

[TABLE]

It is the system for determining equilibrium points of the model (1). Thus, the proposition is proved.

4 Stability analysis

In this section we give sufficient conditions for the scheme (3) to preserve the stability properties of equilibrium points of the model (1). For the purpose of easy tracking we recall the following results of the necessary and sufficient conditions for an equilibrium point to be locally asymptotically stable [1, Theorem 2.10], [18].

Lemma 1

Assume the functions $f(x,y)$ and $g(x,y)$ have continuous first-order partial derivatives in $x$ and $y$ on some open set in $\mathbb{R}^{2}$ that contains the point $(x^{*},y^{*})$ . Then the equilibrium point $(x^{*},y^{*})$ of the nonlinear system

[TABLE]

is locally asymptotically stable if the eigenvalues of the Jacobian matrix J evaluated at the equilibrium satisfy $|\lambda_{i}|<1$ iff

[TABLE]

The equilibrium is unstable if some $|\lambda_{i}|>1$ , that is, if any one of three inequalities is satisfied,

[TABLE]

4.1 The extinction equilibrium point $P_{0}^{*}=(0,0)$

Proposition 3

For the case $m_{1}>r(0)$ and $m_{2}>s(0)$ consider the difference scheme (3) under the assumptions of Proposition 1. Then the extinction equilibrium point $P_{0}^{*}=(x_{0}^{*},y_{0}^{*})=(0,0)$ is locally asymptotically stable if $m_{1}>r(0)$ and $m_{2}>s(0)$ , and unstable otherwise.

Proof

Computing the Jacobian matrix of system (3), evaluated at the extinction equilibrium point $P_{0}^{*}=(x_{0}^{*},y_{0}^{*})=(0,0)$ , one obtains

[TABLE]

In this case, the eigenvalues are

[TABLE]

Since $\alpha_{j}$ and $\beta_{j}$ satisfy Proposition 1 we have $\lambda_{1}>0,\,\,\lambda_{2}>0$ . On the other hand

[TABLE]

Therefore, $|\lambda_{1}|,|\lambda_{2}|<1$ iff $m_{1}>r(0)$ and $m_{2}>s(0)$ . By Lemma 1 Proposition is proved.

Next, consider scheme (3) under the assumptions of Proposition 1 in the case $m_{1}\geq r(0)$ and $m_{2}\geq s(0)$ . Notice that when $m_{1}=r(0)$ or $m_{2}=s(0)$ the equilibrium point $P_{0}^{*}=(0,0)$ becomes non-hyperbolic one. Therefore, it is impossible to use Lemma 1 for proving the locally asymptotical stability. Moreover, the equilibrium point $P_{0}^{*}$ , as shown in [21, Theorem 1], is globally asymptotically stable. As for the continuous system we shall also use Lyapunov’s stability theorem to prove that the point $P_{0}^{*}=(0,0)$ is a globally asymptotically stable equilibrium point of the scheme (3). For this purpose we consider a family of functions

[TABLE]

where $\alpha,\beta,\gamma,\delta>0$ are parameters, which are selected so that the function $V(x_{k},y_{k})$ satisfies the conditions of the Lyapunov’s stability theorem [15, Theorem 4.20].

Obviously, the function $V(x_{k},y_{k})$ defined by (6) is continuous on $\mathbb{R}_{+}^{2}$ and $V(x_{k},y_{k})\to\infty$ as $||(x_{y},y_{k})||\to\infty$ . Moreover, $V(P_{0}^{*})=0$ and $V(x_{k},y_{k})>0$ for any $(x_{k},y_{k})\in\mathbb{R}_{+}^{2}$ , $(x_{k},y_{k})\neq(0,0)$ . Therefore, in order to show that $V(x_{k},y_{k})$ satisfies [15, Theorem 4.20] it suffices to determine conditions for

[TABLE]

We have

[TABLE]

Now we rewrite (5) in the form

[TABLE]

According to the estimates in Theorem 1 in [21] we have: if $m_{1}\geq r(0)$ and $m_{2}\geq s(0)$ , then $r(x)-m_{1}\leq 0$ for all $x\geq 0$ and $s(y)-m_{2}\leq 0$ for all $y\geq 0$ . Therefore, from (8) we see that if $m_{1}\geq r(0)$ , then $x_{k+1}\leq x_{k}$ for all $k\geq 0$ . Consequently, from (7) it implies that if $y_{k+1}-y_{k}<0$ then $\Delta V(x_{k},y_{k})<0$ . Hence, it is sufficient to consider only the case $y_{k+1}\geq y_{k}$ .

Using the mean value theorem for two variables function $u(x_{k},y_{k})=x_{k}y_{k}$ we have

[TABLE]

where $\xi_{x_{k}}\in(x_{k+1},x_{k})$ and $\xi_{y_{k}}\in(y_{k},y_{k+1})$ . On the other hand, since $x_{k+1}\leq x_{k}$ we have

[TABLE]

From (9), (10), (7) and (8) we obtain the estimate

[TABLE]

where $Q(x_{k},y_{k})$ is defined by

[TABLE]

From here it follows that if $Q(x_{k},y_{k})<0$ for all $(x_{k},y_{k})\in\mathbb{R}_{+}^{2}\backslash\{(0,0)\}$ then $\Delta V(x_{k},y_{k})<0$ for all $(x_{k},y_{k})\in\mathbb{R}_{+}^{2}\backslash\{(0,0)\}$ . On the other hand we have

[TABLE]

where

[TABLE]

Notice that $Q_{1}(x_{k},y_{k})>0$ for all $(x_{k},y_{k})\in\mathbb{R}_{+}^{2}\backslash\{(0,0)\}$ . Therefore, $Q(x_{k},y_{k})<0$ for all $(x_{k},y_{k})\in\mathbb{R}_{+}^{2}\backslash\{(0,0)\}$ if $\tau_{i}<0$ ( $i=\overline{1,6}$ ). Taking into account (2) we have $0<r(x_{k})\leq r(0)$ and $0<\phi(x_{k})\leq\phi(0)$ for all $x_{k}\geq 0$ . Therefore from (12) it follows

$\tau_{1}<0$ if $c<\dfrac{\beta}{\alpha}$ , 2. 2.

$\tau_{2}=-\alpha_{2}c\alpha r(x_{k})+\alpha_{6}c\alpha m_{1}-\beta\beta_{6}m_{2}$ $\leq$ $-\alpha_{2}c\alpha r(0)+\alpha_{6}c\alpha m_{1}-\beta\beta_{6}m_{2}<0$

if $\dfrac{-\alpha_{2}cr(0)+\alpha_{6}cm_{1}}{\beta_{6}m_{2}}<\dfrac{\beta}{\alpha}$ , 3. 3.

$\tau_{3}<0$ if $\alpha_{4}+\beta_{4}<0$ , 4. 4.

$\tau_{4}<0$ if $c<\dfrac{\gamma}{\delta}$ , 5. 5.

$\tau_{5}=-\alpha_{2}c\delta r(x_{k})+c\delta\alpha_{6}m_{1}-\beta_{6}\gamma m_{2}\leq-\alpha_{2}c\delta r(0)+c\delta\alpha_{6}m_{1}-\beta_{6}\gamma m_{2}<0$

if $\dfrac{-\alpha_{2}cr(0)+c\alpha_{6}m_{1}}{\beta_{6}m_{2}}<\dfrac{\gamma}{\delta}$ , 6. 6.

$\tau_{6}=\alpha_{4}c\delta\phi(x_{k})-\alpha\beta_{6}m_{2}\leq\alpha_{4}c\delta\phi(0)-\alpha\beta_{6}m_{2}<0$

if $\dfrac{c\alpha_{4}\phi(0)}{\beta_{6}m_{2}}<\dfrac{\alpha}{\delta}$ .

In summary the function $V(x_{k},y_{k})$ defined by (6) satisfies $\Delta V(x_{k},y_{k})<0$ for all $(x_{k},y_{k})\in\mathbb{R}_{+}^{2}\backslash\{(0,0)\}$ if

[TABLE]

Once, the scheme (3) is fixed the selection of the parameters $\alpha,\beta,\gamma,\delta>0$ satisfying the above relations is completely possible. Thus, we obtain the following theorem of the global stability of the equilibrium point $P_{0}^{*}=(0,0)$ .

Theorem 3

For the case $m_{1}\geq r(0)$ and $m_{2}\geq s(0)$ consider the difference scheme (3) under the assumptions of Proposition 1. If additionally assume that

[TABLE]

then the extinction equilibrium point $P_{0}^{*}=(0,0)$ is globally asymptotically stable.

4.2 The equilibrium point $P_{1}^{*}=(K,0)$ (the equilibrium point of extinction of the predator species)

Proposition 4

For the case $m_{1}<r(0)$ and $m_{2}>s(0)+cK\phi(K)$ consider the difference scheme (3) under the assumptions of Proposition 1. If additionally assume that

[TABLE]

then the equilibrium point of the form $P_{1}^{*}=(K,0)$ is locally asymptotically stable if $m_{1}<r(0)$ and $m_{2}>s(0)+cK\phi(K)$ , and unstable otherwise.

Proof

Recall that $m_{1}<r(0)$ is necessary and sufficient condition for the existence of the equilibrium point $P_{1}^{*}$ (see Theorem 1). Computing the Jacobian matrix of system (3), evaluated at the extinction equilibrium point $P_{1}^{*}=(K,0)$ , one obtains

[TABLE]

In this case, the eigenvalues are

[TABLE]

Therefore, $\lambda_{2}<1$ iff $m_{2}>s(0)+cK\phi(K)$ . Besides, due to $r^{\prime}(K)<0$ there holds $\lambda_{1}<1$ . On the other hand we have

[TABLE]

Therefore, if (15) is fulfilled then $\lambda_{1}>-1$ , $\lambda_{2}>-1$ . In result we have $|\lambda_{1}|<1$ and $|\lambda_{2}|<1$ . Hence, by Lemma 1 the proposition is proved.

4.3 The equilibrium point $P_{2}^{*}=(0,M)$ (the equilibrium point of extinction of the prey species)

Proposition 5

For the case $m_{1}>r(0)-M\phi(0)$ and $m_{2}<s(0)$ consider the difference scheme (3) under the assumptions of Proposition 1. If additionally assume that

[TABLE]

then the equilibrium point of the form $P_{2}^{*}=(0,M)$ is locally asymptotically stable if $m_{1}>r(0)-M\phi(0)$ and $m_{2}<s(0)$ , and unstable otherwise.

Proof

Notice that according to Theorem 1, $m_{2}<s(0)$ is necessary and sufficient condition for the existence of the equilibrium point $P_{2}^{*}$ . Computing the Jacobian matrix of system (3), evaluated at the extinction equilibrium point $P_{2}^{*}=(0,M)$ , one obtains

[TABLE]

In this case, the eigenvalues are

[TABLE]

Therefore, $\lambda_{1}<1$ iff $m_{1}>r(0)-M\phi(0)$ . Besides, due to $s^{\prime}(M)<0$ there holds $\lambda_{2}<1$ . On the other hand we have

[TABLE]

Hence, if (16) is fulfilled then $\lambda_{2}>-1$ and $\lambda_{1}>-1$ . In result we have $|\lambda_{1}|<1$ and $|\lambda_{2}|<1$ . By Lemma 1 the proposition is proved.

4.4 The equilibrium point $P_{3}^{}=(x^{},y^{*})$ (the ecological stability equilibrium)

Proposition 6

Consider the difference scheme (3) under the assumptions of Proposition 1. If additionally assume that

[TABLE]

then, if the equilibrium point of the form $P_{3}^{*}=(x^{*},y^{*})$ belong to $\Omega$ , then it is locally asymptotically stable.

Proof

Notice that the conditions for the existence of the equilibrium point $P_{3}^{*}$ are given in Theorem 1. Computing the Jacobian matrix of system (3), evaluated at the extinction equilibrium point $P_{3}^{*}=(x^{*},y^{*})$ , one obtains

[TABLE]

In this case, we have

[TABLE]

where

[TABLE]

According to Theorem 2 in [21] we have $s^{\prime}(y^{*})<0,\,r^{\prime}(x^{*})-y^{*}\phi^{\prime}(x^{*})<0,$ consequently

[TABLE]

On the other hand $\det(J(P_{3}^{*}))<1$ iff

[TABLE]

This equivalent to

[TABLE]

The first term of the left-hand side of the above inequality is positive, therefore, if $T_{5}>0$ then the above inequality is valid. It implies $\det(J(P_{3}^{*}))<1$ .

Finally, we see that $1+Tr(J(P_{3}^{*}))+\det(J(P_{3}^{*}))>0$ iff

[TABLE]

The third term of the above sum is always positive, so if

[TABLE]

then $1+Tr(J(P_{3}^{*}))+\det(J(P_{3}^{*}))>0$ . It is easy to verify that if $T_{6},T_{7}>0$ then (18) holds. Thus, we have proved that if (17) satisfies then the three conditions of Lemma 1 are satisfied. Therefore, the point $P_{3}^{*}$ is locally asymptotically stable. The proposition is proved.

Remark 1

The system of conditions (14)-(17) for $\alpha_{j}$ , $\beta_{j}$ $(j=\overline{1,6})$ has many solutions. Simply, we can choose $\alpha_{4},\alpha_{6},\beta_{6}>0$ sufficiently large, $\alpha_{2},\beta_{2},\beta_{4}<0$ sufficiently small and $\alpha_{4}+\beta_{4}<0$ .

Now we summarize the obtained results above in the following theorem on the NSFD schemes preserving the dynamical properties of the model (1).

Theorem 4

The NSFD scheme (3) is dynamically consistent with (1) if the parameters $\alpha_{j}$ , $\beta_{j}$ $(j=\overline{1,6})$ satisfy the conditions listed in Table 1, where the columns list sufficient conditions for the scheme (3) preserve corresponding properties of the model (1) for different cases of the parameters. The symbol ${}^{\prime\prime}*^{\prime\prime}$ means that the set of equilibrium points of (1) is always preserved by the scheme (3).

Remark 2

There are infinitely many ways for selecting the parameters $\alpha_{j}$ , $\beta_{j}$ $(j=\overline{1,6})$ satisfying the conditions listed in Table 1 ( see Remark 1). This shows the existence of NSFD schemes dynamically consistent with the system (1).

5 Numerical Simulations

Numerical examples presented in this section show that the obtained theoretical results of the NSFD preserving the properties of the general predator-prey system are valid.

Let us consider the example of the predator-prey model presented in [21]. In this example

[TABLE]

for $6$ cases of the parameters $(m_{1},m_{2})$ in Corollary 1 in [21]. Namely,

[TABLE]

Many numerical simulations, such as [2, 3, 8, 10, 11, 12, 13, 31, 32, 35] …agree that standard difference schemes do not preserve dynamical properties of continuous models for large step sizes, i.e., are not dynamically consistent with continuous systems. This confirms the advantages of NSFD schemes. For supporting this confirmation in this example we use the explicit Euler scheme and four stage Runge-Kutta (RK4) scheme compared with the constructed NSFD schemes for the system (1). The numerical solutions obtained for these schemes in Case (i) of $m_{1},\ m_{2}$ are depicted in Figures 1-4. From the figures we see that the property of positivity and stability of the system are destroyed. The numerical experiments in other cases of the parameters are analogous.

The numerical solutions obtained by these NSFD schemes are depicted in Figures 5-10, respectively. Comparing these results with the numerical simulations in Section 4 in [21] we see that all properties of the continuous model are preserved.

6 Conclusion

In this paper we have used NSFD schemes for converting a general predator-prey model to a dynamically consistent discrete system. It should be emphasized that in the continuous as in the discrete system there is a non-hyperbolic equilibrium point, whose global asymptotical stability was proved by means of the Lyapunov stability theorem. The numerical simulations for the model considered in [21] for various collections of parameters confirm the validity of the obtained theoretical results. In the future we shall develop the techniques used in this paper for constructing and investigating NSFD schemes for other dynamical models including ones having non-hyperbolic equilibrium points.

Acknowledgments

This work is supported by Vietnam National Foundation for Science and Technology Development (NAFOSTED) under the grant number 102.01-2014.20.

References

Bibliography35

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1[1] L. J. S. Allen, An Introduction to Mathematical Biology, Prentice Hall, New Jersey, 2007.
2[2] R. Anguelov, J. M. -S Lubuma, Nonstandard finite difference method by nonlocal approximations, Mathematics and Computers in Simulation, 61 (2003) 465-475.
3[3] R. Anguelov, Y. Dumont, J.M.-S. Lubuma, M. Shillor, Dynamically consistent nonstandard finite difference schemes for epidemiological models, Journal of Computational and Applied Mathematics, 255 (2014) 161-182.
4[4] N. Bairagi, M. Biswas, A predator-prey model with Beddington-De Angelis functional response: a non-standard finitedifference method, Journal of Difference Equations and Applications, 22 (4) (2016) 581-593.
5[5] A. A. Berryman, The origins and evolution of predator-prey theory, Ecology, 73 (1992) 1530-1535.
6[6] F. Brauer, C. Castillo - Chavez, Mathematical Models in Population Biology and Epidemiology, Springer New York, 2001.
7[7] J. Cresson, F. Pierret, Non standard finite difference scheme preserving dynamical properties, Journal of Computational and Applied Mathematics, 303 (2016) 15-30.
8[8] Q. A. Dang, M. T. Hoang, Dynamically consistent discrete metapopulation model, Journal of Difference Equations and Applications, (2016) DOI: 10.1080/10236198.2016.1197213.

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Videos

Nonstandard finite difference schemes

Abstract

keywords:

1 Introduction

2 Mathematical model

Theorem 1

Theorem 2

3 Construction of NSFD scheme

Proposition 1

Proposition 2

4 Stability analysis

Lemma 1

4.1 The extinction equilibrium point P0∗=(0,0)P_{0}^{*}=(0,0)P0∗​=(0,0)

Proposition 3

Theorem 3

4.2 The equilibrium point P1∗=(K,0)P_{1}^{*}=(K,0)P1∗​=(K,0) (the equilibrium point of extinction of the predator species)

Proposition 4

4.3 The equilibrium point P2∗=(0,M)P_{2}^{*}=(0,M)P2∗​=(0,M) (the equilibrium point of extinction of the prey species)

Proposition 5

4.4 The equilibrium point P3∗=(x∗,y∗)P_{3}^{*}=(x^{*},y^{*})P3∗​=(x∗,y∗) (the ecological stability equilibrium)

Proposition 6

Remark 1

Theorem 4

Remark 2

5 Numerical Simulations

6 Conclusion

Acknowledgments

4.1 The extinction equilibrium point $P_{0}^{*}=(0,0)$

4.2 The equilibrium point $P_{1}^{*}=(K,0)$ (the equilibrium point of extinction of the predator species)

4.3 The equilibrium point $P_{2}^{*}=(0,M)$ (the equilibrium point of extinction of the prey species)

4.4 The equilibrium point $P_{3}^{}=(x^{},y^{*})$ (the ecological stability equilibrium)