- Email: [email protected]
Temporal localization of limit cycles in a noise-driven chemical oscillator
Physica A 307 (2002) 331 – 353
www.elsevier.com/locate/physa
Temporal localization of limit cycles in a noise-driven chemical oscillator Mazen Al-Ghoul Chemistry Department and Center for Advanced Mathematical Sciences, American University of Beirut, Beirut, Lebanon Received 26 July 2001
Abstract We study the in*uence of multiplicative random external perturbation on the Selkov model exhibiting limit cycles. The dynamics are described by stochastic di0erential equations. The numerically computed ensemble averages of the concentrations exhibit temporal localization of the limit cycles for small amplitudes of the noise. The localization time is smaller when the amplitude of the noise is higher. To understand this behavior, a stochastic complex Ginzburg–Landau equation is derived using normal form transformations on the stochastic evolution equations. In the case of weak noise, the equations for the radial and the phase components of the amplitude separate and can be solved exactly. The solution is used to compute the damping factors, which were found to be Gaussian. We also show that asymptotically, the radial part *uctuates near the average value and the phase is found to be a Gaussian random variable with time-dependent mean and variance and may be responsible for the temporal localization of the limit cycle. c 2002 Elsevier Science B.V. All rights reserved. PACS: 05.10.−a; 05.40.−a; 05.45.−a; 82.20.−w Keywords: Multiplicative noise; Chemical kinetics
1. Introduction The in*uence of internal and external noise on bifurcating dynamical systems has been the subject of several theoretical and experimental studies [1–8]. Environmental *uctuations such as temperature on these systems may alter their behavior compared to those systems without external noise, giving rise, for example, to noise-induced transitions [9], stochastic resonance [10 –12], nonequilibrium @rst-order phase transitions [13], and noise-induced oscillations [14]. This kind of external noise may enter multiplicatively in the macroscopic equations in contrast to internal thermal noise which appears c 2002 Elsevier Science B.V. All rights reserved. 0378-4371/02/$ - see front matter PII: S 0 3 7 8 - 4 3 7 1 ( 0 1 ) 0 0 6 1 2 - 4
332
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
additively in the governing equations that describe the dynamics of the system. Those equations are, in general, nonlinear multidimensional stochastic di0erential equations. Recently, Vlad et al. considered the e0ect of multiplicative noise on a reaction–di0usion system [15] and showed, using the characteristic functional approach, that *uctuations induce a decay or a temporal localization in the averages of the oscillations which can be considered as a temporal analog of Anderson localization of the wavefunction in a disordered system [16]. Vlad later extended his formalism to study localization of limit cycles in the presence of fractal noise [17]. He de@ned localization times, derived damping factors of the oscillations, and found that in the case of white noise, those factors are Gaussian. In this paper, we use a di0erent approach, that of stochastic differential equations, from Vlad’s characteristic functionals to study and characterize the temporal localization of limit cycles in the presence of multiplicative white noise. With this approach, we will be able to obtain similar expressions for the localization time and for the damping factors. Our method is parallel to the study and characterization of the local behavior of deterministic multidimensional nonlinear systems, where one may use the center manifold theorem to reduce the dimensionality of the system, and the normal form method to eliminate nonlinearities by means of successive nonlinear transformations [18]. The obtained normal form equations have a universal character and they are derived in several di0erent areas of physics [19]. The extension of the normal form reduction scheme to stochastic di0erential equations requires the restrictive physical assumptions of weak noise. Several studies have implemented the “stochastic” center manifold to reduce stochastic systems to a more simpli@ed form. Using this approach, one can obtain a reduced Fokker–Planck equation [20,21] that can be written in a normal form by means of local nonlinear transformations. On the other hand, one may apply the normal form transformations directly to the stochastic di0erential equations [22]. For a Hopf bifurcation with external random parameters, Coullet et al. [23] derived a general normal form which contains new “stochastic resonance terms” due to the stochastic nature of the equations. Both approaches yield similar results but the former is relatively easier and can be used directly to interpret the numerical solutions of the stochastic di0erential equations as will be shown in this paper. In this work, we study the dynamics in a model of a limit cycle with imposed multiplicative random forcing. We use the modi@ed Selkov model [24], which contains four adjustable parameters, and we vary randomly one dimensionless parameter, ; which appears multiplicatively in the corresponding chemical kinetic equations. The resulting stochastic di0erential equations are solved numerically in the case of weak noise; the solutions and their ensemble averages for di0erent realizations of the noise are computed and it was found that these ensemble averages exhibit damping of the limit cycle. To understand this damping, we derive a stochastic Landau–Ginzburg (LG) equations using nonlinear transformations of coordinates. Under the white weak noise assumption, the stochastic LG equations separate into two independent equations for the radial and phase components of the complex amplitude. This decay in the oscillations may due to the stochastic nature of the phase which we are going to characterize in this paper by calculating the damping factors in real and Fourier spaces. In Section 2, we present and discuss the numerical solutions of the stochastic di0erential equations. In Section 3, we use normal form transformations to derive a stochastic
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
333
Ginzburg–Landau equation which is solved exactly in the case of weak noise. Then localization times and damping factors are easily de@ned. In the conclusion, we discuss possible extensions of this work to colored noise. In Appendix A, we perform a linear stability analysis of the deterministic equations and de@ne di0erent quantities and expressions. Appendix B contains the numerical scheme used to integrate the stochastic di0erential equations. 2. Weak random forcing The modi@ed Selkov model is de@ned through the following kinetic equations (see Appendix A for de@nitions of various quantities): dt X = B − X + X 2 Y − KX 3 ; dt Y = A − Y − X 2 Y + KX 3 ;
(1)
which can be put in the following more compact form: dt Z = R ; where
Z=
and
(3)
Y
R=
X
(2)
B − X + X 2 Y − KX 3 A − Y − X 2 Y + KX 3
:
(4)
We now assume that the external variations are given in terms of the variation of the parameter : This parameter is assumed to vary randomly according to the following equation: = L + dt W (t) ;
(5)
where L is an average value, dt = d=dt; and dt W (t) is a Gaussian white noise term which we assume to be the time derivative of the Wiener process W (t); is the amplitude of the random noise. Therefore, the properties of the noise are given by dt W (t) = 0 ;
(6)
dt W (t1 )dt W (t2 ) = 2 (t2 − t1 ) :
(7)
If we substitute Eq. (5) into Eq. (2), we obtain the following Ito stochastic di0erential equation: dZ = R dt + Z dW (t) ; where is given by 0 0 = : 0 1
(8)
(9)
334
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
In the deterministic case, i.e., = 0; the system described by Eq. (2) undergoes a Hopf bifurcation when A = 0:5; B = 0:09; and L 6 0:1241: In this paper, the numerical results were obtained by taking L = 0:1: In the case of weak noise,
0:1; we show in the next section that the structure of the phase diagram is not altered and Eq. (8) still admits a set of stationary points very close to the deterministic stationary points. This phase diagram comprises an oscillations region and a stable region separated by the Hopf bifurcation line. The noise is chosen such that con@nes the system inside the oscillations region and crossing to the stable region is an extremely rare event. We integrated numerically Eq. (8) for several small noise amplitudes ranging from
≈ 0:0001 to 0:03. In every case, we compute and we record the times its value crosses back and forth the oscillations region. For in the interval [0:0001; 0:03] no crossing to the stable regions was observed. We describe the numerical scheme in Appendix B. For every realization of the noise we obtained a trajectory. We performed at least 50 realizations for every case we studied in order to perform an ensemble average; this number was suMcient in most of the cases studied. All cases in Fig. 1a–f lie entirely in the oscillations region; we see that for an amplitude lower than 0:01, the noise does not a0ect much a given realization, and the trajectories resemble those of the deterministic case [25]. It is important to note here, that the noise did not delay the bifurcation and a stable limit cycle continue to exist around stationary points that are very close to stationary point of the deterministic system (X0 , and Y0 from Eq. (64)) [25]. However, as increases we can see that the behavior begins to change towards what may be a chaotic regime. As for the average over noise realizations of the concentrations, for any small value we take for ; we can see in Fig. 2a–f that the oscillations decay toward a region near the deterministic values of the steady state. The higher the amplitude of the noise ; the faster the decay. The computed spectrum in Fig. 3, computed for the average shows that an increase in the noise amplitude results in reduction of the number of harmonics present initially in the deterministic system. In addition, the amplitude of the principal frequency decreases and broadens towards the @rst harmonic. Above = 0:03; the peaks are broad enough to demonstrate that the ensemble average has a chaotic character: 3. Normal forms and localization In this section, we present an approximate analysis of the kinetic equations of the stochastically driven Selkov model, where we use an appropriate linearization scheme of the stochastic di0erential equations and then invoke nonlinear coordinate transformations, similar to normal form methods in the theory of deterministic dynamical systems, to bring the equations into an easier solvable form. We show that near the Hopf bifurcation and for small noise amplitudes, the obtained stochastic normal form equation can be used to describe qualitatively the behavior of the system. Because we need to make a series of operations and coordinate transformation, it is better to interpret the stochastic di0erential equation (8) as a Stratonovitch di0erential equation, where we can make use of the easier rules of the ordinary calculus. To do
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
335
Fig. 1. The variation of the concentration of species X in time for one realization of the noise with di0erent intensities : (a) = 0:0001, (b) = 0:005, (c) = 0:009, (d) = 0:01, (e) = 0:02, (f) = 0:03:
336
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
Fig. 2. The variation of the ensemble average of concentration of species X in time for di0erent noise intensities : (a) = 0:0001, (b) = 0:005, (c) = 0:009, (d) = 0:01, (e) = 0:02, (f) = 0:03:
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
337
Fig. 3. Power spectra of oscillations in the presence of noise of di0erent amplitudes: = 1:0 × 10−4 dashed line, = 1:0 × 10−2 solid line.
so, we correct the drift coeMcient R; which is calculated in Appendix A, as follows: ˜ = R − 1 2 (Z) @ (Z) ; R 2 @Z T 0 0
2 0 ; =R− 2 Y 0 1
2 =R− 2 =R−
0
0
0
1
X
Y
2 (Z) : 2
; (10)
The Stratonovitch stochastic di0erential equations corresponding to the drift coeMcient ˜ is given by R ˜ + (Z) ◦ dW (t) ; dt Z = R
(11)
where ◦dW indicates a Stratonovitch type of measure. For the weak noise limit, i.e.,
is small, we can linearize Eq. (11) around the state Z0 = (X0 ; Y0 ) which is an unstable ˜ steady state [26,27] of the nonstochastic system (2). The nonlinear drift coeMcient R is linearized by replacing it by the following Jacobian matrix D evaluated at (X0 ; Y0 ) : 0 X02 −1 + 2X0 Y0 − 3X02 ˜ i
2 0 @R : (D)ij = − = @Zj 2 0 1 −2X0 Y0 + 3X02 −X02 − L (X0 ; Y0 )
(12)
338
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
The di0usion term is already linear and preserves its form. The linearized system assumes the following form: dt Zl = DZl + (Zl ) ◦ dW :
(13)
Since the equation is already linear in the noise term, and because the normal form depends only on the number and nature of the eigenvalues of the linearized problem, we suggest to reduce Eq. (11) to a normal form by a series of nonlinear transformations of coordinates exactly like the deterministic case [28] and leave the noise term intact. The resulting equation corrects for the linearized equation (13). Let = ( ) + i!( ) and L = ( ) − i!( ) be two complex eigenvalues of the matrix D near the Hopf bifurcation (( ) and !( ) are real), and let ’; ’L be the corresponding complex eigenvectors. For small ; ( ) and !( ) are given by ( ) = 0 − 14 2 + O( 4 ) ; !( ) = !0 +
(14)
0 2
+ O( 4 ) ; 4!0
(15)
where 0 and !0 are, respectively, the real and imaginary part of the eigenvalues of the problem without noise. It is clear from Eqs. (14) and (15) that noise can delay the bifurcation by shifting 0 to ( ); however, in the limit of weak noise, ( ) ≈ 0 ; and !( ) ≈ !0 ; and therefore the structure of the phase diagram is not altered near the Hopf bifurcation point in the limit of weak noise. The next step is to make a coordinates transformation from the real and imaginary part of one of the eigenvectors to bring Eq. (11) into its corresponding normal form. More precisely, if for example, ’ = (a1 + ib1 ; a2 + ib2 ) we can construct the following transition matrix T a b b2 a2 = ; (16) T= 1 0 b1 a1 which acts on Z as J = T−1 Z
(17)
to give a new set of coordinates. J is the vector (U; V )T . In this new coordinate system, Eq. (11) has the following form: dJ = [T−1 DTJ + T−1 N(TJ)] dt + ( T−1 TJ) ◦ dW L ◦ dW ; L + N) L dt + J = (DJ where
N=
N1 (U 2 ; UV; U 2 V; U 3 ) N2 (U 2 ; UV; U 2 V; U 3 )
(18)
;
(19)
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
L = T−1 DT = D L = T
−1
T =
−!
!
1
0
− ba
0
339
;
(20)
(21)
and NL = T−1 N(TJ): Using the normal form theorem [28], a couple of standard nonlinear transformations (near-identity transformations) eliminate all quadratic and cubic nonlinearities in NL except for a couple of cubic nonresonant terms. We @nally end up with the following normal form equation in complex form to cubic order a dz = [( + i!)z + (cr + ici )|z|2 z] dt + 1 − Re(z) + O( 3 ) ◦ dW ; b (22) where z =U +iV . cr and ci can be computed analytically [29]. For a supercritical Hopf bifurcation, cr should be negative. For the case of weak noise, we will drop the higher order terms and the nonlinear terms of O( 3 ): Also, the term a=b; when computed numerically, depends weakly on so that in the weak noise limit, it is assumed that it is constant. Therefore, Eq. (22) is written as dz = [( + i!)z + (cr + ici )|z|2 z] dt + Re(z) ◦ dW ;
(23)
where = (1−a=b): Now we will represent Eq. (23) in polar form. Let z =R exp(i ); then Eq. (23) is decomposed into the following radial and phase components: dR = (R + cr R3 ) dt + R cos2 ( ) ◦ dW ;
(24)
d = (! + ci R2 ) dt − 12 sin(2 ) ◦ dW :
(25)
Motivated by results of Fig. 3, that in the weak noise limit the perturbed limit cycle oscillates with a frequency that is close to the frequency of the unperturbed limit cycle, !0 ; and in order to make some analytical progress we propose the following mean-@eld approximation to simplify further Eqs. (24) and (25): in the weak noise limit,
0:1, the phase appearing in the right-hand sides of those equations is replaced by an average quantity → ≈ !0 t + 0 , which is the expected behavior in the limit of zero noise. This approximation will later provide a substantial agreement with the damping factor of Vlad et al. [15,17], and with the numerical solutions of the initial stochastic di0erential equations (8). Therefore, Eqs. (24) and (25) are approximated by dR = (R + cr R3 ) dt + R cos2 (!0 t + 0 ) ◦ dW ;
(26)
d = (!0 + ci R2 ) dt − 12 sin(2!0 t + 2 0 ) ◦ dW :
(27)
To test the approximation, we solved Eqs. (24) and (25) and Eqs. (26) and (27) numerically for the same parameters and same realization of the Wiener process. We
340
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
Fig. 4. Solution of the radial part, R(t); of the amplitude equations from Eqs. (24) and (25) [higher curve], and Eqs. (26) and (27) [lower curve] for the same realization of the Wiener process Wt . = 0:0005: Notice the convergence of the two curves after some time.
also used the numerical scheme described in Appendix B. In Fig. 4, we see how the two solutions for the radial component of the amplitude of the two versions of the amplitude equations converge to each other after some time in the weak noise limit. As the noise increases to about = 0:05, the solutions do not match very well. 3.1. The radial part of the amplitude Eq. (26) is a nonlinear stochastic di0erential equation. The change of variable Y = R−2 transform Eq. (26) to a linear stochastic di0erential equation which can then be solved explicitly to give t exp( t + 0 cos2 (!0 s + 0 ) dWs ) R(t) = ; (28) t t 1 + 2R20 0 exp(2 s + 2 0 cos2 (!0 s + 0 ) dWs ) ds t where = =|cr |: The stochastic integral I(t) = 0 cos2 (!0 s + 0 ) dWs is a Gaussian stochastic process with zero mean and variance given by [27]
t Var(I) = cos4 (!0 s + 0 ) ds 0
3 2 sin(!0 t + 0 ) cos3 (!0 t + 0 ) + 3 sin[2(!0 t + 0 )] = t+ : 8 8!0
(29)
In long times, the variance in Eq. (29) is approximated as Var(I) = 38 t because the second term on the RHS of Eq. (29) is bounded. Therefore, by Doob’s theorem [30], asymptotically the stochastic process I(t) resembles the Wiener process Wt . If we
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
341
make the transformation ! = 38 t or rescale the amplitude of the noise ; then we can write R(t) in Eq. (28) as ˜ t) R0 exp( t + W ; t ˜ s ) ds 1 + 2R20 0 exp(2 s + 2 W
R(t) =
(30)
where = =|cr |; ˜ = 83 ; Wt is a Wiener process and R0 is an initial condition at t = 0. In the case of zero noise, i.e., = 0; R(t) is given by R0 exp( t)
R(t) =
1 + (R20 =)[exp(2 t) − 1]
→
as t → ∞
(31)
as expected. We anticipate that the asymptotic behavior of the stochastic process R(t) for small to exhibit small *uctuations around a certain average. To prove this, we need to study the statistics of the stochastic process R(t) in Eq. t (30). Therefore, we ˜ s ) ds @rst need to understand the time dependence of the integral 0 exp(2 s + 2 W which appears in the denominator of Eq. (30). The idea is to show that the process R(t) very closely to an average value; this is the case only if the integral t *uctuates ˜ s ) vary at the same rate so that ratio ˜ s ) ds and the exp(2 s + 2 W exp(2 s + 2
W 0 (30) *uctuates very close to a constant value. Evidently, the stochastic process R(t) is not Gaussian. However, if the moments higher than the second of this process are small, then we can safely approximate it with a Gaussian. We notice that the process R(t) in Eq. (30) is a solution of the following Ito stochastic di0erential equation: ˜ dW ; dR = ( R − R3 ) d! + R
(32)
2 where = − 12 ˜ and ! = |cr |t: The ensemble average R(t) is easily obtained from averaging over the following stationary distribution function P(R) derived from the Fokker–Planck equation with natural boundary conditions: 2
P(R) = KR" exp(−R2 = ˜ ) ;
(33)
where K is a normalization constant given by −"−1
2 ˜ K= #("=2 + 12 )
(34)
and the exponent " is given by "( ) =
2 3( ) −2= 2 −3: 2 4 cr ( )(1 − a( )=b( ))2
˜
(35)
For P(R) to be a probability density [31], it has to be integrable over R+ , i.e., "+1 ¿ 0. For A = 0:5; B = 0:09; L = 0:1; (" + 1) is positive at a noise amplitude 6 0:6642: Therefore, it is safe to assume that a probability density of the form P(R) in Eq. (33) exists whenever the noise is weak. Hence, all the moments of the stochastic process
342
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
Fig. 5. The process R(t) and its ensemble average R(t) (thicker line) for di0erent noise amplitudes: (a)
= 0:01, (b) = 0:05, (c) = 0:1, (d) = 0:3:
R can be computed from the following Riemann integral: Rm =
0
∞
+ "=2 + 12 ) : #("=2 + 12 )
m #(m=2
Rm P(R) dR = ˜
(36)
In particular, the average and second moment are given respectively by #("=2 + 1) R = ˜ ; #("=2 + 12 )
(37)
2 R2 = 12 ˜ (" + 1) :
(38)
In the weak noise limit, as → 0, the moments in Eq. (36) tend to zero very quickly. Therefore, it is a good approximation to assume that asymptotically R(t) is a Gaussian random variable characterized by Eqs. (37) and (38). Moreover, the second moment (38) is small for weak noise. The numerical evaluation of R(t) in Eq. (30) is shown in Fig. 5a and b where the amplitude of the noise is considered small. In this case,
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
343
the realization of R(t) *uctuates closely to the average R(t) in long times. As becomes large, the *uctuations grow around the ensemble average [Fig. 5c and d], in a similar way to that in Fig. 5a and b, where the average clearly lies in between the realizations of the process R(t): When takes large values, no stationary solution of the corresponding Fokker–Planck exists, and high amplitude *uctuations in the realizations of R(t) starts to appear [Fig. 5d]. 3.2. The phase di/usion and damping functions The equation for phase (25) will then be given by d = (!0 + ci R2 ) dt − 12 sin(2!0 t + 2 0 ) ◦ dW ; = !˜ dt − 12 sin(2!0 t + 2 0 ) dW ;
(39)
where !˜ = !0 + ci R2 = !0 + ci 2 ("=2 + 12 ); and 0 is the initial phase at t = 0. In this case, ◦dW is replaced by dW: This is a linear stochastic di0erential equation in the phase which can be solved exactly; it has the following solution:
1 t (t) = 0 + !t ˜ −
sin(2!0 s + 2 0 ) dWs : (40) 2 0 Therefore, to characterize (t); we need to unfold the stochastic process given by the t integral Z(t) = 0 sin(2!0 s + 0 ) dWs : Let 0 = 0: This integral is a Gaussian process with zero mean and variance given by
t Var(Z(t)) = sin2 (2!0 s + 2( 0 = 0)) ds 0
sin(2!0 t) 1 1 = t− ∼ t !0 2 2
in the long time limit :
(41)
The phase in Eq. (40) can then be evaluated asymptotically in the limit of long time t as Wiener process of the form
(t) = !t ˜ − √ W (t) : 2 2 Using the property Wt lim = 0 with probability1 ; t→∞ t it is clear that
(t) = !˜ : lim t→∞ t
(42)
(43)
(44)
Because of the linearity of Eq. (39), the mean and variance of the phase can be computed directly using the properties of the Ito integral. The average is computed as d (t) = ! ˜ dt − 12 sin(2!0 t + 0 ) dW = ! ˜ dt = !˜ dt ;
(45)
344
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
where the expectation over the Ito integral sin(2!0 t + 0 ) dW always vanishes [27]. The mean of the phase is then given by (t) = !t ˜ + 0 : The variance of the phase is given by the following equation:
2 2 2 sin (2!0 t + 2 0 ) dt ; d = 2 (t)!˜ + 4
(46)
(47)
hence for 0 = 0 and 2 0 = 0; 2 = !˜ 2 t 2 +
sin(4!0 t)
2 : t − 2 8 32!0
(48)
So to a weak noise approximation, the phase is a stochastic Gaussian process with mean and variance given in Eqs. (46) and (48), respectively. The average amplitude vector z(t) can then be approximated in the asymptotic limit of long time as z(t) = R(t) exp(i ) ∼ R(t)exp(i )
2 1 2 2 #("=2 + 1) 2 sin(4!0 t) =
!˜ t + exp (i!t) ˜ − t−
2 8 32!0 #("=2 + 12 ) 2
˜ 1 2 2 #("=2 + 1) !˜ t + t ∼
exp − exp(i!t) ˜ 2 8 #("=2 + 12 ) ˜ exp(i!t) = 'F(!; ˜ ) ˜ ;
(49)
(50) (51)
where ' is given by ' =
#("=2 + 1) #("=2 + 12 )
˜ is the damping function given by and F(!; ˜ ) 2 1
˜ 2 2 ˜ !˜ t + t F(!; ˜ ; t) = exp − : 2 8
(52)
This damping function is Gaussian and similar to the damping function derived by Vlad et al. using the characteristic functional approach [17]. The localization time, !, can be de@ned as the time needed for the damping factor to decrease from the initial ˜ !) = 1=e; which means that 1 (!˜ 2 !2 + ( ˜2 =8)!) = 1; which gives the value to F(!; ˜ ; 2 following in the limit of small
√
2 1 4 2 2 2 − != (
+ 512 ! ˜ −
) ≈ ; (53) !˜ 16!˜ 2 16!˜ 2
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
345
which means that the localization time decreases as the square of the noise amplitude as con@rmed by the numerical solutions. The power spectrum of the average z(t) is given by ∞ 1=2 ˜ P(!; ) = 'F(!; ˜ ; t) exp[i(!˜ − !)t] dt −∞
√
2'+
4 (!˜ − !)2 : (54) = − exp − !˜ 2!˜ 2 512!˜ 2 The principal frequency, !p ; is the frequency for which P(!; ) is maximum, i.e., ˜ which upon substitution in Eq. dP(!; )=d!|!p = 0; which is solved to give !p = !; (54), gives the dependence of the power of the principal peak on the amplitude of the noise √ 2'+
4 : (55) P(!p ; ) = exp − !˜ 512!˜ 2 Eq. (55) is the localization in Fourier space of the principal peak due to the presence of noise and is shown in Fig. 6a. It shows that there is a strong localization of the principal frequency as the amplitude of the noise, , increases. In Fig 3, the power spectra of oscillations in the presence of noise of di0erent amplitudes also exhibit stronger localization as noise is increased. It is also clear from that @gure that the amplitudes of the higher harmonics are also subject to localization. Fig. 6b is a plot of the amplitude of the principal peak of the power of computed ensemble averages as a function of the amplitude of the noise. The trend is similar from that predicted in Eq. (55) and Fig. 6a. 4. Conclusion In this paper, we study numerically and theoretically the temporal localization of the oscillations exhibited by the Selkov model due to the presence of multiplicative white noise. The numerical approach consists of integrating and analysing the stochastic di0erential equations (SDEs) that describe the dynamics. It was found that the computed ensemble averages exhibit localization of the limit cycles for small amplitudes of the noise (0:1–10% of the average value ). L The decays is faster when the amplitude of the noise is higher. To understand the numerical results, our theoretical approach relies on stochastic normal form transformations performed on SDEs and lead to a stochastic complex Landau–Ginzburg equation, which is then used to derive the damping factors that characterize the temporal localization as a function of the amplitude of the noise. This method provides an alternative to the method by Vlad et al. [15] who used characteristic functionals to calculate damping factors. It is found that the decay of the limit cycle is due to the stochastic nature of the phase in the presence noise. The phase is found to have a Gaussian random variable with time-dependent mean and variance. This behavior is similar to the Kubo–Anderson oscillator [32–34] which is a simple linear oscillator with a random frequency. It is worth noting here the di0erence between the e0ects of internal and external *uctuations
346
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
Fig. 6. Localization in Fourier space. (a) Plot of the maximum of principal peak, P(!p ; ) versus the amplitude of the noise as derived from Eq. (55). (b) Localization of the principle amplitude computed numerically from ensemble averages at di0erent noise amplitudes. As the amplitude of noise increases, the amplitude corresponding to the principal frequency falls rapidly. The *uctuations shown can be smoothed out by including more trajectory in the computation of the average.
regarding limit cycles in chemical reactions. Vance and Ross [35] investigated internal *uctuations near a limit cycle using a master equation approach. They found that the variance in phase grows linearly in time. Their work can be repeated using the Langevin equation approach in which the internal noise is additive. In this paper, it was found that weak noise induced a stronger quadratic dependence on time for the variance of the phase. This is an e0ect of multiplicative noise. In the future, we will be looking at the e0ect of colored noise on the Selkov model using the same approach used in this paper. Preliminary studies show that in the presence of weak noise, the averages of concentrations decay completely towards the
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
347
stationary states when the system is driven by Ornstein–Uhlenbeck noise, whereas only a limited decay takes place when the system is driven by noise derived from an exponential distribution function. Acknowledgements This work was supported by a University Research Board grant, American University of Beirut. Appendix A In the study of control and dissipation in chemical engines, Richter et al. modi@ed the Selkov model [24] which consists of three elementary reaction steps k1
AS ; k−1
k2
S + 2P 3P ; k−2
k3
PB:
(A.1)
k−3
In these reactions A and B are kept at @xed concentrations and the intermediates S and P change in time. The kinetic equations describing the time evolution of chemical species in the Selkov model (A.1) are given by dt /p = k2 /s /2p − k−2 /3p − k3 /p + k−3 /b ;
(A.2)
dt /s = k1 /a − k−1 /s − k2 /s /2p + k−2 /3p ;
(A.3)
where dt = d=dt; ki and k−i represent the reaction rate constants for the forward and reverse reactions in the ith step, respectively, and /a ; /b ; /s ; /p are the concentration of species A; B; S; P, respectively. Since it is convenient to work with dimensionless equations, the following reduced variables are introduced t → k3 t ;
(A.4)
X = /p =(k3 =k2 )1=2 ;
Y = /s =(k3 =k2 )1=2 ;
A = (k1 =k3 )(k2 =k3 )1=2 /a ; K = k−2 =k2 ;
B = (k−3 =k3 )(k2 =k3 )1=2 /b ;
= k−1 =k3 :
(A.5) (A.6) (A.7)
With these reduced variables, the evolution equations become Eqs. (1) – (4). Richter et al. [24] studied responses of the Selkov model to periodic variations of the
348
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
experimentally adjusted parameters A and B which appear additively in the kinetic equation (1). In this paper, we take K = 1: The steady states in the deterministic Selkov model are given through the following de@nition: X0 =0 (A.8) d t Z0 = d t Y0 or R(X0 ; Y0 ) = 0 :
(A.9)
Written out explicitly, this equation yields a pair of the algebraic equations B − X0 + X02 Y0 − X03 = 0 ;
(A.10)
A − Y0 − X02 Y0 + X03 = 0 :
(A.11)
There can be one or three real roots of the algebraic equations, depending on the values of the parameters. If B = 0:09 and A = 0:5; one steady state is available. To analyze the stability of the homogeneous steady states in the time domain, we de@ne the following *uctuations of the solution from the steady state (X0 ; Y0 ) and therefore linearize Eq. (1) around them. Thus, with the de@nitions x = X − X0 ;
y = Y − Y0 ;
(A.12)
we obtain dt Zl = R0 Zl :
(A.13)
Here various symbols are de@ned as follows: Xl ; Zl = Yl R0 =
−1 + 2X0 Y0 − 3X02
X02
−2X0 Y0 + 3X02
−X02 −
(A.14) :
(A.15)
The characteristic equation of the matrix R0 is given by !02 + (1 − 2X0 Y0 + 4X02 + )!0 + (1 − 2X0 Y0 + 3X02 ) + X02 = 0 :
(A.16)
For A=0:5 and B=0:09; the product of the roots of Eq. (A.16) is always positive for all positive ; which means they are of the same sign. The sum of the roots of Eq. (A.16) changes sign at = 0:12406: Therefore, the sign of the sum indicates linear stability of Eq. (2). For ¡ 0:1241; the system is in an unstable state where the eigenvalues are complex with positive real parts (). For = 0:1241; we have pure imaginary roots ±i!0 (!0 is a positive real number). For 0:12406 ¡ ¡ 0:1844 the system has a stable focus and for ¿ 0:1844 the system has a node. After the coordinate translation of
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
349
Eq. (A.12), the nonlinear reaction part in Eq. (10) is given by 2X0 xy − (3X0 − Y0 ) + x2 y − x3 ; R= −[2X0 xy − (3X0 − Y0 ) + x2 y − x3 ]
(A.17)
where X0 and Y0 are reference values. They are chosen to be the equilibrium values of the deterministic kinetic equation of the Selkov model. The inverse of the simplest transition matrix (16) which transform the linear part [Eq. (12)] have the following entries: 0 1 : (A.18) T−1 = 1 a − b b The transformed reaction part in the new coordinate system is F(U; V ) L = N G(U; V ) =
( ab + 1)V 3 + (3X0 − Y0 + 2 )V 3 −( a +ab+a+b b2
+
[ (ab+b) b 2 Y0
2 2a b X0 )V
−
(A.19)
− ( 1b )UV 2 − ( 2b X0 )UV
(2a2 +3ab+2a+3b) X0 ]V 2 b2
+
2 ( a+1 b2 )UV
:
(A.20) The stability condition for a limit cycle in a supercritical Hopf bifurcation is computed and it is given by cr =
1 [!0 (Gxxy + Fxxx +xyy +Gyyy ) + Fxx Fxy + Fyy Gyy 16!0 − Gxy Gyy + Fyy Fxy − Fxx Gxx − Gxy Gxx ] :
(A.21)
When cr is negative, the limit cycle is stable, and this is the case of a supercritical Hopf bifurcation: ci =
1 2 2 2 [ − 2Gyy − 2Fxx + 5Fxy Gyy + Fxy Gxx − 5Fyy − 5Gyy Gxx 48!0 2 2 2 − 5Gxx − 2Fxy + 5Gxy Fxx + Fyy Gxy − 5Fyy Fxx − 2Gxy
+ 3!0 (Gxxx + Gxyy − Fyyy − Fxxy )] ;
(A.22)
350
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
where the subscripts on the functions F(U; V ) and G(U; V ) indicate multiple derivatives with respect to the subscript Gxxx = Gxxy = Gxx = Gxy = Fxxx = Fxxy = Fxx = 0 ; 2X0 a 2 Fxyy = − ; Fxy = − ; Fyy :=6X0 − 2Y0 + 4 ; b b b a (a2 + ab + a + b) ; + 1 ; Gyyy = −6 Fyyy = 6 b b2 Gyy :=
2(2a2 + 3ab + 2a + 3b)X0 2(ab + b)Y0 − : b2 b2
In this paper, the adjustable parameters A; B; and L are chosen to be A = 0:5;
B = 0:09
and
L = 0:1 ;
therefore, the deterministic steady states are calculated to be X0 = 0:3396;
and
Y0 = 2:5038 :
Appendix B Eq. (8) is interpreted as an Ito di0erential equation. It can be solved using appropriate numerical methods since the methods for solving ordinary di0erential equations cannot be used for the Wiener process is, almost surely, nowhere di0erentiable [27] and Ito calculus should be used to derive approximation schemes [36]. Let Q = 4Z and let the uniform step size = tn+1 − tn : Time discretization of Eq. (8) yields k Zn+1 = Znk + Rkn I(0) + L0 Rkn I(0; 0) + L1 Rkn I(1; 0) + Qnk I(1)
+ L0 Qnk I(0; 1) + L1 Qnk I(1; 1) + L1 L1 Qnk I(1; 1; 1) + h:o:t: ;
(B.1)
where the superscript k indicates the kth component of the vector and the subscript n denotes the nth time step. If we replace the stochastic Ito integrals I(i; j; k) by their approximations [36] we get k Zn+1 = Znk + Rkn + 12 L0 Rkn 2 + L1 Rkn SZn + Qnk SWn + L0 Qnk ( SWn − SZn )
+ 12 L1 Qnk ((SWn )2 − ) + 16 L1 L1 Qnk ((SWn )3 − 3 SWn ) ;
(B.2)
where the operators L0 and L1 involve di0erential operators and given by L0 =
2
Rk (Z)
2 k=1
2
(B.3)
k=1 j=1
k=1
L1 =
2
@ @ 1 j k + QQ ; k @Z 2 @Z j @Z k
Qk
@ @Z k
(B.4)
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
351
and the processes SWn and SZn will be de@ned later. The derivatives in the expressions of the operators L0 and L1 can be approximated by replacing them with their corresponding @nite di0erences of order 2:0. Platen et al. [36] suggested the following explicit scheme of strong order 1:5 1 k L − Rk (ZL − )]n SZn = Znk + Qnk SWn + √ [Rk (Z) Zn+1 2 +
k L [R (Z + ) + Rk (ZL − ) + 2Rk (Zn )]n 4
+
1 [Qk (ZL + ) + Qk (ZL − ) − 2Qk ]n ( SWn − SZn ) 2
1 + √ [Qk (ZL + ) − Qk (ZL + )]n ((SWn )2 − ) 4 +
1 [Qk (ZLL + ) − Qk (ZLL − ) − Qk (ZL + ) + Qk (ZL − )]n ((SWn )3 − 3 SWn ) ; 12 (B.5)
with the de@nitions
√ k (ZL ± )n = Znk + Rkn ± Qnk ; √ k k k (ZLL ± )n = (ZL + )n ± Q((ZL + )n ) :
(B.6)
We still need to @nd a way to approximate the stochastic process SZn . The mean of SZn is given by tn+1 s tn+1 s (B.7) dWs2 ds = dWs2 ds = 0 : SZn = tn
tn
tn
tn
The variance of SZn is given by (SZn )2 = ( SWn − I(0; 1) )2 = 2 (SWn )2 − 2 SWn I(0; 1) + (I(0; 1) )2
tn+1
tn+1 = 3 − 2 s ds + s2 ds tn
3
= − 2
1 2
2
tn
3 1 3 1 3 = + = 1−1+ : 3 3 3
(B.8)
In Eq. (85), we have used the following identity from the properties of the Ito integrals [27]
I(f)I(g) = f(s; !0)g(s; !0) ds for any f ∈ L2 over any @xed interval : (B.9)
352
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
The covariance or correlation SZn SWn is also calculated SZn SWn = ( SWn − I(0; 1) )SWn = (SWn )2 − I(0; 1) SWn = · − 12 2 = 12 2 : To generate SWn and SZn from uniform random numbers picked from the [0; 1]; we generate two normally distributed random numbers U1 and U2 with and variance and set √ SWn = U1 ;
1 SZn = SWn + √ U2 : 2 3
(B.10) interval mean 0 (B.11) (B.12)
Numerically, this is precisely the way the processes SWn and SZn are evaluated; in addition, the Levy construction, 1 which consists of re@ning the Wiener process at points halfway between those points at which it has already been generated, is also used. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] 1
H. Haken, Synergetics, An Introduction, 2nd Edition, Springer, New York, 1978. K.H. Ho0mann, Z. Phys. B 49 (1982) 245. J. Olarrea, F.J. de la Rubia, Phys. Rev. Lett. 53 (1996) 268. R. Graham, Phys. Rev. A 25 (1982) 3234; R. Graham, T. Tel, Phys. Rev. A 33 (1985) 1322; R. Graham, T. Tel, Phys. Rev. A 35 (1987) 1328. W. Vance, J. Ross, J. Chem. Phys. 91 (1989) 7654; W. Vance, J. Ross, J. Chem. Phys. 105 (1996) 479; W. Vance, J. Ross, J. Chem. Phys. 108 (1997) 2088. C. Meunier, A.D. Verga, J. Stat. Phys. 50 (1988) 345. S. Kabashima, T. Kawakubo, Phys. Lett. A 70 (1979) 375; S. Kabashima, T. Kawakubo, T. Okada, J. Appl. Phys. 50 (1979) 6296. F. Mauricio, S. Valsco, J. Stat. Phys. 43 (1986) 521. W. Horthemke, R. Lefever, Noise-Induced Transitions, Springer, Berlin, 1984. L. Gammaitoni, P. Hanggi, P. Jung, F. Marchesoni, Rev. Mod. Phys. 70 (1998) 223. W. Hohmann, J.M. Schneider, J. Phys. Chem. 100 (1996) 5388. T. Alarc, M. Rubi, J. Chem. 108 (1998) 7367. A.A. Zaikin, J. Garcia-Ojalvo, L. Schimansky-Geier, Phys. Rev. E 60 (1999) 6275. Z. Hou, H. Xin, Phys. Rev. E 60 (1999) 6329. M.O. Vlad, J. Ross, F.W. Schneider, Phys. Rev. E 57 (1998) 4003. P.W. Anderson, Phys. Rev. 109 (1958) 1492; K. Ishi, Prog. Theor. Phys. Suppl. 53 (1973) 77. M.O. Vlad, F.W. Schneider, J. Ross, submitted for publication. J. Guckenheimer, P. Holmes, Nonlinear Oscillations, Dynamical systems, and Bifurcations of Vector Fields, Springer, New York, 1983. G. Nicolis, Introduction to Nonlinear Science, Cambridge University Press, Cambridge, 1995.
Numerically, that the Levy construction consists of re@ning the Wiener process at points halfway between those points at which it has already been generated: W (t + h=2) = 12 [W (t) + W (t + h)] + N(0; h=4). N(0; h=4) is a normally distributed random number of 0 mean and h=4 variance.
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
353
[20] C. Van den Broeck, M. Malek Mansour, F. Baras, J. Stat. Phys. 28 (1982) 557; F. Baras, M. Malek Mansour, C. Van den Broeck, J. Stat. Phys. 28 (1982) 577. [21] E. Knobloch, K.A. Wiesenfeld, J. Stat. Phys. 33 (1983) 611. [22] M. Schumaker, Phys. Lett. A 122 (1987) 317. [23] P.H. Coullet, C. Elphick, E. Tirapegui, Phys. Lett. A 111 (1985) 277. [24] P. Richter, P. Rehmus, J. Ross, Prog. Theor. Phys. 66 (1981) 385. [25] M. Al-Ghoul, PCCP 2 (2000) 3773. [26] R.Z. Khasminsky, Stochastic Stability of Di0erential Equations, Sijtho0 and Noorho0, Alphen naan den Rijn, 1980. [27] L. Arnold, Stochastic Di0erential Equations, Wiley, New York, 1974. [28] S. Wiggins, Introduction to Applied Nonlinear Dynamical Systems and Chaos, Springer, New York, 1990. [29] M. Ipsen, F. Hynne, P.G. Sorensen, preprint chao-dyn=9803006, 1998. [30] J.L. Doob, Stochastic Processes, Wiley, New York, 1972. [31] M.C. Mackey, A. Longtin, A. Lasota, J. Stat. Phys. 60 (1990) 735. [32] P.W. Anderson, J. Phys. Soc. Japan 9 (1954) 316. [33] R. Kubo, Adv. Chem. Phys. 15 (1969) 100. [34] N.G. Van Kampen, Phys. Rep. C 24 (1976) 171. [35] W. Vance, J. Ross, J. Chem. Phys. 91 (1989) 7654. [36] P. Kloeden, E. Platen, Numerical Solution of Stochastic Di0erential Equations, Springer, Berlin, 1992.
www.elsevier.com/locate/physa
Temporal localization of limit cycles in a noise-driven chemical oscillator Mazen Al-Ghoul Chemistry Department and Center for Advanced Mathematical Sciences, American University of Beirut, Beirut, Lebanon Received 26 July 2001
Abstract We study the in*uence of multiplicative random external perturbation on the Selkov model exhibiting limit cycles. The dynamics are described by stochastic di0erential equations. The numerically computed ensemble averages of the concentrations exhibit temporal localization of the limit cycles for small amplitudes of the noise. The localization time is smaller when the amplitude of the noise is higher. To understand this behavior, a stochastic complex Ginzburg–Landau equation is derived using normal form transformations on the stochastic evolution equations. In the case of weak noise, the equations for the radial and the phase components of the amplitude separate and can be solved exactly. The solution is used to compute the damping factors, which were found to be Gaussian. We also show that asymptotically, the radial part *uctuates near the average value and the phase is found to be a Gaussian random variable with time-dependent mean and variance and may be responsible for the temporal localization of the limit cycle. c 2002 Elsevier Science B.V. All rights reserved. PACS: 05.10.−a; 05.40.−a; 05.45.−a; 82.20.−w Keywords: Multiplicative noise; Chemical kinetics
1. Introduction The in*uence of internal and external noise on bifurcating dynamical systems has been the subject of several theoretical and experimental studies [1–8]. Environmental *uctuations such as temperature on these systems may alter their behavior compared to those systems without external noise, giving rise, for example, to noise-induced transitions [9], stochastic resonance [10 –12], nonequilibrium @rst-order phase transitions [13], and noise-induced oscillations [14]. This kind of external noise may enter multiplicatively in the macroscopic equations in contrast to internal thermal noise which appears c 2002 Elsevier Science B.V. All rights reserved. 0378-4371/02/$ - see front matter PII: S 0 3 7 8 - 4 3 7 1 ( 0 1 ) 0 0 6 1 2 - 4
332
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
additively in the governing equations that describe the dynamics of the system. Those equations are, in general, nonlinear multidimensional stochastic di0erential equations. Recently, Vlad et al. considered the e0ect of multiplicative noise on a reaction–di0usion system [15] and showed, using the characteristic functional approach, that *uctuations induce a decay or a temporal localization in the averages of the oscillations which can be considered as a temporal analog of Anderson localization of the wavefunction in a disordered system [16]. Vlad later extended his formalism to study localization of limit cycles in the presence of fractal noise [17]. He de@ned localization times, derived damping factors of the oscillations, and found that in the case of white noise, those factors are Gaussian. In this paper, we use a di0erent approach, that of stochastic differential equations, from Vlad’s characteristic functionals to study and characterize the temporal localization of limit cycles in the presence of multiplicative white noise. With this approach, we will be able to obtain similar expressions for the localization time and for the damping factors. Our method is parallel to the study and characterization of the local behavior of deterministic multidimensional nonlinear systems, where one may use the center manifold theorem to reduce the dimensionality of the system, and the normal form method to eliminate nonlinearities by means of successive nonlinear transformations [18]. The obtained normal form equations have a universal character and they are derived in several di0erent areas of physics [19]. The extension of the normal form reduction scheme to stochastic di0erential equations requires the restrictive physical assumptions of weak noise. Several studies have implemented the “stochastic” center manifold to reduce stochastic systems to a more simpli@ed form. Using this approach, one can obtain a reduced Fokker–Planck equation [20,21] that can be written in a normal form by means of local nonlinear transformations. On the other hand, one may apply the normal form transformations directly to the stochastic di0erential equations [22]. For a Hopf bifurcation with external random parameters, Coullet et al. [23] derived a general normal form which contains new “stochastic resonance terms” due to the stochastic nature of the equations. Both approaches yield similar results but the former is relatively easier and can be used directly to interpret the numerical solutions of the stochastic di0erential equations as will be shown in this paper. In this work, we study the dynamics in a model of a limit cycle with imposed multiplicative random forcing. We use the modi@ed Selkov model [24], which contains four adjustable parameters, and we vary randomly one dimensionless parameter, ; which appears multiplicatively in the corresponding chemical kinetic equations. The resulting stochastic di0erential equations are solved numerically in the case of weak noise; the solutions and their ensemble averages for di0erent realizations of the noise are computed and it was found that these ensemble averages exhibit damping of the limit cycle. To understand this damping, we derive a stochastic Landau–Ginzburg (LG) equations using nonlinear transformations of coordinates. Under the white weak noise assumption, the stochastic LG equations separate into two independent equations for the radial and phase components of the complex amplitude. This decay in the oscillations may due to the stochastic nature of the phase which we are going to characterize in this paper by calculating the damping factors in real and Fourier spaces. In Section 2, we present and discuss the numerical solutions of the stochastic di0erential equations. In Section 3, we use normal form transformations to derive a stochastic
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
333
Ginzburg–Landau equation which is solved exactly in the case of weak noise. Then localization times and damping factors are easily de@ned. In the conclusion, we discuss possible extensions of this work to colored noise. In Appendix A, we perform a linear stability analysis of the deterministic equations and de@ne di0erent quantities and expressions. Appendix B contains the numerical scheme used to integrate the stochastic di0erential equations. 2. Weak random forcing The modi@ed Selkov model is de@ned through the following kinetic equations (see Appendix A for de@nitions of various quantities): dt X = B − X + X 2 Y − KX 3 ; dt Y = A − Y − X 2 Y + KX 3 ;
(1)
which can be put in the following more compact form: dt Z = R ; where
Z=
and
(3)
Y
R=
X
(2)
B − X + X 2 Y − KX 3 A − Y − X 2 Y + KX 3
:
(4)
We now assume that the external variations are given in terms of the variation of the parameter : This parameter is assumed to vary randomly according to the following equation: = L + dt W (t) ;
(5)
where L is an average value, dt = d=dt; and dt W (t) is a Gaussian white noise term which we assume to be the time derivative of the Wiener process W (t); is the amplitude of the random noise. Therefore, the properties of the noise are given by dt W (t) = 0 ;
(6)
dt W (t1 )dt W (t2 ) = 2 (t2 − t1 ) :
(7)
If we substitute Eq. (5) into Eq. (2), we obtain the following Ito stochastic di0erential equation: dZ = R dt + Z dW (t) ; where is given by 0 0 = : 0 1
(8)
(9)
334
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
In the deterministic case, i.e., = 0; the system described by Eq. (2) undergoes a Hopf bifurcation when A = 0:5; B = 0:09; and L 6 0:1241: In this paper, the numerical results were obtained by taking L = 0:1: In the case of weak noise,
0:1; we show in the next section that the structure of the phase diagram is not altered and Eq. (8) still admits a set of stationary points very close to the deterministic stationary points. This phase diagram comprises an oscillations region and a stable region separated by the Hopf bifurcation line. The noise is chosen such that con@nes the system inside the oscillations region and crossing to the stable region is an extremely rare event. We integrated numerically Eq. (8) for several small noise amplitudes ranging from
≈ 0:0001 to 0:03. In every case, we compute and we record the times its value crosses back and forth the oscillations region. For in the interval [0:0001; 0:03] no crossing to the stable regions was observed. We describe the numerical scheme in Appendix B. For every realization of the noise we obtained a trajectory. We performed at least 50 realizations for every case we studied in order to perform an ensemble average; this number was suMcient in most of the cases studied. All cases in Fig. 1a–f lie entirely in the oscillations region; we see that for an amplitude lower than 0:01, the noise does not a0ect much a given realization, and the trajectories resemble those of the deterministic case [25]. It is important to note here, that the noise did not delay the bifurcation and a stable limit cycle continue to exist around stationary points that are very close to stationary point of the deterministic system (X0 , and Y0 from Eq. (64)) [25]. However, as increases we can see that the behavior begins to change towards what may be a chaotic regime. As for the average over noise realizations of the concentrations, for any small value we take for ; we can see in Fig. 2a–f that the oscillations decay toward a region near the deterministic values of the steady state. The higher the amplitude of the noise ; the faster the decay. The computed spectrum in Fig. 3, computed for the average shows that an increase in the noise amplitude results in reduction of the number of harmonics present initially in the deterministic system. In addition, the amplitude of the principal frequency decreases and broadens towards the @rst harmonic. Above = 0:03; the peaks are broad enough to demonstrate that the ensemble average has a chaotic character: 3. Normal forms and localization In this section, we present an approximate analysis of the kinetic equations of the stochastically driven Selkov model, where we use an appropriate linearization scheme of the stochastic di0erential equations and then invoke nonlinear coordinate transformations, similar to normal form methods in the theory of deterministic dynamical systems, to bring the equations into an easier solvable form. We show that near the Hopf bifurcation and for small noise amplitudes, the obtained stochastic normal form equation can be used to describe qualitatively the behavior of the system. Because we need to make a series of operations and coordinate transformation, it is better to interpret the stochastic di0erential equation (8) as a Stratonovitch di0erential equation, where we can make use of the easier rules of the ordinary calculus. To do
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
335
Fig. 1. The variation of the concentration of species X in time for one realization of the noise with di0erent intensities : (a) = 0:0001, (b) = 0:005, (c) = 0:009, (d) = 0:01, (e) = 0:02, (f) = 0:03:
336
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
Fig. 2. The variation of the ensemble average of concentration of species X in time for di0erent noise intensities : (a) = 0:0001, (b) = 0:005, (c) = 0:009, (d) = 0:01, (e) = 0:02, (f) = 0:03:
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
337
Fig. 3. Power spectra of oscillations in the presence of noise of di0erent amplitudes: = 1:0 × 10−4 dashed line, = 1:0 × 10−2 solid line.
so, we correct the drift coeMcient R; which is calculated in Appendix A, as follows: ˜ = R − 1 2 (Z) @ (Z) ; R 2 @Z T 0 0
2 0 ; =R− 2 Y 0 1
2 =R− 2 =R−
0
0
0
1
X
Y
2 (Z) : 2
; (10)
The Stratonovitch stochastic di0erential equations corresponding to the drift coeMcient ˜ is given by R ˜ + (Z) ◦ dW (t) ; dt Z = R
(11)
where ◦dW indicates a Stratonovitch type of measure. For the weak noise limit, i.e.,
is small, we can linearize Eq. (11) around the state Z0 = (X0 ; Y0 ) which is an unstable ˜ steady state [26,27] of the nonstochastic system (2). The nonlinear drift coeMcient R is linearized by replacing it by the following Jacobian matrix D evaluated at (X0 ; Y0 ) : 0 X02 −1 + 2X0 Y0 − 3X02 ˜ i
2 0 @R : (D)ij = − = @Zj 2 0 1 −2X0 Y0 + 3X02 −X02 − L (X0 ; Y0 )
(12)
338
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
The di0usion term is already linear and preserves its form. The linearized system assumes the following form: dt Zl = DZl + (Zl ) ◦ dW :
(13)
Since the equation is already linear in the noise term, and because the normal form depends only on the number and nature of the eigenvalues of the linearized problem, we suggest to reduce Eq. (11) to a normal form by a series of nonlinear transformations of coordinates exactly like the deterministic case [28] and leave the noise term intact. The resulting equation corrects for the linearized equation (13). Let = ( ) + i!( ) and L = ( ) − i!( ) be two complex eigenvalues of the matrix D near the Hopf bifurcation (( ) and !( ) are real), and let ’; ’L be the corresponding complex eigenvectors. For small ; ( ) and !( ) are given by ( ) = 0 − 14 2 + O( 4 ) ; !( ) = !0 +
(14)
0 2
+ O( 4 ) ; 4!0
(15)
where 0 and !0 are, respectively, the real and imaginary part of the eigenvalues of the problem without noise. It is clear from Eqs. (14) and (15) that noise can delay the bifurcation by shifting 0 to ( ); however, in the limit of weak noise, ( ) ≈ 0 ; and !( ) ≈ !0 ; and therefore the structure of the phase diagram is not altered near the Hopf bifurcation point in the limit of weak noise. The next step is to make a coordinates transformation from the real and imaginary part of one of the eigenvectors to bring Eq. (11) into its corresponding normal form. More precisely, if for example, ’ = (a1 + ib1 ; a2 + ib2 ) we can construct the following transition matrix T a b b2 a2 = ; (16) T= 1 0 b1 a1 which acts on Z as J = T−1 Z
(17)
to give a new set of coordinates. J is the vector (U; V )T . In this new coordinate system, Eq. (11) has the following form: dJ = [T−1 DTJ + T−1 N(TJ)] dt + ( T−1 TJ) ◦ dW L ◦ dW ; L + N) L dt + J = (DJ where
N=
N1 (U 2 ; UV; U 2 V; U 3 ) N2 (U 2 ; UV; U 2 V; U 3 )
(18)
;
(19)
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
L = T−1 DT = D L = T
−1
T =
−!
!
1
0
− ba
0
339
;
(20)
(21)
and NL = T−1 N(TJ): Using the normal form theorem [28], a couple of standard nonlinear transformations (near-identity transformations) eliminate all quadratic and cubic nonlinearities in NL except for a couple of cubic nonresonant terms. We @nally end up with the following normal form equation in complex form to cubic order a dz = [( + i!)z + (cr + ici )|z|2 z] dt + 1 − Re(z) + O( 3 ) ◦ dW ; b (22) where z =U +iV . cr and ci can be computed analytically [29]. For a supercritical Hopf bifurcation, cr should be negative. For the case of weak noise, we will drop the higher order terms and the nonlinear terms of O( 3 ): Also, the term a=b; when computed numerically, depends weakly on so that in the weak noise limit, it is assumed that it is constant. Therefore, Eq. (22) is written as dz = [( + i!)z + (cr + ici )|z|2 z] dt + Re(z) ◦ dW ;
(23)
where = (1−a=b): Now we will represent Eq. (23) in polar form. Let z =R exp(i ); then Eq. (23) is decomposed into the following radial and phase components: dR = (R + cr R3 ) dt + R cos2 ( ) ◦ dW ;
(24)
d = (! + ci R2 ) dt − 12 sin(2 ) ◦ dW :
(25)
Motivated by results of Fig. 3, that in the weak noise limit the perturbed limit cycle oscillates with a frequency that is close to the frequency of the unperturbed limit cycle, !0 ; and in order to make some analytical progress we propose the following mean-@eld approximation to simplify further Eqs. (24) and (25): in the weak noise limit,
0:1, the phase appearing in the right-hand sides of those equations is replaced by an average quantity → ≈ !0 t + 0 , which is the expected behavior in the limit of zero noise. This approximation will later provide a substantial agreement with the damping factor of Vlad et al. [15,17], and with the numerical solutions of the initial stochastic di0erential equations (8). Therefore, Eqs. (24) and (25) are approximated by dR = (R + cr R3 ) dt + R cos2 (!0 t + 0 ) ◦ dW ;
(26)
d = (!0 + ci R2 ) dt − 12 sin(2!0 t + 2 0 ) ◦ dW :
(27)
To test the approximation, we solved Eqs. (24) and (25) and Eqs. (26) and (27) numerically for the same parameters and same realization of the Wiener process. We
340
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
Fig. 4. Solution of the radial part, R(t); of the amplitude equations from Eqs. (24) and (25) [higher curve], and Eqs. (26) and (27) [lower curve] for the same realization of the Wiener process Wt . = 0:0005: Notice the convergence of the two curves after some time.
also used the numerical scheme described in Appendix B. In Fig. 4, we see how the two solutions for the radial component of the amplitude of the two versions of the amplitude equations converge to each other after some time in the weak noise limit. As the noise increases to about = 0:05, the solutions do not match very well. 3.1. The radial part of the amplitude Eq. (26) is a nonlinear stochastic di0erential equation. The change of variable Y = R−2 transform Eq. (26) to a linear stochastic di0erential equation which can then be solved explicitly to give t exp( t + 0 cos2 (!0 s + 0 ) dWs ) R(t) = ; (28) t t 1 + 2R20 0 exp(2 s + 2 0 cos2 (!0 s + 0 ) dWs ) ds t where = =|cr |: The stochastic integral I(t) = 0 cos2 (!0 s + 0 ) dWs is a Gaussian stochastic process with zero mean and variance given by [27]
t Var(I) = cos4 (!0 s + 0 ) ds 0
3 2 sin(!0 t + 0 ) cos3 (!0 t + 0 ) + 3 sin[2(!0 t + 0 )] = t+ : 8 8!0
(29)
In long times, the variance in Eq. (29) is approximated as Var(I) = 38 t because the second term on the RHS of Eq. (29) is bounded. Therefore, by Doob’s theorem [30], asymptotically the stochastic process I(t) resembles the Wiener process Wt . If we
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
341
make the transformation ! = 38 t or rescale the amplitude of the noise ; then we can write R(t) in Eq. (28) as ˜ t) R0 exp( t + W ; t ˜ s ) ds 1 + 2R20 0 exp(2 s + 2 W
R(t) =
(30)
where = =|cr |; ˜ = 83 ; Wt is a Wiener process and R0 is an initial condition at t = 0. In the case of zero noise, i.e., = 0; R(t) is given by R0 exp( t)
R(t) =
1 + (R20 =)[exp(2 t) − 1]
→
as t → ∞
(31)
as expected. We anticipate that the asymptotic behavior of the stochastic process R(t) for small to exhibit small *uctuations around a certain average. To prove this, we need to study the statistics of the stochastic process R(t) in Eq. t (30). Therefore, we ˜ s ) ds @rst need to understand the time dependence of the integral 0 exp(2 s + 2 W which appears in the denominator of Eq. (30). The idea is to show that the process R(t) very closely to an average value; this is the case only if the integral t *uctuates ˜ s ) vary at the same rate so that ratio ˜ s ) ds and the exp(2 s + 2 W exp(2 s + 2
W 0 (30) *uctuates very close to a constant value. Evidently, the stochastic process R(t) is not Gaussian. However, if the moments higher than the second of this process are small, then we can safely approximate it with a Gaussian. We notice that the process R(t) in Eq. (30) is a solution of the following Ito stochastic di0erential equation: ˜ dW ; dR = ( R − R3 ) d! + R
(32)
2 where = − 12 ˜ and ! = |cr |t: The ensemble average R(t) is easily obtained from averaging over the following stationary distribution function P(R) derived from the Fokker–Planck equation with natural boundary conditions: 2
P(R) = KR" exp(−R2 = ˜ ) ;
(33)
where K is a normalization constant given by −"−1
2 ˜ K= #("=2 + 12 )
(34)
and the exponent " is given by "( ) =
2 3( ) −2= 2 −3: 2 4 cr ( )(1 − a( )=b( ))2
˜
(35)
For P(R) to be a probability density [31], it has to be integrable over R+ , i.e., "+1 ¿ 0. For A = 0:5; B = 0:09; L = 0:1; (" + 1) is positive at a noise amplitude 6 0:6642: Therefore, it is safe to assume that a probability density of the form P(R) in Eq. (33) exists whenever the noise is weak. Hence, all the moments of the stochastic process
342
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
Fig. 5. The process R(t) and its ensemble average R(t) (thicker line) for di0erent noise amplitudes: (a)
= 0:01, (b) = 0:05, (c) = 0:1, (d) = 0:3:
R can be computed from the following Riemann integral: Rm =
0
∞
+ "=2 + 12 ) : #("=2 + 12 )
m #(m=2
Rm P(R) dR = ˜
(36)
In particular, the average and second moment are given respectively by #("=2 + 1) R = ˜ ; #("=2 + 12 )
(37)
2 R2 = 12 ˜ (" + 1) :
(38)
In the weak noise limit, as → 0, the moments in Eq. (36) tend to zero very quickly. Therefore, it is a good approximation to assume that asymptotically R(t) is a Gaussian random variable characterized by Eqs. (37) and (38). Moreover, the second moment (38) is small for weak noise. The numerical evaluation of R(t) in Eq. (30) is shown in Fig. 5a and b where the amplitude of the noise is considered small. In this case,
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
343
the realization of R(t) *uctuates closely to the average R(t) in long times. As becomes large, the *uctuations grow around the ensemble average [Fig. 5c and d], in a similar way to that in Fig. 5a and b, where the average clearly lies in between the realizations of the process R(t): When takes large values, no stationary solution of the corresponding Fokker–Planck exists, and high amplitude *uctuations in the realizations of R(t) starts to appear [Fig. 5d]. 3.2. The phase di/usion and damping functions The equation for phase (25) will then be given by d = (!0 + ci R2 ) dt − 12 sin(2!0 t + 2 0 ) ◦ dW ; = !˜ dt − 12 sin(2!0 t + 2 0 ) dW ;
(39)
where !˜ = !0 + ci R2 = !0 + ci 2 ("=2 + 12 ); and 0 is the initial phase at t = 0. In this case, ◦dW is replaced by dW: This is a linear stochastic di0erential equation in the phase which can be solved exactly; it has the following solution:
1 t (t) = 0 + !t ˜ −
sin(2!0 s + 2 0 ) dWs : (40) 2 0 Therefore, to characterize (t); we need to unfold the stochastic process given by the t integral Z(t) = 0 sin(2!0 s + 0 ) dWs : Let 0 = 0: This integral is a Gaussian process with zero mean and variance given by
t Var(Z(t)) = sin2 (2!0 s + 2( 0 = 0)) ds 0
sin(2!0 t) 1 1 = t− ∼ t !0 2 2
in the long time limit :
(41)
The phase in Eq. (40) can then be evaluated asymptotically in the limit of long time t as Wiener process of the form
(t) = !t ˜ − √ W (t) : 2 2 Using the property Wt lim = 0 with probability1 ; t→∞ t it is clear that
(t) = !˜ : lim t→∞ t
(42)
(43)
(44)
Because of the linearity of Eq. (39), the mean and variance of the phase can be computed directly using the properties of the Ito integral. The average is computed as d (t) = ! ˜ dt − 12 sin(2!0 t + 0 ) dW = ! ˜ dt = !˜ dt ;
(45)
344
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
where the expectation over the Ito integral sin(2!0 t + 0 ) dW always vanishes [27]. The mean of the phase is then given by (t) = !t ˜ + 0 : The variance of the phase is given by the following equation:
2 2 2 sin (2!0 t + 2 0 ) dt ; d = 2 (t)!˜ + 4
(46)
(47)
hence for 0 = 0 and 2 0 = 0; 2 = !˜ 2 t 2 +
sin(4!0 t)
2 : t − 2 8 32!0
(48)
So to a weak noise approximation, the phase is a stochastic Gaussian process with mean and variance given in Eqs. (46) and (48), respectively. The average amplitude vector z(t) can then be approximated in the asymptotic limit of long time as z(t) = R(t) exp(i ) ∼ R(t)exp(i )
2 1 2 2 #("=2 + 1) 2 sin(4!0 t) =
!˜ t + exp (i!t) ˜ − t−
2 8 32!0 #("=2 + 12 ) 2
˜ 1 2 2 #("=2 + 1) !˜ t + t ∼
exp − exp(i!t) ˜ 2 8 #("=2 + 12 ) ˜ exp(i!t) = 'F(!; ˜ ) ˜ ;
(49)
(50) (51)
where ' is given by ' =
#("=2 + 1) #("=2 + 12 )
˜ is the damping function given by and F(!; ˜ ) 2 1
˜ 2 2 ˜ !˜ t + t F(!; ˜ ; t) = exp − : 2 8
(52)
This damping function is Gaussian and similar to the damping function derived by Vlad et al. using the characteristic functional approach [17]. The localization time, !, can be de@ned as the time needed for the damping factor to decrease from the initial ˜ !) = 1=e; which means that 1 (!˜ 2 !2 + ( ˜2 =8)!) = 1; which gives the value to F(!; ˜ ; 2 following in the limit of small
√
2 1 4 2 2 2 − != (
+ 512 ! ˜ −
) ≈ ; (53) !˜ 16!˜ 2 16!˜ 2
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
345
which means that the localization time decreases as the square of the noise amplitude as con@rmed by the numerical solutions. The power spectrum of the average z(t) is given by ∞ 1=2 ˜ P(!; ) = 'F(!; ˜ ; t) exp[i(!˜ − !)t] dt −∞
√
2'+
4 (!˜ − !)2 : (54) = − exp − !˜ 2!˜ 2 512!˜ 2 The principal frequency, !p ; is the frequency for which P(!; ) is maximum, i.e., ˜ which upon substitution in Eq. dP(!; )=d!|!p = 0; which is solved to give !p = !; (54), gives the dependence of the power of the principal peak on the amplitude of the noise √ 2'+
4 : (55) P(!p ; ) = exp − !˜ 512!˜ 2 Eq. (55) is the localization in Fourier space of the principal peak due to the presence of noise and is shown in Fig. 6a. It shows that there is a strong localization of the principal frequency as the amplitude of the noise, , increases. In Fig 3, the power spectra of oscillations in the presence of noise of di0erent amplitudes also exhibit stronger localization as noise is increased. It is also clear from that @gure that the amplitudes of the higher harmonics are also subject to localization. Fig. 6b is a plot of the amplitude of the principal peak of the power of computed ensemble averages as a function of the amplitude of the noise. The trend is similar from that predicted in Eq. (55) and Fig. 6a. 4. Conclusion In this paper, we study numerically and theoretically the temporal localization of the oscillations exhibited by the Selkov model due to the presence of multiplicative white noise. The numerical approach consists of integrating and analysing the stochastic di0erential equations (SDEs) that describe the dynamics. It was found that the computed ensemble averages exhibit localization of the limit cycles for small amplitudes of the noise (0:1–10% of the average value ). L The decays is faster when the amplitude of the noise is higher. To understand the numerical results, our theoretical approach relies on stochastic normal form transformations performed on SDEs and lead to a stochastic complex Landau–Ginzburg equation, which is then used to derive the damping factors that characterize the temporal localization as a function of the amplitude of the noise. This method provides an alternative to the method by Vlad et al. [15] who used characteristic functionals to calculate damping factors. It is found that the decay of the limit cycle is due to the stochastic nature of the phase in the presence noise. The phase is found to have a Gaussian random variable with time-dependent mean and variance. This behavior is similar to the Kubo–Anderson oscillator [32–34] which is a simple linear oscillator with a random frequency. It is worth noting here the di0erence between the e0ects of internal and external *uctuations
346
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
Fig. 6. Localization in Fourier space. (a) Plot of the maximum of principal peak, P(!p ; ) versus the amplitude of the noise as derived from Eq. (55). (b) Localization of the principle amplitude computed numerically from ensemble averages at di0erent noise amplitudes. As the amplitude of noise increases, the amplitude corresponding to the principal frequency falls rapidly. The *uctuations shown can be smoothed out by including more trajectory in the computation of the average.
regarding limit cycles in chemical reactions. Vance and Ross [35] investigated internal *uctuations near a limit cycle using a master equation approach. They found that the variance in phase grows linearly in time. Their work can be repeated using the Langevin equation approach in which the internal noise is additive. In this paper, it was found that weak noise induced a stronger quadratic dependence on time for the variance of the phase. This is an e0ect of multiplicative noise. In the future, we will be looking at the e0ect of colored noise on the Selkov model using the same approach used in this paper. Preliminary studies show that in the presence of weak noise, the averages of concentrations decay completely towards the
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
347
stationary states when the system is driven by Ornstein–Uhlenbeck noise, whereas only a limited decay takes place when the system is driven by noise derived from an exponential distribution function. Acknowledgements This work was supported by a University Research Board grant, American University of Beirut. Appendix A In the study of control and dissipation in chemical engines, Richter et al. modi@ed the Selkov model [24] which consists of three elementary reaction steps k1
AS ; k−1
k2
S + 2P 3P ; k−2
k3
PB:
(A.1)
k−3
In these reactions A and B are kept at @xed concentrations and the intermediates S and P change in time. The kinetic equations describing the time evolution of chemical species in the Selkov model (A.1) are given by dt /p = k2 /s /2p − k−2 /3p − k3 /p + k−3 /b ;
(A.2)
dt /s = k1 /a − k−1 /s − k2 /s /2p + k−2 /3p ;
(A.3)
where dt = d=dt; ki and k−i represent the reaction rate constants for the forward and reverse reactions in the ith step, respectively, and /a ; /b ; /s ; /p are the concentration of species A; B; S; P, respectively. Since it is convenient to work with dimensionless equations, the following reduced variables are introduced t → k3 t ;
(A.4)
X = /p =(k3 =k2 )1=2 ;
Y = /s =(k3 =k2 )1=2 ;
A = (k1 =k3 )(k2 =k3 )1=2 /a ; K = k−2 =k2 ;
B = (k−3 =k3 )(k2 =k3 )1=2 /b ;
= k−1 =k3 :
(A.5) (A.6) (A.7)
With these reduced variables, the evolution equations become Eqs. (1) – (4). Richter et al. [24] studied responses of the Selkov model to periodic variations of the
348
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
experimentally adjusted parameters A and B which appear additively in the kinetic equation (1). In this paper, we take K = 1: The steady states in the deterministic Selkov model are given through the following de@nition: X0 =0 (A.8) d t Z0 = d t Y0 or R(X0 ; Y0 ) = 0 :
(A.9)
Written out explicitly, this equation yields a pair of the algebraic equations B − X0 + X02 Y0 − X03 = 0 ;
(A.10)
A − Y0 − X02 Y0 + X03 = 0 :
(A.11)
There can be one or three real roots of the algebraic equations, depending on the values of the parameters. If B = 0:09 and A = 0:5; one steady state is available. To analyze the stability of the homogeneous steady states in the time domain, we de@ne the following *uctuations of the solution from the steady state (X0 ; Y0 ) and therefore linearize Eq. (1) around them. Thus, with the de@nitions x = X − X0 ;
y = Y − Y0 ;
(A.12)
we obtain dt Zl = R0 Zl :
(A.13)
Here various symbols are de@ned as follows: Xl ; Zl = Yl R0 =
−1 + 2X0 Y0 − 3X02
X02
−2X0 Y0 + 3X02
−X02 −
(A.14) :
(A.15)
The characteristic equation of the matrix R0 is given by !02 + (1 − 2X0 Y0 + 4X02 + )!0 + (1 − 2X0 Y0 + 3X02 ) + X02 = 0 :
(A.16)
For A=0:5 and B=0:09; the product of the roots of Eq. (A.16) is always positive for all positive ; which means they are of the same sign. The sum of the roots of Eq. (A.16) changes sign at = 0:12406: Therefore, the sign of the sum indicates linear stability of Eq. (2). For ¡ 0:1241; the system is in an unstable state where the eigenvalues are complex with positive real parts (). For = 0:1241; we have pure imaginary roots ±i!0 (!0 is a positive real number). For 0:12406 ¡ ¡ 0:1844 the system has a stable focus and for ¿ 0:1844 the system has a node. After the coordinate translation of
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
349
Eq. (A.12), the nonlinear reaction part in Eq. (10) is given by 2X0 xy − (3X0 − Y0 ) + x2 y − x3 ; R= −[2X0 xy − (3X0 − Y0 ) + x2 y − x3 ]
(A.17)
where X0 and Y0 are reference values. They are chosen to be the equilibrium values of the deterministic kinetic equation of the Selkov model. The inverse of the simplest transition matrix (16) which transform the linear part [Eq. (12)] have the following entries: 0 1 : (A.18) T−1 = 1 a − b b The transformed reaction part in the new coordinate system is F(U; V ) L = N G(U; V ) =
( ab + 1)V 3 + (3X0 − Y0 + 2 )V 3 −( a +ab+a+b b2
+
[ (ab+b) b 2 Y0
2 2a b X0 )V
−
(A.19)
− ( 1b )UV 2 − ( 2b X0 )UV
(2a2 +3ab+2a+3b) X0 ]V 2 b2
+
2 ( a+1 b2 )UV
:
(A.20) The stability condition for a limit cycle in a supercritical Hopf bifurcation is computed and it is given by cr =
1 [!0 (Gxxy + Fxxx +xyy +Gyyy ) + Fxx Fxy + Fyy Gyy 16!0 − Gxy Gyy + Fyy Fxy − Fxx Gxx − Gxy Gxx ] :
(A.21)
When cr is negative, the limit cycle is stable, and this is the case of a supercritical Hopf bifurcation: ci =
1 2 2 2 [ − 2Gyy − 2Fxx + 5Fxy Gyy + Fxy Gxx − 5Fyy − 5Gyy Gxx 48!0 2 2 2 − 5Gxx − 2Fxy + 5Gxy Fxx + Fyy Gxy − 5Fyy Fxx − 2Gxy
+ 3!0 (Gxxx + Gxyy − Fyyy − Fxxy )] ;
(A.22)
350
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
where the subscripts on the functions F(U; V ) and G(U; V ) indicate multiple derivatives with respect to the subscript Gxxx = Gxxy = Gxx = Gxy = Fxxx = Fxxy = Fxx = 0 ; 2X0 a 2 Fxyy = − ; Fxy = − ; Fyy :=6X0 − 2Y0 + 4 ; b b b a (a2 + ab + a + b) ; + 1 ; Gyyy = −6 Fyyy = 6 b b2 Gyy :=
2(2a2 + 3ab + 2a + 3b)X0 2(ab + b)Y0 − : b2 b2
In this paper, the adjustable parameters A; B; and L are chosen to be A = 0:5;
B = 0:09
and
L = 0:1 ;
therefore, the deterministic steady states are calculated to be X0 = 0:3396;
and
Y0 = 2:5038 :
Appendix B Eq. (8) is interpreted as an Ito di0erential equation. It can be solved using appropriate numerical methods since the methods for solving ordinary di0erential equations cannot be used for the Wiener process is, almost surely, nowhere di0erentiable [27] and Ito calculus should be used to derive approximation schemes [36]. Let Q = 4Z and let the uniform step size = tn+1 − tn : Time discretization of Eq. (8) yields k Zn+1 = Znk + Rkn I(0) + L0 Rkn I(0; 0) + L1 Rkn I(1; 0) + Qnk I(1)
+ L0 Qnk I(0; 1) + L1 Qnk I(1; 1) + L1 L1 Qnk I(1; 1; 1) + h:o:t: ;
(B.1)
where the superscript k indicates the kth component of the vector and the subscript n denotes the nth time step. If we replace the stochastic Ito integrals I(i; j; k) by their approximations [36] we get k Zn+1 = Znk + Rkn + 12 L0 Rkn 2 + L1 Rkn SZn + Qnk SWn + L0 Qnk ( SWn − SZn )
+ 12 L1 Qnk ((SWn )2 − ) + 16 L1 L1 Qnk ((SWn )3 − 3 SWn ) ;
(B.2)
where the operators L0 and L1 involve di0erential operators and given by L0 =
2
Rk (Z)
2 k=1
2
(B.3)
k=1 j=1
k=1
L1 =
2
@ @ 1 j k + QQ ; k @Z 2 @Z j @Z k
Qk
@ @Z k
(B.4)
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
351
and the processes SWn and SZn will be de@ned later. The derivatives in the expressions of the operators L0 and L1 can be approximated by replacing them with their corresponding @nite di0erences of order 2:0. Platen et al. [36] suggested the following explicit scheme of strong order 1:5 1 k L − Rk (ZL − )]n SZn = Znk + Qnk SWn + √ [Rk (Z) Zn+1 2 +
k L [R (Z + ) + Rk (ZL − ) + 2Rk (Zn )]n 4
+
1 [Qk (ZL + ) + Qk (ZL − ) − 2Qk ]n ( SWn − SZn ) 2
1 + √ [Qk (ZL + ) − Qk (ZL + )]n ((SWn )2 − ) 4 +
1 [Qk (ZLL + ) − Qk (ZLL − ) − Qk (ZL + ) + Qk (ZL − )]n ((SWn )3 − 3 SWn ) ; 12 (B.5)
with the de@nitions
√ k (ZL ± )n = Znk + Rkn ± Qnk ; √ k k k (ZLL ± )n = (ZL + )n ± Q((ZL + )n ) :
(B.6)
We still need to @nd a way to approximate the stochastic process SZn . The mean of SZn is given by tn+1 s tn+1 s (B.7) dWs2 ds = dWs2 ds = 0 : SZn = tn
tn
tn
tn
The variance of SZn is given by (SZn )2 = ( SWn − I(0; 1) )2 = 2 (SWn )2 − 2 SWn I(0; 1) + (I(0; 1) )2
tn+1
tn+1 = 3 − 2 s ds + s2 ds tn
3
= − 2
1 2
2
tn
3 1 3 1 3 = + = 1−1+ : 3 3 3
(B.8)
In Eq. (85), we have used the following identity from the properties of the Ito integrals [27]
I(f)I(g) = f(s; !0)g(s; !0) ds for any f ∈ L2 over any @xed interval : (B.9)
352
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
The covariance or correlation SZn SWn is also calculated SZn SWn = ( SWn − I(0; 1) )SWn = (SWn )2 − I(0; 1) SWn = · − 12 2 = 12 2 : To generate SWn and SZn from uniform random numbers picked from the [0; 1]; we generate two normally distributed random numbers U1 and U2 with and variance and set √ SWn = U1 ;
1 SZn = SWn + √ U2 : 2 3
(B.10) interval mean 0 (B.11) (B.12)
Numerically, this is precisely the way the processes SWn and SZn are evaluated; in addition, the Levy construction, 1 which consists of re@ning the Wiener process at points halfway between those points at which it has already been generated, is also used. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] 1
H. Haken, Synergetics, An Introduction, 2nd Edition, Springer, New York, 1978. K.H. Ho0mann, Z. Phys. B 49 (1982) 245. J. Olarrea, F.J. de la Rubia, Phys. Rev. Lett. 53 (1996) 268. R. Graham, Phys. Rev. A 25 (1982) 3234; R. Graham, T. Tel, Phys. Rev. A 33 (1985) 1322; R. Graham, T. Tel, Phys. Rev. A 35 (1987) 1328. W. Vance, J. Ross, J. Chem. Phys. 91 (1989) 7654; W. Vance, J. Ross, J. Chem. Phys. 105 (1996) 479; W. Vance, J. Ross, J. Chem. Phys. 108 (1997) 2088. C. Meunier, A.D. Verga, J. Stat. Phys. 50 (1988) 345. S. Kabashima, T. Kawakubo, Phys. Lett. A 70 (1979) 375; S. Kabashima, T. Kawakubo, T. Okada, J. Appl. Phys. 50 (1979) 6296. F. Mauricio, S. Valsco, J. Stat. Phys. 43 (1986) 521. W. Horthemke, R. Lefever, Noise-Induced Transitions, Springer, Berlin, 1984. L. Gammaitoni, P. Hanggi, P. Jung, F. Marchesoni, Rev. Mod. Phys. 70 (1998) 223. W. Hohmann, J.M. Schneider, J. Phys. Chem. 100 (1996) 5388. T. Alarc, M. Rubi, J. Chem. 108 (1998) 7367. A.A. Zaikin, J. Garcia-Ojalvo, L. Schimansky-Geier, Phys. Rev. E 60 (1999) 6275. Z. Hou, H. Xin, Phys. Rev. E 60 (1999) 6329. M.O. Vlad, J. Ross, F.W. Schneider, Phys. Rev. E 57 (1998) 4003. P.W. Anderson, Phys. Rev. 109 (1958) 1492; K. Ishi, Prog. Theor. Phys. Suppl. 53 (1973) 77. M.O. Vlad, F.W. Schneider, J. Ross, submitted for publication. J. Guckenheimer, P. Holmes, Nonlinear Oscillations, Dynamical systems, and Bifurcations of Vector Fields, Springer, New York, 1983. G. Nicolis, Introduction to Nonlinear Science, Cambridge University Press, Cambridge, 1995.
Numerically, that the Levy construction consists of re@ning the Wiener process at points halfway between those points at which it has already been generated: W (t + h=2) = 12 [W (t) + W (t + h)] + N(0; h=4). N(0; h=4) is a normally distributed random number of 0 mean and h=4 variance.
M. Al-Ghoul / Physica A 307 (2002) 331 – 353
353
[20] C. Van den Broeck, M. Malek Mansour, F. Baras, J. Stat. Phys. 28 (1982) 557; F. Baras, M. Malek Mansour, C. Van den Broeck, J. Stat. Phys. 28 (1982) 577. [21] E. Knobloch, K.A. Wiesenfeld, J. Stat. Phys. 33 (1983) 611. [22] M. Schumaker, Phys. Lett. A 122 (1987) 317. [23] P.H. Coullet, C. Elphick, E. Tirapegui, Phys. Lett. A 111 (1985) 277. [24] P. Richter, P. Rehmus, J. Ross, Prog. Theor. Phys. 66 (1981) 385. [25] M. Al-Ghoul, PCCP 2 (2000) 3773. [26] R.Z. Khasminsky, Stochastic Stability of Di0erential Equations, Sijtho0 and Noorho0, Alphen naan den Rijn, 1980. [27] L. Arnold, Stochastic Di0erential Equations, Wiley, New York, 1974. [28] S. Wiggins, Introduction to Applied Nonlinear Dynamical Systems and Chaos, Springer, New York, 1990. [29] M. Ipsen, F. Hynne, P.G. Sorensen, preprint chao-dyn=9803006, 1998. [30] J.L. Doob, Stochastic Processes, Wiley, New York, 1972. [31] M.C. Mackey, A. Longtin, A. Lasota, J. Stat. Phys. 60 (1990) 735. [32] P.W. Anderson, J. Phys. Soc. Japan 9 (1954) 316. [33] R. Kubo, Adv. Chem. Phys. 15 (1969) 100. [34] N.G. Van Kampen, Phys. Rep. C 24 (1976) 171. [35] W. Vance, J. Ross, J. Chem. Phys. 91 (1989) 7654. [36] P. Kloeden, E. Platen, Numerical Solution of Stochastic Di0erential Equations, Springer, Berlin, 1992.