Minimizing Aliasing in Multiple Frequency Harmonic Balance Computations

The harmonic balance method has emerged as an efficient and accurate approach for computing periodic, as well as almost periodic, solutions to nonlinear ordinary differential equations. The accuracy of the harmonic balance method can however be negatively impacted by aliasing. Aliasing occurs because Fourier coefficients of nonlinear terms in the governing equations are approximated by a discrete Fourier transform (DFT). Understanding how aliasing occurs when the DFT is applied is therefore essential in improving the accuracy of the harmonic balance method. In this work, a new operator that describe the fold-back, i.e. aliasing, of unresolved frequencies onto the resolved ones is developed. The norm of this operator is then used as a metric for investigating how the time sampling should be performed to minimize aliasing. It is found that a time sampling which minimizes the condition number of the DFT matrix is the best choice in this regard, both for single and multiple frequency problems. These findings are also verified for the Duffing oscillator. Finally, a strategy for oversampling multiple frequency harmonic balance computations is developed and tested.


Mathematics Subject Classification 65L70 1 Introduction
Many processes found in engineering and scientific applications are time-periodic.Some examples are fluid flows inside turbomachines, vibrations of structures, or voltages within AC circuits.In order to understand these processes, they may be simulated by integrating their governing equations in time until a periodic solution is obtained.In cases when the initial condition is far from being periodic and/or transient phenomena decay slowly, it may however take a long time until the solution reaches a periodic state.In terms of a numerical simulation that integrates between discrete time steps, this translates into a large number of iterations, and consequently, a high computational cost [10,17].
An alternative approach that can be more computationally efficient is to directly seek solutions that belong to some finite dimensional vector space spanned by a set of periodic basis functions.This paper considers the special case when these basis functions are sinusoids with different frequencies.In this case, the solution thus takes the form of a truncated Fourier series, for which the unknowns become the Fourier coefficients.There exist several methods in the literature for determining these Fourier coefficients.One approach is to integrate the equations obtained from substituting the Fourier series into the governing equations against each of the basis functions that span the vector space.This represents a Galerkin method, and will yield one equation per basis function.In cases when the governing equations are linear, the equations resulting from the integration uncouple with respect to the Fourier coefficients.If the governing equations on the other hand are nonlinear, each Fourier coefficient will in general be present in each of the final equations.When the Galerkin method is applied to linear problems it is sometimes referred to as the Linear Frequency Domain method [7,15].Some applications of the aforementioned approach to nonlinear problems can be found in e.g.[31,39].
Unfortunately, the Galerkin integral can become very complicated to solve analytically in the nonlinear case [16,31].This is especially true when the dimension of the vector space is large and/or when the governing equations are complex in nature.Other methods have therefore been developed for nonlinear problems.One common approach is to require that the equations resulting from the substitution only are satisfied at a discrete set of time instances, rather than in a variational sense as in the Galerkin method.This leads to a class of methods which here will be referred to as harmonic balance methods.Harmonic balance methods can be formulated in both the time and the frequency domain, with the main difference being that the solution variables are time samples in the former, and Fourier coefficients in the latter [10,26].Within the electronics community, harmonic balance is most often formulated in the frequency domain [10], and the modern version of harmonic balance is acredited to Nakhla and Valch [34].Within the fluid dynamics community, the time domain harmonic balance method was first introduced by Hall et al. [16] and the frequency domain formulation by McMullen [33].A good review of the history and theory of harmonic balance can be found in a two part review paper by Gilmore and Steer [10,11].
The harmonic balance method is based on applying a discrete Fourier transform (DFT) to either calculate the time derivative term (time-domain formulation), or the Fourier transform of the nonlinear terms (frequency domain formulation) [10,26].Applying a DFT to nonlinear problems can however introduce aliasing, which is known to reduce the accuracy of the harmonic balance method [24,29,32], and sometimes even lead to nonphysical solutions [28].Several methods have therefore been developed to address these problems.Huang and Ekici propose to add a time-spectral viscosity operator to the harmonic balance equations [19].This operator acts by damping the sinusoids which it is applied to, and has been shown by Huang and Ekici to improve convergence in cases when aliasing otherwise prevents it [19].Another approach investigated by La Bryer and Attair [28] is to filter the solution between each nonlinear iteration.The main idea of the filtering approach is to limit the amplitude of the sinusoids that have the highest frequencies, since these in general are the ones most contaminated by aliasing.In cases when the nonlinearity of the governing equations is known, exact frequency domain filters which keep the lower frequency sinusoids perfectly free from aliasing can in fact be constructed by using a sharp cut-off frequency based on Orszag's rule [28,36].It should also be noted that a frequency domain filter with a sharp cut-off frequency is equivalent to the oversampling strategy employed by Frey et al. in their harmonic balance solver [8].
In addition to filtering, a filtered inverse discrete Fourier transform has also been used to improve the convergence and stability of the harmonic balance method [3,18].This approach is particularly useful in cases when discontinuities and/or strong gradients lead to slow harmonic convergence (Gibbs phenomenon).The use of a filtered reconstruction operator is thouroughly discussed by Gottlieb and Shu [12].For its application to Fourier spectral methods alongside so-called reprojection approaches the reader to referred to [9].
Harmonic balance methods were originally developed for problems where the solution contains one fundamental frequency and its harmonics.Several problems in nature can however be expected to have solutions which contain a more general set of frequencies.In these cases, the solution belongs to a wider class of functions known as almost periodic functions [1,26].A fundamental problem that arises when the harmonic balance method is applied to almost periodic functions is to construct the DFT.In the original harmonic balance method, the standard discrete Fourier transform based on M uniformly distributed points between 0 and the reciprocal of the lowest resolved frequency is almost exclusively employed.Here, M ≥ 2K + 1, where K is the number of frequencies included in the Fourier series expansion.The reason for choosing this sampling is that it satisfies the assumptions of the Whittaker-Kotel'nikov-Shannon sampling theorem [23,38,40].A uniform time sampling that satisfies the sampling theorem can also be constructed for almost periodic signals as long as the frequencies considered are commensurable.Such a sampling may however require a very large number of sampling points, which makes this approach unfavourable from a computational perspective [21,26].Other approaches have therefore been developed for cases when the solution is an almost periodic function.Some of these are based on introducing a basis for the frequency set, such that each frequency can be written as an integer linear combination of a finite set of fundamental frequencies [26].The method developed by Frey et al. [8] takes into account only frequencies which are multiples of one the fundamental frequencies and applies the so-called harmonic set approach.The advantages of the harmonic set approach are that its computational cost scales linearly with the number of frequencies considered and that the standard discrete Fourier transform can be used within each harmonic set.The harmonic set approach does however neglect nonlinear coupling between different harmonic sets, except through common frequencies such as the zeroth frequency.These coupling terms may be accounted for by introducing several time variables, one for each fundamental frequency in the basis [14,22].This approach once again relies on the standard discrete Fourier transform with uniform time samples, but requires a substantially larger computational cost than the harmonic set approach.
A third approach, which is the main topic of this paper, is to employ the almost periodic Fourier transform (APFT) introduced by Kundert et al. [25,26] in the harmonic balance computation.The APFT differs from the standard discrete Fourier transform in two ways.First, it considers an arbitrary set of frequencies and is therefore well suited for cases when the solution is an almost periodic function.Secondly, it does not require that the time samples are distributed uniformly like in the standard discrete Fourier transform.Instead, Kundert et al. suggest that the time sampling is chosen so that it minimizes the condition number of the DFT matrix used in the APFT [25,26].Finding a time sampling that satisfies this criterion for an arbitrary set of frequencies is however a very complex problem, that to the authors' knowledge has not yet been solved analytically.In the original work on harmonic balance methods for arbitrary frequency sets, Chua and Ushida [2] employ a uniform discretization and use oversampling to avoid ill-conditioning.The same strategy has also been used in recent years by Ekici and Hall [5,6].Kundert et al. [25,26] were the first to obtain a well conditioned DFT matrix based on 2K + 1 time samples.This was done using their Near-Orthogonal Selection Algorithm.Since the work of Kundert et al., several other algorithms have been developed to compute a time sampling that minimizes the condition number of the DFT matrix, see e.g.[13,21,27,35,37].
The rationale of choosing a time sampling that minimizes the condition number of the DFT matrix is that it limits the effect that a perturbation of the sampled signal can have on the resulting DFT.Since all unresolved frequencies manifest themselves as perturbations of the sampling, this implies that a low condition number can limit the amount of aliasing produced by the DFT [26].In [26], it is also noted that the condition number does not say anything about how much the amplitude of a particular unresolved sinusoid affects the resolved ones.This implies that the condition number in itself can not be used to identify an alias free DFT.In order to do this, the mapping of the unresolved sinusoids onto the resolved ones must be defined [26].In this work, an operator that describes this mapping has been developed.The benefit of this operator, hereinafter referred to as the alias operator, is that it directly defines how a particular time sampling affects aliasing.As will be shown later, the norm of this operator also provides an a-priori bound on the amount of aliasing that a particular time sampling gives rise to.
In this paper, the relation between the alias operator and the condition number of the DFT matrix employed in the APFT will be investigated.After this, the relation between the norm of the alias operator, the condition number of the DFT matrix and the actual alias error obtained in several harmonic balance computations will be investigated.Finally, the benefits of employing oversampling in multiple frequency harmonic balance computations based on the APFT will be investigated.

