Analysis of a cubic-quintic nonlinear oscillator using modified harmonic balance method

Abstract

This study develops an enhanced version of the Modified Harmonic Balance Method (MHBM) for accurately analyzing strongly nonlinear, externally forced, and heavily damped oscillatory systems with cubic–quintic nonlinearities. The proposed formulation addresses the reduced accuracy and convergence limitations of conventional harmonic balance techniques under strong nonlinear conditions. The enhanced MHBM incorporates a systematic power-series expansion of the Fourier coefficients within the harmonic balance framework, improving convergence while maintaining computational efficiency. A representative cubic–quintic nonlinear oscillator is analyzed under different damping ratios, nonlinearity strengths, forcing amplitudes, and initial conditions. The analytical predictions are validated through high-resolution numerical simulations. Results: The proposed method exhibits excellent agreement with numerical solutions over a broad range of system parameters. Compared with the classical MHBM, the enhanced formulation provides improved analytical accuracy, faster convergence, and a more faithful representation of nonlinear dynamic behaviors, including amplitude and phase modulation, particularly in strongly nonlinear regimes. The enhanced MHBM significantly extends the applicability of the classical harmonic balance framework without increasing computational complexity. It provides a reliable, accurate, and computationally efficient semi-analytical approach for investigating strongly nonlinear dynamical systems encountered in engineering and applied physics.

Keywords

nonlinear damped oscillator standard harmonic balance method modified harmonic balance method

1. Introduction

Most physical and technical problems, such as solid mechanics, chemical reactions, electrical circuits, mechanical oscillations, biological systems, and plasma physics, involve nonlinear differential equations (NDEs). The movement of a damped nonlinear mechanical oscillator under the influence of a periodic external force could be represented by the formula. Biological systems where periodic external impacts are important, including population dynamics or biochemical processes, are frequently modeled using nonlinear differential equations. Therefore, it is clear that applied mathematicians, physicists, and engineers must be able to analyze, solve, and understand differential equations. These solutions show the characteristics and actions of the systems. These nonlinear oscillators rarely have suitable solutions. As a result, numerous researchers and scientists have concentrated on creating both analytical and numerical techniques. Numerical approaches use incremental improvement to determine true values at discrete places. Having appropriate initial estimate values is essential for the successful application of numerical approaches.

Although these methods are usually simple, they can require a large amount of computing work and suitable approximations in order to produce the necessary results. Furthermore, nonlinear dynamical systems cannot be fully understood using numerical methods. However, scientists, physicists, engineers, and applied mathematicians have paid close attention to approximation approaches due to their analytical formulations and applicability for parametric analysis. Many methods, such as the averaging approach, the multi-scale method,^1–3 the KBM method,^4,5 and others, can be used to solve nonlinear differential equations with nonlinear stiffness and damping.

Several approximation techniques have been studied to address nonlinear oscillators, including the perturbation approach,^1–7 where perturbation methods are widely used for weakly nonlinear oscillators. Jones improved the breadth and accuracy of the conventional perturbation technique for both small and large parameters,⁸ while Cheung et al. (2010) changed the Lindstedt-Poincaré technique based on Jones concepts.⁹ Also, a modified Lindstedt-Poincaré method for controlling oscillators with significant nonlinearities was shown by Alam et al. (2011).¹⁰ Additional to perturbation method, we have homotopy analysis method,^11,12 homotopy perturbation method,^13–16 variational iteration method,^17,18 harmonic balance method (HBM),^18–23 modified multi-level residue harmonic balance method,^24–26 as well as the modified harmonic balance method (HBM).^27–30

For calculating periodic solutions of nonlinear oscillators, the HBM and MHBM are also useful techniques. These techniques entail choosing a solution that is a shortened Fourier series. In classical HBM, the values of the unknown coefficients are found by numerically solving a set of nonlinear algebraic equations. This approach was refined by several authors.^18–23 For example, Wagner and Lentz and Rahman et al.^19,20 examined the HBM to solve the nonlinear oscillators and the Van der Pol equation. HBM and MHBM were created by a number of writers^21–32 to address nonlinear physical systems.

In recent years there has been substantial growth in applying both analytical techniques and dynamic systems to describe, understand and solve nonlinear and fractional wave equation problems. In this regard, Wang and Wei³³ recently examined the paraxial wave equation as an example of a nonlinear optics problem and used various new optical soliton type solutions to gain greater insight into localized wave behavior in nonlinear optical materials. On top of their work, Wang³⁴ expanded the application of analytical techniques by studying the fractional order version of the Drinfeld-Sokolov system. The study included using conformable fractional derivatives, bifurcation theory, and variational principles to identify a wide range of dynamical properties of the system (i.e. periodic motion, chaotic structure, and bright solitons). In addition, Wang³⁵ further developed the analytical tools he had developed previously by examining the fractionally ordered version of the Kaup-Newell system. The research included developing new solitary wave forms, sensitivity analyses, and bifurcations in relation to each other; it was also found that changes in fractional orders are key factors in determining how waves evolve/stabilize over time. Overall, Wang’s research provides compelling evidence for the development of useful models to examine and analyze complex physical phenomenon related to fractional and nonlinear wave systems.

