\documentclass[reqno]{amsart} \usepackage{graphicx} \AtBeginDocument{{\noindent\small 2004 Conference on Diff. Eqns. and Appl. in Math. Biology, Nanaimo, BC, Canada.\newline {\em Electronic Journal of Differential Equations}, Conference 12, 2005, pp. 9--19.\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 2005 Texas State University - San Marcos.} \vspace{9mm}} \setcounter{page}{9} \begin{document} \title[\hfilneg EJDE/Conf/12 \hfil Hopf bifurcations for non-standard numerical schemes] {$\mathcal{O}(\ell)$ shift in Hopf bifurcations for a class of non-standard numerical schemes} \author[M. E. Alexander, S. M. Moghadas \hfil EJDE/Conf/12 \hfilneg] {Murray E. Alexander, Seyed M. Moghadas} % in alphabetical order \address{Murray E. Alexander \hfill\break Institute for Biodiagnostics, National Research Council Canada, Winnipeg, Manitoba, Canada R3B 1Y6} \email{Murray.Alexander@nrc-cnrc.gc.ca} \address{Seyed M. Moghadas \hfill\break Institute for Biodiagnostics, National Research Council Canada, Winnipeg, Manitoba, Canada R3B 1Y6} \email{Seyed.Moghadas@nrc-cnrc.gc.ca} \date{} \thanks{Published April 20, 2005.} \thanks{Partially supported by the Natural Sciences and Engineering Research Council of Canada} \subjclass[2000]{65P30, 65C20, 37M05} \keywords{Finite-difference methods, Codimension-zero bifurcations, \hfill\break\indent Quantitative analysis, Scheme failure} \begin{abstract} Quantitative aspects of models describing the dynamics of biological phenomena have been mostly restricted to results of numerical simulations, often by employing standard numerical methods. However, several studies have shown that these methods may fail to reproduce the actual dynamical behavior of the underlying continuous model when the integration time-step, model parameters, or initial conditions vary in their respective ranges. In this paper, a non-standard numerical scheme is constructed for a general class of positivity-preserving system of ordinary differential equations. A connection between the dynamics of the system and that of the scheme is established in terms of codimension-zero bifurcations. It is shown that when the continuous model undergoes a bifurcation with a simple eigenvalue passing through zero (pitchfork, transcritical or saddle-node bifurcation), the scheme exhibits a corresponding bifurcation at the same bifurcation parameter value. On the other hand, for a Hopf bifurcation there is in general an $\mathcal{O}(\ell)$ shift in the bifurcation parameter value for the numerical scheme, where $\ell$ is the time-step. Partial results for the bifurcations of codimension-1 and higher are also discussed. Finally, the results are detailed for two examples: predator-prey system of Gause-type and the Brusselator system representing an autocatalytic oscillating chemical reaction. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{thm}{Theorem}[section] \newtheorem{rem}[thm]{Remark} \section{Introduction} Dynamical systems theory provides powerful techniques to qualitatively and quantitatively analyze models describing the dynamics of biological phenomena. During the last three decades, much focus has been on the qualitative aspects of these models, such as the stability of the associated equilibria and their bifurcation behaviors. However, the analysis of the quantitative aspects of these models has received less attention. Most studies emphasize the illustrations of the qualitative behaviors, frequently through numerical simulations, rather than analyzing the actual dynamics of quantitative features. The study of the quantitative behavior of dynamical systems requires some fundamental elements, such as methods which discretize the underlying continuous models. There are several standard discrete methods which have been widely used in the literature for the purpose of numerical simulations \cite{L}. However, a number of studies have shown that these methods may fail to reproduce the actual dynamical behavior of the underlying continuous models, as the time-step, model parameters, or initial conditions vary in their respective ranges. The failure to capture the real dynamics of the continuous models by standard methods, which we will call `scheme-failures', has motivated the introduction of non-standard methods \cite{M1,M2,MAC,MACG}. One of the most common scheme-failures by standard methods occurs when the continuous model undergoes a Hopf bifurcation leading to its oscillatory behavior \cite{MAC,MACG}. In this case, the numerical method may diverge, or converge to a solution that does not correspond to any possible solutions of the model. This type of scheme-failure has been explored as a common characteristic of standard methods, as they are generally affected by the length of time-step (step-size), model parameters and initial values of the model variables. The most probable cause of scheme-failure has been the use of insufficiently small step-sizes. However, we have shown in a recent study \cite{MACG} that even with adaptive step-size, the method may still fail to illustrate the qualitative behavior of the model. In this paper, we shall concentrate on the quantitative behavior of a positivity-preserving system of ordinary differential equations, by constructing a non-standard numerical scheme. The goal is to establish a connection between the dynamics of the continuous system (CS) and the numerical scheme (NS), in particular in terms of their bifurcation behavior. The system is given by \begin{equation}\label{cs} \dot x= g_i(x_1,\dots,x_m)-x_ih_i(x_1,\dots,x_m)\equiv f_i,\quad i=1,\dots,m \end{equation} where $g_i,h_i:\mathbb{R}^m\to\mathbb{R}_+$ are non-negative real valued functions. Many models of biological phenomena can be represented in the form of (\ref{cs}); for example, predator prey systems in ecological modelling \cite{F}, and the Brusselator system representing an autocatalytic oscillating chemical reaction \cite{T}. This paper is organized as follows. In Section 2, a non-standard finite-difference scheme is constructed to approximate (\ref{cs}). The scheme is analyzed for its bifurcations in Section 3. To illustrate the results, two examples are given in Section 4: a predator-prey system of Gause-type, and the Brusselator system. The paper ends with a brief discussion of the results. \section{Non-standard numerical scheme} The construction of the numerical scheme is based on the positivity invariance of the CS. Note that on the hyperplane $x_i=0$, we have \begin{equation*} \textbf{f}\cdot \textbf{e}_i=(f_1,\dots,f_m)\cdot(0,\dots,\underset{\overset{\downarrow}{i \text{th}}}{1},\dots,0)=g_i\geq0,\quad i=1,\dots,m \end{equation*} which implies that the region $\mathbb{R}^m_+=\{\mathbf{x}=(x_1,\dots,x_m): x_i\geq0\}$ is positively invariant for system (\ref{cs}). To preserve this qualitative property of the CS in the NS, the right hand side of (\ref{cs}) is approximated by the advanced stage $x_1^{n+1},\dots,x^{n+1}_i$ for the function $g_i$ and $h_i$, and the Gauss-Seidel method \cite{V} is used to compute $x^{n+1}_{i}$ for $i=2,\dots,m$. The approximations are given by \begin{align*} \frac{x_1^{n+1}-x_1^n}{\ell} &= g_1(x_1^n,\dots,x_m^n)-x_1^{n+1}h_1(x_1^n,\dots,x_m^n),\\ \frac{x_i^{n+1}-x_i^n}{\ell} &= g_i(x_1^{n+1},\dots,x_{i-1}^{n+1},x_i^n,\dots,x_m^n) -x_i^{n+1}h_1(x_1^{n+1},\dots,x_{i-1}^{n+1},x_i^n,\dots,x_m^n), \end{align*} which leads to the NS \begin{align} x^{n+1}_1&= \frac{x^n_1+\ell g_1(x^n_1,\dots,x^n_m)}{1+\ell h_1(x^n_1,\dots,x^n_m)}\equiv F_1(x_1^n,\dots,x_m^n),\label{x1}\\ x^{n+1}_i&= \frac{x^n_i+\ell g_i(x^{n+1}_1,\dots,x^{n+1}_{i-1},x^n_i,\dots,x^n_m)}{1+\ell h_i(x^{n+1}_1,\dots,x^{n+1}_{i-1},x^n_i,\dots,x^n_m)}\nonumber\\ &\equiv F_i(x_i^{n+1},\dots,x_{i-1}^{n+1},x_i^n,\dots,x_m^n),\label{xi} \end{align} for $i=2,\dots,m$, where $\ell>0$ is the step-size. It can be easily shown (by induction) that the scheme remains in $\mathbb{R}^m_+$ as long as the initial value $\mathbf{x}^0=(x^0_1,\dots,x^0_m)$ is chosen in $\mathbb{R}^m_+$. Thus, the scheme preserves the positivity property of (\ref{cs}). Here, we provide some preliminarily results which are needed to analyze the scheme for its bifurcation behavior in the next section. Let $\mathbf{x}^{n+1}=\mathbf{F}(\mathbf{x}^n,\mathbf{x}^{n+1})$ denote the scheme (\ref{x1})-(\ref{xi}), where \begin{equation*} \mathbf{F}=(F_1,\dots,F_m):\mathbb{R}^m_+\times \mathbb{R}^m_+\to \mathbb{R}^m_+. \end{equation*} It is easy to check that if $\mathbf{x}_*\in\mathbb{R}^m_+$ is a critical point of system (\ref{cs}), then $\mathbf{x}_*$ is a fixed point of the scheme, that is, $\mathbf{x}_*=\mathbf{F}(\mathbf{x}_*,\mathbf{x}_*)$. The stability of the fixed point $\mathbf{x}_*$ is determined by the absolute values of the eigenvalues $z$ of $J_\mathbf{F}$, where $z$ is a solution of the equation \begin{equation}\label{jns} \det{\big[J_{(\mathbf{F},n)}-z(\mathbb{I}-J_{(\mathbf{F},n+1)})\big]}=0, \end{equation} where $\mathbb{I}$ is the identity matrix and \begin{equation*} \big(J_{(\mathbf{F},n)}\big)_{ij}\equiv\frac{\partial F_i}{\partial x^n_j}\Big|_{\mathbf{x}_*},\ \ \ \ \ \ \big(J_{(\mathbf{F},n+1)}\big)_{ij}\equiv\frac{\partial F_i}{\partial x^{n+1}_j}\Big|_{\mathbf{x}_*}, \ \ \text{for}\ \ i,j=1,\dots,m. \end{equation*} According to the Fixed-Point theorem \cite{Smith}, the fixed point $\mathbf{x}_*$ of the NS is stable if $|z|<1$, and unstable if $|z|>1$. For more detailed analysis and derivation of (\ref{jns}), the reader may consult \cite{MAC}. \section{Bifurcations in numerical scheme} To establish a connection between bifurcation behaviors of the NS (\ref{x1})-(\ref{xi}) and those of the CS (\ref{cs}), we define $\lambda$ by $$ z=\frac{1+\ell\lambda/2}{1-\ell\lambda/2}, $$ which maps the left-half of the $\lambda$-plane to the interior of the unit circle in $z$-plane. Let $J_\mathbf{f}=J^U_\mathbf{f}+J^L_\mathbf{f}$, where \begin{equation*} J^U_\mathbf{f}=\begin{pmatrix} J_{11} & J_{12} & \dots & J_{1m} \\ 0 & J_{22} & \dots & J_{2m} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & J_{mm} \\ \end{pmatrix},\quad J^L_\mathbf{f}=\begin{pmatrix} 0 & 0 & \dots & 0 \\ J_{21} & 0 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ J_{m1} & \dots & J_{m\ m-1} & 0\\ \end{pmatrix}. \end{equation*} and \begin{equation*} J_{ij}\equiv\frac{\partial f_i}{\partial x_j}\Big|_{\mathbf{x}_*}, \end{equation*} for $i,j=1,\dots,m$. Then, by defining \begin{equation*} \mathcal{D}=\mathop{\rm diag}\big(h_1,\dots,h_m\big)\big|_{\mathbf{x}_*}, \end{equation*} and noting that \begin{gather*} J_{(\mathbf{F},n)}\big|_{\mathbf{x}_*}= \mathbb{I}+\ell(\mathbb{I}+\ell\mathcal{D})^{-1}J^U_\mathbf{f},\\ J_{(\mathbf{F},n+1)}\big|_{\mathbf{x}_*}= \ell(\mathbb{I}+\ell\mathcal{D})^{-1}J^L_\mathbf{f}, \end{gather*} it follows that \begin{equation*} J_{(\mathbf{F},n)}\big|_{\mathbf{x}_*}+J_{(\mathbf{F},n+1)}\big|_{\mathbf{x}_*}=\mathbb{I}+\ell(\mathbb{I}+\ell\mathcal{D})^{-1}J_\mathbf{f}. \end{equation*} Therefore, from (\ref{jns}), we have \begin{align*} 0&= \det{\big[J_{(\mathbf{F},n)}-z(\mathbb{I}-J_{(\mathbf{F},n+1)})\big]}\\ &= \det{\Big[J_{(\mathbf{F},n)}-\frac{1+\ell\lambda/2}{1-\ell\lambda/2}(\mathbb{I}-J_{(\mathbf{F},n+1)})\Big]}\\ &= \frac{1}{(1-\ell\lambda/2)^m}\det{\Big[(1-\ell\lambda/2)J_{(\mathbf{F},n)}-(1+\ell\lambda/2)(\mathbb{I}-J_{(\mathbf{F},n+1)})\Big]}\\ &= \frac{1}{(1-\ell\lambda/2)^m}\det{\Big[\ell(\mathbb{I}+\ell\mathcal{D})^{-1}J_\mathbf{f}-\frac{\ell\lambda\mathbb{I}}{2} -\frac{\ell\lambda\big\{\mathbb{I}+\ell(\mathbb{I}+\ell\mathcal{D})^{-1}(J^U_\mathbf{f}-J^L_\mathbf{f})\big\}}{2}\Big]}\\ &= \frac{\ell^m}{(1-\ell\lambda/2)^m\det{(\mathbb{I}+\ell\mathcal{D})}}\det{\big[J_\mathbf{f}-\lambda\big\{\mathbb{I}+\ell\mathcal{D}+\ell(J^U_\mathbf{f}-J^L_\mathbf{f})/2\big\}\big]} \end{align*} Thus, equation (\ref{jns}) can be written in terms of the Jacobian of (\ref{cs}), $J_\mathbf{f}$, at the critical point $\mathbf{x}_*$ as \begin{equation}\label{jcs} \det{\big[J_\mathbf{f}-\lambda(\mathbb{I}+\ell V)\big]}=0, \end{equation} where \begin{equation*} V\equiv\mathcal{D}+\frac{1}{2}\big(J^U_\mathbf{f}-J^L_\mathbf{f}\big). \end{equation*} is independent of the step-size $\ell$. As an immediate consequence of equation (\ref{jcs}), it can be seen that as $\ell\to0$, $\lambda$ approaches the eigenvalues of $J_\mathbf{f}$. Therefore, the eigenvalues and eigenvectors of the scheme (\ref{x1})-(\ref{xi}) are $\mathcal{O}(\ell)$ perturbations of those of the continuous system (\ref{cs}). Suppose that (\ref{cs}) undergoes a bifurcation at $\mathbf{x}_*$ for which $\alpha_0=0$ is a simple eigenvalue, such as transcritical, pitchfork, or saddle-node bifurcation. Then, at the bifurcation point, the terms including the step-size, $\ell$, disappear from the characteristic equation (\ref{jcs}). Thus, the bifurcation of the NS is consistent with that of the CS (\ref{cs}), and we have the following theorem. \begin{thm}\label{thm1} If the CS (\ref{cs}) undergoes a bifurcation at $\mathbf{x}_*$ for which $\alpha_0=0$ is a simple eigenvalue, then the bifurcation of the NS (\ref{x1})-(\ref{xi}) occurs at $\mathbf{x}_*$ with the same bifurcation parameter value. \end{thm} Now suppose that system (\ref{cs}) undergoes a Hopf bifurcation at $\mathbf{x}_*$. Let us first consider this case for a two-dimensional system (i.e., $m=2$). In this case, the equation (\ref{jcs}) reduces to \begin{equation}\label{eq2d} \begin{aligned} &\Big\{1+ \frac{1}{2}\ell\big[\text{tr}J_\mathbf{f}+2(h_1+h_2)\big]+\frac{1}{4}\ell^2\big[ J_{11}J_{22}+J_{12}J_{21}+2(h_1J_{22}+h_2J_{11})\\ &\quad +4h_1h_2\big]\Big\}\lambda^2 -\Big\{\text{tr}J_\mathbf{f}+\ell(J_{11}J_{22}+h_1J_{22}+h_2J_{11})\Big\}\lambda +\det J_\mathbf{f}=0. \end{aligned} \end{equation} Taking into account that $\det J_\mathbf{f}=J_{11}J_{22}-J_{12}J_{21}>0$ and $\text{tr}J_\mathbf{f}=J_{11}+J_{22}=0$, it can be seen that the coefficient of $\lambda$ in (\ref{eq2d}) is of $\mathcal{O}(\ell)$. Thus, we have established the following theorem. \begin{thm}\label{thm2} If the CS (\ref{cs}) with $m=2$ undergoes a Hopf bifurcation at $\mathbf{x}_*$, then \begin{equation}\label{hopf2d} Hopf_{NS}\big|_{\mathbf{x}_*}=Hopf_{CS}\big|_{\mathbf{x}_*}+\mathcal{O}(\ell). \end{equation} \end{thm} \begin{rem}\rm The equation (\ref{hopf2d}) is interpreted to mean that the bifurcation parameter of the NS and the corresponding eigenvalues and eigenvectors are each subjected to an $\mathcal{O}(\ell)$ shift with respect to those of the CS. \end{rem} Noting that $\text{tr}J_\mathbf{f}=0$ for the Hopf bifurcation of the CS (\ref{cs}) with $m=2$, it follows from equation (\ref{eq2d}) that no $\mathcal{O}(\ell)$ shift occurs for Hopf (Neimark-Sacker) bifurcation of the scheme (\ref{x1})-(\ref{xi}) if one of the following holds: \begin{itemize} \item[(i)] $J_{11}=0$ or $J_{22}=0$; \item[(ii)] $J_{11}J_{22}+h_1J_{22}+h_2J_{11}=0$. \end{itemize} Thus, in these cases the Hopf bifurcation of the NS is consistent with that of the CS. In the following, we consider the problem for a general $m$-dimensional system by applying first-order perturbation theory \cite{Wilk} with step-size $\ell$ as the small parameter and $J_\mathbf{f}$ defining the zeroth-order problem. Let $\{\alpha_1,\dots,\alpha_m\}$ be the discrete spectrum of $J_\mathbf{f}$ at $\mathbf{x}_*$ with $\{\mathbf{u}_1,\dots,\mathbf{u}_m\}$ and $\{\mathbf{w}_1,\dots,\mathbf{w}_m\}$ as the corresponding right and left eigenvectors, respectively, normalized such that \begin{equation}\label{normal} \mathbf{w}_i^T\mathbf{u}_j=\delta_{ij}=\begin{cases} 1 &\mbox{if } i=j,\\ 0 &\mbox{if } i\neq j \end{cases} \end{equation} Let $\{\lambda_1,\dots,\lambda_m\}$ denote the spectrum of $J_\mathbf{f}-\lambda(\mathbb{I}+\ell V)$ at $\mathbf{x}_*$ with $\{\tilde{\mathbf{u}}_1,\dots,\tilde{\mathbf{u}}_m\}$ and $\{\tilde{\mathbf{w}}_1,\dots,\tilde{\mathbf{w}}_m\}$ as the corresponding right and left eigenvectors, respectively, such that $\tilde{\mathbf{w}}_i^T\tilde{\mathbf{u}}_j=\delta_{ij}$. Then, expanding $\lambda_i$, $\tilde{\mathbf{u}}_i$, and $\tilde{\mathbf{w}}_i$, in terms of $\alpha_i$, $\mathbf{u}_i$, and $\mathbf{w}_i$, respectively, gives \begin{align} \lambda_i&= \alpha_i+\ell\alpha_i^{(1)}+\ell^2\alpha_i^{(2)}+\dots,\label{e1}\\ \tilde{\mathbf{u}}_i&= \mathbf{u}_i+\ell\mathbf{u}_i^{(1)}+\ell^2\mathbf{u}_i^{(2)}+\dots,\label{e2}\\ \tilde{\mathbf{w}}_i&= \mathbf{w}_i+\ell\mathbf{w}_i^{(1)}+\ell^2\mathbf{w}_i^{(2)}+\dots,\label{e3} \end{align} where $\alpha_i^{(j)}$, $\mathbf{u}_i^{(j)}$ and $\mathbf{w}_i^{(j)}$ are yet to be determined. Using (\ref{e1}) and (\ref{e2}), it can be seen that for $\ell>0$ \begin{equation} \begin{aligned} 0&= \big[J_\mathbf{f}-\lambda_i(\mathbb{I}+\ell V)\big]\tilde{\mathbf{u}}_i = J_\mathbf{f}(\mathbf{u}_i+\ell\mathbf{u}_i^{(1)}+\ell^2\mathbf{u}_i^{(2)}+\dots)\\ &\quad -(\alpha_i+\ell\alpha_i^{(1)}+\ell^2\alpha_i^{(2)}+\dots)(\mathbb{I}+\ell V)(\mathbf{u}_i+\ell\mathbf{u}_i^{(1)}+\ell^2\mathbf{u}_i^{(2)}+\dots) \end{aligned}\label{order} \end{equation} Grouping terms according to powers of $\ell$ in (\ref{order}) gives: \begin{align} \ell^0-\text{order}:&\quad J_\mathbf{f}\mathbf{u}_i-\alpha_i\mathbf{u}_i=0,\label{order0}\\ \ell^1-\text{order}:&\quad J_\mathbf{f}\mathbf{u}_i^{(1)}-\alpha_i(V\mathbf{u}_i+\mathbf{u}_i^{(1)})-\alpha_i^{(1)}\mathbf{u}_i=0.\label{order1} \end{align} Let $\mathbf{u}_i^{(1)}=\sum_ka_{ik}^{(1)}\mathbf{u}_k$, and substitute into (\ref{order1}) to obtain \begin{equation}\label{coeff} \sum_ka_{ik}^{(1)}(\alpha_k-\alpha_i)\mathbf{u}_k-\alpha_iV\mathbf{u}_i-\alpha_i^{(1)}\mathbf{u}_i=0. \end{equation} Multiplying (\ref{coeff}) by $\mathbf{w}_j^T$, and using (\ref{normal}), leads to \begin{equation*} a_{ik}^{(1)}(\alpha_k-\alpha_i)-\alpha_i(\mathbf{w}_j^TV\mathbf{u}_i)-\alpha_i^{(1)}\delta_{ij}=0 \end{equation*} Therefore, \begin{equation}\label{ai} a_{ik}^{(1)}=\frac{\alpha_i}{\alpha_k-\alpha_i}\mathbf{w}_k^TV\mathbf{u}_i,\quad i\neq k \end{equation} and \begin{equation}\label{alphai} \alpha_i^{(1)}=-\alpha_i(\mathbf{w}_i^TV\mathbf{u}_i),\quad i=k. \end{equation} It should be noted that $a_{ii}^{(1)}$ can be computed by normalizing of $\tilde{\mathbf{u}}_i$ and $\tilde{\mathbf{w}}_i$, for example, $|\tilde{\mathbf{u}}_i|^2=1$, and using $\tilde{\mathbf{w}}_j^T\tilde{\mathbf{u}}_i=\delta_{ij}$. Similar results can be obtained for the left eigenvector $\tilde{\mathbf{w}}_i$. Thus, if the CS undergoes a Hopf bifurcation at $\mathbf{x}_*$ ($\alpha_i\neq0$, $i=1,\dots,m$), then $\mathbf{u}_i^{(1)}$ and $\mathbf{w}_i^{(1)}$ are non-zero vectors, and consequently, the Hopf bifurcation in the NS will be shifted by $\mathcal{O}(\ell)$, unless $\mathbf{w}_i^TV\mathbf{u}_i=0$. Therefore, we have established the following theorem. \begin{thm}\label{thm3} If the CS (\ref{cs}) undergoes a Hopf bifurcation at $\mathbf{x}_*$ with the eigenvalues $\alpha_i={\rm i}\omega$ ($\omega\neq0$), and $\mathbf{w}_i^TV\mathbf{u}_i\neq0$ where $\mathbf{u}_i$ and $\mathbf{w}_i$ are the corresponding right and left eigenvectors, respectively, then \begin{equation}\label{hopfmd} \mathop{\rm Hopf}{}_{NS}\big|_{\mathbf{x}_*}=\mathop{\rm Hopf}{}_{CS}\big|_{\mathbf{x}_*}+\mathcal{O}(\ell). \end{equation} \end{thm} \begin{rem}\rm It is worth noting from equations (\ref{ai}) and (\ref{alphai}), that the eigenvalues and eigenvectors at the center eigenspace of the NS remain unshifted to $\mathcal{O}(\ell)$ for codimension-zero bifurcations with simple zero eigenvalue, with respect to those of the CS. This is consistent with the conclusion of Theorem 3.1 for transcritical, pitchfork, and saddle-node bifurcations. \end{rem} \section{Examples} In this section, we shall detail our results for two examples of biological phenomena in ecological modelling and an oscillation chemical reaction. \subsection{Predator-prey system of Gause-type} A well-known model of species interaction in an ecosystem is the predator-prey model of Gause-type given by \cite{F}: \begin{align} \dot{x}&= rx(1-x/k)-y\phi(x),\label{m1}\\ \dot{y}&= y(\mu\psi(x)-D).\label{m2} \end{align} where $x$ and $y$ are the prey and the predator population size, respectively, and the dot `` $\dot\ $ " represents the operator $d/dt$. The parameters $r$, $\mu$, $D$, and $k$ are positive and denote, respectively, the prey intrinsic growth rate, conversion rate of prey to predator, the predator death rate, and carrying capacity of the prey. The function $\phi(x)$ is called the functional response of predator to prey, and satisfies the following assumptions: \begin{itemize} \item[(A1)] $\phi(0)=0$; \item[(A2)] $\phi'(x)>0$ for $x\geq0$; \item[(A3)] $\phi''(x)<0$, for $x\geq0$; \item[(A4)] $\lim_{x\to\infty}\phi(x)<\infty$. \end{itemize} The predator growth depends on the presence of prey and is proportional to the number of prey. Thus, $\psi(x)$ may be interpreted as the proportion of prey which is eaten by the predators, and satisfies the following assumptions \begin{itemize} \item[(A5)] $\psi(0)=0$; \item[(A6)] $\psi'(x)\geq0$ for $x\geq0$. \end{itemize} There are a number of functional responses satisfying (A1)-(A6), such as Ivlev-type \cite{KZ,S}, Rosenzweig \cite{MDH}, sigmoid \cite{HM}, and Holling-types II and III \cite{HM1,SKM} functional responses. The literature on the dynamics of (\ref{m1})-(\ref{m2}) is vast and we refer the reader to \cite{F,HM1,MA} for general references citing related studies. Using the assumptions (A1)-(A6), the model (\ref{m1})-(\ref{m2}) can be written in the form of (\ref{cs}), where \begin{gather*} g_1(x,y)=2rx, \quad h_1(x,y)=r(1+x/k)+y\tilde\phi(x),\\ g_2(x,y)=2\mu y\psi(x), \quad h_2(x,y)=\mu\psi(x)+D, \end{gather*} with $\phi(x)\equiv x\tilde\phi(x)$. In order to illustrate the predictions by Theorems \ref{thm1} and \ref{thm2}, we use the results derived in our previous work \cite{MA,MAC}. It is easy to see that $E_0=(0,0)$ and $E_k=(k,0)$ are two equilibria on the boundary of the positively invariant region of the model. From the linearized system, it can be seen that $E_0$ is a saddle point with the stable manifold on the $y$-axis and unstable manifold on the $x$-axis. It is also easy to check that $E_k$ is locally asymptotically stable if $\mu\psi(k)D$. In fact, the model undergoes a transcritical bifurcation at $E_k$ when the bifurcation parameter $k$ passes through the critical value $k_0$ for which $\mu\psi(k_0)=D$ (see \cite{MA}). Using the center manifold theorem, it can be seen \cite{MA} that the equation for $\xi\equiv-y/r$ on the center manifold is \begin{equation}\label{center1} \dot\xi=\mu\psi'(k_0)\big[\phi(k_0)\xi+2(k-k_0)\big]\xi+\mathcal{O}(3). \end{equation} It is important to note that the associated numerical scheme also undergoes a transcritical bifurcation at the fixed point $E_k$ when $k$ passes through $k_0$, and the equation for $\theta^n\equiv-\xi^n/r$ on the center manifold is given by \begin{equation}\label{center2} \theta^{n+1}=\theta^n+\frac{\ell\mu\psi'(k_0)}{1+2\ell D}\big[\phi(k_0)\theta^n+2(k-k_0)\big]\theta^n+\mathcal{O}(3). \end{equation} Equations (\ref{center1}) and (\ref{center2}) together imply that the transcritical bifurcation of the numerical scheme at $E_k$ is consistent with that of the continuous model (\ref{m1})-(\ref{m2}), as predicted by Theorem \ref{thm1}. For more details about the derivation of equations (\ref{center1}) and (\ref{center2}), the reader may consult \cite{MA}. In previous work \cite{MAC}, we have shown that the model (\ref{m1})-(\ref{m2}) undergoes a Hopf bifurcation at the positive critical point $E_*=(x_*,y_*)$ in the first quadrant when $k$ passes through a critical value $k_c$, where \begin{equation*}\label{E*}\psi(x_*)=\frac{D}{\mu}, \ \ y_*=\frac{rx_*(1-x_*/k)}{\phi(x_*)}.\end{equation*} Evaluating the Jacobian of (\ref{m1})-(\ref{m2}) at $E_*$ gives $J_{11}=r(1-2x_*/k)-y_*\phi'(x_*)$ and $J_{22}=0$. Thus, $\text{tr}J_\mathbf{f}=J_{11}$, and \begin{equation*} J_{11}J_{22}+h_1J_{22}+h_2J_{11}=h_2J_{11}, \end{equation*} which is proportional to $\text{tr}J_\mathbf{f}$. From the discussion following Theorem 3.2, it follows that the Hopf bifurcation of the NS at $E_*$ is consistent with that of the system (\ref{m1})-(\ref{m2}). It can also be seen that (\ref{m1})-(\ref{m2}) undergoes a Hopf bifurcation at $E_*$ if $k=k_c$ for which $J_{11}(k_c)=0$ (see \cite{MAC}). In summary, the scheme reproduces the qualitative behavior of the continuous system such as bifurcations (and consequently the system's asymptotic behavior as well), regardless of the value of the step-size $\ell$. \subsection{Brusselator system} The Brusselator is a system of differential equations which describes an autocatalytic oscillating chemical system in which a species acts to increase the rate of producing reaction \cite{Fogler}. The reaction equations for the Brusselator system are given by \begin{equation*} A\longrightarrow X\quad B+X\longrightarrow Y+C\quad 2X+Y\longrightarrow 3X\quad X\longrightarrow D \end{equation*} where $A$ and $B$ are input chemical species, $C$ and $D$ are reaction products, and $X$ and $Y$ are two autocatalytic species. This mechanism can be mathematically expressed by the following system of differential equations \cite{NP,T}: \begin{align} \dot X&= B + X^2Y - (1+A)X,\label{x}\\ \dot Y&= AX - X^2Y,\label{y} \end{align} where $A, B>0$. More details about this system can be found in \cite{T}. The system (\ref{x})-(\ref{y}) can be expressed as the form (\ref{cs}), where \begin{gather*} g_1(X,Y)=B+2X^2Y \quad h_1(X,Y)=XY+(1+A)\\ g_2(X,Y)=AX \quad h_2(X,Y)=X^2 \end{gather*} Since (\ref{x})-(\ref{y}) are of the form given in (\ref{cs}), it follows that the first quadrant is a positively invariant region. Also, $E_c=(B,A/B)$ is the unique equilibrium point of the system in $\mathbb{R}^2_+$. In order to analyze the bifurcation of the numerical scheme, we evaluate the Jacobian ($J_c$) of the system at $E_c$. A simple calculation yields that $\det J_c=B^2>0$ and $\text{tr}J_c=A-1-B^2$. Thus, the system undergoes a Hopf bifurcation at $E_c$ if $\text{tr}J_c=0$. Evaluating the characteristic equation (\ref{eq2d}) for the NS gives \begin{equation} \begin{aligned} &\Big\{1+\frac{1}{2}\ell\big[1+5A+B^2\{1+\ell(1/2+2A)\}\big]\Big\}\lambda^2 \\ &\quad +\Big\{\ell(1+2A)B^2-\text{tr}J_c\Big\}\lambda+B^2=0 \end{aligned} \end{equation} It is clear that the term $(1+2A)B^2$ in the coefficient of $\lambda$ is not proportional to $\text{tr}J_c$, and does not vanish as long as $B>0$. Thus, the numerical scheme undergoes a Hopf bifurcation at $E_c$ if $\ell=\ell_c$ in the equation \begin{equation}\label{shift} B^2=\frac{A-1}{1+\ell_c(1+2A)}=A-1+\mathcal{O}(\ell_c). \end{equation} This equation implies that there is an $\mathcal{O}(\ell)$ shift in Hopf bifurcation of the numerical scheme as predicted by Theorem 3.2. Thus, the scheme can reproduce the bifurcation (and consequently asymptotic) behaviors of the Brusselator system whenever $\ell<\ell_c$. To illustrate the results for this system, numerical simulations were performed using parameter values $A=2.5$ and $B=1$. Figure \ref{fig1} illustrates the profiles of some solutions for $\ell=0.08<\ell_c$. Note that, in this case, $\ell_c=0.0833333$, and the Brusselator system has a stable limit cycle surrounding $E_c$. By increasing the step-size to $\ell=0.9>\ell_c$, the scheme fails to reproduce the Hopf bifurcation of the system, and therefore it converges to $E_c$ which, in the CS, is an unstable equilibrium. Hence, the scheme-failures occur for capturing the actual dynamics of the continuous system including bifurcation and asymptotic behaviors. This case is illustrated in Figure \ref{fig2}. It is important to note that, although the NS reproduces the Hopf bifurcation in the CS for $\ell<\ell_c$, the actual limit cycle arising from the Hopf bifurcation in the CS will be captured only as $\ell\to0$ (See Figure \ref{fig3}). \begin{figure}[htb] \begin{center} \includegraphics[width=0.7\textwidth]{fig1} \end{center} \caption{Profiles of some solutions of the Brusselator system for $\ell=0.08<\ell_c$ approaching the stable limit cycle created by the Hopf bifurcation.}\label{fig1} \end{figure} \begin{figure}[htb] \begin{center} \includegraphics[width=0.7\textwidth]{fig2} \end{center} \caption{Profiles of some solutions of the Brusselator system for $\ell=0.09>\ell_c$ approaching the unstable equilibrium point $E_c$ of the CS.}\label{fig2} \end{figure} \begin{figure}[htb] \begin{center} \includegraphics[width=0.7\textwidth]{fig3} \end{center} \caption{Phase portraits showing the limit cycles for different values of $\ell<\ell_c$ and $A=1.5$, $B=1$. The sequence of limit cycles approach the actual limit cycle created by the Hopf bifurcation in the CS as $\ell\to0$.}\label{fig3} \end{figure} \section{Discussion} In this paper, we considered the quantitative dynamics of the $m$-dimensional continuous system (\ref{cs}), by constructing the non-standard finite-difference scheme (\ref{x1})-(\ref{xi}). The scheme was analyzed in order to establish a connection between the bifurcation behaviors of the CS and the NS. We showed that, if the CS undergoes a transcritical, pitchfork, or saddle-node bifurcation, then the NS reproduces the bifurcations with the same parameter values, regardless of the value of the step-size. On the other hand, except for well-defined special cases, the bifurcation parameter value for the NS is in general subject to $\mathcal{O}(\ell)$ shift, if the CS undergoes a Hopf bifurcation. The analysis of the NS reflects the bifurcation behavior of the CS for which the corresponding Jacobian has simple eigenvalues. The results can be extended to examine the bifurcations of the NS with repeated eigenvalues. For example, if the CS undergoes a Bogdanov-Takens bifurcation (as a result of merging a saddle-node and a Hopf bifurcation), then equation (\ref{hopfmd}) shows that the corresponding bifurcation in the NS occurs in general with an $\mathcal{O}(\ell)$ shift in the bifurcation parameter. Bogdanov-Takens is a codimension-1 bifurcation with double zero eigenvalue. However, if $m>2$, for bifurcations with repeated eigenvalues, such as degenerate Hopf and others of codimension-1 and higher, it is necessary first to identify the two-dimensional center eigenspace and center manifold in which the bifurcation occurs. On the center manifold, the normal form of the NS can be explicitly derived from (\ref{x1})-(\ref{xi}) and compared with the corresponding normal form of the CS \cite{MA}. This comparison allows us to determine the relationship between the bifurcation behavior of the CS and that of the NS, in terms of the step-size $\ell$. This can be achieved by an extension of the perturbation techniques employed in bifurcation analysis of the NS in Section 3, and is currently under investigation. \begin{thebibliography}{99} \bibitem{Fogler} H. S. Fogler; \emph{Elements of Chemical Reaction Engineering}, Princton Hall, 1999. \bibitem{F} H. I. Freedman; \emph{Deterministic Mathematical Models in Population Ecology}, 2nd ed. Edmonton: HIFR Conf., 1980. \bibitem{HM} M. Hesaaraki, S. M. Moghadas; Nonexistence of limit cycles in a predator-prey system with a sigmoid functional response, {\it Canad. Appl. Math. Quart}., {\bf 7}(4) (1999), 401-408. \bibitem{HM1} M. Hesaaraki, S. M. Moghadas; Existence of limit cycles for predator prey systems with a class of functional responses, {\it Ecol. Modelling} {\bf 142}(1-2) (2001), 1-9. \bibitem{KZ} R. E. Kooij, A. Zegeling; A predator-prey model with Ivlev's functional response, {\it J. Math. Anal. Appl}. {\bf 188} (1996), 473-489. \bibitem{L} J. D. Lambert; \emph{Computational Methods in Ordinary Differential Equations}, Wiley, New York, 1973. \bibitem{M1} R. E. Mickens; Discretizations of nonlinear differential equations using explicit nonstandard methods, {\it J. Comput. Appl. Math.} {\bf 110} (1999), 181-185. \bibitem{M2} R. E. Mickens; \emph{Nonstandard Finite Difference Models of Differential Equations}, World Scientific, 1994. \bibitem{MA} S. M. Moghadas, M. E. Alexander; Bifurcation and numerical analysis of a generalized Gause-type predator-prey model, Dynamics of Continuous, Discrete, and Impulsive Systems Series B, in press. \bibitem{MAC} S. M. Moghadas, M. E. Alexander, B. D. Corbett; A non-standard numerical scheme for a generalized Gause-type predator-prey model, {\it Phys. D} {\bf 188} (2004) 134-151. \bibitem{MACG} S. M. Moghadas, M. E. Alexander, B. D. Corbett, A. B. Gumel; A positivity preserving Mickens-type discretization of an epidemic model, {\it J. Difference Equ. Appl.}, {\bf 9}(11) (2003), 1037-1051. \bibitem{MDH} M. R. Myerscough, M. J. Darwen, W. L. Hogarth; Stability, persistence and structural stability in a classical predator-prey model, {\it Ecol. Modelling} {\bf 89} (1996), 31-42. \bibitem{NP} G. Nicolis, I. Prigogine; Self-organization in Non-equilibrium Systems, Wiley, New York, 1977. \bibitem{Smith} G. D. Smith; \emph{Numerical Solution of Partial Differential Equations: Finite Difference Methods}, Clarendon Press, Oxford, 1985. \bibitem{S} J. Sugie; Two-parameter bifurcation in a predator-prey system of Ivlev's type, {\it J. Math. Anal. Appl}. {\bf 217} (1998), 349-371. \bibitem{SKM} J. Sugie, R. Kohno, R. Miyazaki; On a predator-prey system of Holling type, {\it Proc. Amer. Math. Soc.}, {\bf 125} (1997), 2041-2050. \bibitem{T} J. J. Tyson; Some further studies of non-linear oscillations in chemical systems, {\it J. Chem. Phys.}, {\bf 58} (1973), 3919-3930. \bibitem{V} R. S. Varga; Matrix iterative Analysis, Second Edition, Springer-Verlag, Berlin Heidelberg, 2000. \bibitem{Wilk} J. H. Wilkinson; The algebraic eigenvalue problem, Chapter 2, Clarendon Press, Oxford, 1965. \end{thebibliography} \end{document}