- Email: [email protected]
Hopf bifurcation of the stochastic model on HAB nonlinear stochastic dynamics
Chaos, Solitons and Fractals 27 (2006) 1072–1079 www.elsevier.com/locate/chaos
Hopf bifurcation of the stochastic model on HAB nonlinear stochastic dynamics Dongwei Huang b
a,b,*
, Hongli Wang a, Jianfeng Feng a, Zhi-wen Zhu
a
a School of Mechanical Engineering, Tianjin University, 300072, PR China School of Science, Tianjin Polytechnic University, No. 63, ChengLinZhuang Road, HeDong District, 300160 Tianjin, PR China
Accepted 12 April 2005
Abstract A nonlinear stochastic dynamical model on a typical HAB algae diatom and dianoflagellate densities was created and presented in this paper. Simplifying the model through a stochastic averaging method, we obtained a two-dimensional diffusion process of averaged amplitude and phase. The singular boundary theory of diffusion process and the invariant measure theory were applied in analyzing the bifurcation of stability and the Hopf bifurcation of the stochastic system. The critical value of the stochastic Hopf bifurcation parameter was obtained and the conclusion that the position of Hopf bifurcation drifting with the parameter increase is presented as a result. Ó 2005 Elsevier Ltd. All rights reserved.
1. Introduction Marine eutrophication is one of the most important reasons for HAB (harmful algal blooms). There are intricate interactions among harmful algae. Undoubtedly probing into the evolving course of the HAB ecosystem is one of the most important ways of exploring the mechanism of HAB. From the beginning of the 1990s in the last century there have been proposed more and more stochastic dynamical models [1–11] on infectious diseases, population competition, etc. There are fewer stochastic models on HAB ecosystems [12–15]. There also exist many stochastic factors effecting and disturbing the densities of harmful algae which can cause HAB in the realistic marine environment. It is necessary to establish stochastic dynamical models on HAB ecosystem. In this paper, we present a nonlinear stochastic dynamical model on typical HAB algae diatoms and dianoflagellates and analyze its characters of stability and Hopf bifurcation.
2. Nonlinear stochastic dynamic modeling In the paper, we select typical algae diatom and dianoflagellates as objects for establishing a deterministic model as follows: * Corresponding author. Address: School of Science, Tianjin Polytechnic University, No. 63, ChengLinZhuang Road, HeDong District, 300160 Tianjin, PR China. Tel.: +86 22 81329643; fax: +86 22 24528289. E-mail addresses: [email protected], [email protected] (D. Huang).
0960-0779/$ - see front matter Ó 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2005.04.086
D. Huang et al. / Chaos, Solitons and Fractals 27 (2006) 1072–1079
8 dP 1 > > ¼ eP 1 aP 21 aP 1 P 2 ; < dt > > dP 2 ¼ eP bP P bP 2 ; : 2 1 2 2 dt
1073
ð1Þ
where P1, P2 represent the densities of the diatom and the dianoflagellate respectively, e is the internal increasing ratio, a, b represent the restrict coefficients of the algae densities, respectively and a, b represent coefficients of the interaction of the two kinds of algae, respectively. There also exist many stochastic factors effecting and disturbing the realistic marine environment considering the change of densities of diatoms and dianoflagellates suffering sunshine, water temperature, tides, typhoons, tsunamis, human activities and many other stochastic factors. We think it is reasonable and necessary to add random terms in the above system. Thus we rewrite system (1) as follows: 8 dP 1 > > ¼ eP 1 aP 21 aP 1 P 2 þ a1 P 1 nðtÞ þ b1 gðtÞ; < dt ð2Þ > dP > : 2 ¼ eP 2 bP 1 P 2 bP 22 þ a2 P 2 nðtÞ þ b2 gðtÞ. dt We assume that n(t) is the multiplicative random excitation that the derivate of algal densities of diatoms and dianoflagellates suffered, which is related to the environmental and their internal factors. In addition, the derivate of algal densities also suffers the external random excitation g(t) directly. We assume that n(t) and g(t) are independent, in possession of zero mean value and standard variance Gauss white noises. i.e. E [n(t)] = E [g(t)] = 0, E [n(t)n(t + s)] = d(s), E [g(t)g(t + s)] = d(s), E [n(t)g(t + s)] = 0. Eq. (4) is a Stratonovich stochastic differential dynamical system under physical meaning, where d(x) is the Dirac function, all the coefficients of the terms in Eq. (2) are non-negative.
3. Stochastic averaging system of a nonlinear stochastic dynamical model There are four equilibrium points of a deterministic system e e ðb aÞe ða bÞe ; . Q0 ð0; 0Þ; Q1 ; 0 ; Q2 0; ; Q3 a b ab ab ab ab It is obvious that the coordinates of the points listed above are non-negative under a realistic meaning. When a = b, Q3 and Q1 ðae ; 0Þ will be coincident and when b = a, Q3 and Q2 ð0; be Þ will be coincident. If a < b and b > a, there exists a unique stable equilibrium point Q1 ðae ; 0Þ. If a > b and b < a, there exists a unique stable equilibrium point Q2 ð0; be Þ; If a = b and b = a, P2 = cP1 can be obtained from the system, it can be regarded as a single population system. Thus whether ba and ab are equal to 1 at the same time can be regarded as the critical value of stability of the equilibrium points. As for nonlinear stochastic dynamical system (2) the random excitation will cause a change in the stability of the equilibrium points. The critical value will diffuse, so the stochastic system has different stable solution and bifurcation from the deterministic system. According to the above analysis, we analyze stochastic stability and stochastic bifurcation of system (2) as follows. Assume that a > b, b > a, we should deduce the perturbation equation. Let y1 ¼ P 1
ðb aÞe ; ab ab
y2 ¼ P 2
ða bÞe ; ab ab
Y ¼ ½y 1 ; y 2 T .
Then by substituting the corresponding variables in Eq. (2), we can obtain Y_ ¼ AY þ f ðY ; nðtÞ; gðtÞÞ.
ð3Þ
So discussing the stability of system (3) at equilibrium point O(0, 0) is equivalent to discussing the stability of system (2) at equilibrium point Q3. Let 2 3
b a a x1 Y ¼ PX ; X ¼ ; P ¼ 4 a b b 5. x2 1 1
1074
D. Huang et al. / Chaos, Solitons and Fractals 27 (2006) 1072–1079
Then by substituting the corresponding variables in the equations, we obtain X_ ¼ P 1 APX þ P 1 f ðPX ; nðtÞ; gðtÞÞ, i.e. ( x_ 1 ¼ a1 x1 a11 x21 a12 x22 þ ðk 10 þ k 11 x1 þ k 12 x2 ÞnðtÞ þ r1 gðtÞ; ð4Þ x_ 2 ¼ a2 x2 a21 x1 x2 a22 x22 þ ðk 20 þ k 21 x1 þ k 22 x2 ÞnðtÞ þ r2 gðtÞ; where the coefficients are denoted as follows: a1 ¼ e;
a11 ¼ ða bÞ1 D;
a12 ¼ aða bÞb1 ;
a21 ¼ ða bÞ1 ð2ab aa bbÞ;
a2 ¼ ðb aÞða bÞeD1 ;
a22 ¼ ðbb aaÞb1 ;
k 10 ¼ ða bÞ½ðb aÞba1 þ ða bÞaa2 eD1 U 1 ;
U ¼ aa þ bb 2ab;
k 11 ¼ ½ðb aÞba1 þ ða bÞaa2 U 1 ;
k 12 ¼ aða bÞða2 a1 ÞU 1 ;
k 20 ¼ ðb aÞða bÞða2 a1 ÞU 1 D1 be;
k 21 ¼ ðb aÞða2 a1 ÞU 1 b;
k 22 ¼ ½ðb aÞba2 þ ða bÞaa1 U 1 ;
1
r1 ¼ ðbb1 þ ab2 Þða bÞU ;
r2 ¼ ððb aÞb2 ða bÞb1 ÞbU 1 .
We set the coordinate transformation x1 ¼ a cos h; x2 ¼ a sin h, and by substituting the variables in Eq. (4), we obtain 8 da ¼ a2 ½a11 cos3 h þ ða12 þ a21 Þ cos hsin2 h þ a22 sin3 h þ aða1 cos2 h þ a2 sin2 hÞ > dt > > > < þ ½k 10 cos h þ k 20 sin h þ aðk 11 cos2 h þ ðk 12 þ k 21 Þ cos h sin h þ k 22 sin2 hÞnðtÞ þ ðr1 cos h þ r2 sin hÞgðtÞ; dh > dt ¼ ða2 a1 Þ cos h sin h þ a½a12 sin3 h a22 cos hsin2 h þ ða11 a21 Þcos2 h sin h > > > : þ ½k 21 cos2 h þ ðk 22 k 11 Þ cos h sin h k 12 sin2 h þ 1a ðk 20 cos h k 10 sin hÞnðtÞ þ 1a ðr2 cos h r1 sin hÞgðtÞ. ð5Þ It is difficult to calculate the exact solution for system (5) today. According to the Khasminskii limit theorem, when the intensities of the white noises (a1, a2, b1, b2) is small enough, the response process {a(t), h(t)} weakly converged to the two-dimensional Markov diffusion process [16,17]. Through the stochastic averaging method, we obtained the Ito stochastic differential equation (6) the process satisfied da ¼ ma dt þ r11 dW a þ r12 dW h ; ð6Þ dh ¼ mh dt þ r21 dW a þ r22 dW h ; where Wa(t) and Wh(t) are the independent and standard Wiener processes. As for the two-dimensional diffusion process, it is necessary to calculate its two-dimensional transition probability density. There is no general and right method for the calculation. As for the concrete system, we could finish the calculation with some special techniques. We set the parameters as follows: 1 l1 ¼ ða1 þ a2 Þ; l2 ¼ 5k 211 þ 5k 222 þ 3k 212 þ 3k 221 þ 6k 12 k 21 2k 11 k 22 ; 2 1 l3 ¼ ðk 210 þ k 220 þ r21 þ r22 Þ; l4 ¼ 3k 211 þ 3k 222 þ k 212 þ k 221 þ 2k 12 k 21 þ 2k 11 k 22 ; 2 1 l5 ¼ ðk 11 þ k 22 Þðk 21 k 12 Þ; l6 ¼ k 211 þ k 222 þ 3k 212 þ 3k 221 2k 12 k 21 2k 11 k 22 . 4 Under the condition r212 ¼ r221 6¼ 0, we rewrote system (6) as follows: 8 < da ¼ l þ l2 a þ l3 dt þ l þ l4 a2 12 dW a þ ðal Þ12 dW h ; 1 3 5 8 a 8 : dh ¼ ðal Þ12 dW þ l3 þ l6 12 dW . 5
a
a2
16
ð7Þ
h
From the diffusion matrix, we can find that the averaging amplitude a(t) is a one-dimensional Markov diffusing process when r212 ¼ r221 ¼ 0, i.e. k11 + k22 = 0 or k21 k12 = 0. Thus we have the equation as follows: h l li l 12 ð8Þ da ¼ l1 þ 2 a þ 3 dt þ l3 þ 4 a2 dW a . 8 a 8 This is an efficient method to obtain the critical point of stochastic bifurcation through analyzing the change of stability of the averaging amplitude a(t) in the meaning of probability.
D. Huang et al. / Chaos, Solitons and Fractals 27 (2006) 1072–1079
1075
4. Stability of the stochastic averaging system with the singular boundary theory of diffusion process 4.1. The stability at l3 = 0 When l3 = 0 the system can be rewritten as follows: da ¼
h
l1 þ
l 12 l2 i a dt þ 4 a2 dW a . 8 8
ð9Þ
Thus a = 0 is the first kind of singular boundary of system (9) when r11 = 0. When a = +1, we can find ma = 1; thus a = +1 is the second kind of singular boundary of system (9). According to the singular boundary theory, we can calculate the diffusion exponent, drifting exponent and characteristic value of boundary a = 0 and the results are as follows: al ¼ 2;
bl ¼ 1;
2 l1 þ l82 a2 2ð8l1 þ l2 Þ 2ma ða 0Þal bl ¼ limþ ¼ . cl ¼ limþ l4 2 l4 a211 a a!0 a!0 8
ð10Þ
8l1 þl2 > 12, the boundary a = 0 is exclusively natural. l4 8l1 þl2 1 < 2, the boundary a = 0 is attractively natural. l4 8l1 þl2 ¼ 12, the boundary a = 0 is strictly natural. l4
So if cl > 1, i.e. If cl < 1, i.e. If cl = 1, i.e.
We can also calculate the diffusion exponent, drifting exponent and characteristic value of boundary a = +1, and the results are as follows: ar ¼ 2;
br ¼ 1;
cr ¼ lim
a!þ1
2 l1 þ l82 a2 2ma aar br 2ð8l1 þ l2 Þ ¼ lim ¼ . l4 2 a!þ1 l4 a a211 8
ð11Þ
8l1 þl2 < 12, the boundary a = +1 is exclusively natural. l4 8l1 þl2 > 12, the boundary a = +1 is attractively natural. l4 8l1 þl2 ¼ 12, the boundary a = +1 is strictly natural. l4
So if cr > 1, i.e. If cr < 1, i.e. If cr = 1, i.e.
We can draw a conclusion as follows. 2 If 8l1lþl < 12, the boundary a = 0 is attractively natural and the boundary a = +1 is exclusively natural. Thus the 4 trivial solution a = 0 is stable, i.e. the stochastic system is stable at the equilibrium point Q3. This demonstrates that the deterministic system is steady at its equilibrium point, it may also be steady in the meaning of probability at its equilibrium point under random excitations. 2 If 8l1lþl > 12 the boundary a = 0 is exclusively natural and the boundary a = +1 is attractively natural. So the trivial 4 solution a = 0 is unstable, i.e. the stochastic system is unstable at the equilibrium point Q3. This demonstrates that although the deterministic system is steady at its equilibrium point, the stochastic system may be unstable in the meaning of probability at its equilibrium under random excitations. 2 If 8l1lþl ¼ 12, both of the boundaries a = 0 and a = +1 are strictly natural. This case is just the line of demarcation of 4 stability. Whether 16l1 + 2l2 l4 is equal to zero or not can be regarded as the critical condition of bifurcation at the equilibrium point. 4.2. The stability at l3 5 0 Analyzing system (8) if l3 5 0, one can find r11 5 0 at a = 0, so a = 0 is a nonsingular boundary of the system. Through some calculations we can find that a = 0 is a regular boundary (reachable). The other result is ma = 1 when a = +1, so a = +1 is the second singular boundary of system (8). The details are presented as follows: ar ¼ 2;
br ¼ 1;
cr ¼ lim
a!þ1
2 l1 þ l82 a þ la3 a 2ma aar br 2ð8l1 þ l2 Þ ¼ lim : ¼ a!þ1 l4 r211 l3 þ l84 a2
ð12Þ
1076
D. Huang et al. / Chaos, Solitons and Fractals 27 (2006) 1072–1079
So if cr > 1, i.e. If cr < 1, i.e. If cr = 1, i.e.
8l1 þl2 < 12, the boundary a = +1 is exclusively natural. l4 8l1 þl2 > 12, the boundary a = +1 is attractively natural. l4 8l1 þl2 ¼ 12, the boundary a = +1 is strictly natural. l4
Thus we can draw the conclusion that the trivial solution a = 0 is unstable, i.e. the stochastic system is unstable at the equilibrium point Q3 no matter whether the deterministic system is stable at equilibrium point Q3 or not. 4.3. Conclusion by the singular boundary theory 2 > 12, the boundary a = +1 is attractively Summarizing the two kinds of discussions above, one can find that if 8l1lþl 4 natural and the system is unstable in the meaning of probability. So it is almost impossible for the Hopf bifurcation to 2 occur near the equilibrium point Q3. If 8l1lþl < 12, the boundary a = +1, is exclusively natural. What is more if l3 = 0 4 the boundary a = 0 is attractively natural and the system is stable at the equilibrium point Q3 in the meaning of probability. And if l3 5 0 the boundary a = 0 is regular and the system is unstable at the equilibrium point Q3 in the meaning of probability, whether the Hopf bifurcation could occur or not is what we will discuss in the next section.
5. Position of stochastic Hopf bifurcation of a stochastic system Calculating extreme values of the invariant measure is one of the most popular efficient methods in studying the bifurcation of a nonlinear dynamical system. The invariant measure is an important characteristic value of stochastic bifurcation. According to the Ito equation of amplitude a(t), we obtain its FPK equation as follows: op o nh l l i o 1 o2 h l i ¼ l1 þ 2 a þ 3 p þ l3 þ 4 a2 p ; 2 ot oa 2 oa 8 a 8
ð13Þ
with the initial value condition pða; tja0 ; t0 Þ ! dða a0 Þ;
t ! t0 ;
ð14Þ
where p(a, t j a0, t0) is the transition probability density of diffusion process a(t). The invariant measure of a(t) is the steady-state probability density pst(a) which is the solution of the degenerate system as follows: 0¼
o nh l l i o 1 o2 h l i l1 þ 2 a þ 3 p þ l3 þ 4 a2 p . 2 oa 2 oa 8 a 8
ð15Þ
Through calculation, we can obtain rffiffiffi 1 3 2 3v 2v l4 2 1 2 l3 pst ðaÞ ¼ 8 Cð2 vÞ C v a2 ðl4 a2 þ 8l3 Þv2 ; ð16Þ p 2 l3 R 1 x1 t where v ¼ ð8l1 þ l2 Þl1 e dt. 4 ; CðxÞ ¼ 0 t According to NamachivayaÕs theory, the extreme value of an invariant measure contains the most important essence of the nonlinear stochastic system. In other words, the invariant measure can uncover the characteristic information of the steady state. When the intension of the noise tends to zero, the extreme values of pst(a) approximate to show the behavior of the deterministic system. If the process a(t) is ergodic then pst(a) can be regarded as the time measurement for staying in the neighborhood of a(t) according to Oseledec ergodic theorem. From the analysis above, we know that the parameters l1 < 0, l3 > 0, l2 > l4 > 0. If pst(a) has a maximum value at a*, the sample trajectory will stay for a longer time in the neighborhood of a*, i.e. a* is stable in the meaning of probability (with a bigger probability). If pst(a) has a minimum value (zero), it is just the opposite. We now calculate the most possible amplitude a* of system (9), i.e. pst(a) has a maximum value at a*. So we have dpst ðaÞ d2 pst ðaÞ ¼ 0; <0 da a¼a da2 a¼a qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 3 2 as 8l1lþl and the solution a = 0 or a ¼ e a ¼ 8l18l < . The probabilities and the positions of the Hopf bifurcaþl2 l4 2 4 tion occurrence with different parameter are listed in Table 1, and the corresponding results can be seen in Fig. 1 as well.
D. Huang et al. / Chaos, Solitons and Fractals 27 (2006) 1072–1079
1077
Table 1 The probabilities and the positions of the Hopf bifurcation occurrence Condition
Parameter
Cond1 Cond2 Cond3 Cond4
l1 = 0.4, l1 = 0.4, l1 = 0.4, l1 = 0.4,
l2 = 3, l2 = 3, l2 = 3, l2 = 3,
l3 = 0.05, l4 = 2.4 l3 = 0.1, l4 = 2.4 l3 = 0.15, l4 = 2.4 l3 = 0.2, l4 = 2.4
a ¼ ~a
pst ðe aÞ
0.392232 0.554700 0.679366 0.784465
0.887301 0.627417 0.512284 0.443651
Fig. 1. The steady-state probability density pst(a) and the position of stochastic Hopf bifurcation at l1 = 0.4, l2 = 3, l4 = 2.4, l3 = 0.05, 0.1, 0.15, 0.2.
d2 pst ðaÞ 1 2þð8l þl 4l Þl1 ¼ 27þ3ð8l1 þl2 4l4 Þl4 l3 1 2 4 4 > 0; da2 a¼0 8l1lþl2 4 3 l4 ð8l1 þ l2 l4 Þ3 8l3 8l18l þl2 l4 d2 pst ðaÞ ¼ < 0. da2 a¼~a 16l23 ð8l1 þ l2 2l4 Þ3 a . In the meantime, pst(a) is 0 (minimum) at a = 0. This means that the system subjected to Thus what we need is a ¼ e random excitations is almost unsteady at the equilibrium point (a = 0) in the meaning of probability. The conclusion is to go all the way with what has been obtained by the singular boundary theory. The original nonlinear stochastic system has a stochastic Hopf bifurcation at a ¼ e a. x21 þ x22 ¼
8l3 ; 8l1 þ l2 l4
ði:e: a ¼ ~aÞ:
ð17Þ
We now choose some values of the parameters in the equations, draw the graphics of pst(a). The curves in the graph belonging to the cond1, 2 3 4 in turn are shown in Fig. 1. It is worth putting forward that calculating the Hopf
Fig. 2. The steady-state probability density pst(a) and the position of stochastic Hopf bifurcation at l1 = 0.476923, l2 = 2.88, l3 = 0.442147, l4 = 2.88, l5 = l6 = 0.
1078
D. Huang et al. / Chaos, Solitons and Fractals 27 (2006) 1072–1079
bifurcation with the parameters in the original system is necessary. If we now have values of the original parameters in system (2), that e = 0.8, a = 0.6, b = 0.5, a = 0.4, b = 0.1, a1 = 0.6, b1 = 0.3, a2 = 0.6, b2 = 0.1. After further calculations we obtain l1 ¼ 0.476923;
l2 ¼ 2.88;
8l þ l2 1 v¼ 1 ¼ 0.324789 < ; l4 2
l3 ¼ 0.442147; pðaÞ ¼
l4 ¼ 2.88;
l5 ¼ l6 ¼ 0;
32.5567a2 ð3.53717 þ 2.88a2 Þ2.32479
.
What is more is that e a ¼ 0.960531 where pst(a) has the maximum value (see Fig. 2).
6. Conclusion To sum up all the analysis given above, we can draw the conclusion as follows: v can be regarded as the threshold value of stability, and l3 can be regarded as the bifurcation parameter of stochastic Hopf bifurcation. The position where the stochastic Hopf bifurcation occurs will become bigger as l3 increases; on the contrary the probability of Hopf bifurcation will decrease. The significance of the result in the realistic problem can be explained as follows. As the stochastic disturbance is inevitable, it is reasonable to study the Hopf bifurcation of the stochastic dynamical system at the equilibrium point more than the stability of the deterministic system at the equilibrium point, and the most possibility with which the trajectory will stay (occur) in the neighborhood of the limit ring. In this case, the densities of the two kinds of harmful algae arrive at a dynamical balance. The position where stochastic Hopf bifurcation occurs will become bigger (drifting) as l3 increases, and when it reaches the threshold value of HAB, it can cause HAB in the meaning of probability in which the Hopf bifurcation happened. On the other hand, if the stochastic dynamical system has no stochastic Hopf bifurcation, and the boundary a = +1 is attractively natural, it means that the densities of the two algae, at least one of them will increase sharply in a short time interval. This case will cause HAB. According to the expression of l3, it is obvious that as the intensities of the random excitations increase l3 becomes bigger, so it means that harmful algae blooms may occur with a large possibility.
Acknowledgement The work reported in this paper is supported by the National Natural Science Foundation of China (Nos. 10472077, 50375107).
References [1] Edwards AM. Adding detritus to a nutrient-phytoplankton-zooplankton model: a dynamical-systems approach. J Plankton Res 2001;23(4):389–413. [2] Allen LJS, Allen EJ. A comparison of three different stochastic population models with regard to persistence time. Theor Popul Biol 2003;64:439–49. [3] Edwards AM, Bees MA. Generic dynamics of a simple plankton population model with a non-integer exponent of closure. Chaos, Solitons & Fractals 2001;12:289–300. [4] Spagnolo B, La Barbera A. Role of the noise on the transient dynamics of an ecosystem of interacting species. Physica A 2002;315:114–24. [5] Anishchenko VS, Vadivasova TE, Strelkova GI, Okrokvertskhov GA. Statistical properties of dynamical chaos. Math Biosci Eng 2004;1(1):161–84. [6] Webb GF, Blaser MJ, Zhu H, Ardal S, Wu J. Critical role of nosocomial transmission in the toronto SARS outbreak. Math Biosci Eng 2004;1(1):1–13. [7] Gakkhar S, Naji RK. Chaos in seasonally perturbed ratio-dependent prey–predator system. Chaos, Solitons & Fractals 2003;15:107–18. [8] Hensonl SM, King AA, Costantino RF, Cushing JM. Explaining and predicting patterns in stochastic population systems. The Royal Society 2003.6.20, Online. [9] Upadhyaya RK, Raib V, Iyengar SRK. How do ecosystems respond to external perturbations? Chaos, Solitons & Fractals 2000;11:1963–82. [10] Rudnickia R. Long-time behaviour of a stochastic prey–predator model. Stoch Proc Appl 2003;108:93–107. [11] Spagnolo B, Valenti D, Fiasconaro A. Noise in ecosystems: a short review. Math Biosci Eng 2004;1(1):185–211.
D. Huang et al. / Chaos, Solitons and Fractals 27 (2006) 1072–1079
1079
[12] Wang S, Liu Z. Mathematical simulation of red tide blooms process in a closed environment. Oceanol etLimnologia Sinica 1998;29(2):163–8. [13] Wang S, Feng G, Duan M, Liu Z. A model on nutrient dynamics for noctiluca scintillans red tide in Dapeng Bay, Tropic. Oceanology 1997;16(1). [14] Wang H-l, Feng J, Li C, Shen F. Research on nonlinear dynamics of multi-species HABs model. Tianjin University, 2003;36(4). [15] Wang H-l, Zhang Q-c. Nonlinear dynamics theory and application. Tianjin Science and Technology Press; 2002. [16] Zhu W-q. Nonlinear dynamics and control—Hamilton theoretical system frame. Beijing: Science Press; 2003. [17] Cao Q-J. Theoretical study on stochastic vibration and bifurcation of nonlinear system. Doctorate dissertation, Tianjin University, 1991.
Hopf bifurcation of the stochastic model on HAB nonlinear stochastic dynamics Dongwei Huang b
a,b,*
, Hongli Wang a, Jianfeng Feng a, Zhi-wen Zhu
a
a School of Mechanical Engineering, Tianjin University, 300072, PR China School of Science, Tianjin Polytechnic University, No. 63, ChengLinZhuang Road, HeDong District, 300160 Tianjin, PR China
Accepted 12 April 2005
Abstract A nonlinear stochastic dynamical model on a typical HAB algae diatom and dianoflagellate densities was created and presented in this paper. Simplifying the model through a stochastic averaging method, we obtained a two-dimensional diffusion process of averaged amplitude and phase. The singular boundary theory of diffusion process and the invariant measure theory were applied in analyzing the bifurcation of stability and the Hopf bifurcation of the stochastic system. The critical value of the stochastic Hopf bifurcation parameter was obtained and the conclusion that the position of Hopf bifurcation drifting with the parameter increase is presented as a result. Ó 2005 Elsevier Ltd. All rights reserved.
1. Introduction Marine eutrophication is one of the most important reasons for HAB (harmful algal blooms). There are intricate interactions among harmful algae. Undoubtedly probing into the evolving course of the HAB ecosystem is one of the most important ways of exploring the mechanism of HAB. From the beginning of the 1990s in the last century there have been proposed more and more stochastic dynamical models [1–11] on infectious diseases, population competition, etc. There are fewer stochastic models on HAB ecosystems [12–15]. There also exist many stochastic factors effecting and disturbing the densities of harmful algae which can cause HAB in the realistic marine environment. It is necessary to establish stochastic dynamical models on HAB ecosystem. In this paper, we present a nonlinear stochastic dynamical model on typical HAB algae diatoms and dianoflagellates and analyze its characters of stability and Hopf bifurcation.
2. Nonlinear stochastic dynamic modeling In the paper, we select typical algae diatom and dianoflagellates as objects for establishing a deterministic model as follows: * Corresponding author. Address: School of Science, Tianjin Polytechnic University, No. 63, ChengLinZhuang Road, HeDong District, 300160 Tianjin, PR China. Tel.: +86 22 81329643; fax: +86 22 24528289. E-mail addresses: [email protected], [email protected] (D. Huang).
0960-0779/$ - see front matter Ó 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2005.04.086
D. Huang et al. / Chaos, Solitons and Fractals 27 (2006) 1072–1079
8 dP 1 > > ¼ eP 1 aP 21 aP 1 P 2 ; < dt > > dP 2 ¼ eP bP P bP 2 ; : 2 1 2 2 dt
1073
ð1Þ
where P1, P2 represent the densities of the diatom and the dianoflagellate respectively, e is the internal increasing ratio, a, b represent the restrict coefficients of the algae densities, respectively and a, b represent coefficients of the interaction of the two kinds of algae, respectively. There also exist many stochastic factors effecting and disturbing the realistic marine environment considering the change of densities of diatoms and dianoflagellates suffering sunshine, water temperature, tides, typhoons, tsunamis, human activities and many other stochastic factors. We think it is reasonable and necessary to add random terms in the above system. Thus we rewrite system (1) as follows: 8 dP 1 > > ¼ eP 1 aP 21 aP 1 P 2 þ a1 P 1 nðtÞ þ b1 gðtÞ; < dt ð2Þ > dP > : 2 ¼ eP 2 bP 1 P 2 bP 22 þ a2 P 2 nðtÞ þ b2 gðtÞ. dt We assume that n(t) is the multiplicative random excitation that the derivate of algal densities of diatoms and dianoflagellates suffered, which is related to the environmental and their internal factors. In addition, the derivate of algal densities also suffers the external random excitation g(t) directly. We assume that n(t) and g(t) are independent, in possession of zero mean value and standard variance Gauss white noises. i.e. E [n(t)] = E [g(t)] = 0, E [n(t)n(t + s)] = d(s), E [g(t)g(t + s)] = d(s), E [n(t)g(t + s)] = 0. Eq. (4) is a Stratonovich stochastic differential dynamical system under physical meaning, where d(x) is the Dirac function, all the coefficients of the terms in Eq. (2) are non-negative.
3. Stochastic averaging system of a nonlinear stochastic dynamical model There are four equilibrium points of a deterministic system e e ðb aÞe ða bÞe ; . Q0 ð0; 0Þ; Q1 ; 0 ; Q2 0; ; Q3 a b ab ab ab ab It is obvious that the coordinates of the points listed above are non-negative under a realistic meaning. When a = b, Q3 and Q1 ðae ; 0Þ will be coincident and when b = a, Q3 and Q2 ð0; be Þ will be coincident. If a < b and b > a, there exists a unique stable equilibrium point Q1 ðae ; 0Þ. If a > b and b < a, there exists a unique stable equilibrium point Q2 ð0; be Þ; If a = b and b = a, P2 = cP1 can be obtained from the system, it can be regarded as a single population system. Thus whether ba and ab are equal to 1 at the same time can be regarded as the critical value of stability of the equilibrium points. As for nonlinear stochastic dynamical system (2) the random excitation will cause a change in the stability of the equilibrium points. The critical value will diffuse, so the stochastic system has different stable solution and bifurcation from the deterministic system. According to the above analysis, we analyze stochastic stability and stochastic bifurcation of system (2) as follows. Assume that a > b, b > a, we should deduce the perturbation equation. Let y1 ¼ P 1
ðb aÞe ; ab ab
y2 ¼ P 2
ða bÞe ; ab ab
Y ¼ ½y 1 ; y 2 T .
Then by substituting the corresponding variables in Eq. (2), we can obtain Y_ ¼ AY þ f ðY ; nðtÞ; gðtÞÞ.
ð3Þ
So discussing the stability of system (3) at equilibrium point O(0, 0) is equivalent to discussing the stability of system (2) at equilibrium point Q3. Let 2 3
b a a x1 Y ¼ PX ; X ¼ ; P ¼ 4 a b b 5. x2 1 1
1074
D. Huang et al. / Chaos, Solitons and Fractals 27 (2006) 1072–1079
Then by substituting the corresponding variables in the equations, we obtain X_ ¼ P 1 APX þ P 1 f ðPX ; nðtÞ; gðtÞÞ, i.e. ( x_ 1 ¼ a1 x1 a11 x21 a12 x22 þ ðk 10 þ k 11 x1 þ k 12 x2 ÞnðtÞ þ r1 gðtÞ; ð4Þ x_ 2 ¼ a2 x2 a21 x1 x2 a22 x22 þ ðk 20 þ k 21 x1 þ k 22 x2 ÞnðtÞ þ r2 gðtÞ; where the coefficients are denoted as follows: a1 ¼ e;
a11 ¼ ða bÞ1 D;
a12 ¼ aða bÞb1 ;
a21 ¼ ða bÞ1 ð2ab aa bbÞ;
a2 ¼ ðb aÞða bÞeD1 ;
a22 ¼ ðbb aaÞb1 ;
k 10 ¼ ða bÞ½ðb aÞba1 þ ða bÞaa2 eD1 U 1 ;
U ¼ aa þ bb 2ab;
k 11 ¼ ½ðb aÞba1 þ ða bÞaa2 U 1 ;
k 12 ¼ aða bÞða2 a1 ÞU 1 ;
k 20 ¼ ðb aÞða bÞða2 a1 ÞU 1 D1 be;
k 21 ¼ ðb aÞða2 a1 ÞU 1 b;
k 22 ¼ ½ðb aÞba2 þ ða bÞaa1 U 1 ;
1
r1 ¼ ðbb1 þ ab2 Þða bÞU ;
r2 ¼ ððb aÞb2 ða bÞb1 ÞbU 1 .
We set the coordinate transformation x1 ¼ a cos h; x2 ¼ a sin h, and by substituting the variables in Eq. (4), we obtain 8 da ¼ a2 ½a11 cos3 h þ ða12 þ a21 Þ cos hsin2 h þ a22 sin3 h þ aða1 cos2 h þ a2 sin2 hÞ > dt > > > < þ ½k 10 cos h þ k 20 sin h þ aðk 11 cos2 h þ ðk 12 þ k 21 Þ cos h sin h þ k 22 sin2 hÞnðtÞ þ ðr1 cos h þ r2 sin hÞgðtÞ; dh > dt ¼ ða2 a1 Þ cos h sin h þ a½a12 sin3 h a22 cos hsin2 h þ ða11 a21 Þcos2 h sin h > > > : þ ½k 21 cos2 h þ ðk 22 k 11 Þ cos h sin h k 12 sin2 h þ 1a ðk 20 cos h k 10 sin hÞnðtÞ þ 1a ðr2 cos h r1 sin hÞgðtÞ. ð5Þ It is difficult to calculate the exact solution for system (5) today. According to the Khasminskii limit theorem, when the intensities of the white noises (a1, a2, b1, b2) is small enough, the response process {a(t), h(t)} weakly converged to the two-dimensional Markov diffusion process [16,17]. Through the stochastic averaging method, we obtained the Ito stochastic differential equation (6) the process satisfied da ¼ ma dt þ r11 dW a þ r12 dW h ; ð6Þ dh ¼ mh dt þ r21 dW a þ r22 dW h ; where Wa(t) and Wh(t) are the independent and standard Wiener processes. As for the two-dimensional diffusion process, it is necessary to calculate its two-dimensional transition probability density. There is no general and right method for the calculation. As for the concrete system, we could finish the calculation with some special techniques. We set the parameters as follows: 1 l1 ¼ ða1 þ a2 Þ; l2 ¼ 5k 211 þ 5k 222 þ 3k 212 þ 3k 221 þ 6k 12 k 21 2k 11 k 22 ; 2 1 l3 ¼ ðk 210 þ k 220 þ r21 þ r22 Þ; l4 ¼ 3k 211 þ 3k 222 þ k 212 þ k 221 þ 2k 12 k 21 þ 2k 11 k 22 ; 2 1 l5 ¼ ðk 11 þ k 22 Þðk 21 k 12 Þ; l6 ¼ k 211 þ k 222 þ 3k 212 þ 3k 221 2k 12 k 21 2k 11 k 22 . 4 Under the condition r212 ¼ r221 6¼ 0, we rewrote system (6) as follows: 8 < da ¼ l þ l2 a þ l3 dt þ l þ l4 a2 12 dW a þ ðal Þ12 dW h ; 1 3 5 8 a 8 : dh ¼ ðal Þ12 dW þ l3 þ l6 12 dW . 5
a
a2
16
ð7Þ
h
From the diffusion matrix, we can find that the averaging amplitude a(t) is a one-dimensional Markov diffusing process when r212 ¼ r221 ¼ 0, i.e. k11 + k22 = 0 or k21 k12 = 0. Thus we have the equation as follows: h l li l 12 ð8Þ da ¼ l1 þ 2 a þ 3 dt þ l3 þ 4 a2 dW a . 8 a 8 This is an efficient method to obtain the critical point of stochastic bifurcation through analyzing the change of stability of the averaging amplitude a(t) in the meaning of probability.
D. Huang et al. / Chaos, Solitons and Fractals 27 (2006) 1072–1079
1075
4. Stability of the stochastic averaging system with the singular boundary theory of diffusion process 4.1. The stability at l3 = 0 When l3 = 0 the system can be rewritten as follows: da ¼
h
l1 þ
l 12 l2 i a dt þ 4 a2 dW a . 8 8
ð9Þ
Thus a = 0 is the first kind of singular boundary of system (9) when r11 = 0. When a = +1, we can find ma = 1; thus a = +1 is the second kind of singular boundary of system (9). According to the singular boundary theory, we can calculate the diffusion exponent, drifting exponent and characteristic value of boundary a = 0 and the results are as follows: al ¼ 2;
bl ¼ 1;
2 l1 þ l82 a2 2ð8l1 þ l2 Þ 2ma ða 0Þal bl ¼ limþ ¼ . cl ¼ limþ l4 2 l4 a211 a a!0 a!0 8
ð10Þ
8l1 þl2 > 12, the boundary a = 0 is exclusively natural. l4 8l1 þl2 1 < 2, the boundary a = 0 is attractively natural. l4 8l1 þl2 ¼ 12, the boundary a = 0 is strictly natural. l4
So if cl > 1, i.e. If cl < 1, i.e. If cl = 1, i.e.
We can also calculate the diffusion exponent, drifting exponent and characteristic value of boundary a = +1, and the results are as follows: ar ¼ 2;
br ¼ 1;
cr ¼ lim
a!þ1
2 l1 þ l82 a2 2ma aar br 2ð8l1 þ l2 Þ ¼ lim ¼ . l4 2 a!þ1 l4 a a211 8
ð11Þ
8l1 þl2 < 12, the boundary a = +1 is exclusively natural. l4 8l1 þl2 > 12, the boundary a = +1 is attractively natural. l4 8l1 þl2 ¼ 12, the boundary a = +1 is strictly natural. l4
So if cr > 1, i.e. If cr < 1, i.e. If cr = 1, i.e.
We can draw a conclusion as follows. 2 If 8l1lþl < 12, the boundary a = 0 is attractively natural and the boundary a = +1 is exclusively natural. Thus the 4 trivial solution a = 0 is stable, i.e. the stochastic system is stable at the equilibrium point Q3. This demonstrates that the deterministic system is steady at its equilibrium point, it may also be steady in the meaning of probability at its equilibrium point under random excitations. 2 If 8l1lþl > 12 the boundary a = 0 is exclusively natural and the boundary a = +1 is attractively natural. So the trivial 4 solution a = 0 is unstable, i.e. the stochastic system is unstable at the equilibrium point Q3. This demonstrates that although the deterministic system is steady at its equilibrium point, the stochastic system may be unstable in the meaning of probability at its equilibrium under random excitations. 2 If 8l1lþl ¼ 12, both of the boundaries a = 0 and a = +1 are strictly natural. This case is just the line of demarcation of 4 stability. Whether 16l1 + 2l2 l4 is equal to zero or not can be regarded as the critical condition of bifurcation at the equilibrium point. 4.2. The stability at l3 5 0 Analyzing system (8) if l3 5 0, one can find r11 5 0 at a = 0, so a = 0 is a nonsingular boundary of the system. Through some calculations we can find that a = 0 is a regular boundary (reachable). The other result is ma = 1 when a = +1, so a = +1 is the second singular boundary of system (8). The details are presented as follows: ar ¼ 2;
br ¼ 1;
cr ¼ lim
a!þ1
2 l1 þ l82 a þ la3 a 2ma aar br 2ð8l1 þ l2 Þ ¼ lim : ¼ a!þ1 l4 r211 l3 þ l84 a2
ð12Þ
1076
D. Huang et al. / Chaos, Solitons and Fractals 27 (2006) 1072–1079
So if cr > 1, i.e. If cr < 1, i.e. If cr = 1, i.e.
8l1 þl2 < 12, the boundary a = +1 is exclusively natural. l4 8l1 þl2 > 12, the boundary a = +1 is attractively natural. l4 8l1 þl2 ¼ 12, the boundary a = +1 is strictly natural. l4
Thus we can draw the conclusion that the trivial solution a = 0 is unstable, i.e. the stochastic system is unstable at the equilibrium point Q3 no matter whether the deterministic system is stable at equilibrium point Q3 or not. 4.3. Conclusion by the singular boundary theory 2 > 12, the boundary a = +1 is attractively Summarizing the two kinds of discussions above, one can find that if 8l1lþl 4 natural and the system is unstable in the meaning of probability. So it is almost impossible for the Hopf bifurcation to 2 occur near the equilibrium point Q3. If 8l1lþl < 12, the boundary a = +1, is exclusively natural. What is more if l3 = 0 4 the boundary a = 0 is attractively natural and the system is stable at the equilibrium point Q3 in the meaning of probability. And if l3 5 0 the boundary a = 0 is regular and the system is unstable at the equilibrium point Q3 in the meaning of probability, whether the Hopf bifurcation could occur or not is what we will discuss in the next section.
5. Position of stochastic Hopf bifurcation of a stochastic system Calculating extreme values of the invariant measure is one of the most popular efficient methods in studying the bifurcation of a nonlinear dynamical system. The invariant measure is an important characteristic value of stochastic bifurcation. According to the Ito equation of amplitude a(t), we obtain its FPK equation as follows: op o nh l l i o 1 o2 h l i ¼ l1 þ 2 a þ 3 p þ l3 þ 4 a2 p ; 2 ot oa 2 oa 8 a 8
ð13Þ
with the initial value condition pða; tja0 ; t0 Þ ! dða a0 Þ;
t ! t0 ;
ð14Þ
where p(a, t j a0, t0) is the transition probability density of diffusion process a(t). The invariant measure of a(t) is the steady-state probability density pst(a) which is the solution of the degenerate system as follows: 0¼
o nh l l i o 1 o2 h l i l1 þ 2 a þ 3 p þ l3 þ 4 a2 p . 2 oa 2 oa 8 a 8
ð15Þ
Through calculation, we can obtain rffiffiffi 1 3 2 3v 2v l4 2 1 2 l3 pst ðaÞ ¼ 8 Cð2 vÞ C v a2 ðl4 a2 þ 8l3 Þv2 ; ð16Þ p 2 l3 R 1 x1 t where v ¼ ð8l1 þ l2 Þl1 e dt. 4 ; CðxÞ ¼ 0 t According to NamachivayaÕs theory, the extreme value of an invariant measure contains the most important essence of the nonlinear stochastic system. In other words, the invariant measure can uncover the characteristic information of the steady state. When the intension of the noise tends to zero, the extreme values of pst(a) approximate to show the behavior of the deterministic system. If the process a(t) is ergodic then pst(a) can be regarded as the time measurement for staying in the neighborhood of a(t) according to Oseledec ergodic theorem. From the analysis above, we know that the parameters l1 < 0, l3 > 0, l2 > l4 > 0. If pst(a) has a maximum value at a*, the sample trajectory will stay for a longer time in the neighborhood of a*, i.e. a* is stable in the meaning of probability (with a bigger probability). If pst(a) has a minimum value (zero), it is just the opposite. We now calculate the most possible amplitude a* of system (9), i.e. pst(a) has a maximum value at a*. So we have dpst ðaÞ d2 pst ðaÞ ¼ 0; <0 da a¼a da2 a¼a qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 3 2 as 8l1lþl and the solution a = 0 or a ¼ e a ¼ 8l18l < . The probabilities and the positions of the Hopf bifurcaþl2 l4 2 4 tion occurrence with different parameter are listed in Table 1, and the corresponding results can be seen in Fig. 1 as well.
D. Huang et al. / Chaos, Solitons and Fractals 27 (2006) 1072–1079
1077
Table 1 The probabilities and the positions of the Hopf bifurcation occurrence Condition
Parameter
Cond1 Cond2 Cond3 Cond4
l1 = 0.4, l1 = 0.4, l1 = 0.4, l1 = 0.4,
l2 = 3, l2 = 3, l2 = 3, l2 = 3,
l3 = 0.05, l4 = 2.4 l3 = 0.1, l4 = 2.4 l3 = 0.15, l4 = 2.4 l3 = 0.2, l4 = 2.4
a ¼ ~a
pst ðe aÞ
0.392232 0.554700 0.679366 0.784465
0.887301 0.627417 0.512284 0.443651
Fig. 1. The steady-state probability density pst(a) and the position of stochastic Hopf bifurcation at l1 = 0.4, l2 = 3, l4 = 2.4, l3 = 0.05, 0.1, 0.15, 0.2.
d2 pst ðaÞ 1 2þð8l þl 4l Þl1 ¼ 27þ3ð8l1 þl2 4l4 Þl4 l3 1 2 4 4 > 0; da2 a¼0 8l1lþl2 4 3 l4 ð8l1 þ l2 l4 Þ3 8l3 8l18l þl2 l4 d2 pst ðaÞ ¼ < 0. da2 a¼~a 16l23 ð8l1 þ l2 2l4 Þ3 a . In the meantime, pst(a) is 0 (minimum) at a = 0. This means that the system subjected to Thus what we need is a ¼ e random excitations is almost unsteady at the equilibrium point (a = 0) in the meaning of probability. The conclusion is to go all the way with what has been obtained by the singular boundary theory. The original nonlinear stochastic system has a stochastic Hopf bifurcation at a ¼ e a. x21 þ x22 ¼
8l3 ; 8l1 þ l2 l4
ði:e: a ¼ ~aÞ:
ð17Þ
We now choose some values of the parameters in the equations, draw the graphics of pst(a). The curves in the graph belonging to the cond1, 2 3 4 in turn are shown in Fig. 1. It is worth putting forward that calculating the Hopf
Fig. 2. The steady-state probability density pst(a) and the position of stochastic Hopf bifurcation at l1 = 0.476923, l2 = 2.88, l3 = 0.442147, l4 = 2.88, l5 = l6 = 0.
1078
D. Huang et al. / Chaos, Solitons and Fractals 27 (2006) 1072–1079
bifurcation with the parameters in the original system is necessary. If we now have values of the original parameters in system (2), that e = 0.8, a = 0.6, b = 0.5, a = 0.4, b = 0.1, a1 = 0.6, b1 = 0.3, a2 = 0.6, b2 = 0.1. After further calculations we obtain l1 ¼ 0.476923;
l2 ¼ 2.88;
8l þ l2 1 v¼ 1 ¼ 0.324789 < ; l4 2
l3 ¼ 0.442147; pðaÞ ¼
l4 ¼ 2.88;
l5 ¼ l6 ¼ 0;
32.5567a2 ð3.53717 þ 2.88a2 Þ2.32479
.
What is more is that e a ¼ 0.960531 where pst(a) has the maximum value (see Fig. 2).
6. Conclusion To sum up all the analysis given above, we can draw the conclusion as follows: v can be regarded as the threshold value of stability, and l3 can be regarded as the bifurcation parameter of stochastic Hopf bifurcation. The position where the stochastic Hopf bifurcation occurs will become bigger as l3 increases; on the contrary the probability of Hopf bifurcation will decrease. The significance of the result in the realistic problem can be explained as follows. As the stochastic disturbance is inevitable, it is reasonable to study the Hopf bifurcation of the stochastic dynamical system at the equilibrium point more than the stability of the deterministic system at the equilibrium point, and the most possibility with which the trajectory will stay (occur) in the neighborhood of the limit ring. In this case, the densities of the two kinds of harmful algae arrive at a dynamical balance. The position where stochastic Hopf bifurcation occurs will become bigger (drifting) as l3 increases, and when it reaches the threshold value of HAB, it can cause HAB in the meaning of probability in which the Hopf bifurcation happened. On the other hand, if the stochastic dynamical system has no stochastic Hopf bifurcation, and the boundary a = +1 is attractively natural, it means that the densities of the two algae, at least one of them will increase sharply in a short time interval. This case will cause HAB. According to the expression of l3, it is obvious that as the intensities of the random excitations increase l3 becomes bigger, so it means that harmful algae blooms may occur with a large possibility.
Acknowledgement The work reported in this paper is supported by the National Natural Science Foundation of China (Nos. 10472077, 50375107).
References [1] Edwards AM. Adding detritus to a nutrient-phytoplankton-zooplankton model: a dynamical-systems approach. J Plankton Res 2001;23(4):389–413. [2] Allen LJS, Allen EJ. A comparison of three different stochastic population models with regard to persistence time. Theor Popul Biol 2003;64:439–49. [3] Edwards AM, Bees MA. Generic dynamics of a simple plankton population model with a non-integer exponent of closure. Chaos, Solitons & Fractals 2001;12:289–300. [4] Spagnolo B, La Barbera A. Role of the noise on the transient dynamics of an ecosystem of interacting species. Physica A 2002;315:114–24. [5] Anishchenko VS, Vadivasova TE, Strelkova GI, Okrokvertskhov GA. Statistical properties of dynamical chaos. Math Biosci Eng 2004;1(1):161–84. [6] Webb GF, Blaser MJ, Zhu H, Ardal S, Wu J. Critical role of nosocomial transmission in the toronto SARS outbreak. Math Biosci Eng 2004;1(1):1–13. [7] Gakkhar S, Naji RK. Chaos in seasonally perturbed ratio-dependent prey–predator system. Chaos, Solitons & Fractals 2003;15:107–18. [8] Hensonl SM, King AA, Costantino RF, Cushing JM. Explaining and predicting patterns in stochastic population systems. The Royal Society 2003.6.20, Online. [9] Upadhyaya RK, Raib V, Iyengar SRK. How do ecosystems respond to external perturbations? Chaos, Solitons & Fractals 2000;11:1963–82. [10] Rudnickia R. Long-time behaviour of a stochastic prey–predator model. Stoch Proc Appl 2003;108:93–107. [11] Spagnolo B, Valenti D, Fiasconaro A. Noise in ecosystems: a short review. Math Biosci Eng 2004;1(1):185–211.
D. Huang et al. / Chaos, Solitons and Fractals 27 (2006) 1072–1079
1079
[12] Wang S, Liu Z. Mathematical simulation of red tide blooms process in a closed environment. Oceanol etLimnologia Sinica 1998;29(2):163–8. [13] Wang S, Feng G, Duan M, Liu Z. A model on nutrient dynamics for noctiluca scintillans red tide in Dapeng Bay, Tropic. Oceanology 1997;16(1). [14] Wang H-l, Feng J, Li C, Shen F. Research on nonlinear dynamics of multi-species HABs model. Tianjin University, 2003;36(4). [15] Wang H-l, Zhang Q-c. Nonlinear dynamics theory and application. Tianjin Science and Technology Press; 2002. [16] Zhu W-q. Nonlinear dynamics and control—Hamilton theoretical system frame. Beijing: Science Press; 2003. [17] Cao Q-J. Theoretical study on stochastic vibration and bifurcation of nonlinear system. Doctorate dissertation, Tianjin University, 1991.