Harmonic Balance Method
The purpose of this paper is to study how aliasing occurs in harmonic balance solutions of first order ordinary differential equations of the following form Here, q ∈ R N is a vector that contains the unknown solution and f : R N × R → R N is a nonlinear function that describes the dynamics of the problem.The nonlinearity of f (q, t) implies that an almost periodic solution to Eq. ( 1) can contain an infinite number of sinusoids with different frequencies.Computing the amplitude of all these sinusoids is however not feasible from a numerical perspective.Therefore, approximate solutions that are spanned by a limited number of sinusoids are considered instead Here, Λ represents the set of frequencies included in the series expansion Note that the only requirement put on these frequencies is that ω −k = −ω k .If a Galerkin projection of Eq. ( 1) onto the subspace spanned by the sinusoids in Eq. ( 2) is performed, the following is obtained Here, Ω Λ = Ω Λ ⊗I, where ⊗ is the Kronecker product, I is an identity matrix of size N × N , and The vectors qΛ and fΛ ( qΛ , Λ) in Eq. ( 4) further consist of 2K + 1 sub-vectors, in which the kth sub-vector contains qk and fk respectively.Calculating the Galerkin projection of Eq. ( 1) can be very complicated in cases when the number of sinusoids considered in Eq. ( 2) is large, and/or the function f (q, t) is very complex in nature.In order to overcome this difficulty, the harmonic balance method approximates fΛ ( qΛ , Λ) by a discrete Fourier transform (DFT) of f (q, t) instead The DFT of f (q, t) in Eq. ( 6) is performed in two steps.In the first step, f (q, t) is evaluated at a set of time instances This will in turn require the realization of Eq. ( 2) at these time instances Here, In the second step, the sampling of f (q, t) is transformed back into the frequency domain using the inverse of E −1 Λ (t), here denoted E Λ (t).For this inverse to be well defined, the columns in E −1 Λ (t) must be linearly independent.Due to the structure of E −1 Λ (t), this holds true if M = 2K + 1 time instances can be found such that the columns in E −1 Λ (t) are linearly independent.As it turns out, this is always possible.In order to prove this, note that for an equidistant sampling, the entries of E −1 Λ (t) take the form For an equidistant time sampling, E −1 Λ (t) thus corresponds to a transposed Vandermonde matrix.Its determinant is therefore given by It is easy to verify that this expression is non-zero for almost all choices of Δt.This shows that for any given set of frequencies one can find an equidistant sampling point distribution t such that E −1 Λ (t), and thereby E −1 Λ (t), is invertible.The time sampling is usually chosen to minimize the condition number of E −1 Λ (t).Note that a small condition number automatically ensures that E −1 Λ (t) is invertible.Some authors use a numerical optimization procedure to minimize the condition number.It is important to note that the cost of this optimization procedure is usually negligible compared to solving the harmonic balance system.
In the literature, the discrete Fourier transform defined by the matrices E −1 Λ (t) and E Λ (t) is often referred to as the almost periodic Fourier transform (APFT) [26].The APFT represents a generalization of the standard DFT to arbitrary sets of frequencies and sampling points, and will in fact be equivalent to the standard DFT when the frequencies in Λ are integer multiples of a single base frequency, and the time samples are distributed uniformly between zero and the reciprocal of this base frequency.
Equation ( 6) represents the frequency domain formulation of the harmonic balance method.An equivalent formulation in the time domain may also be expressed as follows This formulation is easily obtained by premultiplying Eq. ( 6) from the left by E −1 Λ (t) and then using Eq. ( 8).The matrix E −1 Λ (t)E Λ (t) in the above equation will further be equal to the identity matrix when M = 2K +1 time samples are employed.If more than M = 2K +1 time samples on the other hand are used, this matrix will represent a projection that corresponds to the application of a modal filter.
The frequency domain formulation of the harmonic balance method presented in Eq. ( 6) will be used for the remainder of this paper.Due to the equivalence of Eq. ( 6) and ( 12), however, the results presented will also apply to the time domain formulation of the harmonic balance method.

