## Abstract

In this manuscript, we consider generalizations of the classical Mathieu's equation to stochastic systems. Unlike previous works, we focus on internal frequencies that vary continuously between periodic and stochastic variables. By numerically integrating the system of equations using a symplectic method, we determine the Lyapunov exponents for a wide range of parameters to quantify how the growth rates vary in parameter space.

In the nearly periodic limit, we recover the same growth rates as the classical Mathieu's equation. As the stochasticity increases, the maximum growth decreases and the unstable region broadens, which signifies that the stochasticity can both stabilize and destabilize a system in relation to the deterministic limit. In addition to the classical parametric modes, there is also an unstable stochastic mode that arises only for moderate stochasticity. The power spectrum of the solutions shows that an increase in stochasticity tends to narrow the width of the subharmonic peak and increase the decay away from this peak.

## 1. Introduction

Mathieu's equation is an excellent paradigm to study parametric (subharmonic) resonance (or instability), henceforth referred to as PR, because it exemplifies the response of a system to the kind of periodic forcing found in many natural systems. In textbooks on differential equations, it is often derived by perturbing a nonlinear pendulum where the pivot point is oscillated periodically in the vertical direction (Stoker 1950; Nayfeh & Mook 1995). Mathieu's equation also arises in the study of celestial mechanics (Binney 1978, 1981).

In fluid dynamics, we can find many examples of waves being excited by PR. Benjamin & Ursell (1954) determined that Mathieu's equation describes the linear stability of Faraday waves that occur when a container of liquid is oscillated periodically in the vertical direction. The interesting patterns that develop at the surface were first observed by Faraday (1831). The linear stability of Kelvin–Helmholtz instability with time-periodic shear is approximately described by Mathieu's equation (Kelly 1967). Other examples can be found in acoustic instabilities in flames (Sammarco *et al*. 1997) and Rayleigh–Bénard convection (Hall & Seminara 1980; Vittori 1998). In geophysical fluid dynamics, Poulin *et al*. (2003) and Pedlosky & Thomson (2003) derived ordinary differential equations of the Hill type to describe the stability of barotropic and baroclinic shear flows, respectively. Even though Mathieu's equation is a single example of the more general Hill's differential equation, it provides a nice framework to study PR while remaining relevant to a plethora of fluid dynamical phenomena.

Mathieu's equation is also important in describing travelling waves in a spatially uniform medium. For example, if you consider the simple linear, non-dissipative, non-dispersive, two-directional wave equation with wave speed of *c* that is possibly a function of time, it is possible to find a series solution using a Galerkin expansion. (The interested reader is directed to Fauve (1998) for more details about how this approach is applied to the nonlinear Sine–Gordon equation.) The temporal equation that determines the evolution of each wave component is of the form of equation (1.2), where *x* is the displacement from the equilibrium position and the coefficient in brackets is precisely the square of the wave speed *c*^{2}. This yields that the generalized Mathieu's equation can be interpreted as governing the evolution of the amplitude of a particular wave moving through a medium that is varying stochastically in time.

The scenario where the wave speed varies in space is more complicated. Deych *et al*. (1998) solved this problem numerically using Monte Carlo simulations and computed the Lyapunov exponents to determine the stability of waves moving through a random medium. The effect of the randomness can diffuse the incoming wave, as we might expect, however, it can also be strengthened, a process known as Anderson localization (del Rio *et al*. 1995). This process was first discovered in quantum mechanical systems but it has since been determined that it is relevant for many other types of waves.

The linear Mathieu's equation with dissipation is(1.1)where in the classical theory . The parameters *ν*, *ω*_{0}, *ϵ* and *Ω* are the damping coefficient, natural frequency of oscillation, amplitude of the oscillation and frequency of the internal oscillation, respectively. As is well known, if we substitute into the equation, the governing equation for *x*(*t*) is the undamped classical Mathieu's equation,(1.2)where is the square of the natural frequency of the new system. This indicates that *ν* has the effect of not only damping the system but also modifying the natural frequency of the system. For the damped system to be unstable, it is necessary that the growth rate in equation (1.2) exceeds *ν*. Stratonovich (1967) explored the case where both the damping coefficient and the natural frequency are random variables. He concludes that the parametric instability can occur as long as the spectrum of the frequency has a subharmonic component that exceeds the mean of the random damping. Since stochastic damping does not change the behaviour of the system qualitatively and non-zero deterministic damping simply reduces the growth rate by *ν*, we restrict our attention to the case of no damping.