Research related to non-linear oscillators has made significant advancements in recent years. As a result, researchers have developed an increasing number of highly accurate analytical methods for studying the behaviour of non-linear oscillators. References 36 and 37 describe harmonic-balance based analysis strategies that are able to predict the frequency-amplitude relationship, the resonant response and stability properties of cubic-quintic type oscillators (and similar Duffing-type) with very good accuracy. The research described in reference 38 expanded the application of such analytical methodologies to other types of non-linear vibration problems while still demonstrating a very good correlation with computed solutions. In addition to this work, the authors in reference 39 provided further improvements to both analytical and computational methodology as it relates to improving solution accuracy and providing more reliable computation for more complicated non-linear systems. Additionally, reference 40 describes new models for non-linear behaviours which include significant non-linearity and complexity due to non-linear coupling and various dynamic phenomena. Overall, all of the above referenced contributions demonstrate the continued development of robust analytical methodologies for accurately describing the dynamics of strongly non-linearly behaving oscillatory systems.

We consider the nonlinear damped oscillator with a periodic driving force and cubic-quintic nonlinearities:

\ddot{x} + μ \dot{x} + x + ε (x^{3} + x^{5}) = p \sin (ω t),

(1)

where.

• x is the displacement of the system from its resting point.

• μ > 0 is the damping coefficient.

• ɛ is the nonlinearity strength.

• p is the forcing amplitude.

• ω is the angular frequency of the periodic driving force.

2. Modified harmonic balance formulation

The solution of equation (1) is the Fourier series:

x (t) = a_{0} + \sum_{n = 1}^{N} [a_{n} \cos (n ω t) + b_{n} \sin (n ω t)],

(2)

where a_n and b_n denote the Fourier harmonic amplitudes of the response.

Because the right side of equation (1) satisfies $\sin (ω t + π) = - \sin (ω t)$ , and the nonlinear restoring force x³ + x⁵ is an odd function, the equation is invariant under the symmetry transformation.

S : t \to t + \frac{π}{ω}, x \to - x

Given that the system contains damping (μ > 0), the steady-state periodic solution is unique, and satisfies the half-wave symmetry condition: $x (t) = - x (t + π / ω)$ Substituting equation (2) into this symmetry condition, we get:

a_{0} + \sum_{n = 1}^{N} [a_{n} \cos (n ω t + n π) + b_{n} \sin (n ω t + n π)] = - a_{0} - \sum_{n = 1}^{N} [a_{n} \cos (n ω t) + b_{n} \sin (n ω t)]

Using the trigonometric identities cos (nωt + nπ) = (−1)ⁿ cos (nωt) and sin (nωt + nπ) = (−1)ⁿ sin (nωt), this relation simplifies to:

a_{0} + \sum_{n = 1}^{N} {(- 1)}^{n} [a_{n} \cos (n ω t) + b_{n} \sin (n ω t)] = - a_{0} - \sum_{n = 1}^{N} [a_{n} \cos (n ω t) + b_{n} \sin (n ω t)]

Now we balance the coefficients of the corresponding harmonic terms on both sides. This can imply that a₀ = −a₀ ⇒ a₀ = 0. Furthermore, for any even harmonic (n = 2k), the factor (−1)^2k = 1 leads to an algebraic contradiction unless the coefficients themselves vanish. Therefore, can conclude that:

a_{0} = a_{2} = a_{4} = \dots = b_{2} = b_{4} = \dots = 0

This symmetry breaks down if the equation contains even-order nonlinear terms or if the excitation contains a constant term or cosine. In these cases, even harmonics and a constant term may occur.

Therefore, here we assume a periodic solution containing just odd harmonics, and we consider the solution of equation (1) to be as follows and note that a₁, b₁, a₃, b₃ are to be determined.

x (t) = a_{1} \cos (ω t) + b_{1} \sin (ω t) + a_{3} \cos (3 ω t) + b_{3} \sin (3 ω t)

(3)

First, we find the powers of x in equation (1) which are x³ and x⁵ (Please see Appendix A), then substituting (3) into (1) and using trigonometric identities and collecting the like harmonics, we achieve the following Harmonic Balance Equations:

\begin{align} (1 - ω^{2}) a_{1} + μ ω b_{1} + ε [\frac{3}{4} a_{1}^{3} + \frac{3}{4} a_{1}^{2} a_{3} + \frac{3}{2} a_{1} a_{3}^{2} + \frac{3}{4} a_{1} b_{1}^{2} + \frac{3}{2} a_{1} b_{1} b_{3} + \frac{3}{2} a_{1} b_{3}^{2} \\ - \frac{3}{4} a_{3} b_{1}^{2} + \frac{5}{8} a_{1}^{5} + \frac{5}{16} a_{1}^{4} a_{3} + \frac{15}{8} a_{1}^{3} a_{3}^{2} + \frac{5}{4} a_{1}^{3} b_{1}^{2} + \frac{15}{4} a_{1}^{3} b_{1} b_{3} + \frac{15}{4} a_{1}^{3} b_{3}^{2} \\ + \frac{15}{8} a_{1}^{2} a_{3}^{3} - \frac{15}{8} a_{1}^{2} a_{3} b_{1}^{2} + \frac{15}{8} a_{1}^{2} a_{3} b_{3}^{2} + \frac{15}{8} a_{1} a_{3}^{4} + \frac{15}{4} a_{1} a_{3}^{2} b_{1}^{2} + \frac{15}{4} a_{1} a_{3}^{2} b_{1} b_{3} \\ + \frac{15}{4} a_{1} a_{3}^{2} b_{3}^{2} + \frac{5}{8} a_{1} b_{1}^{4} + \frac{5}{4} a_{1} b_{1}^{3} b_{3} + \frac{15}{4} a_{1} b_{1}^{2} b_{3}^{2} + \frac{15}{4} a_{1} b_{1} b_{3}^{3} + \frac{15}{8} a_{1} b_{3}^{4} \\ - \frac{15}{8} a_{3}^{3} b_{1}^{2} - \frac{15}{16} a_{3} b_{1}^{4} - \frac{15}{8} a_{3} b_{1}^{2} b_{3}^{2}] = 0 \end{align}

