Analyzing the chaotic and stability behavior of a duffing oscillator excited by a sinusoidal external force

Abstract

This paper delves into the dynamical, chaotic, and stability aspects of the Duffing oscillator (DO) under sinusoidal external excitation, underscoring its relevance in various scientific and engineering applications. The DO, with its complex behavior, poses a significant challenge to our understanding. A unique perturbation technique, the multiple-scales (MS), is harnessed to tackle this challenge and enrich our knowledge. The analysis entails determining third-order expansions, exploring resonant cases, and factoring in the influence of viscous damping. The accuracy of the analytical solution is cross-validated with numerical results using the Runge–Kutta fourth-order (RK4). The paper also employs effective methods to assess the obtained results qualitatively. As a result, Visual figures, such as bifurcation diagrams and Lyapunov exponents’ spectra (LES), have been presented to illustrate diverse system motions and demonstrate Poincaré diagrams. Moreover, stability analysis has been investigated via resonance curves showing the stable and unstable regions. These aids facilitate our comprehension of the system’s complex behavior and its variations under different conditions. The results are elucidated through displayed curves, offering insights into the dynamics of the DO under sinusoidal excitation and damping effects. The influence of excitation and the quadratic parameter on bifurcation diagrams, LES, and Poincaré maps helps us understand the intricate behavior of nonlinear oscillatory systems, such as DO. The presented oscillator is a compelling example of elucidating the nonlinear behavior observed in numerous engineering and physics phenomena.

Introduction

The DO is a fundamental model in the study of nonlinear dynamics, serving as a key element for understanding and predicting the behavior of complex, oscillatory systems. ¹ Furthermore, the DO’s behavior is characterized by its richness and complexity, encompassing periodic, quasi-periodic, and chaotic motion, which is contingent upon the specific values assigned to its parameters. This makes it a crucial model for understanding and predicting the behavior of various real-world systems, such as mechanical systems with nonlinear springs, electrical circuits with nonlinear elements, and even the dynamics of the human brain. The equation featuring a cubic nonlinearity characterizing an oscillator is called the Duffing equation, which describes an oscillator for any equation with a cubic stiffness term. It can be applied in several fields, like mathematics, engineering, physics, and biology. ² For example, studying large amplitude oscillations, and centrifugal speed control systems, beam nonlinear vibration plates.^3,4 The model was introduced by the German electrical engineer Duffing in 1918. Since then, alongside Van der Pol’s equation, it has emerged as one of the most frequently encountered equations for nonlinear oscillations in both research articles and textbooks. In References 5 and 6.^5,6 the authors illustrated a variation of the classic Duffing equation featuring a negative linear stiffness term. These equations are employed to model the dynamics of a buckling beam or plate, particularly when considering only a single type of vibration. A systematic methodology for designing high-performance oscillator systems was adopted in Lowenstein (2012) and Cveticanin (2014).^7,8 This approach is grounded in a main classification of oscillators, categorized according to their internal timing references. The Duffing equation is a useful model for studying nonlinear processes oscillation and chaotic dynamical systems. ⁹ Chaos is inherently nonlinear, and it exists in many natural and artificial systems that are sensitive to initial conditions. ¹⁰ There is a widespread consensus that mechanical systems, particularly vibration dynamics issues as indicated in Amer et al. (2019)^11–15 hold paramount significance in mechanics, primarily owing to their extensive applications, notably in engineering fields.

The DO has been extensively studied in Lo and Chen (1998) ¹⁶ to provide in-depth analysis and insights into its behavior. Chaos theory stands out as one of the foremost accomplishments within the realm of nonlinear science. ¹⁷ The oscillator is frequently linked with chaotic behavior within mathematical contexts, where such chaotic phenomena are observed in various natural systems, including climate and weather patterns. ¹⁸ Chaos theory finds applications across numerous disciplines, spanning meteorology, sociology, environmental science, and more. As outlined in Strogatz, ¹⁹ chaos manifests as a periodic long-term behavior within deterministic systems that display sensitivity to initial conditions.

