🔬 mesoscale physics

HarmonicBalance.jl: A Julia suite for nonlinear dynamics using harmonic balance

HarmonicBalance.jl is a free Julia package that simplifies the analysis of nonlinear periodic systems by combining the harmonic balance method with symbolic algebra and homotopy continuation to numerically determine all steady-state solutions and their fluctuation dynamics, complete with a graphical user interface for exploring solution topologies.

Original authors: Jan Košata, Javier del Pino, Toni L. Heugel, Oded Zilberberg

Published 2026-09-17

📖 4 min read☕ Coffee break read

CC BY 4.0

Original authors: Jan Košata, Javier del Pino, Toni L. Heugel, Oded Zilberberg

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Many physical systems, from the swinging of a bridge in the wind to the flow of electricity through a circuit, are governed by rules that change over time. When these systems are simple, their behavior is predictable and steady. However, when different parts of a system interact in complex, nonlinear ways, they can produce surprising results. They might suddenly jump to a new state, oscillate in unexpected patterns, or settle into one of several possible stable conditions. Scientists have long known that these nonlinear systems are everywhere, appearing in everything from fluid dynamics to biological processes. Yet, predicting exactly how such a system will behave is notoriously difficult. Traditional methods often rely on guessing a starting point and watching how the system evolves, but this approach has a major flaw: it usually finds only one outcome, missing the other possible states that might exist just as validly. If a system can settle into three different stable patterns, a standard simulation might only reveal one, leaving the full picture incomplete and potentially misleading for engineers and researchers trying to understand the system's true capabilities.

To solve this problem, a team of researchers has developed a new software tool called HarmonicBalance.jl, designed to map out every possible steady state of these complex, time-varying systems. Instead of watching a system evolve from a single starting point, the tool treats the problem as a puzzle of finding all possible solutions at once. It works by assuming that the system's behavior can be described as a combination of repeating waves, or harmonics, that change very slowly over time. By translating the complex, moving equations of the system into a set of static algebraic equations, the software can then use a powerful mathematical technique called homotopy continuation. This method acts like a continuous bridge, starting from a simple, solvable version of the problem and slowly morphing it into the complex, real-world version. As it transforms, it tracks every single path that leads to a solution, ensuring that no possible steady state is left behind. This allows the researchers to see the entire landscape of possibilities, including states that might be unstable or difficult to reach in a real experiment.

The team tested this approach on several classic examples of nonlinear physics, such as a driven oscillator that responds to external forces with a twist of nonlinearity. In these tests, the software successfully identified multiple stable states where traditional methods would have found only one. For instance, in a system of two coupled oscillators, the tool revealed a complex landscape of solutions, showing how the system could jump between different states depending on the direction in which a parameter was changed. This phenomenon, known as hysteresis, was mapped out in detail, showing exactly where the system would snap from one state to another. The researchers also demonstrated that the tool could handle increasingly complex systems, such as a chain of five coupled oscillators, finding hundreds of distinct solutions in a matter of minutes on a standard computer. While the number of theoretical solutions can grow exponentially with the size of the system, the software efficiently filters out the impossible ones, focusing on the real, physical states that actually exist.

Beyond simply finding these states, the tool also analyzes their stability, determining which solutions are robust and which would collapse under a small disturbance. It can even simulate how a system would respond to slow changes in its environment, mimicking the way an experimenter might gradually turn a dial to sweep through different conditions. This capability is crucial for understanding phenomena like spontaneous symmetry breaking, where a system chooses one state over another without any external push. The software is written in Julia, a language known for its speed and ease of use, and is freely available to the scientific community. By providing a way to see the full range of possible behaviors in nonlinear systems, this tool offers a new level of clarity for researchers in fields ranging from structural engineering to optics. It transforms the study of these complex systems from a game of chance, where one hopes to stumble upon the right state, into a comprehensive exploration where every possibility is accounted for.

Technical Summary of "HarmonicBalance.jl: A Julia suite for nonlinear dynamics using harmonic balance"

Problem Statement
Nonlinear ordinary differential equations (ODEs) govern the time evolution of diverse physical systems, from fluid dynamics and robotics to optomechanics and oscillator networks. While linear systems typically relax to a unique steady state, nonlinear systems often exhibit multiple steady-state solutions, leading to phenomena such as hysteresis, spontaneous symmetry breaking, and noise-induced switching. A significant challenge in analyzing these systems is that standard numerical ODE solvers (initial value problems) typically find only a single steady state per run, dependent on initial conditions. Consequently, a complete exploration of the solution landscape requires infinite sampling of the initial condition space, which is computationally infeasible. Furthermore, while transforming driven systems into a rotating frame can render drives stationary, finding all roots of the resulting coupled nonlinear algebraic equations remains difficult for multivariate systems, as standard root-finding algorithms (e.g., Newton-Raphson) only locate solutions near an initial guess.