(4)

\begin{align} (1 - ω^{2}) b_{1} - μ ω a_{1} + ε [\frac{3}{4} a_{1}^{2} b_{1} + \frac{3}{4} a_{1}^{2} b_{3} - \frac{3}{2} a_{1} a_{3} b_{1} + \frac{3}{2} a_{3}^{2} b_{1} + \frac{3}{4} b_{1}^{3} - \frac{3}{4} b_{1}^{2} b_{3} \\ + \frac{3}{2} b_{1} b_{3}^{2} + \frac{5}{8} a_{1}^{4} b_{1} + \frac{15}{16} a_{1}^{4} b_{3} - \frac{5}{4} a_{1}^{3} a_{3} b_{1} + \frac{15}{4} a_{1}^{2} a_{3}^{2} b_{1} + \frac{15}{8} a_{1}^{2} a_{3}^{2} b_{3} \\ + \frac{5}{4} a_{1}^{2} b_{1}^{3} + \frac{15}{8} a_{1}^{2} b_{1}^{2} b_{3} + \frac{15}{4} a_{1}^{2} b_{1} b_{3}^{2} + \frac{15}{8} a_{1}^{2} b_{3}^{3} - \frac{15}{4} a_{1} a_{3}^{3} b_{1} \\ - \frac{15}{4} a_{1} a_{3} b_{1}^{3} - \frac{15}{4} a_{1} a_{3} b_{1} b_{3}^{2} + \frac{15}{8} a_{3}^{4} b_{1} + \frac{15}{4} a_{3}^{2} b_{1}^{3} - \frac{15}{8} a_{3}^{2} b_{1}^{2} b_{3} \\ + \frac{15}{4} a_{3}^{2} b_{1} b_{3}^{2} + \frac{5}{8} b_{1}^{5} - \frac{25}{16} b_{1}^{4} b_{3} + \frac{15}{4} b_{1}^{3} b_{3}^{2} - \frac{15}{8} b_{1}^{2} b_{3}^{3} + \frac{15}{8} b_{1} b_{3}^{4}] = p \end{align}

(5)

\begin{align} (1 - 9 ω^{2}) a_{3} + 3 μ ω b_{3} + ε [\frac{1}{4} a_{1}^{3} + \frac{3}{2} a_{1}^{2} a_{3} - \frac{3}{4} a_{1} b_{1}^{2} + \frac{3}{4} a_{3}^{3} + \frac{3}{2} a_{3} b_{1}^{2} + \frac{3}{4} a_{3} b_{3}^{2} \\ + \frac{5}{16} a_{1}^{5} + \frac{15}{8} a_{1}^{4} a_{3} + \frac{15}{8} a_{1}^{3} a_{3}^{2} - \frac{5}{8} a_{1}^{3} b_{1}^{2} + \frac{5}{8} a_{1}^{3} b_{3}^{2} + \frac{15}{4} a_{1}^{2} a_{3}^{3} \\ + \frac{15}{4} a_{1}^{2} a_{3} b_{1}^{2} + \frac{15}{4} a_{1}^{2} a_{3} b_{1} b_{3} + \frac{15}{4} a_{1}^{2} a_{3} b_{3}^{2} - \frac{45}{8} a_{1} a_{3}^{2} b_{1}^{2} - \frac{15}{16} a_{1} b_{1}^{4} \\ - \frac{15}{8} a_{1} b_{1}^{2} b_{3}^{2} + \frac{5}{8} a_{3}^{5} + \frac{15}{4} a_{3}^{3} b_{1}^{2} + \frac{5}{4} a_{3}^{3} b_{3}^{2} + \frac{15}{8} a_{3} b_{1}^{4} - \frac{5}{4} a_{3} b_{1}^{3} b_{3} \\ + \frac{15}{4} a_{3} b_{1}^{2} b_{3}^{2} + \frac{5}{8} a_{3} b_{3}^{4}] = 0 \end{align}

(6)

\begin{align} (1 - 9 ω^{2}) b_{3} - 3 μ ω a_{3} + ε [\frac{3}{4} a_{1}^{2} b_{1} + \frac{3}{2} a_{1}^{2} b_{3} + \frac{3}{4} a_{3}^{2} b_{3} - \frac{1}{4} b_{1}^{3} + \frac{3}{2} b_{1}^{2} b_{3} + \frac{3}{4} b_{3}^{3} \\ + \frac{15}{16} a_{1}^{4} b_{1} + \frac{15}{8} a_{1}^{4} b_{3} + \frac{5}{4} a_{1}^{3} a_{3} b_{3} + \frac{15}{8} a_{1}^{2} a_{3}^{2} b_{1} + \frac{15}{4} a_{1}^{2} a_{3}^{2} b_{3} \\ + \frac{5}{8} a_{1}^{2} b_{1}^{3} + \frac{15}{4} a_{1}^{2} b_{1}^{2} b_{3} + \frac{45}{8} a_{1}^{2} b_{1} b_{3}^{2} + \frac{15}{4} a_{1}^{2} b_{3}^{3} - \frac{15}{4} a_{1} a_{3} b_{1}^{2} b_{3} \\ + \frac{5}{8} a_{3}^{4} b_{3} - \frac{5}{8} a_{3}^{2} b_{1}^{3} + \frac{15}{4} a_{3}^{2} b_{1}^{2} b_{3} + \frac{5}{4} a_{3}^{2} b_{3}^{3} - \frac{5}{16} b_{1}^{5} + \frac{15}{8} b_{1}^{4} b_{3} \\ - \frac{15}{8} b_{1}^{3} b_{3}^{2} + \frac{15}{4} b_{1}^{2} b_{3}^{3} + \frac{5}{8} b_{3}^{5}] = 0 \end{align}