In the case of a constant natural frequency, it is well known that if *δ*>0(<0) the solutions are oscillatory (exponential), which denotes the fact that the pendulum at the bottom (inverted) position is stable (unstable). If the frequency is time dependent then it is possible for the stability of each of these positions to change for certain choices of the forcing frequency and amplitude. This mechanism is referred to as PR. The transition frequencies are , where *n* is any natural number and *Ω* is the frequency of *ξ*(*t*).

A very efficient way to compute the growth rate of any time-periodic system is to use Floquet theory in conjunction with numerical integration. For the classical Mathieu's equation it is easily found that there are different regions of instability, the so-called Arnold instability tongues (Nayfeh & Mook 1995), which are as illustrated in figure 1. Each integer *n* corresponds to an unstable region and (*n*/2)*Ω* determines the frequency of oscillation of the modes in that tongue as the instability grows. As can be shown theoretically (Nayfeh & Mook 1995), for the undamped case all the unstable regions touch the *x*-axis. This does not happen in figure 1 owing to numerical error. The geometric properties of the stability of equation (1.2) are addressed by Broer & Levi (1995).

To understand the basic mechanism of PR consider equation (1.1) with *ν*=0, and the total frequency is piecewise constant with period *T* defined by(1.3)Farrell & Ioannou (1996*a*,*b*) explain PR in terms of transient growth. When the solution begins the first phase, the non-normality (by which we mean that the eigenvectors are not all orthogonal) induces transient growth. If the frequency is left unchanged the growth will subside and decay away. However, if the frequency changes to another constant value before the transient growth has nearly vanished, the second phase can also have transient growth. If there is any net growth at the end of one period, this will give rise to exponential growth since it occurs at every period. It should be appreciated that PR occurs only for a very particular choice of parameters. Farrell & Ioannou (1996*a*,*b*) state that there are two essential requirements for PR: the system varies in time and the background state is non-normal. As the non-normality increases, so do the growth rates due to PR. An alternative way to quantify transient growth is to consider the pseudo-spectra (Trefethen *et al*. 1993).

The example in equation (1.3) is used by Lindenberg & West (1990) to explain PR from a more geometric approach using phase space. When the frequency is constant, there are trajectories in phase space that form closed elliptical orbits. For the first half of the total period, the trajectory is along an orbit associated with the frequency *Ω*_{1} that passes through the initial phase. When the switch occurs, the phase then continues along a different elliptical orbit associated with the other frequency *Ω*_{2}. After one complete period the net displacement can be away from the centre, which signifies that the solution has acquired energy. Since this acquisition of energy occurs every period an instability results.

Periodic time dependency is very common in nature; however, it is only a special example of how a system can change temporally. This is why it is of great interest to understand the stability of Mathieu's equation with respect to non-periodic time variations since they are more typical. Some criteria for boundedness of solutions to the quasi-periodic Mathieu's equation are obtained by Levi & Zehnder (1995).

There are many different extensions of Mathieu's equation to include stochasticity. One of the simplest is by Van Kampen (2001), which considers equation (1.2) with *δ*=0 and *ξ*(*t*) as a Gaussian random variable with zero mean and *ϵ* as the variance of the internal forcing. By using Bourret's approximation, Van Kampen shows that the second moments can grow exponentially in time giving rise to an instability. However, this method does not give any insight into whether there is a dominant frequency for this unstable mode. This is a characteristic that we feel important since it helps to classify the type of growth.

