arXiv:cond-mat/0502389v1 [cond-mat.stat-mech] 16 Feb 2005Stationary distributions of a noisy logistic processP. F. G´ ora†M. Smoluchowski Institute of Physics and Complex Systems Research Center Jagellonian University, Reymonta 4, 30-059 Krak´ ow, Poland(Received )Stationary solutions to a Fokker-Planck equation corresponding to a noisy logistic equation with correlated Gaussian white noises are con- structed. Stationary distributions exist even if the corresponding deter- ministic system displays an unlimited growth.Positive correlations be- tween the noises can lead to a minimum of the variance of the process and to the stochastic resonance if the system is additionally driven by a periodic signal.PACS numbers: 05.40.Ca, 02.50.Ey, 87.23.Cc1. IntroductionThe logistic equation˙ x = ax(1 − x), a > 0, x > 0,(1)is one of the best-known and most popular models in population dynamics. Perturbing this equation by a multiplicative noise is an obvious generaliza- tion of the deterministic theory, aiming at describing populations that livein an ever-changing environment. The logistic equation with a fluctuating growth rate˙ x = (a + pξ(t))x(1 − x),(2)has been first discussed by Leung in Ref. [1] and later by many other authors. Recently in Ref. [2] we have discussed a further generalization of (2) in whichboth the growth rate and the limiting population level fluctuate, and thesefluctuations are correlated in time:˙ x = (a + pξm(t))x − (b + q ξa(t))x2.(3)†e-mail:

