This version is available at https://doi.org/10.14279/depositonce-9832
This work is licensed under a CC BY-NC-ND 4.0 License (Creative
Commons Attribution-NonCommercial-NoDerivatives 4.0 International). For
more information see https://creativecommons.org/licenses/by-nc-nd/4.0/.
Terms of Use
von Wagner, U., & Lentz, L. (2019). On the detection of artifacts in Harmonic Balance solutions of
nonlinear oscillators. Applied Mathematical Modelling, 65, 408–414. https://doi.org/10.1016/j.
apm.2018.08.013
Utz von Wagner, Lukas Lentz
On the detection of artifacts in Harmonic
Balance solutions of nonlinear oscillators
Accepted manuscript (Postprint)Conference paper |
On the Detection of Artifacts in Harmonic Balance
Solutions of Nonlinear Oscillators
U. von Wagnera,∗, L. Lentza
aInstitut f¨ur Mechanik, Technische Universit¨at Berlin, Einsteinufer 5, 10587 Berlin, Germany
Abstract
Harmonic Balance is a very popular semi-analytic method in nonlinear dynamics. It is
easy to apply and is known to produce good results for numerous examples. Adding
an error criterion taking into account the neglected terms allows an evaluation of the
results. Looking on the therefore determined error for increasing ansatz orders, it can
be evaluated whether a solution really exists or is an artifact. For the low-error solu-
tions additionally a stability analysis is performed which allows the classification of
the solutions in three types, namely in large error solutions, low error stable solutions
and low error unstable solution. Examples considered in this paper are the classical
Duffing oscillator and an extended Duffing oscillator with nonlinear damping and ex-
citation. Compared to numerical integration, the proposed procedure offers a faster
calculation of existing multiple solutions and their character.
Declarations of interest: none.
Keywords: Nonlinear dynamics, Harmonic Balance method, Duffing oscillator
1. Introduction
The investigation of nonlinear dynamical problems in Mechanics and Mechatron-
ics is of permanent interest. For investigating basic properties of such systems, often
simple nonlinear oscillators with one degree of freedom are considered. There are sev-
eral semi-analytic methods known for solving such problems which are described in
textbooks e.g. by Hagedorn [1]. Among them, Harmonic Balance is an easy to apply
method, which is very popular even in actual publications, sometimes compared with
other methods like normal form transformation or multiple scales e.g. [2]. Harmonic
Balance is combined in the present paper by a corresponding error criterion which is
an extension of the error criterion introduced by Urabe et al. [3]. Additionally a sta-
bility analysis is performed for the low error solutions. Another actual attempt for
∗Corresponding author. Tel: +49 30 31421169; Fax: +49 30 31421173
(L. Lentz)
Preprint submitted to Elsevier February 18, 2020
improvements or additions to the Harmonic Balance method is e.g. the reduction of
the resulting systems of nonlinear algebraic equations by neglecting higher order terms
[4], [5]. In the Incremental Harmonic Balance method, a variable number of harmonic
functions is used to ensure a defined accuracy of the solution [6], [7]. The identification
of non-physical solutions, later called in this paper artifacts, is not in the focus of these
publications.
Among the examples of academic nonlinear oscillators, the Duffing oscillator is
very popular [8]. Therefore as examples a classical Duffing oscillator with solely cubic
nonlinearity and an extended Duffing oscillator are chosen. Compared with the papers
of the same authors [9], [10] and [11] where the same (and other) examples are studied,
the parameters investigated in this paper have been changed also resulting in other
phenomena and qualitative behavior. Furthermore a convergence analysis has been
added.
2. Classical Duffing oscillator
Our studies focus in this section on the softening Duffing oscillator where solutions
being suspicious for artifacts are obtained by Harmonic Balance [9]. As in [9], [10] and
[11] additionally to the Harmonic Balance, an accompanying error criterion is applied
which is an extension of the error criterion introduced by Urabe et al. [3]. As examples
for the application and further development of the criterion by Urabe et al., beside the
already mentioned work of the authors, [12] and [13] shall be mentioned.
The softening Duffing oscillator is used in the following for the introduction of the
method and to show some additional results compared to [9]. This Duffing oscillator is
given by
m¨x+d˙x+cx +αx3=F0cos Ωt(1)
with mbeing the oscillator mass, dthe damping coefficient, cthe linear stiffness, αthe
coefficient of the nonlinear stiffness, F0the excitation force amplitude and Ωthe cir-
cular excitation frequency. In the case of the softening Duffing oscillator αis negative,
while m,d,cand F0are positive. This equation is transformed with respect to dimen-
sionless time derivatives by introducing the circular frequency of the undamped free
linear vibrations ω0=√c/m, the damping ratio D=d/(2 √cm) and the dimensionless
time τ=ω0tas
x00(τ)+2Dx0(τ)+x(τ)+εx3(τ)=fcos(ητ) (2)
with ()0=d()/dτ,ε=α/(mω2
0), f=F0/(mω2
0) and η= Ω/ω0. This equation (2) shall
now be solved approximately by using Harmonic Balance with the ansatz
x(τ)=
n
X
k=1
(akcos(kητ)+bksin(kητ)).(3)
It should be emphasized that, as long as nothing else is stated, for the examples shown
in this paper, only odd kare considered, as we have in general cubic nonlinearities.
This means that ak=bk=0 for even kin this case.
2
Introducing (3) into (2) results in a set of nonlinear algebraic equations
n
X
k=1
˜akcos(kητ)+˜
bksin(kητ)=fcos(ητ)−
3n
X
k=n+1
˜akcos(kητ)+˜
bksin(kητ).(4)
Herein the coefficients ˜ak,˜
bkare nonlinear functions of the ansatz coefficients ak,bk.
In the method of Harmonic Balance, the higher order frequency terms
3n
X
k=n+1
˜akcos(kητ)+˜
bksin(kητ)(5)
in (4) are neglected, i.e. not considered for the calculation of the coefficients ak,bk.
This allows the calculation of the coefficients ak,bkfrom the thereby modified equation
(4) by balancing the respective harmonics.
For visualizing the results x(τ), ˆxis defined as
ˆx=max
0≤τ≤2π
η
n
X
k=1
(akcos(kητ)+bksin(kητ))
.(6)
Now as mentioned before an error criterion is introduced by considering the neglected
terms in the Harmonic Balance method ([3], [9], [10], [11]). The corresponding coef-
ficients ˜ak,˜
bkof the neglected higher order terms can be calculated using the solution
for ak,bkand the error ˆeis defined from this by
ˆe=max
0≤τ≤2π
η
3n
X
k=n+1
˜akcos(kητ)+˜
bksin(kητ)
.(7)
Finally the relative error ˜eis introduced as
˜e=ˆe
ˆx.(8)
This relative error is the main evidence for considering a solution as an artifact. If ˆe
remains larger than 1% even for increasing numbers of ansatz orders n, the solution
is considered to be an artifact. According to our experience, this is accompanied by
changing of the shape of the respective solutions with increasing ansatz orders and
non-convergence of the coeffcients akand bkwhich can be used as additional attributes.
With respect to the demonstration of some results, the parameters
D=0.08, ε =−0.1 and f=0.2 (9)
are taken. The results are displayed by ˆx(6) . Here and in the following a classification
of the solutions is done in that way, that the solutions are divided in solutions with
relative errors ˜elarger 1% and solutions with relative errors ˜elower 1%. The former are
sketched with red color. Solutions with relative error lower than 1% are investigated
via Floquet theory for stability and corresponding results are sketched in grey in the
case of unstable and blue in the case of asymptotically stable solutions.
3
n=1
00.511.52
0
2
4
η
ˆx
n=3
00.511.52
0
2
4
η
ˆx
×
×
n=15
00.511.52
0
2
4
η
ˆx
n=49
00.511.52
0
2
4
η
ˆx
Figure 1: Solutions ˆxaccording to (6) for the Duffing oscillator (2) with Harmonic Balance ansatz (3) in case
of parameters D=0.08, ε =−0.1 and f=0.2. Solutions with relative error ˜eaccording to (8) larger than
1% are sketched in red color. Solutions with relative error ˜elower than 1% are sketched in grey color in the
unstable case (not present in these figures) and blue color in the asymptotically stable case. The coefficients
akand bkof the solutions marked by ”x” and ”+” (downside left) are listed in tables 1-4.
Corresponding results are shown in Figure 1 for increasing ansatz orders n. There
are two types of solutions. First the well-known resonance curve starting for small
ηwith also small ˆx, having the resonance peak close to η=1 and going to zero for
η→ ∞. Due to the softening characteristic, the resonance peak is turning to the left.
This solution preserves its basic shape for all ansatz orders. With the error analysis it
can be seen, that this solution shows very small relative errors ˜ebeing smaller than 1%
for n≥3. In the chosen parameter set, the damping is too large to have a region of
three coexisting solutions close to resonance.
The other solution is the ”nose-like” one occurring for small ηand large ˆx. For
n=1, this solution can be found sketched in many textbooks on nonlinear dynamics.
It can also occur in case of using perturbation analysis [9]. Increasing the ansatz order,
it changes its shape and vanishes more and more, while the errors of the remaining parts
are still larger than 1%. Therefore we consider this solution to be an artifact solution.
The change of solution shape is also characteristic for these artifacts according to our
4
experience.
It shall be mentioned, that the error considered by this method is not an error of
comparing an exact solution with the approximate solution but an error in fulfilling the
equation of motion. As a second criterion, also the convergence of the cosine and sine
coefficients akand bkcan be considered, as it is done also e.g. in [14]. This is shown
in the following.
In Tables 1 and 2 the coefficients akand bkfor the solutions marked in case n=15
in Figure 1 by ”x” for η=0.89 are shown. In this case the used error criterion gives a
small error as result and in fact convergence of the coefficients can be observed.
This changes in the case of the solution marked in case n=15 in Figure 1 by ”+”
for η=0.06. Here a large relative error ˜eis observed and in fact Tables 3 and 4 show a
non-convergence of the coefficients akand bk.
It shall be additionally mentioned that in [9] it was furthermore shown, that an
additional unstable solution (not being an artifact) with non-zero mean value can be
calculated also for small η. To get this solution as a result of Harmonic Balance, the
ansatz (3) has to be extended by a constant term a0. We will see in one of the upcoming
examples, that also coefficients ak,bkwith even khave to be taken into consideration
despite the fact, that there are only cubic nonlinearities in the equation of motion, as it
was also e.g. observed in [15].
3. Extended Duffing oscillator
As next example we consider an extended Duffing oscillator
x00 +2Dx0+x+ε1x3+ε2x2x0=f1cos(ητ)+f2x2cos(ητ) (10)
with the same denominations and the same methods as introduced in section 2. The
cubic nonlinearity with parameter ε1is complemented by a cubic damping term x2x0
with parameter ε2and the harmonic excitation with intensity f1is combined with a
nonlinear parameter excitation term with constant f2. Such an equation can e.g. be ob-
tained from piezoceramic continua by adding conservative and non-conservative terms
in piezoelectric coupling and elasticity as described in [16]. The parameters here are
chosen arbitrarily to show certain nonlinear phenomena and do not necessarily repre-
sent a real piezoceramic.
Let’s first consider the ”good-natured” parameter set D=0.08, ε1=0.0, ε2=
−0.1,f1=0.1 and f2=0.1. Using again the ansatz (3) with only odd k, ˆxcan be
calculated according to equation (6) and the relative error ˜ecorresponding to equation
(8). The results are shown in Figure 2. As ε1=0.0, the system shows neither a
stiffening nor a softening character. As shown in [10], [11] there is an ”island” solution,
that shall be considered in the following. In fact the analysis shows, that the island is
not an artifact, but consists of unstable solutions. The same characteristics for non-
artifact solutions as in the section before can be observed, that the low-error solutions
keep their shape and that there is convergence with respect to the stability behavior. In
fact there is no visible change between the cases n=5 and n=49 for this parameter
set.
5
n=1
00.511.52
0
1
2
3
4
η
ˆx
n=3
00.511.52
0
1
2
3
4
η
ˆx
n=5
00.511.52
0
1
2
3
4
η
ˆx
n=49
00.511.52
0
1
2
3
4
η
ˆx
Figure 2: Solutions ˆxaccording to (6) for the extended Duffing oscillator (10) with Harmonic Balance ansatz
(3) in case of parameters D=0.08, ε1=0.0, ε2=−0.1,f1=0.1 and f2=0.1. Solutions with relative
error ˜eaccording to (8) larger than 1% are sketched in red color. Solutions with relative error ˜elower than
1% are sketched in grey color in the unstable case and blue color in the asymptotically stable case.
Less ”good-natured” is the parameter set D=0.08, ε1=−0.4, ε2=0,f1=
0.1 and f2=0.1. Corresponding results are displayed in Figure 3. For n=1 this
oscillator shows a similar behavior to the classic Duffing oscillator in case of Harmonic
Balance solution with n=1 as shown in Figure 1. In this earlier case the ”nose”
solution is identified as an artifact while here an increasing ansatz order shows another
picture. In fact the right part of the ”nose” keeps its shape and the relative errors
are moving below the 1% border. Most of this low relative error part shows instability
while for increasing orders of nalso a small asymptotically stable part can be identified.
A similar behavior is also possible for the classic Duffing oscillator as shown in [12]
and [10]. As can be seen in Figure 3 the left part of the ”nose” maintains its large
errors and therefore is considered to be an artifact. Overall also for the non-artifact
parts, higher ansatz orders are necessary to get convergence compared to the examples
before.
In all examples so far only odd kwere considered in the ansatz (3). Based on the
6
n=1
00.511.52
0
1
2
η
ˆx
n=3
00.511.52
0
1
2
η
ˆx
n=15
00.511.52
0
1
2
η
ˆx
n=33
00.511.52
0
1
2
η
ˆx
n=41
00.511.52
0
1
2
η
ˆx
n=49
00.511.52
0
1
2
η
ˆx
Figure 3: Solutions ˆxaccording to (6) for the extended Duffing oscillator (10) with Harmonic Balance ansatz
(3) in case of parameters D=0.08, ε1=−0.4, ε2=0,f1=0.1 and f2=0.1. Solutions with relative error
˜eaccording to (8) larger than 1% are sketched in red color. Solutions with relative error ˜elower than 1% are
sketched in grey color in the unstable case and blue color in the asymptotically stable case.
experience from [9], that even in the case of a softening Duffing oscillator with purely
7
n=9
00.511.52
0
1
2
η
ˆx
n=9
00.511.52
0
1
2
η
ˆx
00.511.52
0
1
2
η
ˆx
Figure 4: Same problem and parameter set as in Figure 3. Harmonic Balance solution taking also into account
a nonzero mean value a0and ak,bkwith even k(left). Result from numerical integration for comparison
(right)
cubic nonlinearity, solutions with nonzero mean value occur, the ansatz is extended as
x(τ)=a0+
n
X
k=1
(akcos(kητ)+bksin(kητ))(11)
taking now into consideration a constant a0but also even k. Corresponding results are
shown in Figure 4 with the results for the proposed Harmonic Balance method (left)
for n=9. First, there is an additional almost horizontal branch of unstable solutions.
They are similar to those described in [9], namely forced vibrations with nonzero mean
value around the unstable static equilibrium displacements resulting from the softening
Duffing characteristic. Compared to the earlier results in Figure 3 there are also two
small branches bifurcating from the stable part of the ”nose”. In fact it can be observed,
that these solutions contain symmetry breaking even parts in the Harmonic Balance.
Both the existence of the stable part of the ”nose” and the bifurcating asymmetric
solutions could be confirmed by numerical integration as can be seen in figure 4 (right).
This demonstrates again, that some a priori knowledge is indispensable when applying
the Harmonic Balance method.
Some additional comments shall be given at the end of this chapter. While in
this paper only softening or neutral cases, i.e. ε < 0 in equation(2) and ε1≤0 in
8
equation(10) have been investigated, a stiffening case was considered in [11] in case
of equation(10), with the parameter set ε1=0.1, ε2=0.25,D=0.1,f1=0.1 and
f2=0.1. In that case the resonance peak is - as can be expected - turning to the right,
and an additional non-artifact unstable island solution like in figure 2 occurs. In future
work, the stiffening case could also be interesting with respect to possible artifact or
asymmetric solutions.
Looking for potential problems of the proposed method it shall be remembered, that
in general in Harmonic Balance a corresponding set of nonlinear algebraic equations
for the coefficients ak,bkhas to be solved. Depending on the considered system and
ansatz order this can assume large proportions. It can be seen in the example from
figure 3 with the (not very small but) moderate ε1=−0.4, that the stronger nonlinearity
results in slower convergence of the method. This makes higher ansatz numbers and
solution of larger equation sets necessary and complicates the applicability.
4. Conclusions
The actual paper shows results of using an extended Harmonic Balance method ap-
plied to the classic and an extended Duffing oscillator. The Harmonic Balance method
is hereby extended by an error criterion evaluating the neglected terms and comparing
them with the approximate solution. All solutions calculated are classified in the figures
in three types in that way, that solutions with relative errors larger 1% are considered
and solutions with relative errors lower 1% are divided in asymptotically stable and
unstable solutions. Solutions showing large relative errors even for high ansatz orders
are then considered to be artifacts. According to our experience, non-artifact solutions
are accompanied by converging coefficients in the expansion. On the other hand in
the examples investigated so far, solutions with high relative errors change their shapes
for increasing ansatz orders and show non-converging coefficients. The low-error solu-
tions are examined for their stability using Floquet’s theory. The examples considered
in this paper showed solutions like ”noses” and ”islands” and some of them could be
identified as artifacts. On the other hand for increasing ansatz order, parts of ”noses”
or even complete ”islands” could be detected to show low errors for increasing ansatz
orders. Some parts are even stable what could be confirmed also by numerical inte-
gration in time. Even for purely cubic nonlinearities, there exist asymmetric solutions
with non-vanishing mean value and coefficients ak,bkfor even kas could be seen in
the last example. Compared to numerical integration, the proposed Harmonic Balance
procedure offers to obtain a faster calculation of existing multiple solutions and their
character.
References
[1] P. Hagedorn, Nichtlineare Schwingungen, Akademische Verlagsgesellschaft,
Wiesbaden, 1978.
[2] T. Hill, S. Neild, D. Wagg, Comparing the direct normal form method with har-
monic balance and the method of multiple scales, Procedia Engineering (199)
(2017) 869 – 874.
9
[3] M. Urabe, A. Reiter, Numerical computation of nonlinear forced oscillations by
Galerkin’s procedure, Journal of Mathematical Analysis and Applications (14)
(1966) 107–140.
[4] M. Chowdhury, M. Hosen, K. Ahmad, M. Ali, A. Ismail, High-order approxi-
mate solutions of strongly nonlinear cubic-quintic Duffing oscillator based on the
harmonic balance method, Results in Physics, (14) (2017) 3962 – 3967.
[5] M. Saifur, A. Hasan, Modified harmonic balance method for the solution of non-
linear jerk equations, Results in Physics (8) (2018) 893 – 897.
[6] G. Liu, Z. R. Lv, J. K. Liu, Y. Chen, Quasi-periodic aeroelastic response anal-
ysis of an airfoil with external store by incremental harmonic balance method,
International Journal of Non-Linear Mechanics (100) (2018) 10 – 19.
[7] H. Xiong, X. Kong, H. Li, Z. Yang, Vibration analysis of nonlinear systems with
the bilinear hysteric oscillator by using incremental harmonic balance method,
Communications in Nonlinear Science and Numerical Simulation (42) (2017) 437
– 450.
[8] I. Kovacic, M. Brennan, The Duffing Equation: Nonlinear Oscillators and their
Behaviour, John Wiley & Sons, Ltd, 2011.
[9] U. von Wagner, L. Lentz, On some aspects of the dynamic behavior of the soften-
ing Duffing oscillator under harmonic excitation, Archive of Applied Mechanics
8 (86) (2016) 1383–1390.
[10] U. von Wagner, L. Lentz, On artifacts in nonlinear dynamics, Vibration, Control
and Stability of Dynamical Systems, Proceedings of DSTA (2017) 525–535.
[11] U. von Wagner, L. Lentz, On artifact solutions of semi-analytic methods in non-
linear dynamics, Archive of Applied Mechanics (2018) doi 10.1007/s00419–018–
1397–3.
[12] van Dooren R., On the transition form regular to chaotic behaviour in the Duffing
oscillator, Journal of Sound and Vibration 2 (123) (1988) 327–339.
[13] A. Ferri, M. J. Leamy, Error estimates for harmonic-balance solutions of non-
linear dynamical systems, Proceedings of 50th AIAA/ASME/ASCE/AHS/ASC
Structures, Structural Dynamics, and Materials Conference (2009) AIAA 2009–
2667.
[14] L. Liu, J. Thomas, E. Dowell, P. Attar, K. Hall, A comparison of classical and
high dimensional harmonic balance approaches for a Duffing oscillator, Journal
of Computational Physics 1 (215) (2006) 298 – 320.
[15] A. H. Nayfeh, H. E. Sanchez, Bifurcations in a forced softening Duffing oscillator,
International Journal of Non-Linear Mechanics 24 (6) (1989) 483 – 497.
[16] S. K. Parashar, U. von Wagner, P. Hagedorn, Nonlinear shear-induced flexural vi-
brations of piezoceramic actuators: experiments and modeling, Journal of Sound
and Vibration (285) (2005) 989–1014.
10
Table 1: coefficients akof the solution marked by ”x” in Figure 1,η=0.89
na1a3a5a7a9a11 a13 a15
1 0.69512 0 0 0 0 0 0 0
3 0.69570 0.00373 0 0 0 0 0 0
5 0.69570 0.00373 -0.00001 0 0 0 0 0
7 0.69570 0.00373 -0.00001 0 0 0 0 0
9 0.69570 0.00373 -0.00001 0 0 0 0 0
11 0.69570 0.00373 -0.00001 0 0 0 0 0
13 0.69570 0.00373 -0.00001 0 0 0 0 0
15 0.69570 0.00373 -0.00001 0 0 0 0 0
Table 2: coefficients bkof the solution marked by ”x” in Figure 1,η=0.89
nb1b3b5b7b9b11 b13 b15
1 0.80206 0 0 0 0 0 0 0
3 0.79688 -0.00284 0 0 0 0 0 0
5 0.79687 -0.00284 -0.00002 0 0 0 0 0
7 0.79687 -0.00284 -0.00002 0 0 0 0 0
9 0.79687 -0.00284 -0.00002 0 0 0 0 0
11 0.79687 -0.00284 -0.00002 0 0 0 0 0
13 0.79687 -0.00284 -0.00002 0 0 0 0 0
15 0.79687 -0.00284 -0.00002 0 0 0 0 0
Table 3: coefficients akof the solution marked by ”+” in Figure 1,η=0.06
na1a3a5a7a9a11 a13 a15
1 -3.67914 0 0 0 0 0 0 0
3 -3.78828 0.59545 0 0 0 0 0 0
5 -3.62553 0.19122 0.40248 0 0 0 0 0
7 -3.34890 -0.33538 0.64898 -0.15432 0 0 0 0
9 -2.99132 -0.82318 0.54249 0.24145 -0.23127 0 0 0
11 -2.56463 -1.16502 0.12736 0.44813 0.08710 -0.17347 0 0
13 -2.06646 -1.28601 -0.39000 0.20482 0.31504 0.12752 -0.06258 0
15 -1.46211 -1.11565 -0.71877 -0.30981 -0.00469 0.14619 0.16174 0.10051
Table 4: coefficients bkof the solution marked by ”+” in Figure 1,,η=0.06
nb1b3b5b7b9b11 b13 b15
1 0.67137 0 0 0 0 0 0 0
3 1.24476 -0.87185 0 0 0 0 0 0
5 1.79309 -1.17016 0.41445 0 0 0 0 0
7 2.29957 -1.19090 0.04369 0.34652 0 0 0 0
9 2.75390 -0.95800 -0.42070 0.34839 0.13896 0 0 0
11 3.15069 -0.52464 -0.69577 -0.02656 0.28563 0.09942 0 0
13 3.48854 0.03696 -0.60558 -0.41711 -0.02654 0.17824 0.13803 0
15 3.76891 0.65178 -0.11999 -0.36018 -0.32780 -0.17952 -0.03047 0.05862
11