Error Sources
It is well known that the harmonic balance method can give rise to two types of errors: aliasing and harmonic truncation [10,11,26].Aliasing occurs because f (q, t) can contain more frequencies that those accounted for by the DFT in Eq. ( 6).Aliasing will on the other hand not occur if Eq. ( 4) is solved.This is because the Galerkin projection in this case corresponds (up to a constant) to the standard formula for calculating Fourier coefficients.Independently on whether the harmonic balance or Galerkin method is used, however, the harmonic truncation error will always be present.This error occurs as a result of the fact that the contribution of the unresolved sinusoids to the solution, including their nonlinear coupling with the resolved sinusoids, is neglected.
Both the alias error and the harmonic truncation error can naturally be reduced by including more frequencies in Λ.This is however not always possible in practical applications since the computational cost of the harmonic balance method scales almost linearly with the size of Λ, here denoted #Λ.It is therefore often hard to avoid aliasing and harmonic truncation errors in the solution.This highlights the importance of understanding how these errors occur, and how they can be limited.In this paper, the focus is put on the aliasing error.

Alias Operator
As noted in the previous section, aliasing occurs when f (q, t) contains more frequencies than those accounted for by the DFT in Eq. ( 6).Let Λ denote the union of Λ and the set of all frequencies that are contained in f (q, t).In cases when f (q, t) is a polynomial in q, with coefficients that are finite Fourier sums in t, the set Λ will be finite.This in turn implies that the signal f (q, t) can be written as Note that Λ Λ .From this, it follows that Eq. ( 13) can be rewritten as Based on this expression, the DFT of f (q, t) may now be expressed as This shows that the Fourier coefficients obtained from the DFT of f (q, t) consist of two terms.The first term represents the correct (alias free) Fourier coefficients corresponding to the resolved frequencies, and the second term the fold-back of the remaining Fourier coefficients onto the resolved ones.The operator E Λ (t)E −1 Λ \Λ (t) that defines the fold-back of the unresolved Fourier coefficients in Eq. ( 15) will be referred to as the alias operator for the remainder of this paper.Based on Eq. ( 15), it is easy to show that the norm of this operator puts the following bound on the alias error Provided that Λ is known, this relation may be used to obtain an a-priori bound on the alias error that a particular time-sampling t can give rise to.As will be shown next, however, determining Λ for a particular problem is not a trivial task.