[email protected](1)2fp2printed on February 2, 2008Here ξm,aare two Gaussian white noises (GWNs) that satisfy hξi(t)i = 0, hξi(t1)ξi(t2)i = δ(t1−t2), i = m,a, hξm(t1)ξa(t2)i = cδ(t1−t2), p, q are theamplitudes of the two noises and the correlation coefficient c ∈ [−1,1]. For the sake of terminology, we will sometimes call the noise ξm(t) “multiplica- tive” and the noise ξa(t) “additive”; see Eq. (5) below for a rationale behind these names. Please note, though, that on the level of Eq. (3) both these noises are coupled multiplicatively to the process x(t). Note also that ifb > 0, the corresponding deterministic equation converges to a stable fixed point; accordingly, we will call a system with a positive b “convergent”. If b 6 0, the corresponding deterministic system displays unlimited growth. We will call a system with a negative b “exploding”. Recently Mao et al. have shown in Ref. [3] that in the absence of ξm, the system (2) remains positive and bounded even in the “exploding” case.In Ref. [4] we havediscussed certain difficulties that may arise in numerical simulations of such a system. This paper generalizes the result of Mao et al. to the case of both noises present. In Ref. [2] we have shown how the dynamics (3) is related to the problem of a linear stochastic resonance. We have mapped the nonlinear equation (3) into a linear Langevin equation with two correlated noises and used the solutions of the latter to heuristically explain the behaviour of the noisylogistic equation. Specifically, the substitutiony =1 x(4)converts Eq. (3) into a linear equation˙ y = −(a + pξm(t))y + b + q ξa(t)(5)which can be solved exactly for realizations of the process y(t). The process (5) has a convergent mean ifa −1 2p2> 0(6a)and a convergent variance if a stronger conditiona − p2> 0(6b)holds. Using the properties of the process y(t), useful prediction can be made about the noisy logistic process. Since in the presence of correlations,for certain values of parameters the variance of y(t) first shrinks and then grows as a function of the “multiplicative” noise strength, q, one expects a similar behaviour for the logistic process x(t) as well. These predictions have been corroborated numerically in Ref. [2]. In particular, if c = ±1,fp2printed on February 2, 20083bp ∓ aq = 0 and the condition (6b) is satisfied, the variance of the process y(t) vanishes. The relation between the processes y(t) and x(t) = 1/y(t) intuitively means that if almost all realizations of the former asymptotically reach the same constant value, so do almost all realizations of the latter. However, as a formal relation between moments of these processes is not trivial, the predictions based on the properties of y(t) have only a heuristic value. In the following we will construct mathematically exact stationary solu- tions to the Fokker-Planck equation corresponding to Eq. (3) and re-examine the above results from the point of view of these stationary solutions. Fur- thermore, we will show numerically that if the parameters undergo periodic (for example, circaannual or seasonal) oscillations, a positive correlation between the noises leads to a stochastic resonance. In the Appendix we extend to the case of two correlated noises the proof originally proposed by Mao et al. in Ref. [3] that solutions to Eq. (3), when started from a positive initial condition, never become negative almost surely.2. The Fokker-Planck equationThe problem of constructing a Fokker-Planck equation correspondingto a process driven by two correlated Gaussian white noises has been first discussed in Ref. [5], where the two noises have been decomposed into two independent processes. The same result has been later re-derived in [6], where the authors have attempted to avoid an explicit decomposition of the noises but eventually resorted to a disguised form of the decomposition. The general Langevin equation˙ x = h(x) + g1(x)ξm(t) + g2(x)ξa(t),(7)where x(t) is a one-dimensional process and ξm,aare as in Eq. (3), leads to the following Fokker-Planck equation in the Ito interpretation:∂P(x,t) ∂t= −∂ ∂xh(x)P(x,t) +1 2∂2 ∂x2B(x)P(x,t),(8a)where B(x) = [g1(x)]2+ cg1(x)g2(x) + [g2(x)]2.(8b)In the case of Eq. (3) the corresponding Fokker-Planck equation there- fore reads∂P(x,t) ∂t= −∂ ∂x[(a − bx)xP(x,t)] +1 2∂2 ∂x2?x2(p2− 2cpqx + q2x2)P(x,t)?.(9)4fp2printed on February 2, 2008It is apparent that the absolute signs of the two noise amplitudes do not in-fluence the solutions to the above equations, only their relative sign, sgn(pq), does. In the following we will assume that sgn(pq) = +1. This comes at no loss to the generality as Eq. (9) is invariant under a simultaneous change ofsigns of pq and the correlation coefficient, c. Stationary solutions to Eq. (9) are the normalizable solutions to [7, 8]x2(p2−2cpqx+q2x2)dPst(x)dx+2x?2q2x2+ (b−3cpq)x − (a−p2)?Pst(x) = 0.(10)A slight modification of the argument presented originally in Ref. [3] shows that all solutions to Eq. (3) that start from a positive initial condition remain positive almost surely; see the Appendix for a proof. Physically speaking, this results from a presence of an absorbing barrier at x = 0 in Eq. (3): should the population suddenly drop to zero, it would stay there forever. The formal result of Ref. [3], extended here to the case of two noises present, ensures that the population never actually becomes nonpositive although in certain cases (see below) it may dynamically cluster in a close proximity of x = 0+. Therefore, we can divide both sides of (10) by x and obtainx(p2−2cpqx+q2x2)dPst(x)dx+2?2q2x2+ (b − 3cpq)x − (a − p2)?Pst(x) = 0,(11) provided Pst(x) is normalizable over the x > 0 semiaxix. One may be tempted to try to immediately solve Eq. (11) by standardmethods, but a word of caution is needed here: should the coefficient at dPst/dt vanish, special care must be taken.2.1. The case of only one noise presentBefore proceeding to the general case, we will discuss the special cases where only one of the amplitudes p, q does not vanish. Purely “multiplicative” noise. If q = 0, the Langevin equation (3) re- duces to˙ x = (a + pξm(t))x − bx2(12)and we find for Pst(x)Pst(x) = Nx2(a−p2)/p2exp? −2bxp2? ,(13)fp2printed on February 2, 20085where N is a normalization constant.This function is normalizable for x > 0 if b > 0, or if the system is convergent, anda −1 2p2> 0,(14a)which determines the behaviour of the distribution around x = 0. Moreover, ifa − p2> 0,(14b)Pst(x) goes to zero as x → 0+and has a maximum at x = (a − p2)/b. Note that the conditions (14a), (14b) coincide with the conditions (6a), (6b) determining the properties of the linear process (5). For p2> a >1 2p2,the distribution is mildly divergent at x = 0 and decreases monotically with an increasing x. Purely “additive” noise. If p = 0, the Langevin equation takes the form˙ x = (a − bx)x + qx2ξa(t)(15)and we obtain for the stationary distributionPst(x) =N x4exp? −a − 2bxq2x2? (16)which is normalizable whenever a > 0. Note that no bounds on b are im- posed: The stationary distribution exists for both convergent and exploding systems. This special case falls into a broad category discussed recently by Mao et al. in Ref. [3], without further generalizations provided by the present paper. The distribution (16) has a maximum at the positive root of2q2x2+ bx − a = 0.(17)2.2. The general caseIf both noises are present and are not maximally correlated, |c| 6= 1, the solution to (11) readsPst(x) =Nx2(a−p2)/p2(p2−2cpqx+q2x2)(a+p2)/p2exp−2(bp−acq)arctan? qx−cp√1−c2p?√1 − c2p2q.(18) Since the exponential term is limited, the convergence (normalization) prop- erties of (18) are determined by those of the fractional term. The denomi- nator is always strictly positive. For x → ∞, Pst(x) ∼ x−4for all possible6fp2printed on February 2, 200801234567890.00.51.01.52.0Pst(x)xq = 0.101234567890.00.51.01.52.0Pst(x)xq = 0.501234567890.00.51.01.52.0Pst(x)xq = 0.601234567890.00.51.01.52.0Pst(x)xq = 5.0Fig.1. Stationary distributions (18) in a strongly correlated case, c = 0.99. Clock- wise, from top-left q = 0.1, q = 0.5, q = 0.6, and q = 5.0. Other parameters, common for all panels, are p = 0.5, a = b = 1.values of parameters. Pst(x) is therefore normalizable if it does not diverge too rapidly at x → 0+, or again if the condition (14a) holds. Either in this case, no bounds on b are imposed. If the condition (14b) holds as well, the distribution (18) approaches zero as x → 0+, but this condition is no longer associated with the presence of a maximum. The maximum of Pst(x) coincides with the positive root of2q2x2+ (b − 3cpq)x − (a − p2) = 0,(19)cf. Eq. (11), provided such a root exists. It certainly does for a − p2> 0, but it can appear also for p2> a >1 2p2, where the distribution is mildlydivergent. Example stationary distributions in a strongly correlated case are presented on Fig. 1. For large values of the additive noise strength, q, the distributions are highly skewed and squeezed against the x = 0 axis. For comparison, on Fig. 2 we show example stationary distributions for the un- correlated and a strongly anticorrelated cases. The distributions presented are skewed and much wider than the distribution from Fig. 1 with the samefp2printed on February 2, 200870.00.51.01.52.00.00.51.01.52.0Pst(x)xc=00.00.51.01.52.00.00.51.01.52.0Pst(x)xc=-0.99Fig.2.Stationary distributions for the uncorrelated (c = 0, left panel) and a strongly anticorrelated (c = −0.99, right panel) cases. Other parameters are a = b = 1, p = q = 0.5.value of q = 0.5. These distributions also get squeezed as the additive noise strength becomes large. Note that the distribution corresponding to the anticorrelated noises is already more squeezed than the distribution for the uncorrelated case. If the distribution (18) is normalizable, it has a convergent mean and a variance. Its higher moments are divergent.2.3. The maximally correlated caseIf the two noises are maximally (anti)correlated, c = ±1, Eq. (11) takes the formx(p ∓ qx)2dPst dx+ 2?(p ∓ qx)2+ q2x2+ (b ∓ pq)x − a?Pst= 0,(20)leading to the following candidate solution:Pst(x) = Nx2(a−p2)/p2(p ∓ qx)2(a+p2)/p2exp? ∓2(bp ∓ aq) pq(p ∓ qx)? .(21)With our sign convention adopted, sgn(pq) = +1, this solution is nor- malizable if c = −1 and the condition (14a) holds. The distribution (21) with the “+” sign can be obtained from Eq. (18) by taking the limit c → −1. This distribution decreases as x−4with x → ∞. The maximum of (21), if it exists, coincides with the positive root of Eq. (19) with c = −1. The case of c = +1 is more challenging. First, if the resonant conditionbp − aq = 0(22)8fp2printed on February 2, 2008holds, Eq. (20) is solved byPst(x) = δ? x −p q? (23)regardless of the value of p. This result is stronger than that reported in Ref. [2] where we could predict a δ-shaped distribution only in the a−p2> 0 case, or when the variance of the corresponding linear system was conver- gent, as otherwise any predictions based on the linear system failed. Notethat with the condition (22) statisfied, Eq. (3) reduces to a rescaled form of Eq. (2). If c = +1 and the condition (22) does not hold, Eq. (20) does not have a normalizable solution. This observation is slightly surprising, but formally speaking, it results from the fact that the double limitlim x→p/qc→1exp−2(bp − acq)arctan?qx−cp√1−c2p?√1 − c2p2q(24)does not exits: Its value depends on which route the singularity is ap- proached. The nonexistence of stationary solutions in the fully correlated, non-resonant case is, therefore, related to the essential singularity of thecomplex exponential at infinity. The fact that with c = +1, bp−aq 6= 0, thedrift and diffusive terms in Eq. (20) both vanish, but at different points, is the physical reason for this apparent oddity: If the population gets locatedaround the point of the vanishing diffusion, it is washed away by the drift,and if it gets located around the point of the vanishing drift, it diffusively leaks from there. Nevertheless, if a−1 2p2> 0, in numerical simulations the cases of c = 1 and c = 1 − ε with 0 - 2q0.000.050.100.150.200.25012345678910- 2qFig.3. The variance?x2?−hxi2determined form the distribution (18) as a function of the additive noise strength, q. Main panel: p = 0.5, the curves, from bottom to top, correspond to c = 0.99, c = 0.90, c = 0.75, c = 0.50, c = 0.25, c = 0, and c = −0.25, respectively. Inset: p = 1.1, the curves correspond, from bottom to top, to c = 0.99, c = 0.98, c = 0.97, c = 0.96, c = 0.95, and c = 0.94, respectively. Other parameters, common for all curves presented, are a = b = 1.correlated and resonant case, the stationary distribution becomes δ-shaped and its variance indeed vanishes, much as predicted by the linear system. Since we now know the mathematically exact stationary distributions, we can test the behaviour of the variance in the general case directly.Recall that the distribution (18) has the two first moments convergent whenever it is normalizable. Unfortunately, analytical expressions for these moments cannot be obtained, mainly due to the presence of the complicated exponential term. Therefore, we have calculated the moments by numeri- cally integrating over the distribution (18). Results are presented on Fig. 3. If the distribution approaches zero as x → 0+, or when the condition (14b)is satisfied, and if the two noises are positively correlated, 0 a >1 2p2,a distinct minimum in the variance of x also appears but it is present only forfairly large (and po