There are more recent works that studied equation (1.2) for the case of coloured Gaussian noise to include the effects of memory (Bobryk & Chrzeszczyk 2002, 2004; Lee *et al.* 2004). There have also been studies of a stochastic Mathieu's equation using non-Gaussian noise. For instance Zhang *et al*. (1998) used narrow-band Gaussian noise to study the parametric oscillator that is relevant to a microgravity environment. As well, Adams & Bloch (2008) studied a Mathieu's equation that changed the intrinsic period and strength after a complete cycle. Most of the studies of a stochastic Mathieu's equation have assumed white or coloured noise with small correlation time, which means that it is nearly white. We instead consider the alternate limit of extending Mathieu's equation to the case where the internal noise is slightly or moderately stochastic in order to understand the transition better from a periodic to a noisy variation in the natural frequency. We do this by studying different stochastic variables for which, in the limit of infinitesimal stochasticity, we recover the classical Mathieu's equation.

The outline of this article is as follows. In §2 we define the stochastic variables and their statistics. In §3 we discuss the Lyapunov exponent and why we choose this measure to compute the stability of a stochastic differential equation. Section 4 explains the symplectic numerical method we used to integrate equation (1.2) and why this is more desirable than higher order non-symplectic methods. In §5 we present the growth rate calculations using the different stochastic variables and discuss the effect an increase in stochasticity has on the stability of the system. In §6 we present the power spectra for several case examples to illustrate the dominant frequencies and their shape in each of the unstable regions and, finally, in §7 we conclude our findings.

## 2. Stochastic processes

We present two different types of stochastic processes that can be incorporated into a generalized Mathieu's equation. They are chosen because they have autocorrelation functions that are oscillatory (but damped) unlike the Gaussian or bounded uniform distributions (Gardiner 2004). We define each type of noise, its autocorrelation and its power spectra where possible.

### (a) Narrow-band Gaussian noise

The narrow-band Gaussian noise used by Zhang *et al*. (1998) was shown to exhibit parametric instability when introduced into Mathieu's equation. It is defined by(2.1)where *Ω* is the natural frequency of the forcing and *S*_{1}(*t*) and *S*_{2}(*t*) are two Ornstein–Uhlenbeck (OU) processes generated by Gaussian noise with variance *σ*^{2} and correlation time *τ*_{c}. In particular, they are solutions to the following equation:(2.2)for *i*=1, 2. In contrast to Zhang *et al*. (1998) who took the initial conditions for *S*_{1} and *S*_{2} to be zero, we take both the initial conditions to be unity. The motivation is that for time scales shorter than the correlation time, the OU process will be nearly constant, and non-zero, which will yield that equation (2.1) is sinusoidal. When the time scales are larger than *τ*_{c}, the internal noise will no longer be sinusoidal in nature and thus equation (1.2) will not approach the classical Mathieu's equation.

This implies that the mean and autocorrelation functions of *S*_{i} for *i*=1, 2 are,(2.3)(2.4)respectively. In the limit where *t* and *t*′ are much longer than the correlation time, we get the approximate result that(2.5)The autocorrelation function of the noise *ξ*(*t*) is(2.6)where we choose *σ*^{2}=1/2 to ensure that the variance of the solution for long time is a half, the same as the cosine function. The corresponding power spectrum is simply the Fourier transform of the autocorrelation function (Gardiner 2004), which in the limit of long correlation time is(2.7)This spectrum has two peaks at *ω*=±*Ω*. When we restrict our attention to positive frequencies, we find that there should be a peak in the internal noise at the harmonic frequency. We also studied the analogue of equation (2.1) where the OU noises were generated by bounded uniform noise instead of Gaussian noise. The noise is very similar and hence we do not show the details.

### (b) Forced dissipative random oscillator

An alternative type of noise can be generated by solving the equation that describes the dynamics of a dissipative random oscillator that is forced by a complex stochastic variable *r*(*t*),(2.8)The solution depends on *τ*_{c}, *Ω* and *σ*^{2}, which are the correlation time, dominant frequency and variances of the noise, respectively. We pick *r*(*t*) to be uniform bounded white noise restricted to the unit circle with the following statistics:(2.9)The weak form of equation (2.8) for the case of zero initial condition is(2.10)which has a zero mean owing to equation (2.9).