Selection of Unresolved Frequencies
The set Λ is completely defined by the function f (q, t) and the set Λ.In the present work, it is assumed that f (q, t) is a polynomial of degree p in q.This assumption allows the set Λ to be calculated by simple addition and subtraction of the frequencies in Λ.Clearly, this assumption is not valid for all problems.Despite this, however, it can still be a good approximation provided that f (q, t) itself approximates an analytic function, i.e., one that is locally given by a convergent power series (which then coincides with its Taylor series).The choice of Λ directly influences how well the norm of the alias operator will represent the actual alias error obtained in a harmonic balance simulation.To see this, note that Eq. ( 16) represents the largest possible alias error that any fΛ \Λ can give rise to.This means that the bound in Eq. ( 16) may not be representative if some elements in fΛ \Λ are negligible for a given problem.It is also important to keep in mind that Eq. ( 16) gives a relative bound on the alias error.As such, if Λ has been selected to contain all relevant frequencies for a given problem, then fΛ \Λ 2 ≈ 0 and consequently, the actual alias error will be small independent of the norm of the alias operator.

Relation to Condition Number
The condition number of E −1 Λ (t) has been successfully used by several authors in the past as a measure for selecting the time sampling in multiple frequency harmonic balance computations [13, 21, 25-27, 35, 37].This raises the question of how the condition number relates to the alias operator defined in the present work.The condition number can be defined based on the l 2 norm as This definition can be used to derive a bound on the norm of the alias operator as follows Both matrices in the last equation have entries with norm 1.From this, it follows that there exists an upper bound on E −1 Λ \Λ (t) 2 , and a lower bound on E −1 Λ (t) 2 , both independent of t.Equation may therefore be rewritten as This equation shows that the condition number in fact puts a bound on the aliasing operator.Note, however, that this does not imply that a time sampling which minimizes the condition number also minimizes the norm of the alias operator.
The complexity of the alias operator, condition number and α Λ, Λ , t in Eq. ( 19) makes it very hard to investigate how they relate analytically.Numerical optimizations were therefore used to search for time samplings which give a minimal condition number, minimal norm of the alias operator or maximal value of α Λ, Λ , t .These results could then be used to investigate whether a time sampling that gives, e.g., a low condition number, also yields a low norm of the alias operator.For a detailed description of the optimization procedure used in the present work, see "Appendix B".

Almost Periodic Fourier Transform with Oversampling
In this work, the impact of oversampling on the alias error obtained in a harmonic balance computation has also been investigated.The motivation for doing this is that aliasing can never be eliminated when M = 2K + 1 time samples are employed.This can be proven by first noting that the columns in E −1 Λ (t) are required to be linearly independent, and therefore form a basis for C M .From this, it follows that there exists a nonzero matrix W such that E −1 Λ \Λ (t) = E −1 Λ (t)W .But then, the alias operator in Eq. ( 15) must be equal to E Λ (t)E −1 Λ \Λ (t) = W = 0, which proves the desired result.When oversampling is employed, E Λ (t) in Eq. ( 6) will represent a left inverse to E −1 Λ (t).Contrary to the standard inverse of a square matrix, however, the left inverse is not uniquely defined.In cases when Λ is finite, it can therefore be tailored such that aliasing is eliminated.In order to prove this, one can start by noting that such a tailored inverse, here denoted E Λ,Λ (t), must satisfy the following two conditions The above conditions can be combined as follows Note that the matrix within the brackets in the left hand side of the above equation corresponds to E −1 Λ (t).From before, it is known that there exists a time sampling with M = #Λ time samples such that this matrix is invertible.From this, it follows that the above relation may be rewritten as Here, π Λ ,Λ is the projection onto the harmonics in Λ.
The argument above shows that a left inverse that eliminates aliasing always can be defined if Λ is finite.Unfortunately, it is more computationally expensive to employ this left inverse since it requires that f (q, t) is evaluated at M > 2K + 1 time instances.This motivates the use of a left inverse which requires less than M time instances, but still has the potential to reduce aliasing compared to not using oversampling at all.In order to define such a left inverse, introduce a new set of frequencies Λ that satisfies There exists a time sampling with M = #Λ time samples such that E −1 Λ (t) is invertible.Therefore, the following left inverse may be defined The idea behind this left inverse is to eliminate aliasing with respect to the frequencies in Λ \Λ.This left inverse would also completely eliminate aliasing if a time sampling can be found such that the unresolved Fourier coefficients only contribute to the amplitude of the Fourier coefficients whose frequencies are in Λ \Λ.
An alternative choice for the left-inverse which has been suggested in the literature is the Moore-Penrose inverse of E −1 Λ (t) [5,6,21].This left inverse represents a least squares projection onto the subspace spanned by the columns in E −1 Λ (t).From this, it follows the Moore-Penrose inverse can only eliminate aliasing if the columns in E −1 Λ \Λ (t) are orthogonal to those in E −1 Λ (t) (with respect to the standard inner product on C M ).In cases when the frequencies in Λ are commensurable, and f (q, t) is a polynomial, it is possible to ensure this by selecting a uniform time sampling that satisfies Orszag's rule [21] M ≥ ( p + 1) Here, ω max is the largest frequency in Λ and ω base is the greatest common divisor of all frequencies in Λ.
It might also be possible to construct a (possibly non-uniform) time sampling such that the Moore-Penrose inverse eliminates aliasing in the general case, but no such algorithm is known to the authors.As shown in "Appendix A", it is however possible to redefine the inner product on C M such that the columns in E −1 Λ \Λ (t) become orthogonal to those in E −1 Λ (t).There, it is also shown that the Moore-Penrose inverse defined with respect to this new inner product is equivalent to the left inverse defined by Eq. ( 23).