Furthermore, the DO is considered one of the prototypes of nonlinear dynamic systems. ²⁰ In mechanics, the basic form of the Duffing equation can be viewed as a one degree-of-freedom for the motion of a mathematical model, where the damping term may be linear or nonlinear in addition to a nonlinear stiffness. For instance, this equation is applicable to a suspended mass system, where the mass is supported by a parallel arrangement of shock absorbers with constant damping and springs exhibiting nonlinear restoring forces subjected to harmonic excitation. ²¹ It can also be used to describe one or more degrees of freedom of a complex system. Moreover, it can be used to govern the general form of the ship roll equation. ²² In Song et al. ²³ using the fast-slow analysis method, the authors aimed to explore the dynamic behaviors of the entire excitation. A three-dimensional visualization technique is suggested to examine the relaxation oscillations triggered by various modulation modes, offering a comprehensive understanding of the system’s dynamics from different angles. In Barba-Franco et al.(2023), ²⁴ the authors studied the behavior of multiple stable, coexisting rotating waves on a one-way ring consisting of interconnected double-well DO with different numbers of oscillators. Using the analysis of time series, bifurcation diagrams, phase diagrams, and basins of attraction, they demonstrated the existence of multi-stability as the system transitions from a stable equilibrium to hyper-chaos as the coupling strength increases. This transition is characterized by a series of bifurcations, including crisis bifurcations, torus bifurcations, and Hopf bifurcations. In Coccola et al.(2023), ²⁵ the nonlinear DO was examined with fractional damping, a characteristic feature in various physical scenarios. The analysis involved investigating the system under conditions of both smaller (underdamped) and larger (overdamped) damping parameter values. A new method for achieving high-precision periodic movements using an implicit mapping technique, which is validated by the conditions of saddle-node bifurcation domains, was presented in Liu et al (2023). ²⁶ Through an analysis of the cycle development curve movements in a DO with hardening properties, it is discovered that the primary resonance can be utilized to improve the characteristic components in the generated signal.

This paper presents an analytical solution to the governing equation of a DO subjected to sinusoidal external excitation employing the method of MS ²⁷ up to the third approximation. This solution is validated by comparing it with the obtained numerical ones using RK4. The study focuses on the external resonance case, followed by establishing modulation equations (ME) based on the analyzed resonance scenarios. System stability is assessed, and multiple areas near resonance ranges are plotted to examine dynamic frequency response curves across different values of the influencing parameters. The obtained solutions are graphed for consistent parameter values. Furthermore, the behavior is depicted through variations in excitation amplitudes, damping coefficients, and certain parameters. Stability criteria are employed to analyze stable and unstable regions within a system and assess the influence of various parameters. Illustrations such as bifurcation diagrams, LES, and Poincaré diagrams demonstrate several types of system motions. In the design and analysis of mechanical systems exposed to vibration, such as machinery or car suspensions, Duffing oscillators are commonly used. Systems modeled by these oscillators are utilized in energy harvesting to capture ambient mechanical vibrations and convert them into electrical energy. Moreover, they are used to model the nonlinear dynamics of base isolation systems in buildings, bridges, and other structures designed to mitigate the effects of earthquakes. Seismic isolators often have nonlinear stiffness, and the DO helps in predicting how the system will respond to earthquake excitations. The stability behavior ensures the structure’s safety under seismic loads, while chaotic behavior must be avoided to prevent structural failure.

In these applications, grasping the chaotic and stability dynamics of Duffing-like systems is crucial for converting ambient mechanical vibrations into electrical energy through vibration energy harvesting. This understanding assists in controlling chaos for improved system efficiency or avoiding it to preserve stability and protect mechanical, electrical, or biological systems.

Problem’s description and the perturbation procedure

In this section, the following second-order ordinary differential equation (ODE) presents the DO model with excitations, which has several applications in various fields, as mentioned before. This equation illustrates the response of a dimensionless system with quadratic nonlinearities to a sinusoidal excitation ²⁷ modeled as

$\ddot{\mathit{u}}+{\mathit{\omega }}^{\mathbf{2}}\mathit{u}+\mathbf{2}\mathit{\mu }\dot{\mathit{u}}+\mathit{\alpha }{\mathit{u}}^{\mathbf{2}}+\mathit{\eta }{\mathit{u}}^{\mathbf{3}}=\mathit{F}\left(\mathit{t}\right).$