The autocorrelation function is(2.11)We substitute the autocorrelation of the complex noise given in equation (2.9) into equation (2.11), compute the integral and obtain that the real part is(2.12)The corresponding power spectrum in the long time limit is(2.13)which is equivalent to the power spectrum for narrow-band Gaussian noise. For shorter times, the noise is not stationary and we cannot write down the corresponding power spectrum. We pick to ensure that the variance for large time asymptotes to a half, as in the classical Mathieu's equation.

### (c) Kubo oscillator generated by an OU process

As before we assume that *η*(*t*) is Gaussian white noise with mean 0 and variance 1. Next, we define an OU process *S*(*t*) that satisfies equation (2.2) with an initial condition of zero. This yields that the mean and autocorrelation functions are(2.14)

(2.15)

The Kubo oscillator (Risken 1984; Lindenberg & West 1990) is a solution to the following stochastic differential equation:(2.16)where *S*(*t*) is the OU process defined by equation (2.2) and *Ω* is a constant reference frequency. If we restrict out attention to real solutions, then the solution can formally be written as(2.17)This could, equivalently, be taken as the definition of the stochastic process and it is this form that we use to generate the internal noise.

The mean and autocorrelation functions of the solution when time is much longer than the correlation time are then(2.18)(2.19)We are unable to obtain a nice closed-form expression for the associated power spectrum since the noise is non-stationary.

### (d) Kubo oscillator generated by a Wiener process

The final type of noise is a Kubo oscillator generated by a Wiener process. To construct it explicitly we assume that *η*(*t*) is defined as before and recall that a Wiener process can be defined in one of two ways:(2.20)where we have assumed a zero initial condition. The corresponding mean and autocorrelation functions are (Risken 1984)(2.21)(2.22)If we replace *S*(*t*) with *W*(*t*) in either equation (2.16) or (2.17), then we get a Kubo oscillator generated by a Wiener process. This could also be obtained by taking the solution of the previous subsection and assuming that *t*/*τ*_{c} is very large. The corresponding mean and autocorrelation functions, for the solution *ξ*(*t*) given by equation (2.17) generated by a Wiener process, are(2.23)(2.24)where and . Again, the non-stationarity of the solution prevents us from obtaining an analytical expression for the power spectrum.

## 3. Lyapunov exponent

When analysing the stability of a system, it is very important to determine what measure we use to indicate stability. Even though the first moment cannot be used as a sufficient condition for instability, it can yield a necessary condition (Van Kampen 2001; Poulin & Scott 2005). A popular measure for the instability of a stochastic system is to use the second moment. Non-zero growth rates in the second moment imply instability in a typical member of the ensemble. One might expect that any even moment should be equivalent; however, Zhang *et al*. (1998) find that higher order moments may diverge even for an energetically stable state. Also, Lindenberg *et al*. (1981) demonstrate that the higher order even moments are less stable than the lower order even moments. Since our stochastic variables are continuous, we can use the classical dynamical systems theory for deterministic differential equations to quantify the growth rates. In particular, we use the Lyapunov exponent (Coddington & Levinson 1955) which Farrell & Ioannou (1996*a*,*b*) argue is a better measure for instability than any moment.

The Lyapunov exponent is defined as(3.1)where *Φ*(*t*) is the propagator of the solution. For an *N*-dimensional system, we pick the *N* columns of the identity matrix and integrate each of those forward in time. At time *t* these *N* solutions form the column vectors of the propagator, the eigenvalues of which determine the net growth of the system.

The governing equation (1.2) can be rewritten as a linear system(3.2)For the case where *ξ*(*t*) is periodic, as in Mathieu's equation, Floquet theory (Coddington & Levinson 1955) states that the growth rates can be determined by finding the eigenvalues of the propagator evaluated at the period of the internal forcing. This is much less expensive computationally and thus can be computed much more accurately. Using the numerical method discussed in §4, we can integrate equation (3.2) in time and thereby compute the propagator at any finite time. Figure 2 compares the growth rates computed by the Floquet method against those obtained by approximating the Lyapunov exponent as well as the error. The Lyapunov exponent was computed using 10^{3} periods, and the time step used in this and all other calculations is 10^{−2}. We see there is a very good agreement and the error is smaller by approximately two orders of magnitude near the maximum.

## 4. Numerical method

Since we cannot find analytical solutions for the classical Mathieu's equation in terms of the transcendental functions, there is even less likelihood that we can find such solutions for the case when the internal forcing is a stochastic process. Consequently, we use numerical methods to integrate the linear system in equation (3.2). Even though this is a very powerful approach, it does not allow us to find a solution as *t*→∞ as is required in the definition of the Lyapunov exponent (equation (3.1)). However, we surmise that if time is sufficiently large we should find a reasonable approximation to the growth rates. The regime of very large correlation times will necessarily have greater errors than in the case of smaller correlation times. We pick d*t*=*π*/100 and the final time of 10^{4} periods, where one period is 2*π*. Without loss of generality we pick *Ω*=1 so that the mean period is 2*π*.

