Separatrix crossing due to multiplicative colored noise
aa r X i v : . [ n li n . C D ] A ug Separatrix crossing due to multiplicative colored noise.
Jean-R´egis Angilella
Normandie Universit´e, UNICAEN, UNIROUEN, ABTE, ESIX Cherbourg, Caen 14000, France.
Abstract
The effect of weak multiplicative colored noise on the dynamics of a Hamiltonian system isstudied by means of asymptotic methods, in the vicinity of homoclinic or heteroclinic trajectories.A general expression for the probability of noise-induced separatrix crossing is obtained, and isillustrated by means of a two-well Duffing oscillator. It shows how weak noise can significantlyaffect the dynamics near separatrices. In addition, the influence of the degree of non-linearityof the noise amplitude is discussed.Keywords: hamiltonian systems ; heteroclinic orbits ; noise.
Preprint submitted to Physica D August 17, 2020 eparatrix crossing due to multiplicative colored noise.
Jean-R´egis Angilella
Normandie Universit´e, UNICAEN, UNIROUEN, ABTE, ESIX Cherbourg, Caen 14000, France.
1. Introduction
Noise, even weak, can have tremendous effects on dynamical systems. This is a well knownphenomenon which has been studied long ago, and in many contexts [19, 16, 15, 14, 18, 17, 26,10]. Physical systems are particularly sensitive to noise in the vicinity of unstable equilibriumpositions. Here, the dominant forces almost balance each other and noise can thrust the systemtowards unpredictable directions. In the presence of homoclinic (or heteroclinic) cycles relatingsuch equilibrium points, weak noise can create homoclinic bifurcations leading to an extremelycomplex dynamics. The occurrence of this complexity can be predicted by means of stochasticMelnikov functions [4, 21, 7, 6, 5, 24], which emerge as soon as one calculates the jump ofundisturbed Hamiltonian for a particle moving near the separatrix. Such Melnikov functionshave been used in the past to derive necessary conditions for noise-induced escapes [24, 23].However, in spite of these active researches, the probabilistic dynamics in Hamiltonian systemsubmitted to noise is still poorly understood, especially in the case of colored noise.In many situations of interest, noise takes the form of a piecewise regular forcing, whichchanges randomly and at random times. These time intervals being finite, noise has a non-zeroautocorrelation length and enters into the category of colored noise. In a recent paper [2], anexpression for the probability of noise-induced separatrix crossing has been derived for any two-dimensional Hamiltonian system submitted to an additive noise of this kind (Kubo-Andersonprocess). The aim of the present paper is to generalize this result to multiplicative noise. Indeed,many physical systems are submitted to such noise, which is by essence non-uniform in the phasespace. Random terms being proportional to some function of the state variables, some portionsof the phase space might be free of noise, whereas other portions are noisy. The transition fromthe former to the latter can be abrupt if the amplitude of the noise is a non-linear functionof the state variables. For example, a pendulum attached to a non-inertial reference frameundergoing piecewise-constant accelerations will enter into this category. In another context,particles suspended in a turbulent flow modelled by means of a stochastic model [11, 8], are alsosubmitted to a random force of this kind.We have therefore investigated the effect of weak, colored and multiplicative noise on thedynamics of a system in the vicinity of homoclinic or heteroclinic trajectories. To be precise,
Preprint submitted to Physica D August 17, 2020 e consider a Hamiltonian system with two degrees-of-freedom r = ( q, p ) and undisturbedHamiltonian H ( r ). In the absence of perturbation (either deterministic or random), trajectoriesin the phase space correspond to H ( r ( t )) = constant . When perturbations are present, theHamiltonian is no longer conserved along trajectories r ( t ), but one can still use H ( r ( t )) todetermine in which portion of the phase space the state-variable r ( t ) is located. The calculationof the variations of H along a perturbed trajectory r ( t ) therefore allowed us to quantify the effectof noise on the dynamical system. Section 2 presents this methodology, as well as the calculationsand the main results. It is followed by a section devoted to the two-well Duffing equation, wheretheoretical results are compared to numerical simulations. It also contains a discussion aboutthe effect of the non-linearity of the noise amplitude. Conclusions and perspectives are drawnin section 5.
2. General considerations
We consider a perturbed Hamiltonian system of the form˙ q = ∂H∂p + Λ U ( r ) + ε f ( r ) ξ ( t ) , (1)˙ p = − ∂H∂q + Λ U ( r ) + ε f ( r ) ξ ( t ) , (2)where r ( t ) = ( q ( t ) , p ( t )) is the state variable, H ( r ) is the Hamiltonian, ( U , U ) and ( f , f ) aredeterministic vector fields and ε i and Λ are positive constants. Terms Λ U i can be thought ofas a deterministic perturbation of the Hamiltonian system (e.g. friction, gravity, etc.), whereas ε i f i ( r ) ξ i ( t ) are random perturbations. The functions ξ i ( t ) are assumed to be random processes,taking constant random values ξ ik during time intervals [ τ k , τ k +1 ], where τ k are random times.We assume that the time intervals δτ k = τ k +1 − τ k have an exponential distribution, to manifestthe fact that noise is memoryless. Ensemble averages will be denoted by h·i . The average timeduration h δτ k i is denoted δτ in the following. Random processes ξ i ( t ) are therefore coloredKubo-Anderson noises [12, 1, 3], with an exponential auto-correlation function. In addition, weassume that the random pulses are centred and normalized, that is h ξ ik i = 0 and h ξ ik i = 1.In the absence of either deterministic or stochastic perturbations (Λ = ε i = 0), trajectoriesin the phase space correspond to H ( r ( t )) = constant . When perturbations are present, we use H ( r ( t )) to determine the location of the state-variable in the phase space. If AB denotes aheteroclinic trajectory of the undisturbed system (Fig. 1), then H ( r ( t )) < H ( A ) means that r ( t ) is on the left-hand size of arc AB , and vice versa. Therefore, by examining the jump ofHamiltonian ∆ H = H ( r ( τ b )) − H ( r ( τ a )) (3)between two arbitrary times τ a and τ b , one can check whether the state variable r ( t ) crossed theseparatrix during the time interval [ τ a , τ b ] (see for example Refs. [9, 13]). Here, the times τ a and3 igure 1: Sketch of a heteroclinic trajectory (separatrix) between two hyperbolic saddle points A and B , in theunperturbed Hamiltonian system. Vector r ( t ) corresponds to a trajectory along the separatrix, and r ( t ) is asolution of the perturbed system. τ b will be defined as the times where r ( t ) passes nearest to A and B respectively. Injecting Eqs.(1) and (2) into (4) leads to [2]∆ H = ∆ H + ε Z τ b τ a H ,q f ( r ( t )) ξ ( t ) dt + ε Z τ b τ a H ,p f ( r ( t )) ξ ( t ) dt (4)where commas indicate partial derivations and∆ H = Λ Z τ b τ a ∇ H · U dt (5)is the jump of Hamiltonian under the effect of the sole deterministic perturbation, which can beapproximated by: ∆ H ≃ Λ Z BA U ( q, p ) dq − Z BA U ( q, p ) dp ! (6)where integrations are done over the arc AB .
3. Calculation of the jump of Hamiltonian
To exploit further Eq. (4), we need to calculate the two time integrals. Because the trajectoryconsidered here is close to separatrix AB , as both Λ and ε i ’s are small, we write [9]: r ( t ) ≃ r ( t − t )where r ( t ) = ( q ( t ) , p ( t )) is a solution of the undisturbed system, that is ˙ q = H ,p and ˙ p = − H ,q with r ( −∞ ) = A and r (+ ∞ ) = B . In addition, we introduced the time t where r ( t )4s nearest to r (0) ∈ ] A, B [. This leads to:∆ H = ∆ H − ε Z τ b τ a ˙ p ( t − t ) f ( r ( t − t )) ξ ( t ) dt + ε Z τ b τ a ˙ q ( t − t ) f ( r ( t − t )) ξ ( t ) dt. (7)Noise being piecewise constant, the first integral can be split into a sum of elementary integralsof ˙ p f over [ τ k , τ k +1 ]. In addition, if ˙ p keeps a constant sign on [ τ k , τ k +1 ], these elementaryintegrals can be calculated by the theorem of the mean value: Z τ k +1 τ k ˙ p ( t − t ) f ( r ( t − t )) dt = f ( r ( τ k +1 / − t )) Z τ k +1 τ k ˙ p ( t − t ) dt = f ( r ( τ k +1 / − t )) δp k , where τ k +1 / ∈ ] τ k , τ k +1 [ and δp k = p ( τ k +1 − t ) − p ( τ k − t ). The same treatment being donefor the second integral, we get:∆ H = ∆ H − ε N − X k =0 ξ k f ( r ( τ k +1 / − t )) δp k + ε N − X k =0 ξ k f ( r ( τ k +1 / − t )) δq k . (8)We introduce the unit vector perpendicular to the separatrix at point r ( τ k − t ), approximatedas: n k = ( n k , n k ) = (cid:16) − δp k , δq k (cid:17) /δs k where δs k is the discrete arc-length element: δs k = δp k + δq k . The jump of Hamiltoniantherefore takes the form of a sum of random variables X k :∆ H = ∆ H + N − X k =0 X k (9)where X k = X i =1 , ε i ξ ik n ik f ik δs k (10)and f ik = f i ( r ( τ k +1 / − t )). Being a sum of a large number of random variables (apart fromthe deterministic term ∆ H ), the jump of Hamiltonian ∆ H will be calculated below by meansof the generalized central limit theorem.The times τ k are decorrelated from the amplitudes ξ ik of the random force. Hence, theaverage of X k is: h X k i = X i =1 , ε i h ξ ik i h n ik f ik δs k i = 0since h ξ ik i = 0. Following the methodology of Ref. [2] we introduce: Z N = 1 S N N − X k =0 X k (11)5here S N = P k h X k i . If X k satisfies the Lyapunov condition, i.e. if there exists d > S dN N − X k =0 h| X k | d i → N → ∞ , then the random variable Z N converges to a centered Gaussian variable Z with unitvariance as N → ∞ . This allows us to approximate the jump of Hamiltonian as a Gaussianvariable: ∆ H ≃ ∆ H + σ Z, (13)where σ is the limit of S N as N → ∞ . ∆ H The sum of the variances S N can be derived in a compact form in the case where the times ofdiscontinuities τ k , their increments δτ k , and the amplitudes of the random forcing term ξ ik areindependent. We follow the same main lines as in Ref. [2]. First, by noticing that h ξ k ξ k i = 0,and using h ξ ik i = 1, we get: h X k i = X i =1 , ε i h n ik f ik δs k i . (14)By definition, the discrete arc-length element δs k is related to the velocity modulus on theseparatrix by δs k = u k δτ k where u k = | ˙r ( τ k − t ) | . Also, the normal vector components n ik , as well as f ik and u k arefunctions of the times τ k , and are decorrelated from δτ k . Therefore, h n ik f ik u k δτ k i = 2 h n ik f ik u k i δτ where δτ = h δτ k i and we made use of h δτ k i = 2 δτ for exponentially distributed time increments.Finally, we have S N = X k h X k i = 2 δτ X i X k ε i h n ik f ik u k i δτ. (15)The sum over k ∈ [0 , N −
1] in this last equation is of the form: h X k g ( τ k ) δτ i = h X k g ( τ k ) δτ k i + h X k g ( τ k )( δτ − δτ k ) i . (16)The first sum on the right-hand-side of (16) is a Riemann sum that converges to R τ b τ a g ( t ) dt if g ( t ) is integrable. The second sum is zero because τ k is decorrelated from δτ k . We thereforeconclude that, as N ≫ S N can be approximated by σ ≃ δτ X i =1 , ε i Z τ b τ a n i ( r ( t − t )) f i ( r ( t − t )) | ˙r ( τ k − t ) | dt. (17)6lternatively, introducing the infinitesimal arc-length ds along arc AB , defined by ds = | ˙r ( τ k − t ) | dt : σ ≃ δτ X i =1 , ε i Z BA n i ( s ) f i ( s ) |∇ H ( s ) | ds. (18)This expression is of interest since it involves the normals to the separatrix, which can lead tosimple expressions for straight or circular separatrices, or in the case of isotropic noise. A usefulalternative expression can be obtained in terms of integrals of ˙ p dp and ˙ q dq , by noticing that n |∇ H | = − ˙ p and n |∇ H | = ˙ q , so that: σ ≃ δτ (cid:16) ε Z BA f ( q, p ) ˙ p dp + ε Z BA f ( q, p ) ˙ q dq (cid:17) . (19) Finally, having shown that the jump of Hamiltonian ∆ H is a random variable that can beapproximated by a Gaussian variable with mean ∆ H (Eqs. (5) or (6)) and variance σ (Eq.(19)), we can determine the probability of noise-induced separatrix crossing in the same manneras for the case of additive noise [2]. Noise-induce crossing corresponds to ∆ H and ∆ H havingopposite signs, i.e. the effect of noise is opposed to the effect of the deterministic perturbation.The probability of this event is: P = 12 erfc | ∆ H | σ √ ! . (20)Eq. (20), together with (19) (or, alternatively (18)), is the main result of this work, and will beillustrated in the following section. It differs from the additive case [2] in that f and f arenow averaged in the calculation of the variance (19). Note that, to obtain this result, we haveassumed that ˙ p and ˙ q kept a constant sign on [ τ k , τ k +1 ]. This assumption is needed to applythe theorem of the mean value to elementary integrals of ˙ p f and ˙ qf over [ τ k , τ k +1 ] (see Eq.(7) and below). This hypothesis is not so restrictive, since the intervals of integration are small.In practice, it is satisfied in almost all intervals [ τ k , τ k +1 ].
4. Application to the Duffing equation
The probability calculated above can be applied to a wide variety of dynamical systems. Herewe chose to illustrate these results by means of an inverted pendulum, which, in the limit of weakamplitudes, can be modelled as a two-well Duffing equation. It has been shown that additiveand multiplicative noises in this family of dynamical systems had rather different effects [25].The pendulum considered here is assumed to be composed of a rigid stick with length l , on thetop of which a tiny object with mass m has been fixed. It rotates around a fixed horizontal axiscoinciding with its lower extremity, and a spring with stiffness k exerts a torque − kθ that forces7 igure 2: Sketch of the phase portrait of the unperturbed inverted pendulum. The dashed line is a realization ofthe perturbed system. the pendulum back to its vertical position, where θ ( t ) denotes the angle between the verticalaxis and the stick. In addition, viscous friction with coefficient b is present, together with amultiplicative random torque proportional to the angle θ (noise with a ”linear amplitude”).The corresponding dynamical equation is ml ¨ θ ( t ) = − k θ ( t ) + mgl sin θ ( t ) − b ˙ θ ( t ) + ε θ ( t ) ξ ( t ) , (21)where ε is a positive constant. Assuming that | θ | is small we make use of the approximationsin θ ≃ θ − θ / | θ | < π/ ml ¨ θ = ( mgl − k ) θ − mgl θ − b ˙ θ + ε θ ξ. (22)It can be re-written as: ¨ θ = α θ − β θ − Λ ˙ θ + ε θ ξ (23)where α = ( mgl − k ) / ( ml ), β = g/ (6 l ), Λ = b/ ( ml ) and ε = ε/ ( ml ). This is a perturbedhamiltonian system of the form (1)-(2) with ( q, p ) = ( θ, ˙ θ ) and H ( q, p ) = − α q + 14 β q + 12 p (24)and with a deterministic perturbation corresponding to U ( q, p ) = 0 and U ( q, p ) = − p . In thiswork, we will assume that k < mgl , so that α >
0. The phase portrait of the unperturbedsystem is obtained from the lines H = constant , and is sketched in Fig. 2: it corresponds to apair of homoclinic trajectories attached to a hyperbolic saddle point located at the origin, andforming loops around elliptic points located at q = ± p α/β and p = 0. Separatrices intersectthe q axis at q = ± p α/β . When friction is present, the noise-free dynamics corresponds totrajectories spiraling towards the elliptic points, which are asymptotically stable in this case.Any trajectory located inside one of the homoclinic loops will remain in this loop for all times.8oise, even weak, might break this picture. The probability that a pendulum, released nearthe origin and inside a homoclinic loop, exits this loop can be non-zero. Such a trajectory issketched in Fig. 2. According to the theory presented above, the calculation of this probabilityis given by Eq. (20). The jump of Hamiltonian in the absence of noise reads (from Eq. (6)):∆ H = − Λ Z BA p dq = − Z CA p dq for symmetry reasons, where C = ( p α/β,
0) is the intersection between the q axis and thehomoclinic loop AB (with B = A here). On eliminating p since H ( p, q ) = 0 along AC we get∆ H = −
43 Λ α / β . (25)The sign of the deterministic jump of Hamiltonian indicates that friction drives the pendulumtowards the interior of the homoclinic loops, as expected. One can check that the Lyapunovcondition is satisfied here (see below). The jump of Hamiltonian can therefore be approximatedby a Gaussian variable, with standard deviation σ given by Eq. (19). Again, using H ( q, p ) = 0to eliminate p , we get: σ = 4 √ √ √ δτ ε α / β . (26)Finally, by injecting Eqs. (25) and (26) into (20) we obtain the probability of noise-inducedseparatrix crossing for the inverted pendulum: P = 12 erfc √
156 Λ α / ε √ δτ ! . (27)Note that it is independent of the coefficient β of the cubic term in the dynamical equation(24). This theoretical probability is plotted versus ε in Fig. 3, in the case α = β = 1 s − , δτ = 0 . s and for various values of the friction coefficient Λ (in units of s − ). In addition,results from numerical solutions of Eq. (23) involving 50000 stochastic trajectories startingat point (0 . , . q < The theory presented above allows to study the effect of non-linear multiplicative noises, i.e.random forcing where either f or f is a non-linear function of r . In the case of the Duffingoscillator studied in the previous section, the random force is weaker near the origin than nearpoint C , and this non-uniformity can be even more pronounced if the amplitude of this force9 ε P (cid:0)(cid:0)(cid:1)(cid:1) p q Figure 3: Probability of noise-induced separatrix crossing for the inverted pendulum (Eq. (21)) (multiplicativenoise with linear amplitude). Solid line: theoretical result (27). Circles: probability obtained from numericalsolutions of Eq. (23) involving 50000 stochastic trajectories.
10s a non-linear function of ( q, p ): this might affect the probability of noise-induced crossing.To quantify this effect we have replaced the term ε θξ in Eq. (23) by ε θ γ ξ , where γ is a realnumber. This corresponds to f ( r ) = q γ .Prior using the general expression of the probability calculated in section 2, we must checkthat the Lyapunov condition (12) is satisfied. To this end, we first note that S N ≃ K √ δτ , where K is independent of N . This is a consequence of Eq. (19). Also, ε = 0 for the Duffing equation,so that: | X k | d = ε d | ξ k | d | n k f k u k | d | δτ k | d (28)where f k = q γ ( τ k − t ) will be denoted as q γk . On averaging, and taking into account the factthat the random variables ξ k , τ k and δτ k are decorrelated, we get: h| X k | d i = ε d h| ξ k | d i h| n k q γk u k | d i h| δτ k | d i . (29)The time increments δτ k having an exponential distribution with mean δτ we have h| δτ k | d i = δτ d Γ( d +3). We assume that the moment h| ξ k | d i is equal to some finite value h| ξ | d i . Also,we assume that, for all k , h| n k q γk u k | d i is bounded by some positive constant m independentof N . This assumption is satisfied if γ >
0, since | q γk | is bounded, but might not systematicallybe true when γ < | q γk | diverges when q k is close to 0. However, probabilities in the case γ < < S dN N − X k =0 h| X k | d i < N × ε d h| ξ | d i m ( K √ δτ ) d δτ d Γ( d + 3) = constant/N d/ (30)since δτ = constant/N . We therefore conclude that P k h| X k | d i /S dN → d >
0, sothat the Lyapunov condition is satisfied.We will then make use of the general expression for the probability (20). The standarddeviation of the jump now depends on γ and will be denoted by σ γ in the following. By makinguse of expression (19) we get, after some algebra: σ γ = σ f ( γ ) (31)where σ is given by Eq. (26) and f ( γ ) = 15 √ π γ + 1)Γ( γ + 5 / (cid:18) αβ (cid:19) γ − (32)provided γ > −