(7)

These equations form a complete set for the four unknowns a₁, b₁, a₃, b₃ with parameters ω, p, μ, ɛ.

To eliminate the frequency parameter ω² from (4), (6) and (7), we algebraically substitute (5) in them, then discard the terms with higher-order nonlinear products of the small coefficients (a₃, b₃) under the assumption that the fundamental harmonic components (a₁, b₁) dominate the response.

\begin{align} p a_{1} - μ ω (a_{1}^{2} + b_{1}^{2}) + ε [\frac{3}{2} a_{1} b_{1} (a_{1}^{2} + b_{1}^{2}) b_{3} - \frac{3}{2} b_{1}^{2} (a_{1}^{2} + b_{1}^{2}) a_{3} \\ + \frac{5}{4} a_{1} b_{1} {(a_{1}^{2} + b_{1}^{2})}^{2} b_{3} - \frac{5}{4} b_{1}^{2} {(a_{1}^{2} + b_{1}^{2})}^{2} a_{3}] = 0 \end{align}

(8)

\begin{align} 9 μ ω b_{1} a_{3} - 8 a_{1} a_{3} - 3 μ ω a_{1} b_{3} + ε [a_{1} (\frac{1}{4} a_{1}^{3} - \frac{3}{4} a_{1} b_{1}^{2} + \frac{5}{16} a_{1}^{5} - \frac{5}{8} a_{1}^{3} b_{1}^{2} + \frac{5}{16} a_{1} b_{1}^{4}) \\ - 9 a_{3} (\frac{3}{4} a_{1}^{3} + \frac{3}{4} a_{1} b_{1}^{2} + \frac{5}{8} a_{1}^{5} + \frac{5}{4} a_{1}^{3} b_{1}^{2} + \frac{5}{8} a_{1} b_{1}^{4})] = 0 \end{align}

(9)

\begin{align} 9 μ ω b_{1} b_{3} - 8 a_{1} b_{3} + 3 μ ω a_{1} a_{3} + ε [a_{1} (\frac{3}{4} a_{1}^{2} b_{1} - \frac{1}{4} b_{1}^{3} + \frac{15}{16} a_{1}^{4} b_{1} - \frac{15}{8} a_{1}^{2} b_{1}^{3} + \frac{5}{16} b_{1}^{5}) \\ - 9 b_{3} (\frac{3}{4} a_{1}^{3} + \frac{3}{4} a_{1} b_{1}^{2} + \frac{5}{8} a_{1}^{5} + \frac{5}{4} a_{1}^{3} b_{1}^{2} + \frac{5}{8} a_{1} b_{1}^{4})] = 0 \end{align}

(10)

Then (8) is algebraically substituted into (9) and (10) to delete ω from equation (9), (10). After that, as in the previous step, we pick up the linear component in a₃ and b₃, discarding the terms with higher-order nonlinear products of them.

\begin{align} [9 p a_{1} b_{1} - 8 a_{1} (a_{1}^{2} + b_{1}^{2}) - 9 ε (a_{1}^{2} + b_{1}^{2}) (\frac{3}{4} a_{1}^{3} + \frac{3}{4} a_{1} b_{1}^{2} + \frac{5}{8} a_{1}^{5} + \frac{5}{4} a_{1}^{3} b_{1}^{2} + \frac{5}{8} a_{1} b_{1}^{4})] a_{3} \\ - 3 p a_{1}^{2} b_{3} \\ + ε (a_{1}^{2} + b_{1}^{2}) a_{1} (\frac{1}{4} a_{1}^{3} - \frac{3}{4} a_{1} b_{1}^{2} + \frac{5}{16} a_{1}^{5} - \frac{5}{8} a_{1}^{3} b_{1}^{2} + \frac{5}{16} a_{1} b_{1}^{4}) \\ + ε [\frac{27}{2} a_{1}^{2} b_{1}^{2} (a_{1}^{2} + b_{1}^{2}) + \frac{45}{4} a_{1}^{2} b_{1}^{2} {(a_{1}^{2} + b_{1}^{2})}^{2}] b_{3} = 0 \end{align}

(11)

\begin{align} [9 p a_{1} b_{1} - 8 a_{1} (a_{1}^{2} + b_{1}^{2}) - 9 ε (a_{1}^{2} + b_{1}^{2}) (\frac{3}{4} a_{1}^{3} + \frac{3}{4} a_{1} b_{1}^{2} + \frac{5}{8} a_{1}^{5} + \frac{5}{4} a_{1}^{3} b_{1}^{2} + \frac{5}{8} a_{1} b_{1}^{4})] b_{3} \\ + 3 p a_{1}^{2} a_{3} \\ + ε (a_{1}^{2} + b_{1}^{2}) a_{1} (\frac{3}{4} a_{1}^{2} b_{1} - \frac{1}{4} b_{1}^{3} + \frac{15}{16} a_{1}^{4} b_{1} - \frac{15}{8} a_{1}^{2} b_{1}^{3} + \frac{5}{16} b_{1}^{5}) \\ + ε [\frac{27}{2} a_{1}^{2} b_{1}^{2} (a_{1}^{2} + b_{1}^{2}) + \frac{45}{4} a_{1}^{2} b_{1}^{2} {(a_{1}^{2} + b_{1}^{2})}^{2}] a_{3} = 0 \end{align}