Numerical solutions to equation (3.2) with particular initial conditions can be obtained from a variety of different methods; a classical example is the fourth-order Runge–Kutta method (Press *et al*. 1992). We must take care to pick a method that is close to the exact solution for long time. Numerical instability or drift could yield growth rates that are unphysical.

If we define we can restate equation (1.2) as(4.1)which describes the dynamics of a linear pendulum where the natural frequency , in which *g* and *L* are gravity and the length of the rod, is a time-varying function. It is well known from classical mechanics (Goldstein *et al*. 2001) that this equation has a Hamiltonian structure. If is time varying the Hamiltonian is non-autonomous and consequently the energy (Hamiltonian) is generally not conserved in time. However, we can define a related Hamiltonian in a larger phase space that is in fact autonomous. This is done through the following definition:(4.2)(4.3)(4.4)(4.5)which has a Hamiltonian of(4.6)By exploiting the symplectic structure we find that the governing equations are (Leimkuhler & Reich 2004)(4.7)(4.8)(4.9)(4.10)The first two equations determine the evolution of the coordinate *q* and the momentum *p*. The third equation states that *Q* is our time variable and the fourth equation is related to how the Hamiltonian of the non-autonomous system changes in time. To determine the solution we can neglect the latter two equations and simply solve the first two. For equation (4.10) to be well defined we need to be a continuously differentiable function, but for the coloured noises we are considering it as only necessarily continuous. For relatively large correlation times the noise will be smooth, but we recognize that for short correlation times this theory may not be strictly applicable.

There are a variety of numerical methods (called symplectic methods) that preserve the phase space (Hairer *et al*. 2004; Leimkuhler & Reich 2004), which ensure that the error in the Hamiltonian is bounded within some function of the time step for all time. When comparing the first-order symplectic Euler method (Hairer *et al*. 2004; Leimkuhler & Reich 2004) to the fourth-order Runge–Kutta method, we found that the solution in the latter had a Hamiltonian that slowly drifted away from the correct value whereas the former did not.

In order to get higher accuracy, we used the implicit trapezoidal method which is a second-order symplectic method in all of the numerical calculations. For the growth rate calculations, we renormalized the solution to unity after a fixed interval of *Ω*/(2*π*) and used this growth to approximate the Lyapunov exponent for that duration. We average over all these values to find an average Lyapunov exponent of the solution.

## 5. Modes of instability

