Skip to content

The method of harmonic balance

HarmonicBalance.jl focuses on systems with equations of motion that are subject to time-dependent terms of harmonic type.

Harmonic generation in oscillating nonlinear systems

Let us take a general nonlinear system of N second-order ODEs with real variables xi(t), i=1,2,,N, and time t as the independent variable,

x¨(t)+F(x(t),t)=0.

The vector x(t)=(x1(t),...,xN(t))T fully describes the state of the system. In physics, the coordinates x(t) encode the amplitudes of oscillators, e.g., mechanical resonators, voltage oscillations in RLC circuits, oscillating electrical dipole moments, or standing modes of an optical cavity.

As an example, let us first solve the harmonic oscillator in frequency space. The equation of motion is

x¨(t)+γx˙(t)+ω02x(t)=Fcos(ωdt),

where γ is the damping coefficient and ω0 the natural frequency. Fourier-transforming both sides of this equation yields

(ω02ω2+iωγ)x~(ω)=F2[δ(ω+ωd)+δ(ωωd)].

Evidently, x~(ω) is only non-vanishing for ω=±ωd. The system thus responds solely at the driving frequency, i.e., the time evolution can be captured by a single harmonic. This illustrates the general point that linear systems are exactly solvable by transforming to Fourier space, where the equations are diagonal.

The situation becomes more complex if nonlinear terms are present, as these cause frequency conversion. Suppose we add a quadratic nonlinearity βx2(t) to the equations of motion; attempting to Fourier-transform, now, leads to

FT[x2](ω)=x2(t)eiωtdt=+x~(ω)x~(ω)δ(ω+ωω)dωdω,

which couples all harmonics ω,ω,ω such that ω+ω+ω=0. To the lowest order, this means that the induced motion at the drive frequency generates a higher harmonic, ωd2ωd. The frequency conversion couples the response at different frequencies and propagates through the spectrum, thus, coupling an infinite number of harmonics. Hence, the system is not easily solvable in Fourier space anymore!

Harmonic ansatz & harmonic equations

Even though we need an infinity of Fourier components to describe our system exactly, some components are more important than others. The strategy of harmonic balance is to describe the motion of any variable xi(t) in a truncated Fourier space

xi(t)=j=1Miui,j(T)cos(ωi,jt)+vi,j(T)sin(ωi,jt).

Within this space, the system is described by a finite-dimensional vector

u(T)=(u1,1(T),v1,1(T),uN,MN(T),vN,MN(T)).

Under the assumption that u(T) evolves at much slower timescales than the oscillatory harmonics, we may neglect all of its higher order time derivatives. Notice that once the ansatz is used in the equations of motion, all terms become oscillatory - each prefactor of cos(ωi,jt) and sin(ωi,jt) thus generates a separate equation. Collecting these, we obtain a set of coupled 1st order nonlinear ODEs,

du(T)dT=F¯(u),

which we call the harmonic equations. The main purpose of HarmonicBalance.jl is to obtain and solve these Harmonic equations. We are primarily interested in steady states u0 defined by F¯(u0)=0.

The process of obtaining the harmonic equations is best shown using the example below:

Example: the Duffing oscillator

Here, we derive the harmonic equations for a single Duffing resonator, governed by the equation

x¨(t)+ω02x(t)+αx3(t)=Fcos(ωdt+θ).

As explained above, for a periodic driving at frequency ωd and a weak nonlinearity α, we expect the response at frequency ωd to dominate, followed by a response at 3ωd due to frequency conversion.

Single-frequency ansatz

We first attempt to describe the steady states of the Duffing equations of motion using only a single harmonic, ωd, in the ansatz, i.e., our starting point is the following harmonic ansatz for x

x(t)=u(T)cos(ωdt)+v(T)sin(ωdt),

with the harmonic variables u and v. The slow time T is, for now, equivalent to t. Substituting this ansatz into the Duffing equation of motion results in

