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

TL;DR

We address computing stability boundaries for Hill's equation, a time-periodic oscillator that exhibits parametric resonance, by deriving an implicit-function method to trace targeted stability contours in parameter space. The core idea is to treat the boundary as an implicit contour f(a,epsilon)=c and integrate an ODE for a(epsilon) given by da/depsilon, leveraging sensitivity equations to compute derivatives of the monodromy map. This approach yields more accurate and efficient boundary traces than traditional grid-based contouring of Tr Theta(2 pi,0;a,epsilon)=±2, and extends to damped Hill's equation via a damping transformation and the criterion |Tr Phi(2 pi,0)| ≤ 2 cosh(2 pi kappa). The work also discusses numerical integration considerations, extension potential to higher-order systems, and limitations when the contour geometry is non-smooth or lacks a scalar stability criterion. The resulting methods provide tight stability maps useful for design problems involving parametric resonance and vibrational stabilization, including Kapitza-like scenarios.

Abstract

Hill's equation is a common model of a time-periodic system that can undergo parametric resonance for certain choices of system parameters. For most kinds of parametric forcing, stable regions in its two-dimensional parameter space need to be identified numerically, typically by applying a matrix trace criterion. By integrating ODEs derived from the stability criterion, we present an alternative, more accurate and computationally efficient numerical method for determining the stability boundaries of Hill's equation in parameter space. This method works similarly to determine stability boundaries for the closely related problem of vibrational stabilization of the linearized Katpiza pendulum. Additionally, we derive a stability criterion for the damped Hill's equation in terms of a matrix trace criterion on an equivalent undamped system. In doing so we generalize the method of this paper to compute stability boundaries for parametric resonance in the presence of damping.

Figures (5)

• Figure 1: \newlabelfig:stability-boundaries-explanation0 A representative stability diagram for Hill's equation illustrating the implicit function method \ref{['eq:ode-d-a-d-eps']} of finding stability boundaries. The time-periodic forcing in this example is sinusoidal with period $2\pi$. The system experiences parametric resonance with paramters in the shaded regions. All stability boundaries are the curves $(\epsilon, a(\epsilon))$ given by the implicit equation $\left| \mathop{\mathrm{Tr}}\nolimits \Theta(2\pi, 0; a(\epsilon), \epsilon) \right| = 2$, starting at $a = n^2/4,\ n \in \mathbb{N}$ and $\epsilon = 0$. These curves are found by integrating \ref{['eq:ode-d-a-d-eps']} starting at these initial conditions.
• Figure 2: \newlabelfig:solution-algorithm-step0 Illustration of the solution method for stability boundary calculation for Hill's equation. The goal is to find the set of points $\left\{ (a, \epsilon) \right\}$ that satisfies the topmost block, \ref{['eq:ode-d-a-d-eps']}, which is a curve passing through known $(a_0, \epsilon_0)$. For Hill's equation this stability boundary is representable as $a(\epsilon)$. At each time step, the current curve coordinates $(a, \epsilon)$ determine the state transition matrix $\Theta(;a, \epsilon)$ that appears as an input elsewhere. We compute the state transition matrix over one forcing period and use it to integrate systems \ref{['eq:ode-d-phi-d-a']} and \ref{['eq:ode-d-phi-d-eps']}. The final value of $\frac{\partial \Theta}{\partial a}$ and $\frac{\partial \Theta}{\partial \epsilon}$ gives us the value of $\frac{\mathrm{d} a}{\mathrm{d} \epsilon}$ for this time step. Advancing the main integrator in time in turn gives us the new coordinates $(a,\epsilon)$ required to continue the solution. (The dependence of $\Theta(t,0; a, \epsilon)$ and $A(t)$ on time has been suppressed in the notation.)
• Figure 3: \newlabelfig:boundaries-cos-contour-comparison0 Comparison of the stability map in $(\epsilon, a)$ parameter space of Hill's equation with sinusoidal parametric forcing, as computed by (left) the contour method and (right) the implicit function method detailed in \ref{['sec:hill-equation-stability-details']} and \ref{['fig:stability-boundaries-explanation']}. The contour method involves plotting the contour of the function $f(a, \epsilon) = \mathop{\mathrm{Tr}}\nolimits\ \Theta(2\pi, 0; a, \epsilon)$ with absolute value $2$. This function is computed at a grid resolution of $(\Delta \epsilon, \Delta a) = (0.02, 0.02)$, and the contour levels plotted are $1.99 < \left| \mathop{\mathrm{Tr}}\nolimits\ \Theta(2\pi, 0; a, \epsilon) \right| < 2.01$. The implicit function method involves integrating one implicit function for each Arnold tongue boundary \ref{['eq:ode-d-a-d-eps']} with a step of $\Delta \epsilon = 0.05$. The implicit function method is more efficient, and provides much better resolution of the Arnold tongues as well, especially for small $\epsilon$. In our tests, it is also faster by a factor of $10-50$, depending on the desired grid resolution.
• Figure 4: \newlabelfig:stability-diagram-sideways-all0 Stability diagrams in the $(\epsilon, a)$ parameter space for Hill's equation, for four different kinds of zero-mean parametric forcing (inset, in blue). Clockwise from the top, they are sinusoidal, square waves with even and uneven duty cycles, and a periodic ramp. The stability boundaries are computed using the implicit function method (\ref{['fig:solution-algorithm-step']}) for the different forcing functions. For $a > 0$, the plots show the first three Arnold tongues, whose locations as $\epsilon \to 0^+$ are independent of the type of forcing. For $a < 0$, the plots show the $(\epsilon, a)$ parameter "window" that causes vibrational stabilization of the nominally exponentially unstable system.
• Figure 5: \newlabelfig:plot-boundaries-explanation0 Illustration of the method for stability boundary calculations for Hill's equation with and without damping (\ref{['eq:hill-ode', 'eq:hill-ode-damped']} respectively). Left: Hill's equation without damping. Each Arnold tongue boundary corresponds to a contour $\left| \mathop{\mathrm{Tr}}\nolimits \left[ \Theta(2\pi, 0) \right] \right| = 2$ and meets the $a$ axis at $n^2/4$, $n \in \mathbb{N}$. Thus ODE \ref{['eq:ode-d-a-d-eps']} can be integrated starting at $(\epsilon_0, a_0) = \left( 0, n^2/4 \right)$. See \ref{['fig:solution-algorithm-step']} for details of the integration process. Right: Hill's equation with damping coefficient $\kappa$. The method is similar to the undamped case, but the initial conditions $(\epsilon_0, a_0)$ for each Arnold tongue have to be found numerically. The ODE \ref{['eq:ode-d-a-d-eps']} is integrated from the left-most point of each tongue, corresponding to $\left. \mathop{\mathrm{Tr}}\nolimits \left[ \Phi(2\pi, 0) \right] \right|_{(\epsilon_0, a_0)} = 2 \cosh(2\pi \kappa)$ and $\left. \frac{\partial }{\partial a} \mathop{\mathrm{Tr}}\nolimits \left[ \Phi(2\pi, 0) \right] \right|_{(\epsilon_0, a_0)} = 0$. See \ref{['sec:hill-equation-damped-stability-details']} for details. In both cases, finding $a$ as a function of $\epsilon$ requires picking one of the two "branches" of the stability boundary passing through $(\epsilon_0, a_0)$. The upper and lower branches of the stability boundary (colored blue and green respectively) are resolved by perturbing $(\epsilon_0, a_0)$ appropriately.

Theorems & Definitions (4)

• Theorem 1
• Proof 1
• Lemma 1
• Proof 2