(12)

Next we Simplify equations (11) and (12), and solve for a₃ and b₃. Then we have:

\begin{align} a_{3} = - \frac{ε (a_{1}^{2} + b_{1}^{2}) a_{1} (\frac{1}{4} a_{1}^{3} - \frac{3}{4} a_{1} b_{1}^{2} + \frac{5}{16} a_{1}^{5} - \frac{5}{8} a_{1}^{3} b_{1}^{2} + \frac{5}{16} a_{1} b_{1}^{4})}{9 p a_{1} b_{1} - 8 a_{1} (a_{1}^{2} + b_{1}^{2}) - 9 ε (a_{1}^{2} + b_{1}^{2}) (\frac{3}{4} a_{1}^{3} + \frac{3}{4} a_{1} b_{1}^{2} + \frac{5}{8} a_{1}^{5} + \frac{5}{4} a_{1}^{3} b_{1}^{2} + \frac{5}{8} a_{1} b_{1}^{4})} \end{align}

(13)

\begin{align} b_{3} = - \frac{ε a_{1} b_{1} (a_{1}^{2} + b_{1}^{2}) (\frac{3}{4} a_{1}^{2} - \frac{1}{4} b_{1}^{2} + \frac{15}{16} a_{1}^{4} - \frac{15}{8} a_{1}^{2} b_{1}^{2} + \frac{5}{16} b_{1}^{4})}{9 p a_{1} b_{1} - 8 a_{1} (a_{1}^{2} + b_{1}^{2}) - 9 ε (a_{1}^{2} + b_{1}^{2}) (\frac{3}{4} a_{1}^{3} + \frac{3}{4} a_{1} b_{1}^{2} + \frac{5}{8} a_{1}^{5} + \frac{5}{4} a_{1}^{3} b_{1}^{2} + \frac{5}{8} a_{1} b_{1}^{4})} \end{align}

(14)

After that we substitute a₃ and b₃ into equation (8) and solving for b₁, we obtain b₁ as a power series:

b_{1} = a_{1}^{2} λ + a_{1}^{4} λ^{3} + 2 a_{1}^{6} λ^{5} + 5 a_{1}^{8} λ^{7} + \dots,

(15)

where

λ = \frac{μ ω}{p} .

Finally, inserting a₃, b₃, and b₁ into (5), we obtain the equation for the fundamental amplitude a₁, and this equation can be solved numerically for a₁ given the system parameters μ, ɛ, p, and ω.

3. Results and discussion

The graphical results presented in Figures 1–8 provide a comprehensive comparison between the proposed Modified Harmonic Balance Method (MHBM), the standard Harmonic Balance Method (HBM), and the numerical Runge–Kutta (RK4) solutions. First, it is clearly observed that the proposed MHBM exhibits excellent agreement with the RK4 numerical results for all considered cases. In contrast, the standard HBM shows noticeable deviations, particularly in the amplitude and phase of oscillations. This confirms the improved accuracy of the proposed method. In Figure 1 (mild nonlinearity and moderate damping), all methods produce similar responses; however, small discrepancies between HBM and RK4 begin to appear. The MHBM generally follows the numerical solution more closely than the standard HBM, although isolated points may exhibit comparable or slightly smaller errors for the standard HBM.

As the system parameters are changed, the difference in behavior is even clearer. In Figure 2 (with low damping and high nonlinearities), the amplitude of oscillations increases as do the phase differences. While the standard harmonic balance method does not properly model this change in behavior, the modified harmonic balance method continues to be very consistent with the results obtained by using the Runge-Kutta fourth order numerical solution.

In the high-damping case (Figure 3), the system shows less oscillation amplitude and stabilizes more quickly. In this case, both analytical methods work pretty well, but the MHBM is still a better approximation of the numerical solution, especially when it comes to the transient response.

In the case of strong nonlinearity (Figure 4), the standard HBM’s limits become very clear. You can clearly see the differences in both amplitude and phase. The MHBM, on the other hand, keeps giving results that are in line with the RK4 solution, showing that it works well for systems that are very nonlinear.

In Figures 5–8 MHBM’s trajectory is essentially identical to RK4; this suggests that MHBM effectively models the behavior of the system. On the other hand, there were clearly observable deviations from the expected trajectory for HBMs.

In addition to observing the differences graphically, it was also possible to quantify them through use of the data shown in Tables 1–4.

Table 1.

Numerical comparison with high precision for case 1.

t	RK4	MHBM	Std_HBM	Abs_Error_RK_MHBM	Abs_Error_RK_HBM
0.0	-0.23	-0.22601674	-0.23510233	0.00398326	0.00510233
1.0	1.13173384	1.14701044	1.16591091	0.0152766	0.03417706
2.0	1.56436213	1.57607459	1.5887325	0.01171246	0.02437037
3.0	0.42932544	0.42397948	0.43827619	0.00534596	0.00895075
4.0	-0.95648133	-0.96691084	-0.98591534	0.01042951	0.02943401
5.0	-1.62856598	-1.63570092	-1.64689784	0.00713494	0.01833186
6.0	-0.62103861	-0.62123839	-0.63933686	0.00019978	0.01829825
7.0	0.77299156	0.77605816	0.79356825	0.0030666	0.02057669
8.0	1.64757081	1.65334365	1.66414401	0.00577284	0.01657319
9.0	0.80852184	0.81576861	0.83597338	0.00724677	0.02745155
10.0	-0.58010239	-0.57964415	-0.59395817	0.00045825	0.01385577