(1)

Here $F\left(t\right)=F^{*}\hspace{0.17em}\mathit{cos}\hspace{0.17em}\Omega t$ represents an external harmonic force of the sinusoidal excitation, $F^{*}$ and $\Omega$ denote their amplitude and angular frequency, respectively. Moreover, $\mu$ is the damping coefficient, $\alpha$ is the coefficient of the quadratic term, $\eta$ is the nonlinear stiffness coefficients, $\omega$ is the internal frequency, and the dot is the derivative regarding time $\tau$ .

It must be mentioned that system (1) can exhibit a range of responses depending on the parameters, including periodic, quasi-periodic, and chaotic motion. As parameters such as the amplitude and frequency of the excitation are varied, the DO can undergo bifurcations, leading to chaotic behavior. This makes it a rich system for studying complex dynamical phenomena. Due to its nonlinear nature, this oscillator can exhibit hysteresis and jump phenomena in its response, particularly when the excitation frequency is varied. The system may show multiple resonance peaks due to the presence of both linear and nonlinear restoring forces, which can lead to complex resonance behaviors.^28,29

We will now proceed to solve the equation as mentioned earlier (1) by applying a perturbation method known by MS ²⁷ to derive a solution up to the third approximation. Thus, we presume that the solution with respect to the small parameter $\epsilon$ is as follows

$u\left(\tau \right)=\epsilon y\left(\tau ;\epsilon \right).$

(2)

In accordance with the MS method, one may construct the function $y$ in powers of $\epsilon$ as follows

$y={\displaystyle \sum }_{k=0}^2{\epsilon }^k{y}_{k+1}\left({\tau }_0,{\tau }_1,{\tau }_2\right)+O\left({\epsilon }^3\right).$

(3)

Here, we present the fast-scale ${\tau }_0=t$ in addition to the lower’s scales ${\tau }_1=\epsilon t$ and ${\tau }_2={\epsilon }^{2}t$ . Therefore, the following chain rule can be utilized to transform the derivatives of $\tau$ to be ${\tau }_n\left(n=0,1,2\right)$ as

$\frac{d}{d\tau }=\frac{\partial }{\partial {\tau }_0}+\epsilon \frac{\partial }{\partial {\tau }_1}+{\epsilon }^2\frac{\partial }{\partial {\tau }_2},$

$\frac{d^2}{d{\tau }^2}=\frac{{\partial }^2}{\partial {\tau }_0^2}+2\epsilon \frac{{\partial }^2}{\partial {\tau }_0\partial {\tau }_1}+{\epsilon }^2\left(\frac{{\partial }^2}{\partial {\tau }_1^2}+2\frac{{\partial }^2}{\partial {\tau }_0\partial {\tau }_2}\right)+O\left({\epsilon }^3\right).$

(4)

Here, the amplitude $F^{*}$ and the damping coefficient $\mu$ are specified as

$F^{*}={\epsilon }^{3}f,\mu ={\epsilon }^2\stackrel{\sim }{\mu }.$

(5)

Inserting equations (3)-(5) into (1) and equating like powers of $\epsilon$ on both sides, one may derive the below three partial differential equations (PDEs)

Order of $\epsilon$

$\frac{{\partial }^2{y}_1}{\partial {\tau }_0^2}+{\omega }^2{y}_1=0.$

(6)

Order of ${\epsilon }^2$

$\frac{{\partial }^2{y}_2}{\partial {\tau }_0^2}+{\omega }^2{y}_2=-2\frac{{\partial }^2{y}_1}{\partial {\tau }_0\partial {\tau }_1}-\stackrel{\sim }{\alpha }{y}_1^2.$

(7)

Order of ${\epsilon }^3$

$\frac{{\partial }^2{y}_3}{\partial {\tau }_0^2}+{\omega }^2{y}_3=f\hspace{0.17em}\mathit{cos}\left(\Omega {\tau }_0\right)-2\frac{{\partial }^2{y}_2}{\partial {\tau }_0\partial {\tau }_1}-\frac{{\partial }^2{y}_1}{\partial {\tau }_1^2}-2\frac{{\partial }^2{y}_1}{\partial {\tau }_0\partial {\tau }_2}-2\stackrel{\sim }{\mu }\frac{\partial y_1}{\partial {\tau }_0}-2\stackrel{\sim }{\alpha }{y}_1{y}_2-\eta y_1^3.$

