pyHB: an open-source automatic-differentiation-enhanced semi-analytical solver for nonlinear dynamics

The Harmonic Balance (HB) method is widely used to compute and analyze the periodic responses of nonlinear systems. However, its application to high-dimensional complex systems is limited by the burden of handling the partial derivatives of the nonlinearities. This work presents pyHB, an open-source, automatic-differentiation-enhanced semi-analytical framework that integrates the complete HB workflow for general user-defined nonlinear systems. The proposed formulation exploits localized nonlinearities and applies PyTorch-based automatic differentiation (AD) only to the reduced nonlinear force, thereby avoiding the need for user-supplied derivatives of the nonlinear force and maintaining controllable GPU memory usage. Weighted arc-length continuation, sparse matrix assembly, a blocked solution strategy for the augmented continuation equations, and Floquet-based stability analysis are incorporated within a modular architecture that separates model definition from reusable numerical procedures. Hence, pyHB can provide a complete landscape of the nonlinear system's periodic response based solely on the user-defined dynamical equations. Four examples, including a quasi-zero-stiffness isolator, a nonlinear piezoelectric energy harvester, a 284 degrees of freedom (DOFs) aeroengine model, and a 2000 DOFs Bernoulli beam, demonstrate the ability of pyHB to trace stable and unstable solution branches and capture subharmonic resonance, combination resonance, and mixed-order electromechanical responses. Notably, in the Bernoulli beam example with 202000 HB unknowns, the AD-enhanced solver requires approximately 0.44s per continuation point, achieving several-hundred-fold speedup compared to the Newmark-$β$ method and remaining 637.8MB of additional RAM and 243.5MB of GPU memory. The proposed pyHB provides a general, one-stop benchmark platform for HB-based nonlinear dynamics analysis.