Table 2.

Numerical comparison with high precision for case 2.

t	RK4	MHBM	Std_HBM	Abs_Error_RK_MHBM	Abs_Error_RK_HBM
0.0	-0.03	-0.03171627	-0.03523229	0.00171627	0.00523229
1.0	0.73406229	0.76755767	0.80528785	0.03349538	0.07122555
2.0	1.34545603	1.35004762	1.35402047	0.00459158	0.00856443
3.0	0.71914675	0.77567495	0.81172276	0.0565282	0.09257601
4.0	-0.00507874	-0.01962554	-0.02189467	0.0145468	0.01681593
5.0	-0.76652517	-0.83316618	-0.86905277	0.06664101	0.1025276
6.0	-1.35051401	-1.3464485	-1.35055925	0.00406551	4.524e-05
7.0	-0.67113504	-0.71032667	-0.74782867	0.03919162	0.07669363
8.0	0.07641271	0.07118643	0.07920144	0.00522628	0.00278873
9.0	0.86181776	0.8978795	0.93138621	0.03606174	0.06956845
10.0	1.31875519	1.33510102	1.33997761	0.01634583	0.02122242

Table 3.

Numerical comparison with high precision for case 3.

t	RK4	MHBM	Std_HBM	Abs_Error_RK_MHBM	Abs_Error_RK_HBM
0.0	-0.37	-0.36146422	-0.3692883	0.00853578	0.0007117
1.0	0.8312529	0.82989986	0.84633026	0.00135303	0.01507737
2.0	1.00668671	1.02855908	1.04132431	0.02187237	0.0346376
3.0	-0.12251059	-0.09232671	-0.10214836	0.03018389	0.02036223
4.0	-1.14477833	-1.14706285	-1.15524055	0.00228451	0.01046222
5.0	-0.61031386	-0.63859474	-0.65305107	0.02828088	0.04273721
6.0	0.56067012	0.55583623	0.57474545	0.00483389	0.01407533
7.0	1.15426071	1.17138127	1.17985553	0.01712056	0.02559483
8.0	0.19022482	0.18841752	0.18974843	0.0018073	0.0004764
9.0	-0.95395756	-0.97707815	-0.99018827	0.02312059	0.03623071
10.0	-0.88801362	-0.89543959	-0.91032592	0.00742596	0.0223123

Table 4.

Numerical comparison with high precision for case 4.

t	RK4	MHBM	Std_HBM	Abs_Error_RK_MHBM	Abs_Error_RK_HBM
0.0	-0.05	-0.05512255	-0.06663321	0.00512255	0.01663321
1.0	0.65767201	0.70469083	0.76080267	0.04701882	0.10313066
2.0	1.16226067	1.22065892	1.21821364	0.05839824	0.05595297
3.0	0.36671118	0.3832154	0.44085244	0.01650422	0.07414126
4.0	-0.14348574	-0.25683516	-0.31019025	0.11334942	0.16670451
5.0	-1.09674971	-1.15417118	-1.16447061	0.05742147	0.0677209
6.0	-0.75273186	-0.84880207	-0.8863268	0.09607021	0.13359495
7.0	-0.00822531	-0.04475731	-0.05360757	0.03653201	0.04538226
8.0	0.68770008	0.72370318	0.77857322	0.0360031	0.09087314
9.0	1.1433846	1.21374631	1.21205279	0.07036171	0.06866819
10.0	0.3760737	0.36783295	0.42491757	0.00824074	0.04884387

These tables showed that MHBM was able to produce error measurements, which were lower than those produced by HBMs. Finally, as seen in Table 5, the percentage relative improvements demonstrated in the performance of MHBM over HBMs were 67.5%, 40.4%, 34.1%, and 37.5% in all four cases, respectively.

Table 5.

Summary comparison of all four cases (cubic-quintic oscillator).

Case	Description	μ	ϵ	ω	p	Mean err (MHBM)	Mean err (HBM)	Relative improv.
1	Mild case	0.10	0.10	1.00	1.00	0.006421	0.019738	67.5%
2	Low damping	0.05	0.30	0.80	1.50	0.025310	0.042478	40.4%
3	High damping	0.20	0.50	1.20	0.80	0.013347	0.020243	34.1%
4	High nonlinearity	0.15	0.80	0.90	2.00	0.049548	0.079241	37.5%

Note. All values for Improvement % are relative to the Standard HBM error.

The overall results show that system parameters like damping, nonlinearity strength, and forcing amplitude have a big effect on how accurate the proposed method is. The MHBM remains accurate in different situations, which makes it a good analytical tool for solving cubic–quintic nonlinear oscillators. Finally, the present results compare well with previously reported harmonic balance formulations for nonlinear oscillators. The proposed MHBM differs from the conventional HBM in that it takes into account the third-harmonic contributions and a systematic power-series treatment of the Fourier coefficients. The resulting agreement with the RK4 numerical solutions is much improved. This method is especially useful in regimes of strong cubic-quintic nonlinearities where classical harmonic balance approaches often exhibit significant amplitude and phase errors. From the application point of view, cubic–quintic nonlinear oscillators are often encountered in vibration isolation systems, nonlinear mechanical structures, MEMS/NEMS devices, nonlinear electrical circuits, and energy harvesting systems. Therefore, the proposed MHBM can be used as a practical analytical scheme for efficient prediction and design of strongly nonlinear dynamical systems.

