\documentclass[reqno]{amsart}
\usepackage{graphicx}
\usepackage{hyperref}
\AtBeginDocument{{\noindent\small
Sixth Mississippi State Conference on Differential Equations and
Computational Simulations,
{\em Electronic Journal of Differential Equations},
Conference 15 (2007), pp. 67--75.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu (login: ftp)}
\thanks{\copyright 2007 Texas State University - San Marcos.}
\vspace{9mm}}
\begin{document} \setcounter{page}{67}
\title[\hfilneg EJDE-2006/Conf/15\hfil Numerical methods for predator-prey models]
{Nonstandard numerical methods for a class of
predator-prey models with predator interference}
\author[D. T. Dimitrov and H. V. Kojouharov\hfil EJDE/Conf/15 \hfilneg]
{Dobromir T. Dimitrov, Hristo V. Kojouharov} % in alphabetical order
\address{Dobromir T. Dimitrov \hfill\break
Department of Ecology and Evolutionary Biology,
University of Tennessee at Knoxville, Knoxville, TN 37996-1610, USA}
\email{ddimitr1@utk.edu}
\address{Hristo V. Kojouharov \hfill\break
Department of Mathematics, University of Texas at Arlington,
Arlington, TX 76019-0408, USA}
\email{hristo@uta.edu}
\thanks{Published February 28, 2007.}
\subjclass[2000]{37M05, 39A11, 65L12, 65L20}
\keywords{Finite-difference; nonstandard; positive; elementary stable;\hfill\break\indent
predator-prey; predator interference}
\begin{abstract}
We analyze a class of predator-prey models with
Beddington-DeAngelis type functional response. The models
incorporate the mutual interference between predators, which
stabilizes predator-prey interactions even when only a linear
intrinsic growth rate of the prey population is considered.
Positive and elementary stable nonstandard (PESN)
finite-difference methods, having the same qualitative features as
the corresponding continuous predator-prey models, are formulated
and analyzed. The proposed numerical techniques are based on a
nonlocal modelling of the growth-rate function and a nonstandard
discretization of the time derivative. This approach leads to
significant qualitative improvements in the behavior of the
numerical solution. Applications of the PESN methods to specific
Beddington-DeAngelis predator-prey systems are also presented.
\end{abstract}
\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{definition}[theorem]{Definition}
\section{Introduction}
The traditional mathematical model of predator-prey interactions consists of
the following system of two differential equations:
\begin{equation} \label{eq:sysp}
\begin{gathered}
\frac{dx}{dt} = P(x)-F(x,y); \quad x(0)=x_0 \ge 0,\\
\frac{dy}{dt} = eF(x,y)-D(y); \quad y(0)=y_0 \ge 0,
\end{gathered}
\end{equation}
where $x$ and $y$ represent the prey and predator population
sizes, respectively, and functions $P(x)$, $D(y)$ describe the
intrinsic growth rate of the prey and the mortality rate of the
predator, respectively. The function $D(y)$ is assumed to be
linear ($D(y)=dy$), while the function $P(x)$ has a linear
($P(x)=ax$) or logistic ($P(x)=ax(1-\frac{x}{K})$) expression. The
function $F(x,y)$ is called ``functional response'' or ``feeding
rate'' and represents the prey consumption per unit time. The
model assumes a linear correspondence between the prey consumption
and the predator production through the positive constant $e$.
Among the most popular functional responses used in the modelling
of predator-prey systems are the Michaelis-Menten type $
F(x,y)=\frac{fxy}{x+c}$ and the ratio-dependent type $
F(x,y)=\frac{fxy}{x+by}$. However, in some situations they predict
unrealistic population dynamics of the predator or the prey. The
Michaelis-Menten type does not account for the mutual competitions
among predators \cite{may}, while the ratio-dependent type allows
unrealistic positive growth rate of the predator at low densities
\cite{arditi, kuang, hsu}. The Beddington-DeAngelis functional
response $F(x,y)=\frac{fxy}{b+wx+y}$ was introduced independently
by Beddington \cite{bedd} and DeAngelis \cite{deang} as a solution
of the observed problems in the classical predator-prey theory. It
has an extra term in the denominator which models mutual
interference between predators and avoids the ``low densities
problem'' of the ratio-dependent type functional response.
The Beddington-DeAngelis predator-prey model with a linear intrinsic
growth of the prey population, analyzed completely in \cite{pap1},
has the following nondimensional form:
\begin{equation} \label{eq:sys_bd}
\begin{gathered}
\frac{dx}{dt} = x-\frac{axy}{1+x+y}; \quad x(0)=x_0 \ge 0,\\
\frac{dy}{dt} = \frac{exy}{1+x+y}-dy; \quad y(0)=y_0 \ge 0,
\end{gathered}
\end{equation}
where $x$ and $y$ represent the prey and predator population
sizes, respectively, and the positive constants $a$, $e$ and $d$
represent the generalized feeding rate, generalized conversion
efficiency and generalized mortality rate of the predator,
respectively.
Numerical simulations of System (\ref{eq:sys_bd}) can be obtained
by methods based on finite difference approximations, such as
Euler, Runge-Kutta, or Adams methods. However, their use raises
questions about the truncation errors, the stability regions and,
from a dynamical point of view, the accuracy at which the dynamics
of the continuous system are represented by the discrete system.
The nonstandard finite difference schemes, developed by Mickens
\cite{mickens1}, have laid the foundation for designing methods
that preserve the dynamical behavior of the approximated
differential system. The techniques are based on a nonlocal
numerical treatment of the right-hand side functions and more
sophisticated discretizations of time derivatives. Nonstandard
numerical techniques, preserving the stability of equilibria, have
been applied to some specific dynamical systems by Lubuma and Roux
\cite {roux} and Dimitrov and Kojouharov \cite{pap3,pap5}, among
others. The resulting so-called elementary stable nonstandard
(ESN) methods solve the numerical stability problems locally for
an arbitrary time step-sizes. However, the ESN methods, as well as
the standard numerical methods, do not guarantee a positive
discrete solution for all positive initial values. The positivity
condition is natural when interspecies interactions are modeled
and approximated numerically. Failing to satisfy it reflects
negatively on the accuracy and efficiency of the numerical
methods. Several classes of positive and elementary stable
nonstandard (PESN) methods have been recently developed for some
phytoplankton-nutrient \cite{pap2} models and for a
Rosenzweig-MacArthur predator-prey \cite{pap6} model, in which a
constant feeding rate per predator has been assumed.
The predator-dependent Beddington-DeAngelis
functional response $F(x,y)$ leads to richer dynamics
of System (\ref{eq:sys_bd}), which makes the design and
analysis of the PESN numerical methods more challenging
that in earlier cases.
The PESN methods, developed here, preserve both the
positivity of the solutions and the stability of the equilibria
of System (\ref{eq:sys_bd}). In addition,
the designed numerical schemes allow for the discrete systems to
be solved explicitly, which increases the efficiency of the methods.
The paper is organized as follows.
Section \ref{sec:pr} provides some definitions and preliminary results.
Results from the mathematical analysis of the Beddington-DeAngelis predator-prey system followed by the development of the corresponding
PESN numerical methods are presented in Section \ref{sec:mr}.
In the last two sections we illustrate our theoretical results by a set of numerical examples and outline some future research directions.
\section{Definitions and Preliminaries}
\label{sec:pr}
A general $n$-dimensional autonomous system has the following form:
\begin{equation}
\label{eq:sys}
\frac{d \bar x}{dt} = f(\bar x);\hspace{6pt} \bar x(t_0)=\bar x_0,
\end{equation}
where $\bar x=[x^1,x^2,\dots,x^n]^T:[t_0,T) \rightarrow {\mathbb
R}^n$, the function $f=[f^1,f^2,\dots,f^n]^T:{\mathbb R}^n
\mapsto {\mathbb R}^n$ is differentiable and $\bar x_0 \in
{\mathbb R}^n$. The equilibrium points of System (\ref{eq:sys})
are defined as the solutions of $f(\bar x) = 0$.
\begin{definition} \label{def2.1} \rm
Let $\bar x^*$ be an equilibrium of System (\ref{eq:sys}), $J(\bar
x^*)$ be the Jacobian of System (\ref{eq:sys}) at $\bar x^*$ and
$\sigma(J(\bar x^*))$ denotes the spectrum of $J(\bar x^*)$. An
equilibrium $\bar x^*$ of System (\ref{eq:sys}) is called linearly
stable if $Re(\lambda)<0$ for all $\lambda \in \sigma(J(\bar
x^*))$ and linearly unstable if $Re(\lambda)>0$ for at least one
$\lambda \in \sigma(J(\bar x^*))$.
\end{definition}
A numerical scheme with a step size $h$, that approximates the solution
$\bar x(t_k)$ of System (\ref{eq:sys}) can be written in the form
\begin{equation}
\label{eq:num}
D_h(\bar x_k) = F_h(f;\bar x_k),
\end{equation}
where $D_h(\bar x_k) \approx \frac{d\bar
x}{dt}\big|_{t=t_k}$, $\bar x_k \approx \bar x(t_k)$,
$F_h(f;\bar x_k)$ approximates the right-hand side of System
(\ref{eq:sys}) and $t_k=t_0+kh$. Throughout this article, we
assume that System (\ref{eq:sys}) has a finite number of
hyperbolic equilibria, i.e., $Re(\lambda) \neq 0$ for all $\lambda
\in \Omega$, where $\Omega = {\cup}_{\bar x^* \in \Gamma}
\sigma(J(\bar x^*))$ and $\Gamma$ represents the set of all
equilibria of System (\ref{eq:sys}).
\begin{definition} \label{def2.2} \rm
Let $\bar x^*$ be a fixed point of the scheme \eqref{eq:num} and the equation of the perturbed solution $\bar x_k= \bar x^*+\bar \epsilon_k$, where $\bar \epsilon_k$ is small, be linearly approximated by
\begin{equation}
\label{eq:pert}
D_h \bar \epsilon_k = J_h \bar \epsilon_k.
\end{equation}
Here the right-hand side of Equation (\ref{eq:pert}) represents the linear term
in $\bar \epsilon_k$ of the Taylor expansion of $F_h(f;\bar x^* +\bar \epsilon_k)$ around $\bar x^*$.
The fixed point $\bar x^*$ is called stable if $\|\bar \epsilon_k \| \rightarrow 0$ as
$k \rightarrow \infty$, and unstable otherwise, where $\bar \epsilon_k$ is the solution of Equation (\ref{eq:pert}).
\end{definition}
Let System \eqref{eq:num} be expressed in the following explicit way:
\begin{equation}\label{eq:diff}
\bar x_{k+1}=G(\bar x_k),
\end{equation}
where the function $G=[G^1,G^2,\dots,G^n]^T:{\mathbb R}^n \mapsto {\mathbb R}^n$ is differentiable.
If $\bar x^*$ is a fixed point of System (\ref{eq:diff}) then the equation for the perturbed solution $\bar \epsilon_k$, around $\bar x^*$, has the form:
\[
\bar \epsilon_{k+1} = J(\bar x^*) \bar \epsilon_k,
\]
where $J(\bar x^*)$ denotes the Jacobian
$\big( \frac {\partial
G^i(\bar x^*)}{\partial x^j} \big)_{1 \le i,j \le n}$. The
solution $\|\bar \epsilon_k\| \to 0$ when $k \to \infty$ if and
only if all eigenvalues of $J(\bar x^*)$ are less than one in
absolute values.
We introduce the next two definitions
based on definitions given by Anguelov and Lubuma in \cite{anguelov}.
\begin{definition} \label{def:el_st} \rm
The finite difference method \eqref{eq:num} is called elementary stable, if, for any value of the step size h, its only fixed points $\bar x^*$ are those of the differential system (\ref{eq:sys}), the linear stability properties of each $\bar x^*$ being the same for both the differential system and the discrete method.
\end{definition}
\begin{definition} \label{def:nonstandard} \rm
The one-step method \eqref{eq:num} is called a nonstandard
finite-difference method if at least one of the following conditions is
satisfied:
\begin{itemize}
\item $D_h(\bar x_k) = \frac{\bar x_{k+1}-\bar x_k}{\varphi(h)}$, where
$\varphi(h)=h+\mathcal{O}(h^2)$ is a nonnegative function;
\item $F_h(f;\bar x_k) = g(\bar x_k,\bar x_{k+1},h)$, where
$g(\bar x_k,\bar x_{k+1},h)$ is a nonlocal approximation of
the right-hand side of System (\ref{eq:sys}).
\end{itemize}
\end{definition}
The form of predator-prey systems (\ref{eq:sys_bd}) guarantees that
any solution with positive initial conditions remains positive in time.
The corresponding requirement for numerical methods is formulated
in the following definition.
\begin{definition} \label{def:pos} \rm
The finite difference method \eqref{eq:num} is called positive,
if, for any value of the step size $h>0$, and
$\bar{x}_0 \in {\mathbb R}^n_+$ its solution remains positive,
i.e., $\bar{x}_k \in {\mathbb R}^n_+$ for $k=1,2,3,\dots$
\end{definition}
The following theorem deals with the properties of
the general nonstandard schemes based on standard finite-difference methods \cite{anguelov}:
\begin{theorem} \label{th:consist}
If the numerical method \eqref{eq:num} represents a standard
finite-difference scheme that is consistent and zero-stable,
then any corresponding nonstandard finite-difference scheme
in Definition \ref{def:nonstandard} is necessarily consistent.
Furthermore, if the nonstandard scheme is constructed according
to the first bullet of Definition \ref{def:nonstandard},
then this scheme is zero-stable provided that the operator $F_h$ satisfies,
for some $M>0$ independent of $h$ and for any bounded sequences
\{$\bar x_k$\} and \{$\tilde x_k$\} in ${\mathbb R}^n$, the
Lipschitz condition
\[
\sup_k \|F_h(\bar x_k)-F_h(\tilde x_k)\| \le M \sup_k \|\bar x_k-\tilde x_k\|.
\]
\end{theorem}
\section{PESN Methods for Predator-Prey Models with
Beddington-DeAngelis Functional Response}
\label{sec:mr}
Depending on the values of the parameters in System (\ref{eq:sys_bd})
the following equilibria exist:
\begin{enumerate}
\item $E_0=(0,0)$;
\item $E^*=(x^*,y^*)=\big(\frac{ad}{ae-e-ad},\frac{e}{ae-e-ad}\big)$.
The equilibrium $E^*$ exists if and only if $ae-e-ad > 0$.
\end{enumerate}
According to the mathematical analysis of System (\ref{eq:sys_bd}) in \cite{pap1}, the following statements about the stability of the equilibria are true:
\begin {enumerate}
\item The equilibrium $E_0$ is always linearly unstable;
\item The equilibrium $E^*$ is linearly stable if $a < e$ and linearly unstable
if $a > e$;
\end {enumerate}
The dynamics of System (\ref{eq:sys_bd}) are summarized
in the following theorem \cite{pap1}.
\begin{theorem} \label{th:old}
For the global dynamics of System (\ref{eq:sys_bd}) the following holds:
\begin{enumerate}
\item If $ae-e-ad \leq 0$ then all trajectories are unbounded.
In addition, $e>d$ yields that $ \lim _{t \to \infty} x(t) =
\infty$ and $\lim _{t \to \infty} y(t) = \infty$ while $e\le d$
yields that $\lim _{t \to \infty} x(t)=\infty$ and $\lim _{t \to
\infty} y(t)=0$.
\item If $ae-e-ad>0$ and $a0$ and $a>e$ then the interior equilibrium $E^*$
is unstable. If, in addition, $e \geq d+1$ then every trajectory
circles indefinitely and counterclockwise around $E^*$ moving away
from it. If $e0$ and $a=e$ then the interior equilibrium $E^*$
is a global center, which means that all trajectories
(except $E^*$ itself) are periodic orbits containing $E^*$
in their interior.
\end{enumerate}
\end{theorem}
The new PESN methods for solving Beddington-DeAngelis predator-prey systems
with linear intrinsic growth of the prey population can be designed,
based on the following theorem:
\begin{theorem} \label{th:pesn_bd}
Let $\phi$ be a real-valued function on $\mathbb R$ that satisfies
the property:
\begin{equation} \label{eq:phi}
\phi(h)=h+O(h^2) \quad \text{and} \quad 0<\phi(h)<1 \quad
\text{for all } h>0.
\end{equation}
Then the following scheme for solving System (\ref{eq:sys_bd})
represents a PESN method:
\begin{equation} \label{eq:stab_bd}
\begin{gathered}
\frac{x_{k+1}-x_k}{\varphi(h)}=x_k-\frac{ax_{k+1}y_k}{1+x_k+y_k};\\
\frac{y_{k+1}-y_k}{\varphi(h)}=\frac{ex_ky_k}{1+x_k+y_k}-dy_{k+1};
\end{gathered}
\end{equation}
where $\varphi(h)$ belongs to the following class of functions:
\begin{enumerate}
\item If $ae-e-ad \le 0$ then $\varphi(h)=\phi(hq)/q$ for all
$q \ge 0$, with $\varphi(h)=h$ for $q=0$.
\item If $ae-e-ad > 0$ and $a \ne e$ then $\varphi(h)=\phi(hq)/q$ for
$q>\frac{e|a-2|}{|e-a|}$.
\end{enumerate}
\end{theorem}
\begin{proof}
The explicit expression of the nonstandard scheme (\ref{eq:stab_bd})
has the form:
\begin{equation}
\label{eq:stab_bd2}
\begin{gathered}
x_{k+1}=\frac{(1+\varphi(h))x_k}{1+\frac{a\varphi(h)y_k}{1+x_k+y_k}}\\
y_{k+1}=\frac{\left(
1+\frac{e\varphi(h)x_k}{1+x_k+y_k}\right)y_k}{1+\varphi(h)d}.
\end{gathered}
\end{equation}
Since the constants $a$, $e$ and $d$ are positive then the
scheme (\ref{eq:stab_bd2}) is unconditionally positive and its fixed
points are exactly the equilibria $E_0$ and $E^*$ of System (\ref{eq:sys_bd}).
The analysis of the Jacobian $J(E_0)$ of Scheme \ref{eq:stab_bd2}
at the equilibrium $E_0$ shows that $E_0$ is always an unstable
fixed point of System (\ref{eq:stab_bd}).
The Jacobian of the Scheme (\ref{eq:stab_bd2}) at the equilibrium
$E^*$ has the following form:
\[
J(E^*)=\begin{pmatrix}
1+\frac{\varphi(h)d}{(1+\varphi(h))e}
& -\frac{\varphi(h)d(a-1)}{(\varphi(h)+1)e} \\
\frac{\varphi(h)(e-d)}{a(1+\varphi(h)d)}
& 1-\frac{\varphi(h)d}{a(1+\varphi(h)d)}
\end{pmatrix}
\]
The eigenvalues $\lambda_1$ and $\lambda_2$ of the Jacobian
$J(E^*)$
are roots of the quadratic equation:
\[\lambda^2-(A+2)\lambda+(A+B+1)=0,\]
where
\[
A=\frac{\varphi(h)d}{(1+\varphi(h))e}-\frac{\varphi(h)d}{a(1+\varphi(h)d)},\quad
B=\frac{\varphi(h)^2d(ae-ad-e)}{(1+\varphi(h))(1+\varphi(h)d)ae}.
\]
For a quadratic equation of the form $\lambda^2+\alpha \lambda +\beta=0$
both roots satisfy $|\lambda_i|<1, i=1,2$ if and
only if the following three conditions are met:
(a) $1+ \alpha + \beta >0$; (b) $1- \alpha + \beta >0$; and
(c) $\beta<1$ \cite[p.82]{chavez}.
Using this fact, we obtain that $E^*$ is a stable fixed point
of Scheme (\ref{eq:stab_bd}) if and only if the following is true:
\begin{enumerate}
\item[(a)] $B >0$;
\item[(b)] $4+2A+B>0$; and
\item[(c)] $A+B<0$;
\end{enumerate}
and that $E^*$ is an unstable fixed point if at least one of the above conditions fails.
The equilibrium $E^*$ exists when
$ae-ad-e>0$, which implies that $a>1$ and $e>d$.
Therefore $B>0$ and the condition (a) is always true.
Simple calculations yield that $A>-1$, which implies that the condition (b) is also
true.
The condition (c) is equivalent to the following inequality:
\begin{equation}
\label{eq:cond3}
\varphi(h)e(a-2)0$ and Inequality
(\ref{eq:cond3}) is true for $a \le 2$. Let $a>2$ and
$\varphi(h)=\phi(hq)/q$, where $q>\frac{e|a-2|}{|e-a|}$. Since
$0<\varphi(h)<\frac{1}{q}$ then $ \varphi(h)<\frac{e-a}{e(a-2)}$
for all $h>0$. Therefore the condition (c) is satisfied and the
interior equilibrium $E^*$ is a stable fixed point of the Scheme
(\ref{eq:stab_bd}) with $\varphi(h)=\phi(hq)/q$.
If the equilibrium $E^*$ is unstable then $e-a<0$ and Inequality
(\ref{eq:cond3}) is not true for $a \ge 2$. This implies that the
condition (c) fails and $E^*$ is an unstable fixed point of Scheme
(\ref{eq:stab_bd}). If $a<2$ and $\varphi(h)=\phi(hq)/q$, where
$q>\frac{e|a-2|}{|e-a|}$, then the condition (c) is not satisfied
and the interior equilibrium $E^*$ is an unstable fixed point of
the Scheme (\ref{eq:stab_bd}) with $\varphi(h)=\phi(hq)/q$.
\end{proof}
The designed nonstandard methods (\ref{eq:stab_bd})
require that the dynamical system (\ref{eq:sys_bd}) has only
hyperbolic equilibria and also that the parameter $q$ is selected
above a specific threshold value.
In that case, the PESN methods replicate the properties $(1)$,
$(2)$, and $(3)$ of System (\ref{eq:sys_bd})
stated in Theorem \ref{th:old}, i.e.,
the methods guarantee that the stabilities of the equilibria $E_0$
and $E^*$ of the differential system are the same as the stabilities
of $E_0$ and $E^*$ as fixed points of the numerical method for
an arbitrary step-size $h$.
\section{Numerical Simulations}
\begin{figure}
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=0.46\textwidth]{figures/fig1a}
& \includegraphics[width=0.46\textwidth]{figures/fig1b} \\
{\footnotesize (a) $h=0.1$, $x_0=5.0$, $y_0=3.5$}
&{\footnotesize (b) $h=1.0$, $x_0=5.0$, $y_0=3.5$}
\\
\includegraphics[width=0.46\textwidth]{figures/fig1c}
& \includegraphics[width=0.46\textwidth]{figures/fig1d} \\
{\footnotesize (c) $h=1.0$, $x_0=0.5$, $y_0=6.5$}
&{\footnotesize (d) $h=0.4$, $x_0=7.5$, $y_0=5.0$}
\\
\includegraphics[width=0.46\textwidth]{figures/fig1e}
& \includegraphics[width=0.46\textwidth]{figures/fig1f} \\
{\footnotesize (e) $h=1.27$, $x_0=2.0$, $y_0=3.5$}
&{\footnotesize (f) $h=2.56$, $x_0=2.0$, $y_0=3.5$}
\end{tabular}
\end{center}
\caption{Numerical approximations of solutions to
System \eqref{eq:sys_bd}.}
\label{fig1}
\end{figure}
To illustrate the performance of the designed new PESN finite-difference methods, we consider the
Beddington-DeAngelis system (\ref{eq:sys_bd}) with $a=3.0$, $d=2.25$ and $e=4.0$.
Mathematical analysis of the predator-prey system shows that
the equilibrium $E_0=(0,0)$ is unstable, while the equilibrium $E^*=(\frac{27}{5},\frac{16}{5})$ is globally asymptotically stable in the interior of the first quadrant \cite{pap1}.
The first set of simulations (Fig.~\ref{fig1} (a)-(c)) compares the performance of the PESN method (\ref{eq:stab_bd}), the explicit Euler method, and the ESN Euler method with $\theta=0$ \cite{pap5}.
For the above choice of system parameters the ESN method is elementary stable for $q_{_{ESN}}>1.25$
\cite[Theorem 1]{pap5}
and the PESN method is elementary stable for $q_{_{PESN}}>4$
(see Theorem \ref{th:pesn_bd}),
so we use $q_{_{ESN}}=1.3$ and $q_{_{PESN}}=4.5$, respectively, in our set of numerical experiments.
Numerical approximations of the solution of System (\ref{eq:sys_bd})
with initial conditions $x_0=5.0$ and $y_0=3.5$ show that for a relatively small step-size (Fig.~\ref{fig1} (a)) all three methods follow the correct dynamics of the original system, but an increase of the step-size
(Fig.~\ref{fig1} (b)) causes the solution of the
explicit Euler method, which is not elementary stable, to evolve away from
the stable equilibrium $E^*$.
Moving the initial condition closer to the $x$-axis,
e.g. $x_0=0.5$ and $y_0=6.5$, causes the numerical solution of the
non positivity-preserving ESN Euler method,
after an initial jump of the solution's
trajectory outside of the first quadrant (Fig.~\ref{fig1} (c)),
to evolve away from the equilibrium $E^*$.
The next simulation compares the approximations of the solution, originating at $(x_0,y_0)=(7.5,5.0)$, obtained by the PESN method (\ref{eq:stab_bd})
and the Patankar Euler method \cite{burchard}.
Even though the Patankar Euler method is positive,
its numerical solution leads to a machine blowup for
a variety of step-sizes (e.g. $h=0.4$ in Fig.~\ref{fig1} (d)),
because the method is not elementary stable.
The last two simulations compare the PESN method with the second- and fourth-order Runge-Kutta methods.
In the selected numerical examples the two Runge-Kutta methods, which are not elementary stable, fail to preserve the correct
asymptotic behavior of System (\ref{eq:sys_bd}) for a variety of step-sizes
(e.g. $h=1.27$ and $h=2.56$ in Fig.~\ref{fig1} (e) and (f), respectively).
\section*{Conclusions}
The designed new PESN methods preserve two of the most important dynamical characteristics of predator-prey system with Beddington-DeAngelis functional response (\ref{eq:sys_bd}), namely the stability of all equilibria and the positivity of all solutions with positive initial conditions.
In all of the numerical experiments (Fig.~\ref{fig1}) the PESN solution
converges toward the asymptotically stable interior equilibrium for
arbitrary time-steps, while all of the other numerical methods fail to
preserve either the stability of the equilibrium or the positivity
of the trajectories for a variety of step-sizes.
The solution of the PESN method follows the predicted by the mathematical analysis
(Theorem \ref{th:old}) asymptotic behavior.
It also preserves the local stability of the equilibria and the positivity of the solutions of System (\ref{eq:sys_bd}) for any step-size.
The explicit form of
the PESN methods makes them a computationally effective tool in the numerical simulations of predator-prey systems.
Future research directions include the design and analysis of PESN methods for the general predator-prey system (\ref{eq:sysp}) and
the construction of similar nonstandard schemes
for biological systems with non-hyperbolic equilibria.
\begin{thebibliography}{00}
\bibitem{anguelov} Anguelov R., Lubuma J. M.-S.,
\emph{Contributions to the mathematics of the nonstandard finite
difference method and applications},
Numer. Methods Partial Diff. Equations \textbf{17:5}, (2001) 518-543.
\bibitem{bedd} Beddington J.R.,
\emph{Mutual interference between parasites or predators and its
effect on searching efficiency}, J. Animal Ecol. \textbf{44}, (1975) 331-340.
\bibitem{arditi} Berezovskaya F., Karev G., Arditi R.,
\emph{Parametric analysis of the ratio-dependent predator-prey model},
J. Math. Biol. \textbf{43}, (2001) 221-246.
\bibitem{chavez} Brauer F., Castillo-Chavez C.,
\emph{Mathematical Models in Population Biology and Epidemiology}
(Springer-Verlag, New York 2001)
\bibitem{burchard} Burchard H., Deleersnijder E., Meister A.,
\emph{A high-order conservative Patankar-type discretization for stiff systems
of production-destruction equations}, Appl. Numer. Math. \textbf{47:1}, (2003) 1-30.
\bibitem{deang} DeAngelis D. L., Goldstein R. A., O'Neill R. V.,
\emph{A model for trophic interaction}, Ecology \textbf{56}, (1975) 881-892.
\bibitem{pap1} Dimitrov D. T., Kojouharov H. V.,
\emph{Complete mathematical analysis of
predator-prey models with linear prey growth and
Beddington-DeAngelis functional response},
Appl. Math. Comput., \textbf{162:2}, (2005) 523-538.
\bibitem{pap2} Dimitrov D. T., Kojouharov H. V.,
\emph{Analysis and numerical simulation of
phytoplankton-nutrient systems with nutrient loss}, Math. Comput. Simulation,
\textbf{70:1}, (2005) 33-43.
\bibitem{pap3} Dimitrov D. T., Kojouharov H. V.,
\emph{Nonstandard finite-difference schemes
for general two-dimensional autonomous dynamical systems},
Appl. Math. Lett., \textbf{18:7}, (2005) 769-774.
\bibitem{pap5} Dimitrov, D. T. and Kojouharov, H. V.,
\emph{Stability-preserving finite-difference methods for general
multi-dimensional autonomous dynamical systems},
Int. J. Numer. Anal. Model., \textbf{4:2,} (2007) 282-292.
\bibitem{pap6} Dimitrov, D. T. and Kojouharov, H. V.,
\emph{Positive and Elementary Stable
Nonstandard Numerical Methods with Applications to Predator-Prey
Models}, J. Comput. Appl. Math., \textbf{189:1-2}, (2006) 98-108.
\bibitem{hsu} Hsu S.-B., Hwang T.-W., Kuang Y.,
\emph{Global analysis of the Michaelis-Menten-type ratio-dependent
predator-prey system}, J. Math. Biol. \textbf{42}, (2001) 489-506.
\bibitem{kuang} Kuang Y. and Beretta E.,
\emph{Global qualitative analysis of a ratio-dependent predator-prey system},
J. Math. Biol. \textbf{36}, (1998) 389-406.
\bibitem{roux} Lubuma J. M.-S. and Roux A.,
\emph{An improved theta-method for systems of ordinary differential equations},
J. Differ. Equations Appl. \textbf{9:11}, (2003) 1023-1035.
\bibitem{may} May, R. M.,
\emph{Stability and Complexity in Model Ecosystem}
(Princeton Univ. Press, Princeton 1974)
\bibitem{mickens1} Mickens, R. E.,
\textit{Nonstandard finite difference model of differential equations}
(World Scientific, Singapore 1994)
\end{thebibliography}
\end{document}