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