Aliasing for M = #Λ
This section considers the case when the number of sampling points corresponds to the theoretical minimum of M = #Λ.To begin with, it is demonstrated how κ(E −1 Λ ) and Λ \Λ 2 relate in this case.This is done for both the case when Λ contains a single as well as multiple fundamental frequencies.After this, it is investigated how κ(E −1 Λ ) and E Λ E −1 Λ \Λ 2 influence the alias error obtained when the Duffing oscillator is simulated with the harmonic balance method.This is done for both the single and multiple frequency case as well.

Relation Between Ä(E
In Sect.2.2.2 it was shown that E Λ E −1 Λ \Λ 2 is bounded by κ(E −1 Λ ) through Eq. ( 19).This fact can be illustrated by a diagonal line in a κ(E −1 Λ )− E Λ E −1 Λ \Λ 2 plane, if the slope of this line is taken to be the maximum value of α Λ, Λ , t .In the same plane, the minimal values of κ(E −1 Λ ) ≥ 1 and E Λ E −1 Λ \Λ 2 ≥ 0 1 can also be represented by a vertical and a horizontal line respectively.The region between these three lines will then contain all possible values of κ(E −1 Λ ) and E Λ E −1 Λ \Λ 2 that a time sampling can give rise to.In Fig. 1a, this region is illustrated for the case when Λ contains a single fundamental frequency and Λ is generated by a second degree polynomial ( p = 2).The case when Λ is generated by a third degree polynomial ( p = 3) is further depicted in Fig. 1b.The limit values of κ and α Λ, Λ , t shown in these figures have all been computed with the OPT algorithm [13], as described in "Appendix B".A large set of random time samples was also generated to investigate how κ(E −1 Λ ) and the region between the three lines in Fig. 1.They also indicate that the maximal norm of the alias operator is proportional to κ(E −1 Λ ), as would be expected from Eq. (19).For the single frequency case shown in Fig. 1, it did in fact turn out that both the minimum value of κ(E −1 Λ ) and the minimum value of E Λ E −1 Λ \Λ 2 were obtained with M = #Λ uniformly distributed time samples between 0 and T .The results in Fig. 1 thus point to the fact that the standard DFT is optimal in terms of aliasing for the single frequency problem.Several other tests with different numbers of harmonics in Λ and different sets Λ have confirmed this statement.
The relation between κ(E −1 Λ ) and E Λ E −1 Λ \Λ 2 for the case when Λ contains two fundamental frequencies and Λ is generated by a second degree polynomial ( p = 2) is depicted in Fig. 2a.The case when Λ is generated by a third degree polynomial ( p = 3) is further depicted in Fig. 2b.The results shown in Fig. 2 very much resemble the single frequency case in Fig. 1.A couple of differences can however be noted.To begin with, the minimal value of E Λ E −1 Λ \Λ 2 can be seen to increase when going from the single frequency to the multiple frequency case.For both cases, it can also be seen that the minimal value of increases when the nonlinearity of the problem increases.Both these trends can be explained by the fact that the number of unresolved frequencies increases when the frequencies in Λ are no longer harmonically related and/or when the nonlinearity of the problem increases.In addition, it can be seen in Fig. 2 that the time sampling which minimizes κ(E −1 Λ ) (blue square marker) and the one that minimizes E Λ E −1 Λ \Λ 2 (yellow circle marker) are not the same in the multiple frequency case.The difference between these two samplings can however be seen to be very small.One possible explanation of this difference could be that the OPT algorithm has not converged.Alternatively, there may be a trade off between κ(E −1 Λ ) and E Λ E −1 Λ \Λ 2 for the multiple frequency case.The trend of the randomly generated samplings shown in Fig. 2 does however suggest that this trade-off is very small.Overall, Figs. 1 and 2 suggest that κ(E −1 Λ ) is a reasonable metric for selecting a time sampling when M = #Λ time samples are used in the single and multiple frequency case.These figures also show that it is important that the optimization algorithm used to minimize κ(E −1 Λ ) is able to come close to the global optimum, since the time sampling will otherwise not minimize the norm of the alias operator.

Duffing Oscillator
The Duffing oscillator is a model for a damped and driven oscillator with a nonlinear stiffening spring.The equation governing the displacement of the weight in the presence of two driving forces may be written as Here, ζ and F i respectively denote the damping coefficient and the amplitude of the ith driving force.Furthermore, a dot represents differentiation with respect to time.Equation ( 27) may be rewritten as a first order ordinary differential equation by introducing the velocity of the weight, y = ẋ, as an additional degree of freedom.This yields the following system of equations where Equation ( 28) is discretized in time using the frequency domain harmonic balance method presented in Eq. ( 6).The resulting nonlinear system of equations have been implemented in the Python programming languange and solved using the fsolve routine available in SciPy [20].The exact Jacobian is also provided to fsolve to avoid problems associated with finite-difference approximations.
The implementation of the frequency domain harmonic balance solver for the Duffing oscillator has been validated against the publically available MATLAB tool NLvib [24].To begin with, only a single driving force is considered (F 2 = 0).The presence of a single driving force admits the standard DFT to be used in the harmonic balance method.Figure 3a shows a comparison between the Python solver and NLvib for ζ = 0.1, F 1 = 1.25 and K = 5.Both solutions were obtained by starting with a driving frequency close to 0 and then gradually increase the frequency.Each new frequency computation is then started from the previous solution.Figure 3a shows that the Python solver and NLvib yield identical results for this case.Note, however, that the hysteresis branch could only be computed with NLvib since no arc length continuation method was implemented in the Python solver.It should also be noted that both solutions were computed with M = 21 time instances to ensure that Orszag's rule (Eq.( 26)) was satisfied for the cubic nonlinearity present in the governing equations.
Figure 3b shows results obtained with the Python solver for a multiple frequency problem.In this case, two incommensurable driving frequencies with a ratio of ω 2 /ω 1 = √ 2 are used.The damping factor and amplitudes of the driving forces are further set to ζ = 0.1 and F 1 = F 2 = 0.1 respectively.No higher harmonics of the driving frequencies are considered in the computation.The presence of a cubic nonlinearity in Eq. ( 29) will however generate a large number of additional frequencies.To avoid aliasing, the DFT defined by Eq. ( 23) is used.This requires the use of M = #Λ = 25 sampling points, which were selected to minimize the condition number of the DFT matrix.Unfortunately, NLvib does not support the APFT.direct comparison between the two solvers was therefore not possible in the multiple frequency case.Instead, two NLvib computations were run at each of two driving frequencies and with M = 5 time samples to obtain alias free solutions.These results are provided for reference in Fig. 3b.This figure clearly shows that the APFT solution and the reference single-frequency solution deviate.To ensure that this deviation is due only to the nonlinear coupling resolved by the APFT, the Python solver was rerun with F 2 = 0 and all other settings intact.The results from this computation can be seen to agree perfectly with the corresponding NLvib solution.
As shown previously, both E Λ E −1 Λ \Λ 2 and κ(E −1 Λ ) put a bound on the amount of aliasing produced by a DFT.These metrics will now be compared to the actual amount of aliasing obtained with the aforementioned frequency domain harmonic balance solver for the Duffing oscillator.The comparison has been performed for the two sets of frequencies that were previously used for the validation.For each of these sets, 90 random time samplings were generated by drawing M = #Λ samples from a uniform distribution and then checking the condition number of the corresponding DFT matrix.If the condition number was below 10, the sampling was saved.Otherwise, another sampling was generated until 90 samplings had been obtained in total.
For both the single and multiple frequency cases the corresponding alias free discrete Fourier transforms that were used for the validation were first used to compute a set of reference solutions using the same frequency stepping procedure as described previously.These reference solutions were then used as initial conditions for the simulations that employed the random samplings.For each frequency, several solutions in which the random samplings were shifted in time were computed.These time shifts were employed to account for the fact that the phase of the solution can affect the amount of aliasing obtained.Once all frequency and time shift combinations had been computed for a given random time sampling,

(b) (a)
Fig. 4 Numerical alias error for the Duffing oscillator plotted against κ(E −1 Λ ) and the following numerical alias error was calculated Here, APFT and OS respectively denote the random sampling and the alias-free oversampled solution.The vector xΛ in this equation further contains all Fourier coefficients of the displacement.
The numerical alias error is plotted against κ(E −1 Λ ) with circular markers for the single frequency case in Fig. 4a.This figure shows that κ(E −1 Λ ) puts a bound on the amount of aliasing obtained in a harmonic balance simulation, which is consistent with the results presented in Fig. 1.The red diamond marker and blue square marker in Fig. 4a further denote the alias free reference solution and a solution obtained with the standard DFT.It is interesting to note that the standard DFT gives the lowest numerical alias error among all samplings with M = #Λ, which is consistent with the results in Fig. 1.From Fig. 4a it can also be noted that κ(E −1 Λ ) = 1 for both the oversampled reference solution as well as the standard DFT based on M = #Λ time samples.This result points to the fact that κ(E −1 Λ ) can not distinguish an alias free solution from an aliased one.This is on the other hand possible when the numerical alias error is plotted against E Λ E −1 Λ \Λ 2 , as shown in Fig. 4b.This figure also shows that Λ \Λ 2 puts a bound on the numerical alias error, which is consistent with Eq. ( 16).In Fig. 5a and b the numerical alias errors obtained for the multiple frequency case have been plotted against κ(E −1 Λ ) and E Λ E −1 Λ \Λ 2 respectively.These figures very much resemble their single frequency counterparts in the sense that both κ(E −1 Λ ) and E Λ E −1 Λ \Λ 2 can be seen to bound the numerical alias error.The square and triangle markers in Fig. 5 correspond to two samplings which give the lowest value of κ(E −1 Λ ) and E Λ E −1 Λ \Λ 2 respectively.These two results can be seen to be very close to each other.This once again indicates that the trade off between κ(E −1 Λ ) and E Λ E −1 Λ \Λ 2 is very small in the multiple frequency case, and thus that it is sufficient to optimize for κ(E −1 Λ ) in order to obtain a DFT that minimizes aliasing.Note however that this statement only holds true when M = #Λ sampling points are employed since the minimal value of κ(E −1 Λ ) can be the same with or without oversampling.

Aliasing for M > #Λ
Two strategies for employing oversampling with the APFT were presented in Sect.2.3.In the first strategy, E Λ (t) is constructed by explicitly accounting for some of the unresolved frequencies (Eq.( 25)), and in the second strategy, E Λ (t) is computed as the Moore-Penrose inverse of E −1 Λ (t).In this section, it will be investigated whether the additional knowledge about Λ built into the first definition of E Λ (t) is beneficial from an aliasing perspective.This will be done by investigating how the alias error obtained when the Duffing oscillator is simulated with the harmonic balance method varies M for both a single and a multiple frequency case.In relation to this, it is also shown how the numerical alias error relates to the norm of the alias operator and the condition number of the DFT matrix by plotting their variation with M as well.

Duffing Oscillator
The first case that has been used to investigate the impact of oversampling with the APFT is the Duffing oscillator with a single driving force.The alias free reference solution for this case was computed with a uniform sampling that satisfies Orszag's rule (Eq.( 26)) for cubic nonlinearities.The two definitions of E Λ (t) presented in Sect.2.3 were evaluated using two different sets of time samples.The first set consists of M uniformly distributed points between 0 and the reciprocal of the lowest frequency.The second set consists of M non-uniformly distributed points created by perturbing each point in the first sampling set by a random value drawn from the uniform distribution U (−0.15T /M , 0.15T /M ).When E Λ (t) was computed based on Eq. ( 25), Λ was selected to contain the M lowest frequencies in Λ .
Once all samplings had been generated for all values of M , the numerical alias error was calculated according to the procedure outlined in Sect.3.1.2.The results from this calculation are presented in Fig. 6a.This figure shows that the two definitions of E Λ (t) give equivalent results for the uniform sampling.This is because both these discrete Fourier transforms are equivalent to the standard DFT in this case.This fact also explains why both definitions of E Λ (t) eliminate aliasing once Orszag's rule is satisfied.From Fig. 6a it can also be seen that aliasing is not eliminated for any value of M when a non-uniform sampling is used and E Λ (t) is selected to be the Moore-Penrose inverse of E −1 Λ (t).This is on the other hand possible when all the information about Λ is built into the definition of E Λ (t) for the non-uniform sampling.
The variation of E Λ E −1 Λ \Λ 2 with M for the single frequency case is further presented in Fig. 6b.This figure shows that the norm of the alias operator follows the same trend as the numerical alias error when the number of sampling points increase.In particular, it can be seen that both E Λ E −1 Λ \Λ 2 and the numerical alias error becomes zero for the same number of sampling points.This highlights the fact that E Λ E −1 Λ \Λ 2 is able to identify an alias free sampling.As noted previously, this is on the other hand not possible if only the condition number is considered.To demonstrate this, κ(E −1 Λ ) is plotted against M in Fig. 6c.In general, Fig. 6 indicates that both the definition of E Λ (t) as well as the choice of time sampling can have a substantial impact on the amount of aliasing obtained in a harmonic balance simulation.This raises the question regarding which combination is best.Based only on the results shown in Fig. 6, it appears that it is beneficial to define E Λ (t) based on Eq. (25).It also appears that the uniform sampling, which gives the lowest value of κ(E −1 Λ ), is a better choice than a non-uniform sampling.This indicates that the best choice in terms of aliasing is to construct E Λ (t) based on Eq. ( 25), and then select a time sampling which minimizes κ(E −1 Λ ).Here, care must however be taken in the definition of the condition number.In this work, it is suggested to compute the condition number based on the square matrix E −1 Λ (t), and not based on Eq. ( 17) in which some columns/rows of E −1 Λ /E Λ (t) are neglected.The first reason for this is that, in practice, E Λ (t) is computed from a numerical inverse of E −1 Λ (t) (see Eq. ( 25)).For this preliminary step to be well-posed it is thus necessary that κ(E −1 Λ ) is small.The second reason is that a minimal value of the condition number has been seen to make the APFT robust against aliasing for the non-oversampled cases, independently on the choice of Λ .In the suggested approach, E Λ (t) is thus computed by first computing a square matrix E Λ (t) that gives as little aliasing as possible, and then explicitly eliminate aliasing with respect to some of the unresolved Fourier coefficients, i.e., those that correspond to frequencies in Λ \Λ.
In order to investigate the usefulness of the suggested approach for computing E Λ (t), the Duffing oscillator with two incommensurable driving frequencies was once again considered.The alias free reference solution for this case was computed using the DFT defined in Eq. ( 23).The time sampling for this case was selected by minimizing κ(E −1 Λ ) using the OPT algorithm [13].For each M , two different discrete Fourier transforms were then constructed.The first one was calculated as the Moore-Penrose inverse of E −1 Λ (t), using a time sampling which minimized κ(E −1 Λ ), and the second one using Eq. ( 25) and a time sampling that minimized κ(E −1 Λ ).These minima were also obtained with the OPT algorithm.The variation of the numerical alias error with M for the multiple frequency case is presented in Fig. 7a.From this figure, it can be seen that for this case the suggested approach effectively eliminates aliasing once M ≥ 17.Calculating E Λ (t) as the Moore-Penrose inverse of E −1 Λ (t) does on the other hand not eliminate aliasing for any M .These results once again indicate that the definition of E Λ (t) can have a significant effect on aliasing.The fact that aliasing also is effectively eliminated despite the fact that not all the unresolved frequencies are accounted for in the construction of E Λ (t) is also interesting, since this suggests that there could exist a generalization of Orszag's rule for multiple frequency problems.That is, it is not necessary to construct E Λ (t) in Eq. ( 25) such that all the unresolved sinusoids can be distinguished, only such that the unresolved ones alias onto those with frequencies in Λ \Λ.
The variation of E Λ E −1 Λ \Λ 2 with M for the multiple frequency case is shown in Fig. 7b.This figure shows that the trend of E Λ E −1 Λ \Λ 2 once again follows that of the numerical alias error as the number of sampling points increase.
A comment on the point corresponding to M = 15 in Fig. 7 should also be made.In this point, it can be seen that the numerical alias error is the largest for the case when E Λ (t) is computed based on Eq. (25).At first, this might seem counter intuitive since more information about the unresolved sinusoids is included in the DFT compared to e.g.M = 5.A possible explanation for this can be found in Fig. 7c.This figure shows that the largest value of κ(E −1 Λ ) is obtained for M = 15.As such, it can be expected that E Λ (t) is the least robust against aliasing for M = 15.Thus, if not all the unresolved frequencies alias onto the frequencies in Λ \Λ, which can not be guaranteed, this sampling may very well generate more aliasing than another one that uses less sampling points and have a lower condition number.It is however not known at this point whether the relatively large value of κ(E −1 Λ ) for M = 15 is due to the fact that the OPT algorithm was unable to find the global optima or not.

Conclusions
It is well known that aliasing may occur when the harmonic balance method [16,33,34] is applied to nonlinear problems.In the present paper, aliasing is studied in detail for the case when the harmonic balance method is used to solve first order ordinary differential equations which are nonlinear functions of the solution variables.It is shown that aliasing occurs as a result of the fact that the harmonic balance method approximates the Fourier coefficients of a nonlinear function by a discrete Fourier transform (DFT).Although this is already well known, the aliasing of each unresolved frequency onto the resolved ones is here given a precise meaning by introducing a new operator that governs aliasing.This operator, referred to as the alias operator, is derived for the general case when the almost periodic Fourier transform (APFT) [26] is used in the harmonic balance method.As such, the alias operator can be used to study aliasing in both single frequency as well as multiple frequency harmonic balance computations.Using the norm of the newly introduced alias operator as a metric, the best sampling strategy for single and multiple frequency harmonic balance computations is investigated.It is found that for single frequency problems, the widely adopted uniform sampling approach is optimal.This sampling is also known to minimize the condition number of the DFT matrix [26].For the general case with multiple fundamental frequencies, the time sampling that minimizes the condition number of the DFT matrix seems to be close, but not identical, to the time sampling that minimizes the norm of the alias operator.A low condition number and norm of the alias operator is also shown to give a smaller alias error in numerical simulations of the Duffing Oscillator.The results presented in this paper therefore demonstrate that the condition number of the DFT matrix is a reasonable metric for selecting a time sampling, as has previously been done by e.g.[13, 21, 25-27, 35, 37].
Another important finding from the present work is that the alias error is only bounded by, but not strongly correlated with, the condition number.This does in turn imply that the time sampling must give a condition number which is close to the global optimum if the alias error is to be minimized.Therefore, in cases when an optimization algorithm is used to find the time sampling, it is important that it converges well.In the present work, it is found that the OPT algorithm introduced in [13] performs very well.
Employing the APFT in combination with oversampling is also considered in this work.In this case, the DFT matrix is defined as a left inverse of the inverse DFT matrix.Previous studies have used the Moore-Penrose inverse for this purpose [5,6,21].In the present work, however, it is shown that this left inverse in general can not eliminate aliasing in the multiple frequency case.Based on this finding, a new left inverse that can eliminate aliasing completely is introduced.This new left inverse is also shown to perform better in general than the Moore-Penrose inverse on the Duffing Oscillator.Based on these observations, it is suggested that the newly introduced left inverse is used in favor of the Moore-Penrose inverse when oversampling is employed in multiple frequency harmonic balance computations based on the APFT.
The above findings for the case of oversampling apply both to the abstract metric defined by the norm of the aliasing operator and to the measured alias error in the numerical experiments.In contrast, the conditon number could not be used as a criterion to distinguish samplings that were optimal with regard to aliasing.It is therefore concluded that the norm of the alias operator contains valuable information about the quality of a given sampling strategy.
It should finally be noted that both the alias error and the harmonic truncation error encountered in the harmonic balance method can be reduced by increasing the number of frequencies in the computation.In many practical applications, this is unfortunately not feasible, and good control over the alias error is thus important to achieve accurate results.In cases when more frequencies can be afforded, however, one must be mindful of the computational complexity of the APFT, which scales as O(M 2 ).To overcome this limitation, the Non Uniform Fast Fourier Transform (NUFFT) of type three could be used instead [4,30].The theory developed in this paper would however need to be adapted for the NUFFT.