An Implicit Function Method for Computing the Stability Boundaries of Hill's Equation Free

Second-order oscillators with time-periodically varying system parameters can experience the related phenomena of parametric resonance and vibrational stabilization. Parametric resonance involves exponential amplitude growth under certain combinations of the oscillator's ostensible “natural frequency” and parametric forcing amplitude. Parametric forcing refers to the effect of a time-periodic system parameter that (typically) multiplies the state, and differs from the more typical additive harmonic forcing in this way. Parametric resonance is observed in a variety of natural and engineered systems, such as in hydrodynamic instabilities [1,2], time-periodic axial loading of columns [3], and in MEMS devices subjected to alternating voltage [4].

A canonical model for linear second-order systems that can undergo parametric resonance is Hill's equation [5], a harmonic oscillator with a time-periodically varying spring constant. Certain combinations of the system's natural frequency and parametric forcing amplitude can cause instability, even at arbitrarily low forcing amplitudes under low or zero damping conditions. These unstable parameter ranges form the well-known Arnold tongue structures [6].

Vibrational stabilization, the flip side of parametric resonance, involves the stabilization via time-periodic parametric forcing of nominally unstable systems, typically a second-order oscillator with a negative spring constant. A common example of this is the Kapitza pendulum, an inverted pendulum stabilized via sinusoidal base motion that acts as sensorless feedback. Although previously considered a high-frequency phenomenon, recent work by Berg and Wickramasinghe [7] shows that vibrational stabilization can be achieved via parametric forcing at lower frequencies with precisely chosen amplitudes.

The stability boundaries in parameter space of both Hill's equation and the linearized Kapitza pendulum are determined by specific contours of an implicit function of the system parameters. This provides a two-dimensional stability map for both kinds of systems. By stability, we refer here to both neutral stability in the Lyapunov sense (boundedness for all time) and to asymptotic stability in the presence of damping. Such maps have been used in the design of quadrupole ion traps for mass spectrometry [8]. While canonical models of both phenomena are typically of second order, higher-order systems can be modeled as collections of Hill oscillators [9] or their unstable counterparts [10] under certain conditions, and stability criteria can be found as the composition or overlap of the individual two-dimensional stability maps of the component systems.

It is possible to derive analytical criteria for parametric resonance in the limit of small parametric forcing [6] using the classical Floquet theory [7] and for vibrational stabilization using higher-order averaging methods [11]. However, these results are valid in limited regimes—respectively, in the limit of low parametric forcing amplitudes and frequencies for parametric resonance and vibrational stabilization. More general estimates require numerical computation. Typically, these boundaries are computed and visualized via contour plots of the system, as computed using a marching squares algorithm [12]. In this work, we propose a more efficient method involving direct numerical calculation and integration of this implicit function. This method leans on the idea that we are only interested in specific contours and not the overall picture, reducing the numerical complexity of the problem.

To illustrate the method, we cover the general idea of integrating an implicit function of two variables in Sec. 2. For Hill's equation, this function does not have a closed form, so we cover some relevant properties in Sec. 3.1, additional solution steps required in Sec. 3.2, and provide some examples of the results of this method. Section 3.3 extends the numerical method to the case of Hill's equation with damping. Finally, we cover some considerations for numerical integration in Sec. 4.

The general idea of the method is as follows: given a smooth function $f(a,ϵ)$ of variables a, $ϵ$⁠, we wish to plot a specific contour $f(a,ϵ)=c$ in $(a,ϵ)$ space passing through $(a0,ϵ0)$⁠, where $f(a0,ϵ0)=c$⁠.

where at least one of the two equations is well-defined in a neighborhood of $(a0,ϵ0)$⁠. We note that neither the ordinary differential equations (ODEs) nor the initial conditions in (2.1) involve c—these ODEs are satisfied by the contour of f that passes through $(a0,ϵ0)$ irrespective of its value.

The latter can be avoided by rescaling a and $ϵ$ appropriately locally around $(a,ϵ)$⁠. In both cases, stitching together the piecewise numerical solutions for the contour $f(a,ϵ)=c$ can lead to a loss of accuracy, which we describe in further detail in Sec. 4.

which represents a harmonic oscillator with a periodically varying spring constant. $a>0$ is the mean value of the spring constant, and $ϵ$ is the parametric forcing amplitude. The parametric forcing $p(t)$ has zero mean. Such systems can experience parametric resonance for specific choices of a and $ϵ$⁠. Allowing a to be negative gives us the linearized model of the Kapitza pendulum, an inverted pendulum with sinusoidal vertical base motion. This nominally unstable system is vibrationally stabilized for specific choices of a and $ϵ$⁠. Our treatment of stability boundaries covers both regimes of behavior.