We compute the Lyapunov exponents of the solutions to equation (3.2) with narrow-band Gaussian noise and forced dissipative random oscillator. The two contour plots in figure 3 depict the growth rates for the range of parameters where *ϵ*=0.1, *δ*∈[0,0.5] and *τ*_{c}∈[10^{−1},10^{4}]. We obtain that for correlation times of 10^{4}, the growth rates are very similar to those predicted by the classical Mathieu's equation. The growth rates are the largest for the narrow-band noise, but both cases yield qualitatively similar results. As the correlation time decreases below 10^{4}, the maximum growth decreases and there is also a widening of the unstable region. This denotes that the effect of the stochasticity is to stabilize the parametrically unstable mode. However, stochasticity also destabilizes particular modes, since the region of instability widens with the decreasing correlation time. The parametric instability continues to decrease and eventually it is completely stabilized.

Near a correlation time of 100 we see that for small values of *δ* there is a new unstable region that is introduced. It achieves larger maximum growth rates than the parametric region and we refer to this as the stochastic mode. It is important to differentiate between these two modes since they arise in very different scenarios and have different spectra, as we shall see in the next section. We are unaware of any other studies that have demonstrated this transition from a parametric mode at large correlation times to a stochastic mode for small correlation times.

To understand the origin of the stochastic mode, consider equation (1.2). If the squared frequency is positive we have oscillatory behaviour, whereas for negative values we have exponential behaviour. When *δ* is substantially large on average, it will dominate over the stochasticity. Conversely, when *δ* is very small the squared frequency is essentially the stochastic component. Since the random variable has mean zero, it will have some phases when the solution is oscillatory and others for which it will be evanescent. In the latter case we can expect the solution to grow exponentially, which is a mechanism by which instability can be generated. When the correlation time gets very small, this effect weakens due to the fact that the noise is even more random and it is harder to maintain the coherency necessary to give rise to instability.

The stochastic mode is similar to the one which occurs in Van Kampen's (2001) simple example of a stochastic linear harmonic oscillator. One important difference is that Van Kampen took the internal noise to be Gaussian whereas we are focused on non-Gaussian noise. Our results seem to indicate that the stochastic mode vanishes for small correlation time, but we do not approach the white noise limit so we are outside the range of applicability of Van Kampen's theory.

The two plots in figure 4 are contour plots of the Lyapunov exponents of the stochastic Mathieu's equation in the case where the noise is the Kubo oscillator for *ϵ*=0.1. In the figure 4*a* the variance is fixed and the correlation time is varied over the same range as in figure 3, whereas in figure 4*b*, we have the Kubo oscillator generated by a Wiener process and it is the variance that is varied over a wide range, *σ*∈[10^{−2},10^{1}]. Decreasing the correlation time in the first three noises and increasing the variance in the fourth noise, all have the effect of increasing the stochasticity in the internal noise. This is clear from the definition of the OU process, equation (2.2), since both act to increase the coefficient of the additive white noise, . An important difference is that the coefficient of the linear term in equation (2.2) is dependent on the correlation time but not on the variance, and so the effects of changing *σ*^{2} and are not exact. Note that the correlation times and the variance are of the noise that generates the internal noise and not of the internal noise itself.

Figure 4*a* illustrates that the stochastic mode is also introduced at correlation times of order 1. The decrease in correlation time decreases the magnitude of the parametric mode; however, there is no widening of the parametric mode and it is not entirely stabilized for the same range of *τ*_{c}. Figure 4*b* indicates that an increase in the variance gives more similar plots to those in figure 3, but inverted in the *y*-direction. The advantage of the Kubo oscillator is that in the limit of small variance we recover exactly the same growth rates as in the classic Mathieu's equation. This is unlike both the narrow-band noise and forced dissipative random oscillator, which yield approximations to the deterministic Mathieu's equation.

Figure 5 is the analogue of figure 4*b* except that *ϵ*=0.5. Figure 5*a*,*b* illustrates the effect of the increase in variance on both the first subharmonic, centred at *δ*=0.25, and first harmonic, centred at *δ*=1, respectively. The increase in *ϵ* by a factor of 5 increases the maximum growth rate by a factor slightly larger than 4. The first subharmonic is much larger and there is still a widening of the mode as *σ* increases. There is still evidence of the stochastic mode but it is very faint since it is dwarfed by the parametric instability. The maximum growth rate of the first harmonic is smaller than the first subharmonic by a factor of 10. It is clear that the stochasticity also stabilizes the first harmonic, and from this we extrapolate that the stochasticity should stabilize all the modes. As well, there appears to be an unstable path connecting the first two modes, which is caused by the spreading of the Arnold instability tongue with increasing stochasticity.

