2004 Conference on Diff. Eqns. and Appl. in Math. Biology, Nanaimo, BC, Canada. Electronic Journal of Differential Equations, Conference 12, 2005, pp. 9–19. ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu ftp ejde.math.txstate.edu (login: ftp)

O(ℓ) SHIFT IN HOPF BIFURCATIONS FOR A CLASS OF NON-STANDARD NUMERICAL SCHEMES

MURRAY E. ALEXANDER, SEYED M. MOGHADAS

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 positivitypreserving 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 O(ℓ) shift in the bifurcation parameter value for the numerical scheme, where ℓ 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.

1. 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

2000 Mathematics Subject Classification. 65P30, 65C20, 37M05. Key words and phrases. Finite-difference methods, Codimension-zero bifurcations, Quantitative analysis, Scheme failure. \begin{matrix}c \\ ◯ \end{matrix}2005 Texas State University - San Marcos. Published April 20, 2005. Partially supported by the Natural Sciences and Engineering Research Council of Canada.

models. There are several standard discrete methods which have been widely used in the literature for the purpose of numerical simulations [6]. 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 [7, 8, 10, 11].

One of the most common scheme-failures by standard methods occurs when the continuous model undergoes a Hopf bifurcation leading to its oscillatory behavior [10, 11]. 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 [11] 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 positivitypreserving 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

ẋ = g_{i}(x_{1}, …, x_{m}) − x_{i}h_{i}(x_{1}, …, x_{m}) ≡ f_{i}, i = 1, …, m  (1.1)

where g_{i}, h_{i} : ℝ^{m} → ℝ_{+} are non-negative real valued functions. Many models of biological phenomena can be represented in the form of (1.1); for example, predator prey systems in ecological modelling [2], and the Brusselator system representing an autocatalytic oscillating chemical reaction [17].

This paper is organized as follows. In Section 2, a non-standard finite-difference scheme is constructed to approximate (1.1). 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.

2. 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

f · e_{i} = (f_{1}, …, f_{m}) · (0, …, 1_{ith^{↓}}, …, 0) = g_{i} ≥ 0, i = 1, …, m

EJDE/CONF/12 HOPF BIFURCATIONS FOR NON-STANDARD NUMERICAL SCHEMES 11

which implies that the region ℝ_{+}^{m} = {x = (x_{1}, …, x_{m}) : x_{i} ≥ 0} is positively invariant for system (1.1). To preserve this qualitative property of the CS in the NS, the right hand side of (1.1) is approximated by the advanced stage x_{1}^{n+1}, …, x_{i}^{n+1} for the function g_{i} and h_{i}, and the Gauss-Seidel method [18] is used to compute x_{i}^{n+1} for i = 2, …, m. The approximations are given by

\frac{x_{1}^{n+1} − x_{1}^{n}}{ℓ} = g_{1}(x_{1}^{n}, …, x_{m}^{n}) − x_{1}^{n+1}h_{1}(x_{1}^{n}, …, x_{m}^{n}),
\frac{x_{i}^{n+1} − x_{i}^{n}}{ℓ} = g_{i}(x_{1}^{n+1}, …, x_{i−1}^{n+1}, x_{i}^{n}, …, x_{m}^{n}) − x_{i}^{n+1}h_{1}(x_{1}^{n+1}, …, x_{i−1}^{n+1}, x_{i}^{n}, …, x_{m}^{n}),

which leads to the NS

x_{1}^{n+1} = \frac{x_{1}^{n} + ℓg_{1}(x_{1}^{n}, …, x_{m}^{n})}{1 + ℓh_{1}(x_{1}^{n}, …, x_{m}^{n})} ≡ F_{1}(x_{1}^{n}, …, x_{m}^{n}),  (2.1)
x_{i}^{n+1} = \frac{x_{i}^{n} + ℓg_{i}(x_{1}^{n+1}, …, x_{i−1}^{n+1}, x_{i}^{n}, …, x_{m}^{n})}{1 + ℓh_{i}(x_{1}^{n+1}, …, x_{i−1}^{n+1}, x_{i}^{n}, …, x_{m}^{n})}
≡ F_{i}(x_{i}^{n+1}, …, x_{i−1}^{n+1}, x_{i}^{n}, …, x_{m}^{n}),  (2.2)

for i = 2, …, m, where ℓ > 0 is the step-size. It can be easily shown (by induction) that the scheme remains in ℝ_{+}^{m} as long as the initial value x^{0} = (x_{1}^{0}, …, x_{m}^{0}) is chosen in ℝ_{+}^{m}. Thus, the scheme preserves the positivity property of (1.1).