(8)

It is worth mentioning that equations (6)–(8) represent a set of PDEs according to the various orders of $\epsilon$ , which can be treated sequentially. To achieve this, let us begin with the subsequent general solution of the homogeneous equation (6)

$y_1=A{e}^{i\omega {\tau }_0}+c.c.$

(9)

Here $A$ is undetermined complex functions and $c.c$ implies the conjugates of the previous terms.

Substituting the first-order solution (9) into the PDE (7) and subsequently eliminating the secular terms (ST), one gets the below second-order solution

$y_2=-\frac{\alpha A^2}{3{\omega }^2}{e}^{2i\omega {\tau }_0}+\frac{2}{{\omega }^2}A\overline{A}+c.c.$

(10)

Using the equations given earlier, we can identify the conditions required to remove ST from the second-order approximations as follows

$2i\frac{\partial A}{\partial {\tau }_1}=0,$

(11)

where the functions $A$ only depends on ${\tau }_2$ .

The third-order solution can be derived by setting equations (9)–(11) into equation (8) and then eliminating the produced ST. The next conditions need to be satisfied to remove such terms

$2i\omega \left(\frac{\partial A}{\partial {\tau }_2}+\stackrel{\sim }{\mu }A\right)+\left(\frac{11\stackrel{\sim }{\alpha }}{3{\omega }^2}+3\eta \right)A^2\overline{A}=0.$

(12)

Moreover, the third-order approximation is solved in the form

$y_3=\frac{f}{2\left({\omega }^2-{\Omega }^2\right)}{e}^{i\Omega {\tau }_0}+\left(\frac{2{\stackrel{\sim }{\alpha }}^2}{24{\omega }^4}+\frac{\eta }{8{\omega }^2}\right)A^3{e}^{3i\omega {\tau }_0}+c.c,$

(13)

with the solutions provided in (9), (10), and (13), it is possible to construct the desired approximate solution up to the third approximation.

It must be mentioned that this solution reduces the complexity of analyzing the DO, making it easier to understand the system’s behavior without solving highly nonlinear differential equations exactly. It provides valuable insights into the qualitative and quantitative behavior of the system, such as identifying stability regions, resonance frequencies, and possible chaotic behavior. Moreover, this solution allows for predictions about the system’s response to various inputs, which is crucial for designing and controlling engineering systems. It enables studying how different parameters affect the system’s dynamics, helping optimize the design and operation of devices modeled by the DO. Given its generality, the approximate solution can be adapted to various applications, including mechanical systems, electrical circuits, and even biological systems, where the DO model is relevant.

The achieved approximate solution is graphed in Figures 1–4 for a various value of the parameters in view of the data: $\mu =0.5,\omega =0.5,\alpha =0.1,f=0.5,$ and $\eta =0.2$ .

Figure 1.

Reveals the waves of the solution $u\left(\tau \right)$ for a various value of (a) $\mu =0.5,1.1,1.5$ , (b) $\alpha =0.1,0.2,0.4$ , (c) $f=0.09,0.2,0.5$ , and (d) $\omega =0.5,0.6,0.7.$

Figure 2.

Demonstrates the comparison between the approximate and numerical solutions.

Figure 3.

Shows the waves of $a\left(\tau \right)$ for various values of: (a) $\mu =0.5,1.1,1.5$ , (b) $\alpha =0.1,0.2,0.4$ , (c) $f=0.09,0.2,0.5$ , and (d) $\omega =0.5,0.6,0.7.$

Figure 4.

Explores the waves of $\theta \left(\tau \right)$ for different values of: (a) $\mu =0.5,1.1,1.5$ , (b) $\alpha =0.1,0.2,0.4$ , (c) $f=0.09,0.2,0.5$ , and (d) $\omega =0.5,0.6,0.7.$