## 6. Power spectra of the modes

In the classical Mathieu's equation, the internal forcing occurs at a single frequency, say *Ω*, whose spectrum is a delta function. When we generalize the problem to non-periodic internal noise *ξ*(*t*), such as those defined in §2, the internal forcing is over a spectrum instead of a single point. In particular, figure 6 illustrates the power spectrum of the internal noise for various variances, *σ*=10^{−2}, 10^{−1}, 1, 10. In figure 6*a–d* we computed the fast Fourier transform of the norm of the solution and then normalized the amplitude and frequency such that the maximum is 1 and the harmonic frequency is 1 (Press *et al*. 1992). Figure 6*a*, which is nearly deterministic, *σ*=10^{−2}, has a spike at the harmonic frequency and a rapid decay away from this peak. In figure 6*b* the internal noise is only slightly non-periodic with *σ*=10^{−1}. The effect of the stochasticity is to broaden the peak and to shallow the slope of the decaying spectrum. The strong stochasticity in figure 6*c* has eliminated everything but a trace of a peak at the harmonic frequency. Finally, figure 6*d* has no evidence of a peak and the uniform spectrum is typical of a random variable.

The result of the harmonic forcing in the traditional Mathieu's equation is that there are a finite number of modes that are excited, and they are all of the form (*n*/2)*Ω*, where *n* is a positive integer and *Ω* is the harmonic frequency. When the integer is odd and even, we refer to the modes as the subharmonics and harmonics, respectively. The most unstable mode is the first subharmonic, *n*=1, and for small values of *ϵ*, the growth rates decay exponentially with *δ*. For each mode, the power spectrum indicates the frequency of the solution in relation to the internal forcing frequency.

Figure 7*a–d* shows the power spectra of the solution to the stochastic Mathieu's equation where the internal noise is the Kubo oscillator generated by a Wiener process, the power spectra of which is illustrated in the previous figure. The two parameters in equation (1.2) are *ϵ*=0.1 and *δ*=0.25, which focus on the first subharmonic Arnold instability tongue in the regime where the instabilities are relatively weak. The plots are normalized in the same fashion as for the internal noise and were also integrated for a thousand periods. Figure 7*a*,*b* with *σ*=10^{−2} and 10^{−1} is almost identical and very similar to the power spectra of the solution to the classical Mathieu's equation for the same choice of parameters. There is a broad peak at the first subharmonic frequency, as should be expected, but there is also evidence of a much smaller peak at the second subharmonic. This signifies that most of the energy in the instability is located at the first subharmonic but there is also a substantial amount of energy located near this critical frequency.

When the variance of the Wiener process is increased to *σ*=1, we observe in figure 7*c* that there is still a peak at the first subharmonic; however, it is not as smooth as the previous case. This irregularity is a classical signature of stochasticity and becomes more prominent with increasing *σ*. The second subharmonic is no longer apparent, the peak frequency is narrower than that in the deterministic limit and the spectrum decays more rapidly away from the peak. In figure 7*d*, where *σ*=10, the differences are even more magnified. This series of calculations indicates that, as stochasticity is introduced into the system, the range of excitable frequencies becomes very narrow but it does not vanish. We have demonstrated this for one particular type of noise but the three other types presented in §2 yield results that are qualitatively the same. The one major difference is that for the other noises the stochasticity is introduced by decreasing the correlation time unlike here where we increase the variance.

Figure 8 is the power spectrum computed for *ϵ*=0.1, *δ*=0.01 and *σ*=1. This choice of parameters focuses on an instability in the stochastically unstable region that is observed in all of our calculations of the Lyapunov exponents for small *δ* and either order 1 correlation time or order 1 variance, depending on the type of noise in question. The plot reveals the following: first, there is a peak much below the first subharmonic frequency; second, the spectrum is bumpy and irregular; and third, the spectrum decays away from the peak in a non-monotonic fashion. There is no evidence of a dramatic peak at the subharmonic or harmonic frequencies since the stochastic mode is not a parametric mode.

