^[8]

Consider the differential equation ${\displaystyle {\ddot {x}}+x^{3}=0}$ . We use the ansatz solution ${\displaystyle x=A\cos(\omega t)}$ , and plugging in, we obtain

Then by matching the ${\displaystyle \cos(\omega t)}$  terms, we have ${\displaystyle \omega ={\sqrt {\frac {3}{4}}}A}$ , which yields approximate period ${\displaystyle T={\frac {2\pi }{\omega }}\approx {\frac {7.2552}{A}}}$ .

For a more exact approximation, we use ansatz solution ${\displaystyle x=A_{1}\cos(\omega t)+A_{3}\cos(3\omega t)}$ . Plugging these in and matching the ${\displaystyle \cos(\omega t)}$ , ${\displaystyle \cos(3\omega t)}$  terms, we obtain after routine algebra:

The cubic equation for ${\displaystyle y}$  has only one real root ${\displaystyle y\approx 0.0448}$ . With that, we obtain an approximate period

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

${\displaystyle 0=F(t,x,{\dot {x}})}$ 

with a sufficiently smooth function ${\displaystyle F:\mathbb {R} \times \mathbb {C} ^{n}\times \mathbb {C} ^{n}\rightarrow \mathbb {C} ^{n}}$  where ${\displaystyle n}$  is the number of equations and ${\displaystyle t,x,{\dot {x}}}$  are placeholders for time, the vector of unknowns, and the vector of time derivatives.

The system is non-autonomous if the function ${\displaystyle t\in \mathbb {R} \mapsto F(t,x,{\dot {x}})}$  is not constant for (some) fixed ${\displaystyle x}$  and ${\displaystyle {\dot {x}}}$ . Nevertheless, we require that there is a known excitation period ${\displaystyle T>0}$  such that ${\displaystyle t\in \mathbb {R} \mapsto F(t,x,{\dot {x}})}$  is ${\displaystyle T}$ -periodic.

A natural candidate set for the ${\displaystyle T}$ -periodic solutions of the system equations is the Sobolev space ${\displaystyle H_{\rm {per}}^{1}((0,T),\mathbb {C} ^{n})}$  of weakly differentiable functions on the interval ${\displaystyle [0,T]}$  with periodic boundary conditions ${\displaystyle x(0)=x(T)}$ . We assume that the smoothness and the structure of ${\displaystyle F}$  ensures that ${\displaystyle F(t,x(t),{\dot {x}}(t))}$  is square-integrable for all ${\displaystyle x\in H_{\rm {per}}^{1}((0,T),\mathbb {C} ^{n})}$ .

The system ${\displaystyle B:=\left\{\psi _{k}\mid k\in \mathbb {Z} \right\}}$  of harmonic functions ${\displaystyle \psi _{k}:=\exp \left(ik{\frac {2\pi t}{T}}\right)}$  is a Schauder basis of ${\displaystyle H_{\rm {per}}^{1}((0,T),\mathbb {C} ^{n})}$  and forms a :Hilbert basis of the Hilbert space ${\displaystyle H:=L^{2}([0,T],\mathbb {C} )}$  of square-integrable functions. Therefore, each solution candidate ${\displaystyle x\in H_{\rm {per}}^{1}((0,T),\mathbb {C} ^{n})}$  can be represented by a Fourier-series ${\displaystyle x(t)=\sum _{k=-\infty }^{\infty }{\hat {x}}_{k}\exp \left(ik{\frac {2\pi t}{T}}\right)}$  with Fourier-coefficients ${\displaystyle {\hat {x}}_{k}:={\frac {1}{T}}\int _{0}^{T}\psi _{k}^{*}(t)\cdot x(t)dt}$  and the system equation is satisfied in the weak sense if for every base function ${\displaystyle \psi \in B}$  the variational equation

${\displaystyle 0=\langle \psi ,F(t,x,{\dot {x}})\rangle _{H}:={\frac {1}{T}}\int _{0}^{T}\psi ^{*}(t)\cdot F(t,x,{\dot {x}})dt}$ 

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 ${\displaystyle \psi }$  in ${\displaystyle B}$ .

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 ${\displaystyle B_{N}:=\{\psi _{k}\mid k\in \mathbb {Z} {\text{ with }}-N\leq k\leq N\}}$ .

This gives the finite-dimensional solution ${\displaystyle x(t)=\sum _{k=-N}^{N}{\hat {x}}_{k}\psi _{k}(t)=\sum _{k=-N}^{N}{\hat {x}}_{k}\exp \left(ik{\frac {2\pi t}{T}}\right)}$  and the finite set of equations

${\displaystyle 0=\langle \psi _{k},F(t,x,{\dot {x}})\rangle \quad {\text{ with }}k=-N,\ldots ,N}$ 

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:

1. Voltages ${\displaystyle V}$  are used to calculate the currents of the linear part, ${\displaystyle I_{\text{linear}}}$  in the frequency domain.
2. Voltages ${\displaystyle V}$  are then used to calculate the currents in the nonlinear part, ${\displaystyle I_{\text{nonlinear}}}$ . Since nonlinear devices are described in the time domain, the frequency-domain voltages ${\displaystyle V}$  are transformed into the time domain, typically using inverse Fast Fourier transforms. The nonlinear devices are then evaluated using the time-domain voltage waveforms to produce their time-domain currents. The currents are then transformed back into the frequency domain.
3. According to Kirchhoff's circuit laws, the sum of the currents must be zero, ${\displaystyle \epsilon =I_{\text{linear}}+I_{\text{nonlinear}}=0}$ . An iterative process, usually Newton iteration, is used to update the network voltages ${\displaystyle V}$  such that the current residual ${\displaystyle \epsilon }$  is reduced. This step requires formulation of the Jacobian ${\displaystyle {\tfrac {d\epsilon }{dV}}}$ .

Convergence is reached when ${\displaystyle \epsilon }$  is acceptably small, at which point all voltages and currents of the steady-state solution are known, most often represented as Fourier coefficients.