Here, we provide some preliminarily results which are needed to analyze the scheme for its bifurcation behavior in the next section. Let x^{n+1} = F(x^{n}, x^{n+1}) denote the scheme (2.1)-(2.2), where

F = (F_{1}, …, F_{m}) : ℝ_{+}^{m} × ℝ_{+}^{m} → ℝ_{+}^{m}.

It is easy to check that if x_{∗} ∈ ℝ_{+}^{m} is a critical point of system (1.1), then x_{∗} is a fixed point of the scheme, that is, x_{∗} = F(x_{∗}, x_{∗}). The stability of the fixed point x_{∗} is determined by the absolute values of the eigenvalues z of J_{F}, where z is a solution of the equation

det [J_{(F,n)} − z(𝕀 − J_{(F,n+1)})] = 0,  (2.3)

where 𝕀 is the identity matrix and

(J_{(F,n)})_{ij} ≡ \frac{∂F_{i}}{∂x_{j}^{n}} |_{x_{∗}}, (J_{(F,n+1)})_{ij} ≡ \frac{∂F_{i}}{∂x_{j}^{n+1}} |_{x_{∗}}, for i, j = 1, …, m.

According to the Fixed-Point theorem [14], the fixed point x_{∗} of the NS is stable if |z| < 1, and unstable if |z| > 1. For more detailed analysis and derivation of (2.3), the reader may consult [10].

3. Bifurcations in numerical scheme

To establish a connection between bifurcation behaviors of the NS (2.1)-(2.2) and those of the CS (1.1), we define λ by

z = \frac{1 + ℓλ/2}{1 − ℓλ/2},

which maps the left-half of the λ-plane to the interior of the unit circle in z-plane. Let J_{f} = J_{f}^{U} + J_{f}^{L}, where

J_{f}^{U} = (\begin{matrix}J_{11} & J_{12} & … & J_{1m} \\ 0 & J_{22} & … & J_{2m} \\ ⋮ & ⋮ & ... & ⋮ \\ 0 & 0 & … & J_{mm} \end{matrix}) , J_{f}^{L} = (\begin{matrix}0 & 0 & … & 0 \\ J_{21} & 0 & … & 0 \\ ⋮ & ⋮ & ... & ⋮ \\ J_{m1} & … & J_{m m−1} & 0 \end{matrix}) .

and

J_{ij} ≡ \frac{∂f_{i}}{∂x_{j}} |_{x_{∗}},

for i, j = 1, …, m. Then, by defining

D = diag (h_{1}, …, h_{m})|_{x_{∗}},

and noting that

J_{(F,n)}|_{x_{∗}} = 𝕀 + ℓ(𝕀 + ℓD)^{−1}J_{f}^{U},
J_{(F,n+1)}|_{x_{∗}} = ℓ(𝕀 + ℓD)^{−1}J_{f}^{L},

it follows that

J_{(F,n)}|_{x_{∗}} + J_{(F,n+1)}|_{x_{∗}} = 𝕀 + ℓ(𝕀 + ℓD)^{−1}J_{f} .

Therefore, from (2.3), we have

0 = det [J_{(F,n)} − z(𝕀 − J_{(F,n+1)})]

= det [J_{(F,n)} − \frac{1 + ℓλ/2}{1 − ℓλ/2}(𝕀 − J_{(F,n+1)})]

= \frac{1}{(1 − ℓλ/2)^{m}} det [(1 − ℓλ/2)J_{(F,n)} − (1 + ℓλ/2)(𝕀 − J_{(F,n+1)})]

= \frac{1}{(1 − ℓλ/2)^{m}} det [ℓ(𝕀 + ℓD)^{−1}J_{f} − \frac{ℓλ𝕀}{2} − \frac{ℓλ\{𝕀 + ℓ(𝕀 + ℓD)^{−1}(J_{f}^{U} − J_{f}^{L})\}}{2}]

= \frac{ℓ^{m}}{(1 − ℓλ/2)^{m} det (𝕀 + ℓD)} det [J_{f} − λ\{𝕀 + ℓD + ℓ(J_{f}^{U} − J_{f}^{L})/2\}]

Thus, equation (2.3) can be written in terms of the Jacobian of (1.1), J_{f} , at the critical point x_{∗} as

det [J_{f} − λ(𝕀 + ℓV)] = 0,  (3.1)

where