Methodology
The authors introduce HarmonicBalance.jl, a Julia package designed to find all steady-state solutions for harmonically-driven nonlinear systems. The methodology integrates three core components:

Harmonic Balance Ansatz: The package assumes that for sufficiently long times, the system's response can be approximated by a finite sum of harmonic terms with slowly varying amplitudes. For a system of $N$ second-order ODEs with harmonic time dependence, the solution $x_i(t)$ is expanded as:
$x_i(t) = \sum_{j=1}^{M} [u_{i,j}(T)\cos(\omega_{i,j}t) + v_{i,j}(T)\sin(\omega_{i,j}t)]$
where $T$ is a coarse-grained timescale. By substituting this ansatz into the original ODEs and filtering out fast oscillations (dropping second-order time derivatives of the amplitudes), the problem is reduced to a set of autonomous first-order ODEs for the harmonic variables $u$ and $v$ . Steady states are found by setting the time derivatives of these variables to zero, resulting in a system of coupled polynomial algebraic equations.
Homotopy Continuation: To solve the resulting polynomial systems, the package employs the homotopy continuation method via the HomotopyContinuation.jl library. Unlike iterative root-finding methods, homotopy continuation deforms an analytically solvable system (with a known number of roots) into the target system. This approach guarantees the discovery of all complex roots. The package utilizes a two-step algorithm to manage computational cost: a "warm-up" phase tracks paths from a generic system to filter out singular paths, followed by a parameter homotopy phase to track the remaining non-singular paths to physical solutions.
Symbolic and Numerical Integration: The package uses Symbolics.jl for symbolic manipulation to derive the harmonic equations and Jacobian matrices automatically. It interfaces with DifferentialEquations.jl for time-dependent simulations, allowing users to verify steady-state results and study adiabatic parameter sweeps (hysteresis).

Key Contributions

Comprehensive Solution Finding: The primary contribution is the ability to numerically determine all steady-state solutions (real and complex) for harmonically-driven nonlinear systems, rather than a single solution dependent on initial conditions.
Unified Framework: The package unifies symbolic algebra, harmonic balance, and homotopy continuation into a single, open-source Julia environment, offering high performance and flexibility.
Stability and Fluctuation Analysis: It provides tools to automatically classify solutions as real or complex and stable or unstable based on the eigenvalues of the Jacobian matrix. It also supports the calculation of linear response spectra for stable fixed points.
Visualization: The suite includes a graphical user interface and plotting routines to visualize 1D and 2D phase diagrams, solution branches, and bifurcation structures, facilitating the exploration of complex solution topologies.

Results and Examples
The paper demonstrates the package's capabilities through several examples:

Driven Duffing Oscillator: The package successfully identifies the multiple steady-state solutions of a single Duffing oscillator, including the superharmonic response ( $3\omega_d$ ) generated by nonlinearity. It reproduces standard resonance curves and phase diagrams showing regions with one, two, or three solutions.
Coupled Duffing Resonators: For a system of two coupled oscillators, the package maps the solution landscape across parameter space, revealing intricate topologies. It contrasts the "complete" solution set found via homotopy with the limited subset accessible via adiabatic time-sweeps, demonstrating hysteresis and jumps between branches.
Parametrically Driven Oscillator: The package analyzes a parametrically driven Duffing oscillator (Mathieu equation), identifying regions of single and double stable solutions and characterizing the pitchfork bifurcation where the number of stable states changes.
Performance Scaling: Benchmarks on a chain of $N$ coupled Duffing oscillators show that while the Bézout bound (theoretical maximum number of solutions) scales exponentially ( $3^{2N}$ ), the number of physical real solutions is significantly smaller. The package can solve systems with $N=5$ (involving hundreds of paths) on a single CPU core within minutes, highlighting the efficiency of the path-tracking strategy.

Significance
The authors position HarmonicBalance.jl as a tool to bridge the gap between theoretical analysis and experimental observation in nonlinear dynamics. By enabling the exploration of the complete solution landscape, the package allows researchers to identify experimentally observable phenomena such as hysteresis and symmetry breaking that are often missed by standard time-evolution simulations. The software is designed to be accessible to the academic community, offering an "ab-initio" approach where arbitrary equations of motion can be defined without specialized GUI constraints. The authors note that the tool is applicable to a wide range of fields, including structural dynamics, nonlinear optics, optomechanics, and Ising machines. While the current version focuses on deterministic dynamics and steady states, the modular design is intended to support future extensions, such as limit cycle detection and higher-order averaging methods.

More like this