Finally, figure 9 gives the power spectrum for the same range of variances but now focusing on the first harmonic mode, *δ*=1. Figure 9*a* lies in the regime of the deterministic Mathieu's equation, the maximum peak is clearly indicated at the harmonic frequencies and the second and third harmonic frequencies are also apparent. As the variance increases we again see that the peak persists but narrows, the spectrum decays more quickly away from the peak and generally the spectrum is more uneven. The second harmonic persists for a wider range of variances than the second subharmonic in figure 7 but they are both altered in the same way by the stochasticity.

While computing the power spectrum for stable parameters, we found that there was a dominant peak with other smaller peaks in the spectrum and that the dominant peak moved to larger frequencies as *δ* increased. As the variance increased, the spectrum narrowed around the peak and decayed more rapidly away from the peak. This indicates that the dominant frequency of oscillation in the deterministic limit is still prominent in the stochastic limit. In the classical Mathieu's equation, there is one dominant peak and several smaller peaks while in the stochastic case there is only one peak for the regime that we explored.

## 7. Conclusions

The classical Mathieu's equation is the simplest equation that exhibits parametric instability, where there are infinite numbers of modes that are unstable and the most unstable is the first subharmonic frequency. There have been many studies of stochastic generalizations of Mathieu's equation. Nearly all of these have used Gaussian noise and focused on either white noise or coloured noise with small correlation times. In this manuscript we considered four different types of internal noise to force Mathieu's equation. The preferred type is the Kubo oscillator forced by a Wiener process since it tends to a sinusoidal function in the limit of zero variance and becomes less periodic as *σ* increases. It is perhaps a more natural choice to extend Mathieu's equation than a Gaussian variable since it cannot have a magnitude larger than 1.

We used a symplectic method to integrate the system of equations numerically to ensure we have a well-behaved system over long times. To quantify the growth rate, we computed the Lyapunov exponents of the system for a wide region of parameters. It is readily observed that for decreasing stochasticity, *σ*, we recover the same growth rates as the classical Mathieu's equation. As the variance increases the maximum growth decreased and the unstable region broadened, which signified that the stochasticity can both stabilize and destabilize a system in comparison with the periodic limit. We have found that when the variance is near unity, there is a stochastic mode that is introduced, which occurs where the mean frequency is zero or nearly zero.

The power spectrum of the solutions revealed the effect the stochasticity has on the stability of the system. The peak in power spectrum narrowed, the decay rate increased in moving away from the peak and the spectrum became irregular. The power spectrum of the stochastic mode appeared to be quite different, which is not surprising owing to the different nature by which this mode is excited. This is similar to the instability that appeared in Van Kampen's (2004) simple example of random harmonic oscillator forced by Gaussian white noise.

It has been stated by Farrell & Ioannou (1996*a*,*b*) that PR can occur in any system that has a non-normal structure and time variation. This manuscript has built a bridge between PR that arises from periodic forcing and PR that occurs in stochastic systems by studying the growth rates for a wide range of parameters. As the stochasticity increases, the unstable parametric mode weakens and broadens and is eventually dwarfed by the stochastic mode, which is only present for very small values of *δ*. It is of great interest to study the relative importance of these two modes in other complex systems as well as their nonlinear evolution. In future work we will extend this analysis to the stability of shear flows to determine whether the same qualitative features are observed and to understand the nature of the nonlinear equilibration.

## Acknowledgments

The authors would like to thank Matthew Scott, Zoran Miskovic and Marek Stastna for their useful discussions about stochastic processes, a reviewer for pointing out the relevance to Anderson localization and NSERC and NSF for financial support during the research of this manuscript.

## Footnotes

- Received January 7, 2008.
- Accepted March 5, 2008.

- © 2008 The Royal Society