Figure 1 examines the temporal evolution of the attained asymptotic solution, considering the selected parameter values. Parts (a) and (b) depict the plotted standing periodic waves, wherein the solution is influenced by the values of $\mu$ and $\alpha$ . As the values of $\mu$ and $\alpha$ increase, the displayed wave’s amplitude decreases. In part (c), however, the wave’s amplitude rises as the value of $f$ increases. On the other hand, the selected values of $\omega$ do not significantly affect the amplitude height, as noted in section (d). The figure illustrates periodic waves that describe the temporal dynamics of this solution. This simulation concluded that the proposed wave solution remains stable throughout the time.

It must be noted that periodic waves play a crucial role in studying dynamical systems, where they help to describe and predict the behavior of systems that exhibit regular, repeating patterns over time. ³⁰ In mechanical systems, periodic waves are essential in understanding oscillations, vibrations, and stability, which are keys to designing structures and machinery that can withstand repetitive stress. In electrical engineering, they are fundamental in analyzing alternating current circuits and signal processing. Periodic waves also enable the study of resonant phenomena, where the natural frequencies of a system can be exploited for efficient energy transfer and amplification. ³¹ Overall, periodic waves provide a framework for understanding and managing the complex behaviors of dynamical systems in various fields of science and engineering.

Figure 2 illustrates the temporal progressions of the approximate outcomes achieved and numerical ones attained across the same interval using the RK4. A more detailed analysis of the figure demonstrates a strong agreement between these solutions, highlighting the efficacy of the perturbation method utilized in generating precise approximate results that closely align with the numerical ones of the primary governing equation.

Based on the discussion mentioned above, we can conclude that the agreement between analytical and numerical solutions validates the accuracy of the mathematical models and assumptions used in the analysis.^32,33 This ensures that the models reliably represent the real-world system. The matching between them confirms the correctness and reliability of the numerical methods and algorithms being used. This is crucial for complex systems where analytical solutions might be difficult or impossible to obtain. Moreover, the consistency between both solutions increases confidence in the results and predictions derived from the system analysis, which is essential for making informed decisions in engineering and scientific applications. Analytical solutions often provide insight into the system’s fundamental behavior and underlying principles, while numerical solutions can handle more complex scenarios. ³⁴

Evaluation of resonance

In the current section, the resonance cases are classified. These scenarios appear when the denominators of the last two approximations are nullified. ³⁵ Accordingly, from the higher approximations (10) and (13), one can obtain the later resonance cases

(a) When $\omega \approx 0,$ the internal resonance appears.

(b) When $\omega \approx \Omega ,$ the external (primary) resonance occurs.

Therefore, if any of these conditions are satisfied, the system’s behavior becomes notably challenging under such circumstances. ³⁶ Moreover, the approximate solution previously obtained will persist if the oscillations do not align with resonance states. On the other side, to investigate the system’s stability, it’s imperative to meet the requirements for solvability and ME for this system. So, we will study the external resonance case when $\omega$ close to $\Omega$

$\Omega =\omega +\sigma ,\sigma =\epsilon \stackrel{\sim }{\sigma }.$

(14)

Now, by substituting (14) into the equations (7) and (8) addition to eliminate the ST, the solvability conditions can be obtained as follows

- Regarding the second-order approximation PDE

$2i\frac{\partial A}{\partial {\tau }_1}=0.$

(15)

- Regarding the third-order approximation PDE

$f{e}^{i{\tau }_1{\stackrel{\sim }{\sigma }}_1}-2i\omega \left(\frac{\partial A}{\partial {\tau }_2}+\stackrel{\sim }{\mu }A\right)-\frac{11\alpha A^2\overline{A}}{3{\omega }^2}-3\eta A^2\overline{A}=0.$

(16)

We may determine the functions $A$ as in the previous solvability requirements (15) and (16). Therefore, $A$ can be represented in a polar form as below

$A=\frac{1}{2}\stackrel{\sim }{a}{e}^{i\psi },a=\epsilon \stackrel{\sim }{a}.$

(17)

Here, $\stackrel{\sim }{a}$ and $\psi$ refer to the amplitude and phase of the solution $y$ , respectively. The phase is adjusted as below

$\theta ={\tau }_1\stackrel{\sim }{\sigma }-\psi .$

(18)