4. Conclusion

The MHBM has been successfully developed to determine the steady-state periodic response of a forced damped cubic–quintic nonlinear oscillator, where the obtained time histories, phase portraits and the comparisons with numerical simulation indicate that the proposed MHBM gives an accurate and computationally efficient framework for the analysis of the steady-state periodic response of cubic–quintic nonlinear oscillators. A novel approximation technique has been introduced by retaining third harmonic components and expressing some of the Fourier series expansion functions as power series. Thus, the MHBM produces an approximate analytical representation of the system’s behavior which can capture the essential nonlinear aspects of the motion much better than the conventional harmonic balance method. The large amplitude motions associated with the cubic-quintic restoring forces lead to a strong harmonic interaction between the primary and higher order harmonics. Due to its inability to include such higher-order harmonic interactions, the standard harmonic balance method cannot capture the nonlinear energy transfer between modes. As a result, the method will contain errors in both magnitude and phase. However, when the third-order harmonic terms are included in the MHBM, they can effectively describe the nonlinear interaction between the modes. Such inclusion becomes more important as the strength of the non-linearity increases. Additionally, since the interaction between different harmonic components leads to a modification of the amplitude of the fundamental component and therefore modifies the natural frequency of the fundamental mode (nonlinear frequency shift), such modification is inherently captured by the MHBM. Finally, consistent with prior formulations, damping is treated directly within the MHBM formulation; thus, energy dissipated over one cycle is modeled and accurate predictions of amplitude decays and phase lags at each point in time are made during steady-state conditions. Therefore, although computationally efficient, the MHBM delivers nearly-numerically precise results and is capable of accounting for harmonic coupling, frequency shifts due to non-linear effects, dissipative effects and stable periodic motions. Hence, it appears that the MHBM can serve as a reliable analytical tool to study nonlinear oscillators with simultaneous cubic and quintic nonlinearities.

Finally, the present framework can be extended for further investigations of the global nonlinear dynamics such as the bifurcation behavior, Poincaré maps and the Largest Lyapunov Exponents in order to characterize the transition from the periodic to the chaotic motion in the forced cubic–quintic nonlinear oscillators.

Footnotes

Acknowledgements

The authors wish to express their gratitude to Palestine Technical University-Kadoorie for providing the financial support and resources necessary to conduct this research.

ORCID iD

Jihad Asad

Funding

The authors received no financial support for the research, authorship, and/or publication of this article.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Appendix

References

Nayfeh

. Introduction to perturbation techniques. John Wylie, 1993.

Nayfeh

. Perturbation method. John Wiley & Sons, 1983.

Nayfeh

Mook

. Nonlinear oscillations. John Wiley & Sons, 1989.

Krylov

Bogolyubov

. Introduction to nonlinear mechanics. Princeton University, 1948.

Bogoliubov

Mitropolskii

. Asymptotic methods in the theory of nonlinear oscillations. Gordan and Breach, 1961.

Alam

. A unified Krylov-Bogoliubov-Mitropolskiimethod for solving nthorder nonlinear systems. J Franklin Inst 2002; 399(2): 239–248. https://doi.org/10.1016/s0016-0032(02)00020-0

Alam

Azad

MAK

Hoque

. A general Struble’s technique for solving an nth-order weakly non-linear differential system with damping. Int J Non Lin Mech 2006; 41: 905–918.

Jones

. Remarks on the perturbation process for certain conservative systems. Int J Non Lin Mech 1988; 13: 125–128. https://doi.org/10.1016/0020-7462(78)90021-5

Cheung

Chen

Lau

. A modified Lindstedt-Poincare method for certain strongly nonlinear oscillators. Int J Non Lin Mech 1991; 26: 367–378. https://doi.org/10.1016/0020-7462(91)90066-3

10.

Alam

Yeasmin

Ahamed

. Generalization of the modified Lindstedt-Poincare method for solving some strong nonlinear oscillators. Ain Shams Eng J 2019; 10(1): 195–201. https://doi.org/10.1016/j.asej.2018.08.007

11.

Fooladi

Abaspour

Kimiaeifar

, et al. On the analytical solution of Kirchhoff simplified Model for beam by using of homotopy analysis method. World Appl Sci J 2009; 6: 298–302.

12.

Liao

. The proposed homotopy analysis technique for the solution of nonlinear problems. Ph.D. Thesis, Shanghai JiaoTong University, 1992.

13.

. Homotopy perturbation method for nonlinear oscillators with coordinate dependent mass. Results Phys 2018; 10: 270–271. https://doi.org/10.1016/j.rinp.2018.06.015

14.

Uddin

Sattar

Alam

. An approximate technique for solving strongly nonlinear differential systems with damping effects. Indian J Math 2011; 53(1): 83–98.

15.

Uddin

Alom

Ullah

. An analytical approximate technique for solving a certain type of fourth order strongly nonlinear oscillatory differential system with small damping. Far East J Math Sci 2012; 68(1): 59–82.

16.

Wazwaz

. The variational iteration method: a reliable analytic tool for solving linear and nonlinear wave equations. Comput Math Appl 2008; 54: 926–932. https://doi.org/10.1016/j.camwa.2006.12.038

