\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}