Utilizing equations (17) and (18) into (15) and (16) and subsequently segregating the imaginary components from the real ones yields the following ME

$\frac{da}{d\tau }=\frac{f}{2\omega }\mathit{sin}\hspace{0.17em}\theta -a\mu ,a\frac{d\theta }{d\tau }=\frac{f}{2\omega }\mathit{cos}\hspace{0.17em}\theta +a\sigma -\left(\frac{11}{24{\omega }^3}+\frac{3\eta }{8\omega }\right)a^3.$

(19)

The system of equation (19) comprises two ODEs that regulate the adjusted phase $\theta$ and the amplitude $a$ . This system can be solved numerically through the plane $a\theta$ with the initial conditions $\left(0\right)=0.001,$ $\theta \left(0\right)=0$ .

Furthermore, to examine the impact of different parameters, including $\mu ,\alpha ,f,$ and $\omega$ on the response of the system’s behavior, the mentioned figures are plotted. Figures 3 and 4 explain the transient response of the amplitude $a$ and modified phase $\theta$ of the system (19) to show the temporal change of the amplitude over time intervals. Notably, the curves depicted in these figures exhibit oscillation at the outset of the time before reaching stability at the remaining time, suggesting a tendency towards stability behavior within the system. Moreover, the decay behavior of these curves in this interval was expected based on the structure of the ME dynamical system (19), where a negative time-dependent exponential term can represent its solutions. In part (a) of Figure 3, the curves representing $a$ are affected by the damping coefficient $\mu$ , as $\mu$ ’s value rises, the displayed waves’ amplitudes diminish. As can be seen in part (b) of Figure 3, altering the quadratic coefficient $\alpha$ affects the system by varying $\mu$ , except for an increase in the oscillation observed at the beginning of the time interval. Conversely, in sections (c) and (d) of Figure 3, the wave of the amplitude $a$ escalates alongside higher values of the excitation $f$ and the natural frequency $\omega$ . Figure 4 shows the variation of the modified phase $\theta$ via time $\tau$ when the parameters $\mu ,\alpha ,f,$ and $\omega$ have various values. It was found that the curves of parts (b) and (d) have a number of oscillations greater than the waves of parts (a) and (c). So, the change in the values $\mu ,$ $f$ have a slight effect on the waves contrarily the parameters $\alpha ,$ and $\omega$ .

Figure 5 illustrates the projection of the amplitude $a$ and phase curve $\theta$ onto a plane. The graphical representation depicts spiraling curves that converge towards a singular point. This convergence may arise from the full alignment between the trajectories’ time of the amplitude and its corresponding phase. Consequently, the system’s dynamic motion remains stable throughout the studied time interval.

Figure 5.

Presents the phase plane $\theta a$ for a various value of: (a) $\mu =0.5,1.1,1.5$ , (b) $\alpha =0.1,0.2,0.4$ , (c) $f=0.09,0.2,0.5$ , and (d) $\omega =0.5,0.6,0.7.$

The chaotic response

This section delves into the intriguing realm of chaotic behavior demonstrated by our dynamical model. Bifurcation diagrams, LES, and Poincaré maps are investigated to provide a measure of chaos and stability. It is crucial to note that even slight parameter adjustments can considerably influence the system’s behavior. Consequently, the system exhibits two distinct modes of motion, each of which will be elucidated below.

To illustrate the distinct types of motions, we simulated the bifurcation diagram of $u$ versus the excitation $f$ as a bifurcation parameter and the LES of ${\lambda }_n\left(n=1,2\right)$ in Figures 6–8 under various value for the parameter $\omega =1.0,1.1,1.2$ , respectively.^37,38

Figure 6.

Bifurcation diagram and LES of the system $\omega =1$ .

Figure 7.

Bifurcation diagram and LES of the system at $\omega =1.1$ .

Figure 8.

Bifurcation diagram and LES of the system at $\omega =1.2$ .