17.

. Some asymptotic methods for strongly nonlinear equations. Int J Mod Phys B 2006; 20: 1141–1199. https://doi.org/10.1142/s0217979206033796

18.

Mickens

. A generalization of the method of harmonic balance. J Sound Vib 1986; 111: 515–518. https://doi.org/10.1016/s0022-460x(86)81410-9

19.

Rahman

Haque

Shanta

. Harmonic balance solution of nonlinear differential equation (non-conservative). Journal of Advances in Vibration Engineering 2010; 9(4): 334–346.

20.

Wagner

Lentz

. On the detection of artifacts in harmonic balance solutions of nonlinear oscillators. Applied Mathematical Modeling 2019; 65: 408–414.

21.

. The harmonic balance method for Yao-Cheng oscillator. J Low Freq Noise Vib Act Control 2019; 38: 1816–1818. https://doi.org/10.1177/1461348419847307

22.

Yeasmin

Sharif

Rahman

, et al. Analytical technique for solving the quadratic nonlinear oscillator. Results Phys 2020; 18: 103303. https://doi.org/10.1016/j.rinp.2020.103303

23.

Rahman

Lee

. New modified multi-level residue harmonic balance method for solving nonlinearly vibrating double-beam problem. J Sound Vib 2018; 406: 295–327. https://doi.org/10.1016/j.jsv.2017.06.017

24.

Rahman

Hasan

ASMZ

Yeasmin

. Modified multi-level residue harmonic balance method for solving nonlinear vibration problem of beam resting on nonlinear elastic foundation. Journal of Applied and Computational Mechanics 2019; 5(4): 628–638.

25.

Hasan

ASMZ

Rahman

Lee

, et al. Multi-level residue harmonic balance solution for the nonlinear natural frequency of axially loaded beams with an internal hinge. Mech Adv Mater Struct 2018; 24(13): 1074–1085. https://doi.org/10.1080/15376494.2016.1227492

26.

Lee

. Analytic solution for nonlinear multimode beam vibration using a modified harmonic balance approach and Vieta’s substitution. Shock Vib 2016; 2016: 6.

27.

Ullah

Uddin

Rahman

. A modified harmonic balance method for solving strongly generalized nonlinear damped forced vibration systems. Nonlinear Dynam Syst Theor 2021; 21(5): 544–552.

28.

Ullah

Uddin

Rahman

. An analytical technique for handling forced Van der Pol vibration equation. J Bangladesh Acad Sci 2021; 45(2): 231–240. https://doi.org/10.3329/jbas.v45i2.56654

29.

Ullah

Rahman

Uddin

. Free vibration analysis of nonlinear axially loaded beams using a modified harmonic balance method. Partial Differential Equations in Applied Mathematics 2022; 6: 100414. https://doi.org/10.1016/j.padiff.2022.100414

30.

Ullah

Uddin

Rahman

. Analytical solution of modified forced Van der Pol vibration equation using modified harmonic balance method. Khulna University Studies. Khulna Univ Stud 2022: 892–903. https://doi.org/10.53808/KUS.2022.ICSTEM4IR.0180-se

31.

Ullah

Rahman

Uddin

. A modified harmonic balance method for solving forced vibration problems with strong nonlinearity. J Low Freq Noise Vib Act Control 2021; 40(2): 1096–1104. https://doi.org/10.1177/1461348420923433

32.

Kandil

Hamed

Awrejcewicz

. Harmonic balance method to analyze the steady-state response of a controlled mass-damper-spring model. Symmetry 2022; 14: 1247. https://doi.org/10.3390/sym1406124

33.

Wang

Wei

. Novel optical soliton solutions to nonlinear paraxial wave model. Modern Physics Letters B 2025; 39(12): 2450469. https://doi.org/10.1142/S0217984924504694

34.

Wang

. Bifurcation, chaotic, variational principle and new soliton solutions for the nonlinear fractional Drinfeld–Sokolov system. Fractals 2026; 34: 2650062. https://doi.org/10.1142/s0218348x26500623

35.

Wang

. Investigation of soliton structures and dynamical behaviors for fractional Kaup–Newell system. Fractals 2026; 2026: 2650081. https://doi.org/10.1142/s0218348x26500817

36.

Amer

Starosta

Elameer

, et al. Analyzing the stability for the motion of an unstretched double pendulum near resonance. Applied Sciences 2021; 11(20): 9520. https://doi.org/10.3390/app11209520

37.

Amer

Bek

Hassan

. The dynamical analysis for the motion of a harmonically two degrees of freedom damped spring pendulum in an elliptic trajectory. Alexandria Engineering Journal 2022; 61(2): 1715–1733. https://doi.org/10.1016/j.aej.2021.06.063

38.

Amer

Starosta

, et al. Resonance in the cart-pendulum system—an asymptotic approach. Applied Sciences 2021; 11(23): 11567. https://doi.org/10.3390/app112311567

39.

Moatimid

Amer

Galal

. Inspection of some extremely nonlinear oscillators using an inventive approach. Journal of Vibration Engineering Technologies 2024; 12(Suppl 2): 1211–1221. https://doi.org/10.1007/s42417-024-01469-y

40.

Moatimid

Amer

Mohamed

. Chaotic and dynamic vibration analysis of a time-delayed nonlinear mathieu oscillator via non-perturbative approach. Scientific Reports 2026; 16(1): 12219. https://doi.org/10.1038/s41598-026-45062-7