Integrable nonlinear oscillators with polynomial invariants: construction, Poincare geometry, and an analytic stability boundary

Starting from the nonlinear ODE $z^{\mathrm{\prime }\mathrm{\prime }}+f(t)z+g(t)z^m=0$ with $m>1$, we show that after a suitable normal-form reduction of any Hill equation one may, without loss of generality, fix the linear part as $f(t)\equiv {\omega }^2$ (with $\omega >0$ constant). For the class $z^{\mathrm{\prime }\mathrm{\prime }}+{\omega }^{2}z+g(t)z^m=0$ with $m>1$, our goal is to compile a catalogue of all possible integrable cases. We restrict attention to integrals that are polynomial in the variables $z$ and $p=z^{\mathrm{\prime }}$. The Hamiltonian does not provide such an integral because it is explicitly time dependent. Instead, we search for invariants that are quadratic in $p=z^{\mathrm{\prime }}$. We show that such invariants exist precisely when ${\alpha }_2(t):=g(t{)}^{-2\mathrm{/}(m+3)}$ satisfies the linear third-order ODE ${\alpha }_2^{\mathrm{\prime }\mathrm{\prime }\mathrm{\prime }}+4{\omega }^2{\alpha }_2^{\mathrm{\prime }}=0$. This yields the three-parameter solution $g(t)=[a_0+a_1\cos (2\omega t)+a_2\sin (2\omega t){]}^{-(m+3)\mathrm{/}2}$. For $m=2$ this reproduces the trigonometric structure with exponent $-5\mathrm{/}2$ found in Hagel–Bouquet (1992). In addition we present a detailed stability analysis based on the invariant using Poincaré sections and find full agreement with numerical simulations.