Figure 6(a) elucidates the bifurcation diagram of $u$ , while Figure 6(b) depicts the LES characterizing distinct types of motion at $\omega =1.0$ . The mentioned presents the bifurcation diagram and the corresponding LES for the system under study, as shown in Figure 6(a) and (b). The system exhibits periodic behavior for excitation amplitude $f$ values less than 0.5 and greater than 1.2, as depicted in Figure 6(a). Within the range of $f$ from 0.51 to 0.56, the system undergoes period doubling, transitioning from a single periodic state to increasingly complex periodic states. As the excitation amplitude $f$ increases further, the system enters a chaotic regime between $f=0.56$ and $f=1.20$ .

The LES in Figure 6(b) corroborates these observations. Negative LES indicate periodic behavior, consistent with the bifurcation diagram for $f<0.5$ and $f>1.2$ . Positive LES, observed in the range $f=0.56$ to $f=1.20$ , confirm the system’s transition into chaos.

A Poincaré map is provided in Figure 9 to illustrate the observed chaotic behavior further. This map offers a visual representation of the system’s dynamics in the chaotic regime, highlighting the complex and aperiodic nature of the state.

Figure 9.

Poincaré map for the approximate solution at $\omega =1$ .

Figure 7(a) shows the bifurcation diagram of $u$ , while Figure 7(b) elucidates the LES characterizing discriminate types of motion at $\omega =1.1$ . Under studying the bifurcation diagram and its corresponding LES, it is shown that the system exhibits periodic behavior of discontinuous interval for excitation amplitude values of $f$ , as depicted in Figure 7(a). As in the range of $f$ from $0.58$ to $0.7$ , the system undergoes period doubling, transitioning from a single periodic state to increasingly complex periodic states. As the excitation amplitude $f$ increases further, the system enters a chaotic regime between $f=0.49$ and $f=1.20$ .

Corroborating these observations in Figure 7(b), zero or negative LES indicate periodic behavior, and positive LES confirm the system’s transition into chaos. Also, the Poincaré map in Figure 10 shows chaotic behavior; the points within the diagram begin to diverge randomly. This divergence serves as a clear indication, as mentioned earlier, of the system’s chaotic response.

Figure 10.

Poincaré map for the approximate solution at $\omega =1.1$ .

Figure 8 shows the same behavior in Figure 7 but with less discontinuous interval, and this occurs when $\omega =1.2$ . Also, the Poincaré map in Figure 11 shows a chaotic response that confirms the system’s complexion.

Figure 11.

Poincaré map for the approximate solution at $\omega =1.2$ .

It is crucial to state that the chaotic response of a dynamical system, LES, and Poincaré maps enable a deeper understanding of the complex dynamics of nonlinear systems, aiding in the prediction, control, and analysis of such systems across various scientific and engineering disciplines.^39,40

Moreover, we were eager to observe the effect of the quadratic term on the studied oscillator. Therefore, we investigated the influence of its coefficient $\alpha$ , as a bifurcation parameter on the system’s chaotic behavior. Figure 12(a) shows the bifurcation diagram of $u$ , while Figure 12(b) explains the LES characterizing different types of motion. The system exhibits periodic behavior for $\alpha$ values that less than 0.282, as included in Figure 12(a). In the range of $\alpha$ from 0.282 to 0.36, the system, in the case of period doubling, transitions from a single periodic state to increasingly complex periodic states. The LES in Figure 12(b) confirms these observations. The negative LES indicates a periodic behavior consistent with the bifurcation diagram for $\alpha <0.282$ . The positive LES observed in the range $\alpha >0.282$ confirm the system’s transition to chaos. To further illustrate the observed chaotic behavior, the Poincaré map is shown in Figure 13. This plot provides a visual representation of the system’s dynamics in the chaotic state and highlights its complexity and aperiodic nature.

Figure 12.

Bifurcation diagram and LES of the system.

Figure 13.

Poincaré map reveals the chaotic behavior of the system at $\alpha =0.35$ .

Analysis of stability and steady-state

In this section, we illustrate the stability analysis and the determination of steady-state solutions based on resonance response curves. Steady-state oscillations emerge once transient processes are mitigated by an effective damping system. ⁴¹ This stabilization results in the time derivatives of the phase and amplitude within the ME becoming zero (19). Subsequently, we transform these equations from ODEs into algebraic equations to simplify the analysis and facilitate the characterization of stable equilibria, as follows