[u¨+2ωdv˙+u(ω02ωd2)+3α(u3+uv2)4+Fcosθ]cos(ωdt)+[v¨2ωdu˙+v(ω02ωd2)+3α(v3+u2v)4Fsinθ]sin(ωdt)+α(u33uv2)4cos(3ωdt)+α(3u2vv3)4sin(3ωdt)=0.

We see that the x3 term has generated terms that oscillate at 3ωd, describing the process of frequency upconversion. We now Fourier-transform both sides of the above equations with respect to ωd to obtain the harmonic equations. This process is equivalent to extracting the respective coefficients of cos(ωdt) and sin(ωdt). Here, the distinction between t and T becomes important: since the evolution of u(T) and v(T) is assumed to be slow, they are treated as constant for the purpose of the Fourier transformation. Since we are interested in steady states, we drop the higher-order derivatives and rearrange the resulting equation to

ddT(uv)=18ωd(4v(ω02ωd2)+3α(v3+u2v)4Fsinθ4u(ωd2ω02)3α(u3+uv2)4Fcosθ).

Steady states can now be found by setting the l.h.s. to zero, i.e., assuming u(T) and v(T) constant and neglecting any transient behaviour. This results in a set of 2 nonlinear polynomial equations of order 3, for which the maximum number of solutions set by Bézout's theorem is 32=9. Depending on the parameters, the number of real solutions is known to be between 1 and 3 [see proof].

Sidenote: perturbative approach

The steady states describe a response that may be recast as x0(t)=X0cos(ωdt+ϕ), where X0=u2+v2 and ϕ=atan(v/u). Frequency conversion from ωd to 3ωd can be found by setting x(t)x0(t)+δx(t) with |δx(t)||x0(t)| and expanding the equations of motion to first-order in δx(t). The resulting equation

δx¨(t)+[ω02+3αX024]δx(t)=αX034cos(3ωdt+3ϕ),

describes a simple harmonic oscillator, which is exactly soluble. Correspondingly, a response of δx(t) at frequency 3ωd is observed. Since this response is obtained 'on top of' each steady state of the equations of motion, no previously-unknown solutions are generated in the process.

Two-frequency ansatz

An approach in the spirit of harmonic balance is to use both harmonics ωd and 3ωd on the same footing, i.e., to insert the ansatz

x(t)=u1(T)cos(ωdt)+v1(T)sin(ωdt)+u2(T)cos(3ωdt)+v2(T)sin(3ωdt),

with u1,u2,v1,v2 being the harmonic variables. As before we substitute the ansatz into the equations of motion, drop second derivatives with respect to T and Fourier-transform both sides. Now, the respective coefficients correspond to cos(ωdt), sin(ωdt), cos(3ωdt) and sin(3ωdt). Rearranging, we obtain

du1dT=12ωd[(ω02ωd2)v1+3α4(v13+u12v1+u12v2v12v2+2u22v1+2v22v12u1u2v1)+Fsinθ],dv1dT=12ωd[(ωd2ω02)u13α4(u13+u12u2+v12u1v12u2+2u22u1+2v22u1+2u1v1v2)Fcosθ],du2dT=16ωd[(ω029ωd2)v2+α4(v13+3v23+3u12v1+6u12v2+3u22v2+6v12v2)],dv2dT=16ωd[(9ωd2ω02)u2α4(u13+3u23+6u12u23v12u1+3v22u2+6v12u2)].

In contrast to the single-frequency ansatz, we now have 4 equations of order 3, allowing up to 34=81 solutions (the number of unique real ones is again generally far smaller). The larger number of solutions is explained by higher harmonics which cannot be captured perturbatively by the single-frequency ansatz. In particular, those where the 3ωd component is significant. Such solutions appear, e.g., for ωdω0/3 where the generated 3ωd harmonic is close to the natural resonant frequency. See the examples for numerical results.