We review a few salient properties of Hill's equation for applying the implicit function method to compute its parametric stability.

Since $detΘ(0,0)=detI=1$⁠, it is 1 for all time.

Since $Θ(2π,0)$ is real, we have $Tr{Θ(2π,0)}∈ℝ$⁠. Since $Θ(2π,0)$ is $2×2$⁠, we have $−2≤Tr{Θ(2π,0)}≤2$ iff $σ¯(Θ(2π,0))≤1$⁠. Thus, system (3.1) is stable iff $|Tr{Θ(2π,0)}|≤2$⁠, with equality occurring when both eigenvalues of $Θ(2π,0)$ are at $±1$⁠.

The condition $TrΘ(2π,0;a,ϵ)=±2$ is a function of two variables whose contours in $(a,ϵ)$ space with values $±2$ form the boundary between regions of stability and regions of parametric resonance. These boundaries in parameter space form the well-known Arnold tongue structures associated with parametric resonance (Fig. 1).

Further, it can be shown [6] that there are no other values of $a>0$ (with $ϵ→0+$⁠) for which the system is unstable.

The stability boundaries are typically computed in the following way:

1. Create a grid in $(a,ϵ)$ space that covers the parameter ranges of interest.

2. Find the monodromy map (and thus its trace) numerically at all points on this grid.

3. Plot the contours of the map that have values $±2$⁠.

This method has three shortcomings. The first is the computational burden of calculating the value of $Tr{Θ(2π,0)}$ at all grid points. Essentially, this involves scanning a two-dimensional space in search of a one-dimensional curve. The second is that the slope of the curve is not known a priori, so the spatial grid resolution cannot be adjusted to take advantage of this. As a result, several runs with progressively finer grids are required to obtain a suitably detailed picture. Finally, there are regions in the parameter space where $Tr{Θ(2π,0)}$ changes rapidly with a, making it difficult to resolve fine details in these areas in a contour map created with a typical marching squares algorithm.

Applying the implicit function method outlined in Sec. 2 avoids all of these shortcomings, but it requires knowledge of at least one point on the contour of interest. Specifically, for each contour such that $Tr{Θ(2π,0;a0,ϵ0)}=±2$⁠, we need to know the corresponding $(a0,ϵ0)$⁠. From Eq. (3.4), this point is $(n2/4,0)$⁠.

where we have suppressed the explicit dependence of $Θ$ and A on a and $ϵ$ in the notation. The linear ODEs (3.7) and (3.8) can be integrated from $t=0$ to $2π$ to give $(∂Θ/∂a)(2π,0;a,ϵ)$ and $(∂Θ/∂ϵ)(2π,0;a,ϵ)$⁠, respectively. Since $Θ(0,0;a,ϵ)=I$ uniformly for any choice of a and $ϵ$⁠, the initial conditions are uniformly zero in the above equations, and they are driven by the state-transition matrix $Θ(t,0;a,ϵ)$ at $(a,ϵ)$⁠, which appears as an additive input. This recipe for finding stability boundaries is illustrated in Figs. 2 and 5.

Figure 3 compares the stability maps for Hill's equation as computed using the contour method and the implicit function method introduced in this paper. The implicit function method provides better resolution with significantly less computational effort.

Figure 4 shows the stability boundaries computed using this method for Hill's equation with a few different parametric forcing functions.

The variable substitution $z(t):=eκtθ(t)$ may be used to represent Eq. (3.9) as an equivalent undamped Hill's equation [7]. Using this transformation, it is possible to derive a stability criterion for Hill's equation with viscous damping.

Proof. See the  Appendix.▪

Since $2cosh(2πκ)≥2$⁠, the Arnold tongues for Hill's equation with damping do not meet on the a axis, and there is no equivalent to criterion (3.4). Integrating Eq. (3.5) along the contour $|Tr[Φ(2π,0;a,ϵ)]|=2cosh(2πκ)$ thus requires finding one point along each “branch” of the stability boundary.

We find these initial conditions $(ϵ0,a0)$ numerically, using the criteria $|Tr[Φ(2π,0;a0,ϵ0)]|=2cosh(2πκ)$ and $(∂/∂a)[Φ(2π,0;a,ϵ)]|(ϵ0,a0)=0$⁠. The rest of the treatment is identical to the undamped case covered in Sec. 3.3. Figure 5 illustrates the method.

In this section, we briefly mention some details of and considerations for the numerical integration of Eq. (3.5).