$\frac{f}{2\omega }\mathit{sin}\hspace{0.17em}\theta -a\mu =0,\frac{f}{2\omega }\mathit{cos}\hspace{0.17em}\theta +a\sigma -\left(\frac{11}{24{\omega }^3}+\frac{3\eta }{8\omega }\right)a^3=0.$

(20)

Furthermore, upon eliminating the variable $\theta$ algebraically from (20) and (21) will transform into a correlation between $a$ and $\sigma$ in the following manner

$f^2=4{\omega }^2\left\{a^2{\mu }^2+{\left[\left(\frac{11}{24{\omega }^3}+\frac{3\eta }{8\omega }\right)a^3-a\sigma \right]}^2\right\}.$

(21)

The stability of oscillations at the steady state is demonstrated by examining the oscillations in proximity to a fixed position. According to Lyapunov’s perspective, fixed points such as stationary solutions exhibit asymptotic stability only if the real parts of all eigenvalues are negative. ⁴² We depict the stability of the fixed point utilizing the resonance response curves derived from the equation’s numerical solution (21). In Figure 14, the frequency response is illustrated in the $\sigma a$ plane across different values of specific parameters, highlighting the notable influence of these parameters. The impact on stability is evident from the depiction. We observe that the stability of the fixed points plays a crucial role in determining the system’s overall stability. In certain regions, fixed points undergo a loss of stability, transitioning into unstable states. It is worth emphasizing that stable fixed points are represented by blue, red, and green lines, each corresponding to different parameter configurations, whereas black dots identify unstable points. It is clear from these plots that, with the increasing values of the variables $\mu ,\alpha$ and $\omega$ , the stability zone increases as shown in parts in Figures 14(a), (b) and (d).

Figure 14.

Demonstrates the frequency response in the plane $\sigma a$ for the values of: (a) $\mu \left(=0.5,1.1,1.5\right)$ , (b) $\alpha \left(=0.1,0.2,0.4\right)$ , (c) $f\left(=0.09,0.2,0.5\right)$ , and (d) $\omega \left(=0.5,0.6,0.7\right)$ .

Finally, we conducted comparisons between the nonlinear oscillators Song et al. (2023), Liu et al. (2023), and Moatimid et al. (2023) and the examined one, as in Table 1:

Table 1.

Shows comparisons between the studied oscillator and other oscillators.

No.	Comparison items	The current study	The examined work in Song et al. 23	The investigated work in Liu et al. 26	The obtained work in Moatimid et al. 29
1	The quadratic term	It is taken into account	-	-	-
2	Approximate solution	It has been achieved up to third approximate	-	-	It was examined up to second approximate
3	Comparison between approximate solution and the numerical one	It has been done	-	-	It was obtained
4	Bifurcation diagrams	Achieved	Achieved	Achieved	-
6	LES	They have been achieved and analyzed	-	-	-
7	Poincaré maps	The Poincaré maps have been simulated	-	-	-

Conclusion

An analytical solution was provided for the DO subjected to sinusoidal external excitation. The solution to the DO’s equation of motion was obtained through Lagrange’s equation and analytically solved up to the third-order approximation utilizing MS. The validity of the obtained solution was confirmed by comparing it with the numerical solution of the governing equation, demonstrating a significant level of matching. External resonance was also examined to obtain the solvability conditions and ME. To elucidate the behavior of the dynamical at a particular moment, we created and discussed a specific time history plot of the obtained solution, considering the precise parameter values utilized. The ME of the system was numerically solved to demonstrate the amplitude and fitted phase of the corresponding data, considering the parameters under examination. The stability criterion is employed to analyze both stable and unstable regions within the system and examine the impact of various parameters. Bifurcation diagrams and LEs spectra are highlighted to demonstrate diverse types of system motions, along with the illustration of the Poincaré diagram. It can be applied in many fields, like physics, mathematics, engineering, and biology, for example, studying oscillations, vibrating strings and membranes, structural dynamics, and electrical circuits. The study of the DO has provided deep insights into nonlinear dynamic systems, particularly its chaotic behavior and the impact of external forces and damping on system stability. In future research, we can add a harvest device with the control system. Recommendation: Incorporating machine learning and advanced computational tools will enhance the predictive capabilities of Duffing models in chaotic systems.