V ≡ D + \frac{1}{2}(J_{f}^{U} − J_{f}^{L}).

is independent of the step-size ℓ. As an immediate consequence of equation (3.1), it can be seen that as ℓ → 0, λ approaches the eigenvalues of J_{f} . Therefore, the eigenvalues and eigenvectors of the scheme (2.1)-(2.2) are O(ℓ) perturbations of those of the continuous system (1.1). Suppose that (1.1) undergoes a bifurcation at x_{∗} for which α_{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, ℓ, disappear from the characteristic equation (3.1). Thus, the bifurcation of the NS is consistent with that of the CS (1.1), and we have the following theorem.

Theorem 3.1. If the CS (1.1) undergoes a bifurcation at x_{∗} for which α_{0} = 0 is a simple eigenvalue, then the bifurcation of the NS (2.1)-(2.2) occurs at x_{∗} with the same bifurcation parameter value.

Now suppose that system (1.1) undergoes a Hopf bifurcation at x_{∗}. Let us first consider this case for a two-dimensional system (i.e., m = 2). In this case, the equation (3.1) reduces to

\{1 + \frac{1}{2}ℓ[trJ_{f} + 2(h_{1} + h_{2})] + \frac{1}{4}ℓ^{2}[J_{11}J_{22} + J_{12}J_{21} + 2(h_{1}J_{22} + h_{2}J_{11})

+ 4h_{1}h_{2}]}λ^{2} − {trJ_{f} + ℓ(J_{11}J_{22} + h_{1}J_{22} + h_{2}J_{11})}λ + det J_{f} = 0.

(3.2)

Taking into account that det J_{f} = J_{11}J_{22} − J_{12}J_{21} > 0 and trJ_{f} = J_{11} + J_{22} = 0, it can be seen that the coefficient of λ in (3.2) is of O(ℓ). Thus, we have established the following theorem.

Theorem 3.2. If the CS (1.1) with m = 2 undergoes a Hopf bifurcation at x_{∗}, then

Hopf_{NS}|_{x_{∗}} = Hopf_{CS}|_{x_{∗}} + O(ℓ).  (3.3)

Remark 3.3. The equation (3.3) is interpreted to mean that the bifurcation parameter of the NS and the corresponding eigenvalues and eigenvectors are each subjected to an O(ℓ) shift with respect to those of the CS.

Noting that trJ_{f} = 0 for the Hopf bifurcation of the CS (1.1) with m = 2, it follows from equation (3.2) that no O(ℓ) shift occurs for Hopf (Neimark-Sacker) bifurcation of the scheme (2.1)-(2.2) if one of the following holds:

(i) J_{11} = 0 or J_{22} = 0;

(ii) J_{11}J_{22} + h_{1}J_{22} + h_{2}J_{11} = 0.

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 [19] with step-size ℓ as the small parameter and J_{f} defining the zeroth-order problem. Let {α_{1}, …, α_{m}} be the discrete spectrum of J_{f} at x_{∗} with {u_{1}, …, u_{m}} and {w_{1}, …, w_{m}} as the corresponding right and left eigenvectors, respectively, normalized such that

w_{i}^{T} u_{j} = δ_{ij} = \begin{cases}1 & if i = j, \\ 0 & if i ≠ j \end{cases}  (3.4)

Let {λ_{1}, …, λ_{m}} denote the spectrum of J_{f} − λ(𝕀 + ℓV) at x_{∗} with {ũ_{1}, …, ũ_{m}} and {w̃_{1}, …, w̃_{m}} as the corresponding right and left eigenvectors, respectively, such that w̃_{i}^{T} ũ_{j} = δ_{ij}. Then, expanding λ_{i}, ũ_{i}, and w̃_{i}, in terms of α_{i}, u_{i}, and w_{i}, respectively, gives

λ_{i} = α_{i} + ℓα_{i}^{(1)} + ℓ^{2}α_{i}^{(2)} + …,  (3.5)
ũ_{i} = u_{i} + ℓu_{i}^{(1)} + ℓ^{2}u_{i}^{(2)} + …,  (3.6)
w̃_{i} = w_{i} + ℓw_{i}^{(1)} + ℓ^{2}w_{i}^{(2)} + …,  (3.7)

where α_{i}^{(j)}, u_{i}^{(j)} and w_{i}^{(j)} are yet to be determined. Using (3.5) and (3.6), it can be seen that for ℓ > 0

0 = [J_{f} − λ_{i}(𝕀 + ℓV)]ũ_{i} = J_{f} (u_{i} + ℓu_{i}^{(1)} + ℓ^{2}u_{i}^{(2)} + …)
− (α_{i} + ℓα_{i}^{(1)} + ℓ^{2}α_{i}^{(2)} + …)(𝕀 + ℓV)(u_{i} + ℓu_{i}^{(1)} + ℓ^{2}u_{i}^{(2)} + …)

(3.8)

Grouping terms according to powers of ℓ in (3.8) gives:

ℓ^{0} − order : J_{f} u_{i} − α_{i}u_{i} = 0,  (3.9)
ℓ^{1} − order : J_{f} u_{i}^{(1)} − α_{i}(V u_{i} + u_{i}^{(1)}) − α_{i}^{(1)}u_{i} = 0.  (3.10)

Let u_{i}^{(1)} = ∑_{k} a_{ik}^{(1)}u_{k}, and substitute into (3.10) to obtain

∑_{k} a_{ik}^{(1)}(α_{k} − α_{i})u_{k} − α_{i}V u_{i} − α_{i}^{(1)}u_{i} = 0.  (3.11)

Multiplying (3.11) by w_{j}^{T}, and using (3.4), leads to

a_{ik}^{(1)}(α_{k} − α_{i}) − α_{i}(w_{j}^{T} V u_{i}) − α_{i}^{(1)}δ_{ij} = 0

Therefore,

a_{ik}^{(1)} = \frac{α_{i}}{α_{k} − α_{i}} w_{k}^{T} V u_{i}, i ≠ k  (3.12)

and

α_{i}^{(1)} = −α_{i}(w_{i}^{T} V u_{i}), i = k.  (3.13)

It should be noted that a_{ii}^{(1)} can be computed by normalizing of ũ_{i} and w̃_{i}, for example, |ũ_{i}|^{2} = 1, and using w̃_{j}^{T} ũ_{i} = δ_{ij}. Similar results can be obtained for the left eigenvector w̃_{i}. Thus, if the CS undergoes a Hopf bifurcation at x_{∗} (α_{i} ≠ 0, i = 1, …, m), then u_{i}^{(1)} and w_{i}^{(1)} are non-zero vectors, and consequently, the Hopf bifurcation in the NS will be shifted by O(ℓ), unless w_{i}^{T} V u_{i} = 0. Therefore, we have established the following theorem.

Theorem 3.4. If the CS (1.1) undergoes a Hopf bifurcation at x_{∗} with the eigenvalues α_{i} = iω (ω ≠ 0), and w_{i}^{T} V u_{i} ≠ 0 where u_{i} and w_{i} are the corresponding right and left eigenvectors, respectively, then

Hopf_{NS}|_{x_{∗}} = Hopf_{CS}|_{x_{∗}} + O(ℓ).  (3.14)

Remark 3.5. It is worth noting from equations (3.12) and (3.13), that the eigenvalues and eigenvectors at the center eigenspace of the NS remain unshifted to O(ℓ) 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.

4. Examples

In this section, we shall detail our results for two examples of biological phenomena in ecological modelling and an oscillation chemical reaction.

4.1. 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 [2]:

ẋ = rx(1 − x/k) − yϕ(x),  (4.1)
ẏ = y(μψ(x) − D).  (4.2)

where x and y are the prey and the predator population size, respectively, and the dot “ ˙ ” represents the operator d/dt. The parameters r, μ, 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 ϕ(x) is called the functional response of predator to prey, and satisfies the following assumptions:

(A1) ϕ(0) = 0;

(A2) ϕ′(x) > 0 for x ≥ 0; (A3) ϕ′′(x) < 0, for x ≥ 0; (A4) lim_{x→∞} ϕ(x) < ∞.

The predator growth depends on the presence of prey and is proportional to the number of prey. Thus, ψ(x) may be interpreted as the proportion of prey which is eaten by the predators, and satisfies the following assumptions (A5) ψ(0) = 0; (A6) ψ′(x) ≥ 0 for x ≥ 0. There are a number of functional responses satisfying (A1)-(A6), such as Ivlev-

type [5, 15], Rosenzweig [12], sigmoid [3], and Holling-types II and III [4, 16] functional responses. The literature on the dynamics of (4.1)-(4.2) is vast and we refer the reader to [2, 4, 9] for general references citing related studies.

Using the assumptions (A1)-(A6), the model (4.1)-(4.2) can be written in the form of (1.1), where

g_{1}(x, y) = 2rx, h_{1}(x, y) = r(1 + x/k) + yϕ̃(x),
g_{2}(x, y) = 2μyψ(x), h_{2}(x, y) = μψ(x) + D,

with ϕ(x) ≡ xϕ̃(x). In order to illustrate the predictions by Theorems 3.1 and 3.2, we use the results derived in our previous work [9, 10]. 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 μψ(k) < D; and unstable (saddle point) if μψ(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 μψ(k_{0}) = D (see [9]). Using the center manifold theorem, it can be seen [9] that the equation for ξ ≡ −y/r on the center manifold is

ξ̇ = μψ′(k_{0})[ϕ(k_{0})ξ + 2(k − k_{0})]ξ + O(3).  (4.3)

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 θ^{n} ≡ −ξ^{n}/r on the center manifold is given by

θ^{n+1} = θ^{n} + \frac{ℓμψ′(k_{0})}{1 + 2ℓD}[ϕ(k_{0})θ^{n} + 2(k − k_{0})]θ^{n} + O(3).  (4.4)

Equations (4.3) and (4.4) together imply that the transcritical bifurcation of the numerical scheme at E_{k} is consistent with that of the continuous model (4.1)-(4.2), as predicted by Theorem 3.1. For more details about the derivation of equations (4.3) and (4.4), the reader may consult [9].

In previous work [10], we have shown that the model (4.1)-(4.2) 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

ψ(x_{∗}) = \frac{D}{μ} , y_{∗} = \frac{rx_{∗}(1 − x_{∗}/k)}{ϕ(x_{∗})}.

Evaluating the Jacobian of (4.1)-(4.2) at E_{∗} gives J_{11} = r(1 − 2x_{∗}/k) − y_{∗}ϕ′(x_{∗}) and J_{22} = 0. Thus, trJ_{f} = J_{11}, and

J_{11}J_{22} + h_{1}J_{22} + h_{2}J_{11} = h_{2}J_{11},

which is proportional to trJ_{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 (4.1)-(4.2). It can also be seen that (4.1)-(4.2) undergoes a Hopf bifurcation at E_{∗} if k = k_{c} for which J_{11}(k_{c}) = 0 (see [10]). 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 ℓ.

4.2. 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 [1]. The reaction equations for the Brusselator system are given by

A ⟶ X B + X ⟶ Y + C 2X + Y ⟶ 3X X ⟶ D

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 [13, 17]:

Ẋ = B + X^{2}Y − (1 + A)X,  (4.5)
Ẏ = AX − X^{2}Y,  (4.6)

where A, B > 0. More details about this system can be found in [17].

The system (4.5)-(4.6) can be expressed as the form (1.1), where

g_{1}(X, Y) = B + 2X^{2}Y h_{1}(X, Y) = XY + (1 + A)
g_{2}(X, Y) = AX h_{2}(X, Y) = X^{2}

Since (4.5)-(4.6) are of the form given in (1.1), 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 ℝ_{+}^{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 trJ_{c} = A − 1 − B^{2}. Thus, the system undergoes a Hopf bifurcation at E_{c} if trJ_{c} = 0. Evaluating the characteristic equation (3.2) for the NS gives

\begin{matrix}\{1 + \frac{1}{2}ℓ[1 + 5A + B^{2}\{1 + ℓ(1/2 + 2A)\}]\}λ^{2} \\ + \{ℓ(1 + 2A)B^{2} − trJ_{c}\}λ + B^{2} = 0 \end{matrix}  (4.7)

It is clear that the term (1 + 2A)B^{2} in the coefficient of λ is not proportional to trJ_{c}, and does not vanish as long as B > 0. Thus, the numerical scheme undergoes a Hopf bifurcation at E_{c} if ℓ = ℓ_{c} in the equation

B^{2} = \frac{A − 1}{1 + ℓ_{c}(1 + 2A)} = A − 1 + O(ℓ_{c}).  (4.8)

This equation implies that there is an O(ℓ) 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 ℓ < ℓ_{c}.

To illustrate the results for this system, numerical simulations were performed using parameter values A = 2.5 and B = 1. Figure 1 illustrates the profiles of some solutions for ℓ = 0.08 < ℓ_{c}. Note that, in this case, ℓ_{c} = 0.0833333, and the Brusselator system has a stable limit cycle surrounding E_{c}. By increasing the step-size to ℓ = 0.9 > ℓ_{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 2.

It is important to note that, although the NS reproduces the Hopf bifurcation in the CS for ℓ < ℓ_{c}, the actual limit cycle arising from the Hopf bifurcation in the CS will be captured only as ℓ → 0 (See Figure 3).

Figure 1. Profiles of some solutions of the Brusselator system for ℓ = 0.08 < ℓ_{c} approaching the stable limit cycle created by the Hopf bifurcation.

5. Discussion

In this paper, we considered the quantitative dynamics of the m-dimensional continuous system (1.1), by constructing the non-standard finite-difference scheme (2.1)-(2.2). 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 O(ℓ) 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 saddlenode and a Hopf bifurcation), then equation (3.14) shows that the corresponding bifurcation in the NS occurs in general with an O(ℓ) 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 Figure 2. Profiles of some solutions of the Brusselator system for ℓ = 0.09 > ℓ_{c} approaching the unstable equilibrium point E_{c} of the CS.

Figure 3. Phase portraits showing the limit cycles for different values of ℓ < ℓ_{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 ℓ → 0.

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 (2.1)-(2.2) and compared with the corresponding normal form of the CS [9]. 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 ℓ. 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.

References

[1] H. S. Fogler; Elements of Chemical Reaction Engineering, Princton Hall, 1999.

[2] H. I. Freedman; Deterministic Mathematical Models in Population Ecology, 2nd ed. Edmonton: HIFR Conf., 1980.

[3] M. Hesaaraki, S. M. Moghadas; Nonexistence of limit cycles in a predator-prey system with a sigmoid functional response, Canad. Appl. Math. Quart., 7(4) (1999), 401-408.

[4] M. Hesaaraki, S. M. Moghadas; Existence of limit cycles for predator prey systems with a class of functional responses, Ecol. Modelling 142(1-2) (2001), 1-9.

[5] R. E. Kooij, A. Zegeling; A predator-prey model with Ivlev’s functional response, J. Math. Anal. Appl. 188 (1996), 473-489.

[6] J. D. Lambert; Computational Methods in Ordinary Differential Equations, Wiley, New York, 1973.

[7] R. E. Mickens; Discretizations of nonlinear differential equations using explicit nonstandard methods, J. Comput. Appl. Math. 110 (1999), 181-185.

[8] R. E. Mickens; Nonstandard Finite Difference Models of Differential Equations, World Scientific, 1994.

[9] S. M. Moghadas, M. E. Alexander; Bifurcation and numerical analysis of a generalized Gausetype predator-prey model, Dynamics of Continuous, Discrete, and Impulsive Systems Series B, in press.

[10] S. M. Moghadas, M. E. Alexander, B. D. Corbett; A non-standard numerical scheme for a generalized Gause-type predator-prey model, Phys. D 188 (2004) 134-151.

[11] S. M. Moghadas, M. E. Alexander, B. D. Corbett, A. B. Gumel; A positivity preserving Mickens-type discretization of an epidemic model, J. Difference Equ. Appl., 9(11) (2003), 1037-1051.

[12] M. R. Myerscough, M. J. Darwen, W. L. Hogarth; Stability, persistence and structural stability in a classical predator-prey model, Ecol. Modelling 89 (1996), 31-42.

[13] G. Nicolis, I. Prigogine; Self-organization in Non-equilibrium Systems, Wiley, New York, 1977.

[14] G. D. Smith; Numerical Solution of Partial Differential Equations: Finite Difference Methods, Clarendon Press, Oxford, 1985.

[15] J. Sugie; Two-parameter bifurcation in a predator-prey system of Ivlev’s type, J. Math. Anal. Appl. 217 (1998), 349-371.

[16] J. Sugie, R. Kohno, R. Miyazaki; On a predator-prey system of Holling type, Proc. Amer. Math. Soc., 125 (1997), 2041-2050.

[17] J. J. Tyson; Some further studies of non-linear oscillations in chemical systems, J. Chem. Phys., 58 (1973), 3919-3930.

[18] R. S. Varga; Matrix iterative Analysis, Second Edition, Springer-Verlag, Berlin Heidelberg, 2000.

[19] J. H. Wilkinson; The algebraic eigenvalue problem, Chapter 2, Clarendon Press, Oxford, 1965.

Murray E. Alexander Institute for Biodiagnostics, National Research Council Canada, Winnipeg, Manitoba, Canada R3B 1Y6

E-mail address: Murray.Alexander@nrc-cnrc.gc.ca

Seyed M. Moghadas Institute for Biodiagnostics, National Research Council Canada, Winnipeg, Manitoba, Canada R3B 1Y6

E-mail address: Seyed.Moghadas@nrc-cnrc.gc.ca