• Since Hill's equation is a Hamiltonian dynamical system, Eq. (3.6) needs to be solved using a symplectic integrator. Depending on the required accuracy, a symplectic Euler method or an “almost symplectic” trapezoidal method can suffice. This work uses the latter, which is symmetric and second order, with adaptive time-stepping and no numerical damping [15].

• The dynamics of $∂Θ/∂a$⁠, $∂Θ/∂ϵ$⁠, and $da/dϵ$ (Eqs. (3.5), (3.7), and (3.8), respectively) have no special structure and are solved for using an explicit Runge–Kutta (Owren–Zennaro optimized interpolation $5/4$ [16]) method.

• The solution method fails at Arnold tongue “crossover points,” where tongue boundaries intersect at nonzero values of $ϵ$⁠, creating stability “pockets.” This is the case, for instance, for square wave forcing in Fig. 4. At these intersection points, $Tr[(∂Θ/∂a)(2π,0;a,ϵ)]→0$ and the integration for Eq. (3.5) can stall. When using an adaptive time-step method as we do here, this stall can be diagnosed by a progressive decrease in the time-step down to machine precision, before which the integration needs to be suspended. This can be addressed in a few different ways:

• Switch conditionally to the solving the reciprocal problem:

  $dϵda=−Tr[∂Θ∂a(2π,0;a,ϵ)]Tr[∂Θ∂ϵ(2π,0;a,ϵ)]$

• Find new initial conditions for integrating Eq. (3.5) by extrapolating the two tongue boundaries past the pocket, while ensuring that $Tr[(∂Θ/∂a)(2π,0;a,ϵ)]=±2$⁠.

  Both methods cause a loss of accuracy in the solutions past the “crossover” point, which is a limitation of this method.

• Finally, computing the stability map for the damped Hill's equation requires finding points $(ϵ0,a0)$ where $Tr[Φ(2π,0;a0,ϵ0)]=2cosh(2πκ)$ and $(∂/∂a)Tr[Φ(2π,0;a,ϵ)]|(ϵ0,a0)=0$⁠. These may be found by integrating Eqs. (3.6) and ((3.7) using an explicit Runge–Kutta method (Owren–Zennaro optimized interpolation 5/4) followed by a binary search above the locations of the tongues of the undamped system, i.e., in the vicinity of $(ϵ0,a0)=(0,n2/4)$ for $n∈ℕ$⁠.

The implicit function method detailed in Sec. 2 can be applied to compute specific contours of any smooth scalar-valued function more efficiently than evaluating the function on a grid.

While the method provides better resolution of the tongues compared to plotting contours of the trace of the monodromy map, it also rests on certain assumptions about the geometry of the Arnold tongues, for instance, that the stability boundaries are smooth and can be represented in the form $a(ϵ)$⁠. It is thus not a true parametric description of the stability boundaries of the form $(a(η),ϵ(η))$ for some parameter $η$ and requires careful handling of crossover points of the Arnold tongues, where two tongues meet. In particular, it is not yet clear how to retain the accuracy of the numerical solution past the crossover points.

It is also reliant on the existence of a scalar criterion for stability, of the form $f(a,ϵ)=constant$⁠. The stability criterion for parametric resonance in higher-order linear systems involves the eigenvalues of the monodromy map that is not reducible to a criterion involving its trace. When the system cannot be diagonalized into a collection of decoupled second-order time-periodic oscillators, neither the contour method nor the implicit function method discussed in this paper can be directly applied to obtain stability maps. Extending the method to generate the stability maps of higher-order systems, possibly using the symplectic properties of the monodromy map, is a promising avenue of future work.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

In this section, we review properties of the monodromy map of the damped Hill's equation that we use in Sec. 3.3.

The nominal spring constant in the transformed coordinates is lower than the original by $κ2$⁠. We denote the state-transition matrix for system (A3) as $Φ(t,0)$⁠, suppressing the dependence on the parameters a, $ϵ$⁠, and $κ$⁠.

System (3.9) is unstable when $Θ(2π,0)$ has an eigenvalue outside the unit circle. Our goal is to obtain an equivalent stability criterion for the damped system (3.9) in terms of the monodromy map $Φ(2π,0)$ of the transformed system. The implicit function method of Sec. 3 can then be applied to the transformed, undamped system to plot the stability map of the original damped system. To do this, we use the following property of the solution of the damped Hill's equation:

Lemma A.1. The monodromy map$Θ(2π,0)$of the damped Hill's equation (3.9)is a scalar multiple of a symplectic matrix. Specifically,$e2πκΘ(2π,0)$is symplectic.

is a product of three symplectic matrices and is also symplectic.

This result is summarized in Table 1. We note that $e2πκ+e−2πκ≥2$⁠, with the lower bound achieved in the undamped case of $κ=0$⁠.