STOCHASTIC GLOBAL BIFURCATION IN PERTURBED
PROCEEDINGS OF THE FOURTH INTERNATIONAL CONFERENCE ON DYNAMICAL SYSTEMS AND DIFFERENTIAL EQUATIONS May 24 – 27, 2002, Wilmington, NC, USA
pp. 122–131
STOCHASTIC GLOBAL BIFURCATION IN PERTURBED HAMILTONIAN SYSTEMS LORA BILLINGS1 , ERIK BOLLT2 , DAVID MORGAN3 , AND IRA B. SCHWARTZ3 1 Department
of Mathematical Sciences, Montclair State University, Upper Montclair, NJ 07043 of Mathematics & Computer Science and Department of Physics, Clarkson University, Potsdam, NY 13699 3 Naval Research Laboratory, Code 6792, Plasma Physics Division, Washington, DC 20375 2 Departments
Abstract. We study two perturbed Hamiltonian systems in which chaos-like dynamics can be induced by stochastic perturbations. We show the similarities of a class of population and laser models, analytically and topologically. Both systems have similar manifold structure that includes bi-instability and partially formed heteroclinic connections. Noise takes advantage of this structure, inducing a global bifurcation and chaotic-like dynamics which exhibits mixed mode behavior of the original bi-stable solutions. We support these claims with numerical approximations of the transport between basins.
1. Introduction. In this article, we describe how two ordinary differential equation systems, from population dynamics and lasers, have the manifold structure sufficient for exciting noise induced chaos. Though coming from very different fields, we find that the single-mode periodically modulated class-B laser [12] and the seasonally forced nonlinear epidemic model [14] can both be reduced to weakly perturbed Hamiltonian systems, having equivalent Hamiltonians. The perturbations result in the existence of saddle-node bifurcations and lead to a bi-instability, or the existence of at least two unstable periodic saddles. When the manifolds of these saddles are close to a homoclinic or heteroclinic topology, periodic stochastic perturbations, or noise, can essentially complete the intersections to form a topological horseshoe, in a metaphorical sense. Since it is not sufficient to compute only Lyapunov exponents to detect or define stochastic chaos [9], we compute transport of a stochastic flux from one basin to another [7, 6]. This is done using the Galerkin approximation of the stochastic Frobenius-Perron operator [4]. Increasing the standard deviation of the noise reduces the stability of the stable periodic attractors and increases the frequency of a trajectory visiting the neighborhood of the induced homoclinic or heteroclinic tangle. It is well known in the noiseless case that homoclinic or heteroclinic tangles create the horseshoe topology that allows chaotic dynamics [15]. The stochastic perturbations enable transport in a similar fashion, and we call this phenomena noise-induced chaos. We can numerically approximate where pseudo-crossings are induced by identifying the transport between the noiseless basins of attraction of two attracting orbits. Therefore, we can conjecture from numerical evidence how both models undergo a smooth transition to chaos by increasing the standard deviation of the noise. 1991 Mathematics Subject Classification. 34F05, 65P20, 34D08, 37M20. Key words and phrases. Stochastic dynamical systems, Hamiltonian, bifurcation, chaos.
122
STOCHASTIC GLOBAL BIFURCATION IN PERTURBED HAMILTONIAN SYSTEMS
123
We begin by presenting the two systems and their reduction to the same weakly perturbed Hamiltonian system. Then, we describe their manifold structures and conjecture how noise induces chaotic behavior. Finally, we present numerical approximations of the Lyapunov exponents and transport, or flux, in the system due to noise. This is similar to the phenomenal bifurcation described in [1] in which the standard deviation of the noise changes the Lyapunov exponent (integrated with respect to the invariant density) from negative to positive. The conclusion is that noise induces enough transport across basin boundaries that global effects contribute to the Lyapunov exponent bifurcation of stochastic chaos. 2. A model from population dynamics. Consider a compartmental model to describe the dynamics of how measles spreads in a large population. The population can be divided into disjoint subgroups which evolve in time. Based on singular perturbation theory as presented in [14, 10], the standard SEIR model is reduced to two equations. We call this the modified SI model (MSI): S (t) I (t)
= µ − µS(t) − β(t)I(t)S(t) α = β(t)I(t)S(t) − (µ + α)I(t). µ+γ
(1)
The people that may contract the disease are susceptible, S(t). Those individuals capable of transmitting the diseases are infective, I(t). This contact rate fluctuates with the seasons and can be approximated by a sinusoidal forcing of the form β(t) = β0 (1 + δ cos 2πt), with 0 ≤ δ < 1. The parameter values are fixed at µ = 0.02 (year)−1 , α = 35.84 (year)−1 , γ = 100 (year)−1 , β0 = 1575 (year)−1 , and we vary δ, the fluctuating contact rate amplitude. These are the standard parameters used in [10]. We will add noise discretely, not continuously, to both S and I once a period with respect to the forcing to represent fluctuations in the population [3]. When δ = 0, there are two steady states: the disease free equilibrium, (S, I) = (1, 0), and the endemic steady state, (µ + α)(µ + γ) µ αβ0 −1 . , (Se , Ie ) = αβ0 β0 (µ + α)(µ + γ) The stability of these steady states depends on the value of Q, where Q=
αβ0 . (µ + α)(µ + γ)
If Q < 1, the disease free equilibrium is stable, without any endemic state. When Q > 1, the disease free state is unstable, and the endemic state is attracting. Since we wish to understand the dynamics when the disease is prevalent, we will assume Q > 1. With a change of variables, we recenter the endemic equilibrium at the origin and the disease free equilibrium at infinity. Using S = Se (1 + x) and I = Ie (1 + y), the new system is x
=
y
=
S = µQ − µ(1 + x) − µ(1 + δ cos 2πt)(Q − 1)(1 + y)(1 + x) Se I = (µ + α)(1 + y)(δ cos 2πt + (1 + δ cos 2πt)x). Ie
(2)
124
LORA BILLINGS, ERIK BOLLT, DAVID MORGAN, AND IRA B. SCHWARTZ
To check the stability of the origin when δ = 0, we find eigenvalues of the Jacobian of Eq. (2). −µQ −µ(Q − 1) J|(0,0) = µ+α 0 µQ 1 λ± = − ± µ2 Q2 − 4µ(Q − 1)(µ + α). 2 2 1 1 Similar to an approach in [16, 14], recall the fact that µ, µ+γ , and µ+α are O(10−2 ). We introduce the small parameter = µ(Q − 1), and µ + α = ∆ . Using the Q substitution η = Q−1 , 4∆ − η 2 2 η λ± = − ± i. 2 2 The eigenvalues are in the √ form λ± = −r ± νi. To the leading order, the frequency of the oscillation is ν = ∆, and r = η/2. This is a weakly stable spiral sink. One last change of variables will transform the system so that when δ = 0, the linear part will be close to its Jordan form. We also assume that there is a small amplitude periodic perturbation, requiring δ/ to be small. Let δ¯ = δ/, x ¯ = ν x − y, and y¯ = y in Eq. (2). The new system is x ¯
= −ν y¯ − (¯ x(η − 1) + (1 + ν)(1 + y¯)(¯ x + ν δ¯ cos 2πt)) − 2 (¯ y (η − 1) 3 ¯ ¯ +(1 + ν)(1 + y¯)(¯ y+x ¯δ cos 2πt)) − (1 + ν)(1 + y¯)(δ cos 2πt)
y¯
= νx ¯(1 + y¯) + ν 2 (1 + y¯)(δ¯ cos 2πt) +ν(1 + y¯)(¯ y+x ¯δ¯ cos 2πt) + 2 ν y¯(1 + y¯)(δ¯ cos 2πt).
(3)
We can rewrite Eq. (3) as x ¯ y¯
¯ = −ν y¯ + f1 (¯ x, y¯, t, , δ) ¯ 2 (1 + y¯)(cos 2πt) + f2 (¯ ¯ = νx ¯(1 + y¯) + δν x, y¯, t, , δ),
(4)
with ¯ = fi (¯ ¯ fi (¯ x, y¯, t, , δ) x, y¯, t + 1, , δ) fi (0, 0, t, 0, 0) = 0, i = 1, 2. ¯ ¯ Now treat and δ as small parameters. Setting = δ = 0 in Eq. (4), we find the following reduced system, x ¯
y¯
= −ν y¯ = νx ¯(1 + y¯).
(5)
This is the same system that was found in earlier work on both the SEIR and SIR systems [14]. The system is conservative and the solution is described by level curves of the integral function E = x2 + 2y − 2 ln |1 + y|. It is interesting to note that the structure of such a first integral results from unfolding a degenerate Hopf bifurcation [8]. Moreover, as a result of the quadratic nonlinearities, the orbits arising the Hamiltonian may be unbounded. We can also recast Eq. (3) as a weakly perturbed Hamiltonian system by making one more change of variables. Let y¯ = −1 + ez¯ and drop terms of order higher than to find x(η − 1) + ez¯(1 + ν)(¯ x + ν δ¯ cos 2πt)) x ¯ = −ν(−1 + ez¯) − (¯ z¯
¯ 2 cos 2πt + ν(−1 + ez¯ + x = νx ¯ + δν ¯δ¯ cos 2πt).
(6)
STOCHASTIC GLOBAL BIFURCATION IN PERTURBED HAMILTONIAN SYSTEMS
125
The system in Eq. (6) can be rewritten as ∂H ¯ + f˜1 (¯ x, z¯, t, , δ) ∂ z¯ ∂H ¯ 2 ¯ + δν cos 2πt + f˜2 (¯ x, z¯, t, , δ), ∂x ¯
x ¯
= −
z¯
=
(7)
with 1 2 ¯ + ez¯ − z¯ − 1). H(¯ x, z¯) = ν( x 2 We will show below the similarity of this Hamiltonian system to the one we find for the following laser model. 3. The laser model example. Uncovered in a study of global mixed-mode chaos [5, 13], it was shown that the modulated CO2 laser model also reduces to a weakly perturbed Hamiltonian system similar to Eq. (6). They used approximate maps, derived by matched asymptotic expansions and subharmonic Melnikov theory, to prove the existence of period doubling and saddle-node bifurcations. Existence of the subharmonic Melnikov function also allowed for additional conclusions concerning the global dynamics and transverse manifold intersections. The goal was to identify the global heteroclinic connections that lead to stochastic chaotic saddles. We use the identification of a stochastic chaotic saddle, in the presence of a bi-instability, to study how noise can induce chaos. First, we recast the derivation of the Hamiltonian system to make explicit the similarities. The scaled equations for a single-mode periodically modulated class-B laser are x y
= −y − x(a + by) = (1 + y)(x − δ cos(ωt)),
(8)
where x is the population inversion, y is the intensity, = 0.001 is the dissipation parameter, a = 56 and b = 3 are related to the pump, ω = 0.9 is the forcing frequency, and δ is the forcing amplitude. We will add noise discretely to this system once a period with respect to the forcing. This effect is important since laser systems are subject to electrical component noise and spontaneous emission noise. Notice that treating and δ as small parameters, and setting = δ = 0 results in the same reduced system as Eq. (5). The weakly perturbed Hamiltonian system can be found by using the substitutions y = −1 + ez and δ = δ1 : x z
= 1 − ez − x(a − b + bez ) = x − δ1 cos(ωt).
(9)
The system in Eq. (9) can be rewritten in its Hamiltonian form: ∂H + g1 (x, z) ∂ z¯ ∂H + g2 (t), ∂x ¯
x
= −
z
=
with H(x, z) =
1 2 x + ez − z − 1. 2
(10)
126
LORA BILLINGS, ERIK BOLLT, DAVID MORGAN, AND IRA B. SCHWARTZ
ln(I)
−10
−15
−20 0
0.05
δ
0.1
0.15
z=ln(1+y)
0
−10
−20
0.8
1
1.2
δ
1.4
1.6
Figure 1. The top graph is the bifurcation diagram for the MSI model. The bottom graph is a bifurcation diagram for the laser model. The orbits are sampled once a period. The black points show the period doubling cascade of the period one orbit. The gray points represent the coexisting attractor from the saddle-node pair.
The similarity of the MSI and laser model systems affords us the unique opportunity to compare well known results from the fields of optics with epidemiology. Bifurcations and interesting dynamics from one field should have similar counterparts in the other. 4. The topology of the models without noise. Understanding the stochastic effects requires a topological description of the underlying deterministic structure of Eqs. (1) and (8). In our study, we find that a combination of a bi-stability and bi-instability creates the manifold structure that is sufficient to induce chaos-like dynamics from noise. In both systems, when δ = 0, there is a period one orbit. As δ increases, this orbit proceeds through a period doubling cascade. The manifolds of the unstable orbits are critical in the formation of horseshoes, which eventually lead to chaotic attractors. Coexisting with this set of stable periodic orbits is a sequence of saddle-node bifurcations. As proved by a theorem in [14], systems that have a reduced Hamiltonian form as in Eq. (5) can have periodic saddle-node orbits excited by a periodic forcing at the correct amplitude. We focus on the first saddle node orbit that occurs as we increase δ from zero. Saddle node orbits of period N , which bifurcate from the orbits of the unperturbed Hamiltonian of period N , are called primary saddle-node bifurcations (PSNB) [11]. They are crucial in understanding how the global picture of the perturbed models evolves as forcing and damping parameters are increased. For the MSI model, the first PSNB is period three. For the laser model, it is period two. Numerical approximations of the bifurcation diagrams for the stable orbits are shown in Fig. 1. Notice how the saddle-node orbits go through
STOCHASTIC GLOBAL BIFURCATION IN PERTURBED HAMILTONIAN SYSTEMS
127
ln(I) −8
−12
−16
−3.0
−2.7
−2.4
ln(S)
Figure 2. A graph of the manifolds for the MSI model with δ = 0.095. Stable periodic orbits are in black. The basin of the period two orbit is white and the basin of the period three orbit is gray. The manifolds of the unstable period one orbit are solid black and the manifolds of the unstable period three orbit are dashed. Note the forward-heteroclinic crossings where the solid and dashed manifolds intersect. z 0
−5
−10
−15
−20 −6
−4
−2
0
2
4
x
Figure 3. A graph of the manifolds for the laser model with δ = 1.3. The stable periodic orbits are black, and the basins of attraction are white and gray. The manifolds of the period one orbit are solid black and the manifolds of the unstable period two orbit are dashed. a faster period doubling cascade. Eventually, these cascades are destroyed by a crisis, as studied in [11], followed by the creation of another saddle-node pair. Next, we consider parameter values in the bi-instability region to add noise. We find that if the unstable manifolds of the flip saddle, resulting from the period one bifurcation, and the stable manifolds of the regular saddle, from the saddle node pair, intersect transversely or in a way that forms the forward part of a heteroclinic tangle, noise can effectively complete the reverse crossing. In both models, we have this structure. In Fig. 2, we show the stable and unstable manifolds for the periodic orbits of the MSI model when δ = 0.095. The unstable manifolds of the period three saddle and the stable manifolds of the period one saddle cross transversally, creating a
128
LORA BILLINGS, ERIK BOLLT, DAVID MORGAN, AND IRA B. SCHWARTZ
0.2
Lyapunov Exponents
0.1
Chaos ↑
0 −0.1 −0.2 −0.3 −0.4 −0.5 −0.6 0
0.02
0.04 0.06 σ = standard deviation
0.08
0.1
Figure 4. The Lyapunov exponents of the MSI model for a random orbit near the stable period two as a function of the standard deviation, σ, for δ = 0.095. These values are the averaged Lyapunov exponents for 20 realizations of the stochastic model. forward-heteroclinic connection. Therefore, the unstable manifold of the period one saddle does not cross the stable manifold of the period three saddle. However, notice how closely it comes toward the top right of the graph. This is where we conjecture that noise completes the heteroclinic-like connection. In Fig. 3, we show the manifolds for the laser model when δ = 1.3. Here, the unstable manifolds of the period two saddle and the stable manifolds of the period one saddle cross transversally, creating a forward-heteroclinic connection. In this case, the unstable manifolds of the period one are not close to stable manifolds of the period two. Instead, we conjecture that the unstable manifold of the period two orbit forms a homoclinic-like connection, crossing its own stable manifold, with the addition of noise. 5. Adding noise to the models. Consider Eqs. (1) and (8) as discrete dynamical systems sampled at the drive periods T . Then, stochastic perturbations are added discretely on that Poincar´e section. That is, consider each stochastic system as a two-dimensional map, F , of a region D into itself (x, z) (t + T ) = Fµ [(x, z) (t)] + η(t),
(11)
where η is a two-dimensional random variable having a normal distribution given √ 2 2 by v(x) = e−(x /2σ ) / 2πσ 2 , applied once each iteration. The random part, η is assumed to be independent of state, i.i.d., which we tacitly presume to be relatively small, so that the deterministic part, F , has primary influence. In Fig. 4, we show the Lyapunov exponents for the MSI model with δ = 0.095 as a function of the standard deviation, σ. The maximum Lyapunov exponent changes sign in the interval 0.02 < σ < 0.0225 and becomes a chaotic attractor. Notice that the transition to a positive Lyapunov exponent is a smooth function of the standard deviation, implying that different parts of the phase space are becoming weighted differently for larger values of the standard deviation. Transport regions reveal how noise facilitates emerging stochastic chaos. We can also think of this as the excitation of the stochastic chaotic saddle, which were the forward heteroclinic crossings before noise was added. In [2], we developed a
STOCHASTIC GLOBAL BIFURCATION IN PERTURBED HAMILTONIAN SYSTEMS
129
ln(I) −8
−12
−16
−3.0
−2.7
−2.4
ln(S)
Figure 5. The PDF flux of the MSI model for the period two basin with σ = 0.03. Overlaid are the periodic orbits and their manifolds. The shading indicates the relative probability varying from greatest (black) to least (white) of the trajectories escaping to the other basin. Galerkin-Ulam approximation of the stochastic Frobenius-Perron operator to measure the transport across basin boundaries in Eq. (1). We continue with a quick overview of this tool. The stochastic Frobenius-Perron operator is (x−F (y))2 1 2σ 2 e− ρ(y)dy. (12) PF [ρ(x)] = √ 2πσ 2 D We then project the infinite dimensional linear space, with discretely indexed basis functions onto a finite dimensional linear subspace generated by a subset of the basis functions. This projection is realized optimally by the Galerkin method, in terms of the inner product, which we choose to be integration. We choose the basis functions to be a family of characteristic functions in a nested refinement of boxes covering D. Therefore, a Galerkin matrix entry represents how mass flows from cell to another, with column sums of one if the space is compact. We calculate the “PDF flux” by multiplying the probability density function (PDF) by the area flux. Figure 5 shows the flux for the period two basin when σ = 0.03 superimposed on the manifold schematic from Fig. 2. The darkest shading denotes the most probable regions for a trajectory to leave the period two basin and visit the period three basin. Of course the manifolds and basins boundaries do not exist in the traditional sense in the presence of noise, but these structures give an idea of the underlying dynamics. The large dark region is where trajectories near the unstable manifold of the period one can escape along the stable trajectory of the period three saddle, stochastically “completing” the heteroclinic connection. There is additional evidence from the PDF that shows the trajectories frequenting the region between the periodic attractors, where the stochastic chaotic saddle exists. In Fig. 6, the Lyapunov exponents for the laser model with δ = 1.3 are also shown as a function of σ. The maximum Lyapunov exponent changes sign near σ = 0.15, indicating the creation of a chaotic-like attractor. Noise has a slightly different effect on laser model topologically. In Fig. 7, the PDF flux for σ = 0.16
130
LORA BILLINGS, ERIK BOLLT, DAVID MORGAN, AND IRA B. SCHWARTZ
Lyapunov Exponents
0.2 Chaos ↑
0
−0.2
−0.4
−0.6 0
0.05
0.1 0.15 0.2 σ = standard deviation
0.25
0.3
Figure 6. The Lyapunov exponents of the laser model for a random orbit near the stable period two orbit as a function of the standard deviation, σ, for δ = 1.3. These values are the averaged Lyapunov exponents for 20 realizations of the stochastic model.
z 0
−5
−10
−15
−20 −6
−4
−2
0
2
4
x
Figure 7. The PDF flux of the laser model for the period two basin from the period doubling cascade with σ = 0.16 . Overlaid are the periodic orbits and their manifolds. The shading indicates the relative probability varying from greatest (black) to least (white) of the trajectories escaping to the other basin.
shows a high probability (darkest regions) that transport occurs from the basin of the period two orbit in the bottom left region. Notice how the unstable manifolds of the period one orbit are far from the stable manifolds of the period two saddles. Instead, the unstable manifold of the period two saddle is close its own stable manifold. Noise completes this homoclinic-like connection to induce chaotic-like orbits, as suggested by the PDF flux. Also, the PDF supports the claim that a stochastic chaotic saddle is excited by the addition of noise. 6. Conclusions. Although the fields of optics and population dynamics are very far apart, the MSI model and the laser model are very similar dynamically. First, they are perturbations off the same Hamiltonian system. Second, they have the
STOCHASTIC GLOBAL BIFURCATION IN PERTURBED HAMILTONIAN SYSTEMS
131
same bifurcation structure and topology. Therefore, adding stochastic perturbations of the appropriate standard deviation can induce chaos in both systems. We have discussed the homoclinic and heteroclinic connections that can be formed, and presented numerical evidence supporting these claims. Perhaps further study of the parallels in these systems will drive the search for interesting phenomena as yet undiscovered. We also remark that the laser models, appropriately modified, offers an opportunity to perform experiments on a platform that the epidemiology community cannot due to ethical barriers. If the laser can be a test bed for MSI dynamics, it could be used to test vaccination strategies that would control the spread of disease. Work on this area is now progressing to the point where global understanding of population models can be manipulated in real experiments. 7. Acknowledgments. This work was supported by the Office of Naval Research (ONR). L.B. was supported by ONR grant N00173-01-1-G911. E.B. was supported by the National Science Foundation under grant DMS-0071314. D.M. was supported by an NRC postdoctoral fellowship. REFERENCES [1] L. Arnold, Random Dynamical Systems, (Springer-Verlag, New York, 1998). [2] L. Billings, E.M. Bollt, and I.B. Schwartz, Phase-space transport of stochastic chaos in population dynamics of virus spread, Physical Review Letters 88 (2002), art. no. 234101. [3] L. Billings and I. Schwartz, Exciting chaos with noise: unexpected dynamics in epidemic outbreaks, J. Math. Biol. 44 (2002) 31–48. [4] E.M. Bollt, L. Billings, and I.B. Schwartz, A manifold independent approach to understanding transport in stochastic dynamical systems, Physica D 173 2002, 153–177. [5] T.W. Carr, L. Billings, I.B. Schwartz, and I. Triandaf, Bi-instability and the global role of unstable resonant orbits in a driven laser, Physica D 147 (2000), 59–82. [6] R. Courant and D. Hilbert, Methods of Mathematical Physics, Volume 1 (John Wiley & Sons, 1970). [7] M. Dellnitz and O. Junge, On the approximation of complicated dynamical behavior, SIAM J. Numer. Anal. 36 (1999), 491–515. [8] Y. A. Kuznetsov, Elements of Applied Bifurcation Theory, Springer-Verlag, NY 1994. [9] D.A. Rand and H.B. Wilson, Chaotic stochasticity: a ubiquitous source of unpredictability in epidemics, Proc. R. Soc. Lond. B 246 (1991), 179–184. [10] I.B. Schwartz, Multiple stable recurrent outbreaks and predictability in seasonally forced nonlinear epidemic models, J. Math. Biol. 21 (1985), 347–361. [11] I.B. Schwartz, Sequential horseshoe formation in the birth and death of chaotic attractors, Physical Review Letters 60 (1988), 1359-1362. [12] I.B. Schwartz and T. Erneux, Subharmonic hysteresis and period doubling bifurcations for periodically driven laser, SIAM J. Appl. Math. 60 (1994), 1083-1100. [13] I.B. Schwartz and T.W. Carr, Bi-instability as a precursor to global mixed-mode chaos, Phys. Rev. E 59 (1999), no. 6, 6658–6661. [14] I.B. Schwartz and H.L. Smith, Infinite subharmonic bifurcations in an SEIR epidemic model, J. Math. Biol. 18 (1983), 233–253. [15] S. Smale, Differentiable dynamical systems, Bull. Amer. Math. Soc. 73 (1967), 747–817. [16] H.L. Smith, Subharmonnic Bifurcation in an S-I-R Epidemic Model, J. Math. Biol. 17 (1983), 163-177. E-mail address: [email protected]
pp. 122–131
STOCHASTIC GLOBAL BIFURCATION IN PERTURBED HAMILTONIAN SYSTEMS LORA BILLINGS1 , ERIK BOLLT2 , DAVID MORGAN3 , AND IRA B. SCHWARTZ3 1 Department
of Mathematical Sciences, Montclair State University, Upper Montclair, NJ 07043 of Mathematics & Computer Science and Department of Physics, Clarkson University, Potsdam, NY 13699 3 Naval Research Laboratory, Code 6792, Plasma Physics Division, Washington, DC 20375 2 Departments
Abstract. We study two perturbed Hamiltonian systems in which chaos-like dynamics can be induced by stochastic perturbations. We show the similarities of a class of population and laser models, analytically and topologically. Both systems have similar manifold structure that includes bi-instability and partially formed heteroclinic connections. Noise takes advantage of this structure, inducing a global bifurcation and chaotic-like dynamics which exhibits mixed mode behavior of the original bi-stable solutions. We support these claims with numerical approximations of the transport between basins.
1. Introduction. In this article, we describe how two ordinary differential equation systems, from population dynamics and lasers, have the manifold structure sufficient for exciting noise induced chaos. Though coming from very different fields, we find that the single-mode periodically modulated class-B laser [12] and the seasonally forced nonlinear epidemic model [14] can both be reduced to weakly perturbed Hamiltonian systems, having equivalent Hamiltonians. The perturbations result in the existence of saddle-node bifurcations and lead to a bi-instability, or the existence of at least two unstable periodic saddles. When the manifolds of these saddles are close to a homoclinic or heteroclinic topology, periodic stochastic perturbations, or noise, can essentially complete the intersections to form a topological horseshoe, in a metaphorical sense. Since it is not sufficient to compute only Lyapunov exponents to detect or define stochastic chaos [9], we compute transport of a stochastic flux from one basin to another [7, 6]. This is done using the Galerkin approximation of the stochastic Frobenius-Perron operator [4]. Increasing the standard deviation of the noise reduces the stability of the stable periodic attractors and increases the frequency of a trajectory visiting the neighborhood of the induced homoclinic or heteroclinic tangle. It is well known in the noiseless case that homoclinic or heteroclinic tangles create the horseshoe topology that allows chaotic dynamics [15]. The stochastic perturbations enable transport in a similar fashion, and we call this phenomena noise-induced chaos. We can numerically approximate where pseudo-crossings are induced by identifying the transport between the noiseless basins of attraction of two attracting orbits. Therefore, we can conjecture from numerical evidence how both models undergo a smooth transition to chaos by increasing the standard deviation of the noise. 1991 Mathematics Subject Classification. 34F05, 65P20, 34D08, 37M20. Key words and phrases. Stochastic dynamical systems, Hamiltonian, bifurcation, chaos.
122
STOCHASTIC GLOBAL BIFURCATION IN PERTURBED HAMILTONIAN SYSTEMS
123
We begin by presenting the two systems and their reduction to the same weakly perturbed Hamiltonian system. Then, we describe their manifold structures and conjecture how noise induces chaotic behavior. Finally, we present numerical approximations of the Lyapunov exponents and transport, or flux, in the system due to noise. This is similar to the phenomenal bifurcation described in [1] in which the standard deviation of the noise changes the Lyapunov exponent (integrated with respect to the invariant density) from negative to positive. The conclusion is that noise induces enough transport across basin boundaries that global effects contribute to the Lyapunov exponent bifurcation of stochastic chaos. 2. A model from population dynamics. Consider a compartmental model to describe the dynamics of how measles spreads in a large population. The population can be divided into disjoint subgroups which evolve in time. Based on singular perturbation theory as presented in [14, 10], the standard SEIR model is reduced to two equations. We call this the modified SI model (MSI): S (t) I (t)
= µ − µS(t) − β(t)I(t)S(t) α = β(t)I(t)S(t) − (µ + α)I(t). µ+γ
(1)
The people that may contract the disease are susceptible, S(t). Those individuals capable of transmitting the diseases are infective, I(t). This contact rate fluctuates with the seasons and can be approximated by a sinusoidal forcing of the form β(t) = β0 (1 + δ cos 2πt), with 0 ≤ δ < 1. The parameter values are fixed at µ = 0.02 (year)−1 , α = 35.84 (year)−1 , γ = 100 (year)−1 , β0 = 1575 (year)−1 , and we vary δ, the fluctuating contact rate amplitude. These are the standard parameters used in [10]. We will add noise discretely, not continuously, to both S and I once a period with respect to the forcing to represent fluctuations in the population [3]. When δ = 0, there are two steady states: the disease free equilibrium, (S, I) = (1, 0), and the endemic steady state, (µ + α)(µ + γ) µ αβ0 −1 . , (Se , Ie ) = αβ0 β0 (µ + α)(µ + γ) The stability of these steady states depends on the value of Q, where Q=
αβ0 . (µ + α)(µ + γ)
If Q < 1, the disease free equilibrium is stable, without any endemic state. When Q > 1, the disease free state is unstable, and the endemic state is attracting. Since we wish to understand the dynamics when the disease is prevalent, we will assume Q > 1. With a change of variables, we recenter the endemic equilibrium at the origin and the disease free equilibrium at infinity. Using S = Se (1 + x) and I = Ie (1 + y), the new system is x
=
y
=
S = µQ − µ(1 + x) − µ(1 + δ cos 2πt)(Q − 1)(1 + y)(1 + x) Se I = (µ + α)(1 + y)(δ cos 2πt + (1 + δ cos 2πt)x). Ie
(2)
124
LORA BILLINGS, ERIK BOLLT, DAVID MORGAN, AND IRA B. SCHWARTZ
To check the stability of the origin when δ = 0, we find eigenvalues of the Jacobian of Eq. (2). −µQ −µ(Q − 1) J|(0,0) = µ+α 0 µQ 1 λ± = − ± µ2 Q2 − 4µ(Q − 1)(µ + α). 2 2 1 1 Similar to an approach in [16, 14], recall the fact that µ, µ+γ , and µ+α are O(10−2 ). We introduce the small parameter = µ(Q − 1), and µ + α = ∆ . Using the Q substitution η = Q−1 , 4∆ − η 2 2 η λ± = − ± i. 2 2 The eigenvalues are in the √ form λ± = −r ± νi. To the leading order, the frequency of the oscillation is ν = ∆, and r = η/2. This is a weakly stable spiral sink. One last change of variables will transform the system so that when δ = 0, the linear part will be close to its Jordan form. We also assume that there is a small amplitude periodic perturbation, requiring δ/ to be small. Let δ¯ = δ/, x ¯ = ν x − y, and y¯ = y in Eq. (2). The new system is x ¯
= −ν y¯ − (¯ x(η − 1) + (1 + ν)(1 + y¯)(¯ x + ν δ¯ cos 2πt)) − 2 (¯ y (η − 1) 3 ¯ ¯ +(1 + ν)(1 + y¯)(¯ y+x ¯δ cos 2πt)) − (1 + ν)(1 + y¯)(δ cos 2πt)
y¯
= νx ¯(1 + y¯) + ν 2 (1 + y¯)(δ¯ cos 2πt) +ν(1 + y¯)(¯ y+x ¯δ¯ cos 2πt) + 2 ν y¯(1 + y¯)(δ¯ cos 2πt).
(3)
We can rewrite Eq. (3) as x ¯ y¯
¯ = −ν y¯ + f1 (¯ x, y¯, t, , δ) ¯ 2 (1 + y¯)(cos 2πt) + f2 (¯ ¯ = νx ¯(1 + y¯) + δν x, y¯, t, , δ),
(4)
with ¯ = fi (¯ ¯ fi (¯ x, y¯, t, , δ) x, y¯, t + 1, , δ) fi (0, 0, t, 0, 0) = 0, i = 1, 2. ¯ ¯ Now treat and δ as small parameters. Setting = δ = 0 in Eq. (4), we find the following reduced system, x ¯
y¯
= −ν y¯ = νx ¯(1 + y¯).
(5)
This is the same system that was found in earlier work on both the SEIR and SIR systems [14]. The system is conservative and the solution is described by level curves of the integral function E = x2 + 2y − 2 ln |1 + y|. It is interesting to note that the structure of such a first integral results from unfolding a degenerate Hopf bifurcation [8]. Moreover, as a result of the quadratic nonlinearities, the orbits arising the Hamiltonian may be unbounded. We can also recast Eq. (3) as a weakly perturbed Hamiltonian system by making one more change of variables. Let y¯ = −1 + ez¯ and drop terms of order higher than to find x(η − 1) + ez¯(1 + ν)(¯ x + ν δ¯ cos 2πt)) x ¯ = −ν(−1 + ez¯) − (¯ z¯
¯ 2 cos 2πt + ν(−1 + ez¯ + x = νx ¯ + δν ¯δ¯ cos 2πt).
(6)
STOCHASTIC GLOBAL BIFURCATION IN PERTURBED HAMILTONIAN SYSTEMS
125
The system in Eq. (6) can be rewritten as ∂H ¯ + f˜1 (¯ x, z¯, t, , δ) ∂ z¯ ∂H ¯ 2 ¯ + δν cos 2πt + f˜2 (¯ x, z¯, t, , δ), ∂x ¯
x ¯
= −
z¯
=
(7)
with 1 2 ¯ + ez¯ − z¯ − 1). H(¯ x, z¯) = ν( x 2 We will show below the similarity of this Hamiltonian system to the one we find for the following laser model. 3. The laser model example. Uncovered in a study of global mixed-mode chaos [5, 13], it was shown that the modulated CO2 laser model also reduces to a weakly perturbed Hamiltonian system similar to Eq. (6). They used approximate maps, derived by matched asymptotic expansions and subharmonic Melnikov theory, to prove the existence of period doubling and saddle-node bifurcations. Existence of the subharmonic Melnikov function also allowed for additional conclusions concerning the global dynamics and transverse manifold intersections. The goal was to identify the global heteroclinic connections that lead to stochastic chaotic saddles. We use the identification of a stochastic chaotic saddle, in the presence of a bi-instability, to study how noise can induce chaos. First, we recast the derivation of the Hamiltonian system to make explicit the similarities. The scaled equations for a single-mode periodically modulated class-B laser are x y
= −y − x(a + by) = (1 + y)(x − δ cos(ωt)),
(8)
where x is the population inversion, y is the intensity, = 0.001 is the dissipation parameter, a = 56 and b = 3 are related to the pump, ω = 0.9 is the forcing frequency, and δ is the forcing amplitude. We will add noise discretely to this system once a period with respect to the forcing. This effect is important since laser systems are subject to electrical component noise and spontaneous emission noise. Notice that treating and δ as small parameters, and setting = δ = 0 results in the same reduced system as Eq. (5). The weakly perturbed Hamiltonian system can be found by using the substitutions y = −1 + ez and δ = δ1 : x z
= 1 − ez − x(a − b + bez ) = x − δ1 cos(ωt).
(9)
The system in Eq. (9) can be rewritten in its Hamiltonian form: ∂H + g1 (x, z) ∂ z¯ ∂H + g2 (t), ∂x ¯
x
= −
z
=
with H(x, z) =
1 2 x + ez − z − 1. 2
(10)
126
LORA BILLINGS, ERIK BOLLT, DAVID MORGAN, AND IRA B. SCHWARTZ
ln(I)
−10
−15
−20 0
0.05
δ
0.1
0.15
z=ln(1+y)
0
−10
−20
0.8
1
1.2
δ
1.4
1.6
Figure 1. The top graph is the bifurcation diagram for the MSI model. The bottom graph is a bifurcation diagram for the laser model. The orbits are sampled once a period. The black points show the period doubling cascade of the period one orbit. The gray points represent the coexisting attractor from the saddle-node pair.
The similarity of the MSI and laser model systems affords us the unique opportunity to compare well known results from the fields of optics with epidemiology. Bifurcations and interesting dynamics from one field should have similar counterparts in the other. 4. The topology of the models without noise. Understanding the stochastic effects requires a topological description of the underlying deterministic structure of Eqs. (1) and (8). In our study, we find that a combination of a bi-stability and bi-instability creates the manifold structure that is sufficient to induce chaos-like dynamics from noise. In both systems, when δ = 0, there is a period one orbit. As δ increases, this orbit proceeds through a period doubling cascade. The manifolds of the unstable orbits are critical in the formation of horseshoes, which eventually lead to chaotic attractors. Coexisting with this set of stable periodic orbits is a sequence of saddle-node bifurcations. As proved by a theorem in [14], systems that have a reduced Hamiltonian form as in Eq. (5) can have periodic saddle-node orbits excited by a periodic forcing at the correct amplitude. We focus on the first saddle node orbit that occurs as we increase δ from zero. Saddle node orbits of period N , which bifurcate from the orbits of the unperturbed Hamiltonian of period N , are called primary saddle-node bifurcations (PSNB) [11]. They are crucial in understanding how the global picture of the perturbed models evolves as forcing and damping parameters are increased. For the MSI model, the first PSNB is period three. For the laser model, it is period two. Numerical approximations of the bifurcation diagrams for the stable orbits are shown in Fig. 1. Notice how the saddle-node orbits go through
STOCHASTIC GLOBAL BIFURCATION IN PERTURBED HAMILTONIAN SYSTEMS
127
ln(I) −8
−12
−16
−3.0
−2.7
−2.4
ln(S)
Figure 2. A graph of the manifolds for the MSI model with δ = 0.095. Stable periodic orbits are in black. The basin of the period two orbit is white and the basin of the period three orbit is gray. The manifolds of the unstable period one orbit are solid black and the manifolds of the unstable period three orbit are dashed. Note the forward-heteroclinic crossings where the solid and dashed manifolds intersect. z 0
−5
−10
−15
−20 −6
−4
−2
0
2
4
x
Figure 3. A graph of the manifolds for the laser model with δ = 1.3. The stable periodic orbits are black, and the basins of attraction are white and gray. The manifolds of the period one orbit are solid black and the manifolds of the unstable period two orbit are dashed. a faster period doubling cascade. Eventually, these cascades are destroyed by a crisis, as studied in [11], followed by the creation of another saddle-node pair. Next, we consider parameter values in the bi-instability region to add noise. We find that if the unstable manifolds of the flip saddle, resulting from the period one bifurcation, and the stable manifolds of the regular saddle, from the saddle node pair, intersect transversely or in a way that forms the forward part of a heteroclinic tangle, noise can effectively complete the reverse crossing. In both models, we have this structure. In Fig. 2, we show the stable and unstable manifolds for the periodic orbits of the MSI model when δ = 0.095. The unstable manifolds of the period three saddle and the stable manifolds of the period one saddle cross transversally, creating a
128
LORA BILLINGS, ERIK BOLLT, DAVID MORGAN, AND IRA B. SCHWARTZ
0.2
Lyapunov Exponents
0.1
Chaos ↑
0 −0.1 −0.2 −0.3 −0.4 −0.5 −0.6 0
0.02
0.04 0.06 σ = standard deviation
0.08
0.1
Figure 4. The Lyapunov exponents of the MSI model for a random orbit near the stable period two as a function of the standard deviation, σ, for δ = 0.095. These values are the averaged Lyapunov exponents for 20 realizations of the stochastic model. forward-heteroclinic connection. Therefore, the unstable manifold of the period one saddle does not cross the stable manifold of the period three saddle. However, notice how closely it comes toward the top right of the graph. This is where we conjecture that noise completes the heteroclinic-like connection. In Fig. 3, we show the manifolds for the laser model when δ = 1.3. Here, the unstable manifolds of the period two saddle and the stable manifolds of the period one saddle cross transversally, creating a forward-heteroclinic connection. In this case, the unstable manifolds of the period one are not close to stable manifolds of the period two. Instead, we conjecture that the unstable manifold of the period two orbit forms a homoclinic-like connection, crossing its own stable manifold, with the addition of noise. 5. Adding noise to the models. Consider Eqs. (1) and (8) as discrete dynamical systems sampled at the drive periods T . Then, stochastic perturbations are added discretely on that Poincar´e section. That is, consider each stochastic system as a two-dimensional map, F , of a region D into itself (x, z) (t + T ) = Fµ [(x, z) (t)] + η(t),
(11)
where η is a two-dimensional random variable having a normal distribution given √ 2 2 by v(x) = e−(x /2σ ) / 2πσ 2 , applied once each iteration. The random part, η is assumed to be independent of state, i.i.d., which we tacitly presume to be relatively small, so that the deterministic part, F , has primary influence. In Fig. 4, we show the Lyapunov exponents for the MSI model with δ = 0.095 as a function of the standard deviation, σ. The maximum Lyapunov exponent changes sign in the interval 0.02 < σ < 0.0225 and becomes a chaotic attractor. Notice that the transition to a positive Lyapunov exponent is a smooth function of the standard deviation, implying that different parts of the phase space are becoming weighted differently for larger values of the standard deviation. Transport regions reveal how noise facilitates emerging stochastic chaos. We can also think of this as the excitation of the stochastic chaotic saddle, which were the forward heteroclinic crossings before noise was added. In [2], we developed a
STOCHASTIC GLOBAL BIFURCATION IN PERTURBED HAMILTONIAN SYSTEMS
129
ln(I) −8
−12
−16
−3.0
−2.7
−2.4
ln(S)
Figure 5. The PDF flux of the MSI model for the period two basin with σ = 0.03. Overlaid are the periodic orbits and their manifolds. The shading indicates the relative probability varying from greatest (black) to least (white) of the trajectories escaping to the other basin. Galerkin-Ulam approximation of the stochastic Frobenius-Perron operator to measure the transport across basin boundaries in Eq. (1). We continue with a quick overview of this tool. The stochastic Frobenius-Perron operator is (x−F (y))2 1 2σ 2 e− ρ(y)dy. (12) PF [ρ(x)] = √ 2πσ 2 D We then project the infinite dimensional linear space, with discretely indexed basis functions onto a finite dimensional linear subspace generated by a subset of the basis functions. This projection is realized optimally by the Galerkin method, in terms of the inner product, which we choose to be integration. We choose the basis functions to be a family of characteristic functions in a nested refinement of boxes covering D. Therefore, a Galerkin matrix entry represents how mass flows from cell to another, with column sums of one if the space is compact. We calculate the “PDF flux” by multiplying the probability density function (PDF) by the area flux. Figure 5 shows the flux for the period two basin when σ = 0.03 superimposed on the manifold schematic from Fig. 2. The darkest shading denotes the most probable regions for a trajectory to leave the period two basin and visit the period three basin. Of course the manifolds and basins boundaries do not exist in the traditional sense in the presence of noise, but these structures give an idea of the underlying dynamics. The large dark region is where trajectories near the unstable manifold of the period one can escape along the stable trajectory of the period three saddle, stochastically “completing” the heteroclinic connection. There is additional evidence from the PDF that shows the trajectories frequenting the region between the periodic attractors, where the stochastic chaotic saddle exists. In Fig. 6, the Lyapunov exponents for the laser model with δ = 1.3 are also shown as a function of σ. The maximum Lyapunov exponent changes sign near σ = 0.15, indicating the creation of a chaotic-like attractor. Noise has a slightly different effect on laser model topologically. In Fig. 7, the PDF flux for σ = 0.16
130
LORA BILLINGS, ERIK BOLLT, DAVID MORGAN, AND IRA B. SCHWARTZ
Lyapunov Exponents
0.2 Chaos ↑
0
−0.2
−0.4
−0.6 0
0.05
0.1 0.15 0.2 σ = standard deviation
0.25
0.3
Figure 6. The Lyapunov exponents of the laser model for a random orbit near the stable period two orbit as a function of the standard deviation, σ, for δ = 1.3. These values are the averaged Lyapunov exponents for 20 realizations of the stochastic model.
z 0
−5
−10
−15
−20 −6
−4
−2
0
2
4
x
Figure 7. The PDF flux of the laser model for the period two basin from the period doubling cascade with σ = 0.16 . Overlaid are the periodic orbits and their manifolds. The shading indicates the relative probability varying from greatest (black) to least (white) of the trajectories escaping to the other basin.
shows a high probability (darkest regions) that transport occurs from the basin of the period two orbit in the bottom left region. Notice how the unstable manifolds of the period one orbit are far from the stable manifolds of the period two saddles. Instead, the unstable manifold of the period two saddle is close its own stable manifold. Noise completes this homoclinic-like connection to induce chaotic-like orbits, as suggested by the PDF flux. Also, the PDF supports the claim that a stochastic chaotic saddle is excited by the addition of noise. 6. Conclusions. Although the fields of optics and population dynamics are very far apart, the MSI model and the laser model are very similar dynamically. First, they are perturbations off the same Hamiltonian system. Second, they have the
STOCHASTIC GLOBAL BIFURCATION IN PERTURBED HAMILTONIAN SYSTEMS
131
same bifurcation structure and topology. Therefore, adding stochastic perturbations of the appropriate standard deviation can induce chaos in both systems. We have discussed the homoclinic and heteroclinic connections that can be formed, and presented numerical evidence supporting these claims. Perhaps further study of the parallels in these systems will drive the search for interesting phenomena as yet undiscovered. We also remark that the laser models, appropriately modified, offers an opportunity to perform experiments on a platform that the epidemiology community cannot due to ethical barriers. If the laser can be a test bed for MSI dynamics, it could be used to test vaccination strategies that would control the spread of disease. Work on this area is now progressing to the point where global understanding of population models can be manipulated in real experiments. 7. Acknowledgments. This work was supported by the Office of Naval Research (ONR). L.B. was supported by ONR grant N00173-01-1-G911. E.B. was supported by the National Science Foundation under grant DMS-0071314. D.M. was supported by an NRC postdoctoral fellowship. REFERENCES [1] L. Arnold, Random Dynamical Systems, (Springer-Verlag, New York, 1998). [2] L. Billings, E.M. Bollt, and I.B. Schwartz, Phase-space transport of stochastic chaos in population dynamics of virus spread, Physical Review Letters 88 (2002), art. no. 234101. [3] L. Billings and I. Schwartz, Exciting chaos with noise: unexpected dynamics in epidemic outbreaks, J. Math. Biol. 44 (2002) 31–48. [4] E.M. Bollt, L. Billings, and I.B. Schwartz, A manifold independent approach to understanding transport in stochastic dynamical systems, Physica D 173 2002, 153–177. [5] T.W. Carr, L. Billings, I.B. Schwartz, and I. Triandaf, Bi-instability and the global role of unstable resonant orbits in a driven laser, Physica D 147 (2000), 59–82. [6] R. Courant and D. Hilbert, Methods of Mathematical Physics, Volume 1 (John Wiley & Sons, 1970). [7] M. Dellnitz and O. Junge, On the approximation of complicated dynamical behavior, SIAM J. Numer. Anal. 36 (1999), 491–515. [8] Y. A. Kuznetsov, Elements of Applied Bifurcation Theory, Springer-Verlag, NY 1994. [9] D.A. Rand and H.B. Wilson, Chaotic stochasticity: a ubiquitous source of unpredictability in epidemics, Proc. R. Soc. Lond. B 246 (1991), 179–184. [10] I.B. Schwartz, Multiple stable recurrent outbreaks and predictability in seasonally forced nonlinear epidemic models, J. Math. Biol. 21 (1985), 347–361. [11] I.B. Schwartz, Sequential horseshoe formation in the birth and death of chaotic attractors, Physical Review Letters 60 (1988), 1359-1362. [12] I.B. Schwartz and T. Erneux, Subharmonic hysteresis and period doubling bifurcations for periodically driven laser, SIAM J. Appl. Math. 60 (1994), 1083-1100. [13] I.B. Schwartz and T.W. Carr, Bi-instability as a precursor to global mixed-mode chaos, Phys. Rev. E 59 (1999), no. 6, 6658–6661. [14] I.B. Schwartz and H.L. Smith, Infinite subharmonic bifurcations in an SEIR epidemic model, J. Math. Biol. 18 (1983), 233–253. [15] S. Smale, Differentiable dynamical systems, Bull. Amer. Math. Soc. 73 (1967), 747–817. [16] H.L. Smith, Subharmonnic Bifurcation in an S-I-R Epidemic Model, J. Math. Biol. 17 (1983), 163-177. E-mail address: [email protected]
