Harmonic balance

Harmonic balance is a method originally introduced by Michel Nakhla. It is used to calculate the steady-state response of nonlinear differential equations, and is mostly applied to nonlinear electrical circuits
It is a "frequency domain" method for calculating the steady state, as opposed to the various "time-domain" steady state methods. The name "harmonic balance" is descriptive of the method, which starts with Kirchhoff's Current Law written in the frequency domain and a chosen number of harmonics. A sinusoidal signal applied to a nonlinear component in a system will generate harmonics of the fundamental frequency. Effectively the method assumes the solution can be represented by a linear combination of sinusoids, then balances current and voltage sinusoids to satisfy Kirchhoff's law. The method is commonly used to simulate circuits which include nonlinear elements, and is most applicable to systems with feedback in which limit cycles occur.
Microwave circuits were the original application for harmonic balance methods in electrical engineering. Microwave circuits were well-suited because, historically, microwave circuits consist of many linear components which can be directly represented in the frequency domain, plus a few nonlinear components. System sizes were typically small. For more general circuits, the method was considered impractical for all but these very small circuits until the mid-1990s, when Krylov subspace methods were applied to the problem. 
The application of preconditioned Krylov subspace methods allowed much larger systems to be solved, both in size of circuit and in numbers of harmonics. This made practical the present-day use of harmonic balance methods to analyze radio-frequency integrated circuits (RFICs).
Algorithm.
The harmonic balance algorithm is a special version of Galerkin's method. It is used for the calculation of periodic solutions of autonomous and non-autonomous differential-algebraic systems of equations. The treatment of non-autonomous systems is slightly simpler than the treatment of autonomous ones. A non-autonomous DAE system has the representation
with a sufficiently smooth function formula_2
where formula_3 is the number of equations and formula_4 are placeholders for time, the vector of unknowns and the vector of time-derivatives.
The system is non-autonomous if the function formula_5 is not constant for (some) fixed formula_6 and formula_7. Nevertheless, we require that there is a known "excitation period" formula_8 such that formula_5 is formula_10-periodic.
A natural candidate set for the formula_10-periodic solutions of the system equations is the Sobolev space formula_12 of weakly differentiable functions on the interval formula_13 with periodic boundary conditions formula_14.
We assume that the smoothness and the structure of formula_15 ensures that formula_16 is square-integrable for all formula_17.
The system formula_18 of harmonic functions formula_19 is a Schauder basis of formula_12 and forms a :Hilbert basis of the Hilbert space formula_21 of square-integrable functions. Therefore, each solution candidate formula_17 can be represented by a Fourier-series formula_23 with Fourier-coefficients formula_24 and the system equation is satisfied in the weak sense if for every base function formula_25 the variational equation
is fulfilled. This variational equation represents an infinite sequence of scalar equations since it has to be tested for the infinite number of base functions formula_27 in formula_28.
The Galerkin approach to the harmonic balance is to project the candidate set as well as the test space for the variational equation to the finitely dimensional sub-space spanned by the finite base formula_29.
This gives the finite-dimensional solution formula_30 and the finite set of equations
which can be solved numerically.
In the special context of electronics the algorithm starts with Kirchhoff's current law written in the frequency-domain. To increase the efficiency of the procedure, the circuit may be partitioned into its linear and nonlinear parts, since the linear part is readily described and calculated using nodal analysis directly in the frequency domain.
First, an initial guess is made for the solution, then an iterative process continues:
Convergence is reached when formula_39 is acceptably small, at which point all voltages and currents of the steady-state solution are known, most often represented as Fourier coefficients.
Tools.
A harmonic-balance tool, named 
Agile,
for microwave circuits is available for download.
A parallelized version
was also developed, but this version is not available.
Sandia National Labs developed Xyce, a high performance parallel electronic simulator that can perform harmonic balance analysis.
The harmonic balance method is also natively supported for general nonlinear multiphysics finite element simulations in the open source c++ FEM library Sparselizard.