On chaotic behavior,stability analysis,and vibration control of the van der Pol–Mathieu–Duffing oscillator under parametric force and resonance

Abstract

The van der Pol–Mathieu–Duffing oscillator (VMDO) finds applications across diverse fields due to its ability to model complex dynamic behaviors. Oscillators are versatile tools used in fields such as mechanics, biology, and electrical engineering. Its ability to model complex behaviors makes it an important subject of study for researchers and engineers studying nonlinear dynamic systems. A parametric forcing excited the VMDO under feedback control is investigated and discussed in the most severe resonance scenario. The multiple-scales-strategy (MSS) is utilized to obtain the approximate solution (AS). Furthermore, the AS is validated against the numerical solution (NS) obtained using the Runge–Kutta of fourth-order (RK-4) method. A negative-velocity-feedback (NVF) and negative-cubic-velocity-feedback (NCVF) controllers are integrated into the primary system to mitigate unwanted vibrations, which can significantly impact the system’s efficiency, particularly under resonance conditions. The stability analysis is thoroughly examined, and optimal feedback gains are selected to suppress amplitude peaks. A range of response curves is presented to illustrate and compare the effectiveness of the controllers. Bifurcation diagrams and Poincaré maps (PM) are utilized to investigate the system’s diverse motions, providing valuable insights into its complex behavior and the variations it undergoes under different conditions. The results are explained through the displayed curves and provide insights into the system’s dynamics. The VMDO is a versatile model in various scientific and engineering disciplines. Its ability to embody nonlinear dynamics makes it invaluable for studying the stability, resonance, and chaotic behavior of complex systems. As research continues, its applications may expand into new areas of technology.

Introduction

An essential model for examining complicated oscillatory patterns is the Duffing oscillator (DO), a basic but fascinating system in nonlinear dynamics. This oscillator, which dates back to Georg Duffing’s work at the beginning of the twentieth century, is a prime example of how complex, multilayered phenomena like chaos, bifurcations, and multi-stability may appear in even the most basic mechanical systems. In addition, the theoretical physics study has real-world implications across many fields, including engineering and biology, where it helps design and analyze systems subject to nonlinear pressures.

In Ref. 1,2, the authors constructed efficient oscillators and oscillator systems using a methodical approach. The foundation of this method is a basic classification of oscillators based on their inner temporal constraints. In Ref. 3, the authors examined how two periodic excitements with wildly disparate frequencies affect the nonlinear reaction of a bistable VMDO. They demonstrated that a symmetrical fluctuating bistable potential might be transformed into a symmetric fluctuating monostable potential by methodically varying the high-frequency drives and the parametric oscillation’s strengths. Through theoretical and computational analysis, ⁴ explored the vibrational resonance’s potential in a softened DO, looking at both underdamped and overdamped cases. The authors in Ref. 5 investigated the dynamic properties of fractions VMDO under the influence of two periodic excitements and an evenly distributed period delay. The pitchfork bifurcation of the system, which is caused by a low-frequency stimulus from outside and a high-frequency internal excitement, was first examined by the authors. Then, when the two excitement frequencies are the same, they study the fractal behavior of the whole system.

In Ref. 6, the authors examined a forced nonlinear self-sustained oscillator operating near primary resonance, focusing on the suppression of hysteresis effects. The oscillator was subjected to quick forcing to study the inhibition. To identify the influence zone, the quasi-periodic control domain, and the zone of hysteresis area, an analytical technique utilizing perturbations modeling was carried out. In Ref. 7, the authors sought to investigate the dynamic behaviors of the whole excitation utilizing the fast-slow analysis technique. To analyze the relaxation oscillations caused by different modulation modes, a three-dimensional visualization method was proposed, providing a thorough grasp of the system’s dynamics from several perspectives. In Ref. 8, the authors explored the bifurcation behavior in a system of coupled oscillators with periodically varying coupling. They identified the area of coupling strengths where oscillations were effectively regulated, utilizing basin stability analysis and bifurcation diagrams to characterize the system’s behavior. In Ref. 9, a comprehensive investigation of the model, employing standard nonlinear analysis techniques, uncovered intriguing and remarkable dynamic phenomena. These included a tree of eight parallel bifurcations, Hopf ones, various coexisting oscillation modes, and an attractor of eight chaotic spirals, among other features.

In Ref. 10, using both integer and fractional ordering differential operators to investigate nonlinear DO. The fractional-order nonlinear DO was shown to get only one solution under suitable initial circumstances using the fixed-point theory. An analysis of the most dangerous resonance scenario of an externally excited DO under feedback control was conducted and detailed in Ref. 11. In order to obtain the AS for this system, the various time scales approach was used to determine whether the time delay on the control loop exists or not. To lessen the amplitude peak, a suitable stability analysis was also carried out, and appropriate options for the time delay and the feedback gains were identified. To illustrate and contrast controller effects, various response curves are used. Additionally, RK-4 method is used to compare AS with NS. The implementation of a time-delayed controller for a damped nonlinear stimulated DO was the main objective of Ref. 12 To mitigate the system’s nonlinear oscillations, the controller employs position and velocity time delays. A more precise approximation was achieved by applying a sophisticated Homotopy perturbation method. Plotting of this solution’s time development for several parameters was done. The results showed that the increased Homotopy perturbation technique minimizes solution errors and broadens its applicability, making it suitable for various damped nonlinear oscillators. Fractional damping, a trait common to many physical situations, was used to study the nonlinear DO in Ref. 13. Both smaller and larger damping parameter values were used in the investigation to examine the system.

The study in Ref. 14 investigated the dynamics of a Duffing equation with a fifth-order nonlinear element, similar to the developed model in Ref. 5 but with a wider range of bifurcation analysis. In contrast to Ref. 5, this research studies the bifurcation patterns in the system using both analytical and numerical methods. In Ref. 15, a nonlinear quadratic system for the steady-state motion of a controlled superstructure is formulated, including its AS. Response curves and stability charts were employed to investigate bifurcation properties. The complete response of the system was also numerically simulated and demonstrates how well the added controller performs in removing resonant vibrations from the structure and bringing non-resonant unstable motion under control. In Ref. 16, the authors examined how noises affect forced nonlinear dynamics and time-delayed DO, with a focus on their influence on the maximum amplitude of oscillations and the characteristic frequency of the steady state for various time delays and driving force values. The nonlinear vibration properties of a nonlinear Zener model with linear and nonlinear displacement delay feedback control under harmonic excitation were examined in Ref. 17. The MSS was employed to determine the amplitude–frequency curve of the system’s external resonance. The stability theory of Lyapunov was then used to determine the system’s stability requirements. Time delay, displacement feedback gain coefficient, nonlinear term coefficient, damping coefficient, and external excitation were among the parameters whose effects on the behavior of the dynamical system were also examined.

The study in Ref. 18 analyzed the impact of nonlinearity on the exceptional point displacement in a two-coupled-DO system, focusing on the influence of coupling coefficient variations and insertion loss. The study revealed that altering the nonlinearity coefficient reduced the oscillation amplitude and a resonance frequency shift at the exceptional point. The work ¹⁹ studied the DO with a fractional damping term. In Ref. 20, the DO with impact is discretized in order to create sub-implicit mappings. These mapping structures were used to establish distinct impact periodic motions for the DO with effects. The bifurcation trees of impact chatter periodic motions were obtained through a semi-analytical approach. In Ref. 21, a dual-stage approach was proposed for identifying mechanical component faults. The initial step involves detecting faults by studying the trajectory of phase space for the DO. Subsequently, the method of stochastic resonance, based on the Duffing equation, was applied to improve the noise ratio of the original signal. In this step, the optimization of two barrier parameters in this equation, performed with the artificial bee colony algorithm, was essential. In Ref. 22, a study that analyzed the chaotic behavior of the VMDO under external and parametric excitations was presented, using both analytical and numerical methods. Critical curves that distinguish between chaotic and nonchaotic regions were derived.

The resonance response of many damped, linearly coupled DO under an acting sinusoidal force on the first oscillator was explored in Ref. 23. This investigation combined analytical determination of the system’s fixed points with numerical simulations using the RK-4 method. The study in Ref. 24 analyzed two coupled systems, where one functioned as the primary driver and the other as the secondary-driven system. The driver system was characterized as a time-delayed oscillator, while the secondary one had minimal delay. By focusing on the impact of the primary driver system as the sole external force on the secondary system, the study delved into how it affected the amplitude, frequency, and conditions leading to resonance in the output of the secondary system. In Ref. 25, an extensive parametric analysis was performed to evaluate the influence of key design parameters on the system’s overall stability. The findings revealed complex nonlinear dynamics, such as detached resonance curves and transitions from periodic stability to chaotic behavior. In Ref. 26, the authors utilized the global error minimization method to enhance and extend an analytical technique in order to obtain a third-order AS for strongly nonlinear Duffing-harmonic oscillators. A comparison was then made between the results achieved through this method and those obtained from other analytical approaches, as well as the NS.

In Ref. 27, the authors presented the results of synthesizing and further validating previously established mathematical models for two different mechanical oscillators, each with one degree of freedom and subjected to a harmonic excitation force. One oscillator features magnetically adjustable elasticity, leading to a double symmetric potential well, while the other uses springs of linear mechanical with a one-sided motion limiter. In Ref. 28, the authors proposed a new technique for attaining high-precision periodic movements via implicit mapping, which was confirmed by examining the saddle-node bifurcation domains. Through an analysis of the evolution curve of periodic movements in a DO with hardening properties, it was discovered that the primary resonance can be utilized to improve the characteristic components in the generated signal. In Ref. 29, the study examined the master-slave synchronization scheme between a Rayleigh–Duffing system and DO. The authors examined both elastic and dissipative couplings individually, as well as a combination of both. The comparison of results allowed the identification of the optimal coupling mechanism for synchronization between the oscillators. Numerical analysis indicated that, regardless of elastic or dissipative coupling, complete synchronization occurred solely in one state of the slave system. It is widely recognized that dynamical systems, particularly those involving vibrations as explored in Ref. 30–33, play a crucial role in mechanics, given their diverse applications, especially in engineering.

The Variational Iteration Method (VIM) is an iterative technique introduced by Ref. 34 for solving both linear and nonlinear differential equations. It utilizes the Lagrange multiplier method to construct an iterative formula that progressively refines the solution. Unlike conventional numerical methods, VIM does not require linearization, discretization, or restrictive assumptions, making it a more straightforward approach. VIM is widely applied to solving nonlinear and fractional differential equations, as well as engineering problems such as heat conduction, fluid mechanics, elasticity, and vibrations. Additionally, it has been extended to fractional calculus and chaotic systems, offering superior approximations in cases where traditional perturbation methods prove ineffective. Also, the Homotopy Perturbation Method (HPM), introduced in Ref. 35, integrates homotopy theory from topology with perturbation techniques. Unlike traditional perturbation methods, HPM does not require a small perturbation parameter, making it particularly effective for strongly nonlinear problems. It is widely applied in engineering mechanics, reaction-diffusion systems, and nonlinear fluid flow, as well as in computational mechanics and applied physics. HPM has also been extended to fractional differential equations and chaotic dynamics, demonstrating its versatility in handling highly nonlinear systems. Frequency formulation, on the other hand, focuses on directly determining frequency–amplitude relationships for nonlinear oscillators without relying on small perturbation parameters. This approach employs techniques such as the energy balance method, harmonic balance method, and modified Lindstedt–Poincaré technique. It finds applications in aerospace engineering, structural mechanics, and nonlinear oscillators, including the Duffing oscillator, van der Pol oscillator, and predator-prey systems. Moreover, frequency formulation has become increasingly significant in micro-electromechanical systems (MEMS) and nonlinear optics, offering accurate frequency predictions where numerical simulations are computationally expensive. ³⁶ presented a review of the periodic behavior of micro-electro-mechanical systems (MEMS) analyzed using various analytical techniques, including the HPM, variational iteration method, variational theory, He’s frequency formulation, and the Taylor series method. Additionally, it explores fractal MEMS systems and discusses future research directions. The focus of this mini-review is primarily on advancements made in the past decade, and as a result, the references are not intended to be comprehensive. ³⁷ investigated the periodic motion of MEMS, which are influenced by a singularity that complicates obtaining an exact solution and understanding their dynamic behavior. To address this, the study employs frequency formulation to analyze the relationship between frequency and amplitude. The results reveal that when the amplitude surpasses a critical threshold, the periodic motion transitions into pull-in instability. This insight enhances the ability to detect unsafe operating conditions, simplifying the system’s warning mechanism. Furthermore, the frequency-amplitude relationship provides a valuable tool for the precise and reliable optimization of MEMS design.

The study of VMDO plays a pivotal role in the development of various technological and scientific fields. From MEMS devices and power systems to biomedical engineering, aerospace, and robotics, VMDO models provide fundamental insights into the behavior of systems subject to the effects of variable mass and damping. Understanding these dynamics enables improved system performance, enhances stability, and facilitates the development of innovative solutions in engineering and applied science. As VMDO research continues to advance, its applications are expected to expand, contributing to the linkage of theoretical advances with real-world technological progress.

The study in Ref. 38 explores the discrete-time T system, analyzing its fixed points and topological divisions. It establishes the presence of a Neimark–Sacker bifurcation under specific parameters and confirms the Flip-NS bifurcation using explicit criteria. Center manifold theory determines the bifurcation directions, while numerical simulations validate the theoretical findings. The 0-1 test is applied to detect chaos, and a hybrid control method is implemented to suppress chaotic behavior. A novel chaotic system incorporating two flux-controlled memristors, exhibiting symmetric bifurcation and multi-stability within a four-dimensional framework is studied in Ref. 39. Computational simulations reveal complex dynamical behaviors, including symmetric bifurcations, multi-stability, and extreme sensitivity to initial conditions. The work in Ref. 40 investigates the resistive-capacitive shunted Josephson junction with a topologically nontrivial barrier coupled to a linear RLC resonator. The governing rate equations are derived using Kirchhoff’s current and voltage laws. The dynamics and control of the permanent magnet synchronous motor with load torque, identifying three fixed points: one saddle and two saddle-focus points, is examined. ⁴¹ In Ref. 42, the authors analyzed discrete-time predator-prey systems by determining the minimum prey consumption required for predator reproduction and examining the system’s stability and bifurcation behavior.

This paper investigates the chaotic dynamics and stability of the VMDO. NVF and NCVF controllers are implemented to mitigate potentially hazardous vibrations. Designing feedback laws for vibration suppression requires careful consideration of system characteristics, control objectives, and potential disturbances. While conventional controllers are simple, their limited adaptability and robustness necessitate the exploration of advanced strategies. Adaptive and sliding controllers offer promising alternatives, improving performance in various resonance scenarios and enhancing resilience to uncertainties, thus expanding the range of effective vibration suppression methods. Therefore, the control was applied to the system in all resonance cases, demonstrating its effectiveness in each scenario. Consequently, the study was conducted under the worst resonance conditions with this control in place. The MSS method is employed to analyze the AS up to the second approximation, while the RK-4 method is utilized for a graphical determination of the NS, facilitating comparison with the AS to enhance the precision of the perturbation method. Solvability conditions are obtained by eliminating secular terms, allowing for the identification of all resonance cases. Stability is analyzed using the NS of the modulation equations, with resonance curves plotted to show stability and instability regions as parameter values change. Furthermore, the system’s behavior is explored through bifurcation graphs and PM, providing a vivid representation of nonlinear phenomena important for various engineering and physics applications.

System model and mathematical analysis

The VMDO is a complex nonlinear dynamical system that combines features of the VMDO. This combination enables the modeling of a wide range of vibrational phenomena across various scientific and engineering disciplines. The mathematical model for the VMDO subjected to a parametric excitation force, incorporating both a composite NVF and NCVF, is presented in Ref. 43, where the governing equation has the form

$\ddot{u}\left(t\right)+[1-h\hspace{0.17em}\mathit{cos}\left(\omega t\right)]u\left(t\right)-\left[\alpha -\beta x^3\left(t\right)\right]\dot{u}\left(t\right)-\gamma u^3\left(t\right)=q{\Omega }^2\hspace{0.17em}\mathit{cos}\hspace{0.17em}u\hspace{0.17em}\mathit{cos}\left(\Omega t\right)-G_1\dot{u}\left(t\right)-G_2{\dot{u}}^3\left(t\right).$

(1)

The terms in equation (1) correspond to various aspects of the system, representing its physical components as follows: $\alpha$ and $\beta$ are the linear and nonlinear damping coefficients, respectively. Moreover, $u\left(t\right)$ denotes the displacement; $\Omega ,\omega$ are the external frequencies excitations, $h,q$ are the parametric excitation forces, $\gamma$ is the nonlinear stiffness coefficient. Finally, $G_1$ and $G_2$ are feedback gains that are also used, as shown in Figure 1, to avoid increasing the amplitudes during resonance. The schematic diagram in Figure 1 illustrates a nonlinear dynamical system described by a differential equation. It consists of inputs, transformations via gain functions, derivative operators, and nonlinear cubic functions, ultimately producing outputs through sensors. Artificial intelligence (AI) driven problem solving methods can be utilized to analyze and interpret the system’s behavior, particularly for predicting stability, bifurcations, and chaotic dynamics. Additionally, AI can be employed to optimize control parameters, ensuring stability and enhancing system efficiency. These applications are valuable in fields such as vibration analysis, mechanical system health monitoring, and predictive maintenance. By integrating AI-based techniques, the system’s performance can be better understood, optimized, and regulated. AI can extract meaningful insights from sensor data, detect underlying patterns, and support decision-making, making it an essential tool for studying complex dynamical systems like the one depicted in the diagram.

Figure 1.

A schematic diagram for the VMDO with the NVF and NCVF.

Methodology

The goal of this section is to derive the AS of equation (1) by applying MSS. This method is used to analyze nonlinear and weakly perturbed systems accurately by separating time scales, preventing divergence, and providing insights into stability and bifurcation. It is especially useful for vibration systems and energy harvesting models with multiple interacting frequencies. It is applied up to the second approximation, ensuring good accuracy and exploring various resonant conditions. We also use a Taylor expansion to approximate the trigonometric function to the second order, which holds true in the vicinity of static equilibrium. ⁴⁴ Therefore, equation (1) can be rewritten as follows:

$\ddot{u}\left(t\right)+[1-h\hspace{0.17em}\mathit{cos}\left(\omega t\right)]u\left(t\right)-\left[\alpha -\beta u^3\left(t\right)\right]\dot{u}\left(t\right)-\gamma u^3\left(t\right)=q{\Omega }^2\hspace{0.17em}\mathit{cos}\left(\Omega t\right)-\frac{q}{2}{\Omega }^2{u}^2\left(t\right)\mathit{cos}\left(\Omega t\right)-G_1\dot{u}\left(t\right)-G_2{\dot{u}}^3\left(t\right).$

(2)

Applying MSS, then the approximations of $u\left(t\right)$ is expressed in the form of a power series as:

$u\left({\tau }_o,{\tau }_1,\epsilon \right)={\displaystyle \sum }_{i=0}^1{\epsilon }^i{u}_i\left({\tau }_o,{\tau }_1\right)+O\left({\epsilon }^2\right),$

(3)

in which ${\mathit{\tau }}_{\mathit{n}}={\mathit{\epsilon }}^{\mathit{n}}\mathit{\tau }\left(\mathit{n}=\mathbf{0},\mathbf{1}\right)$ are the different time scales, where ${\tau }_0$ represents a fast time scale and ${\tau }_1$ denotes a slow time scale. To determine $u$ as a function of ${\tau }_n$ , we can change $t$ to ${\tau }_n$ according to the chain rule as follows:

$\frac{d}{dt}=D_0+\epsilon D_1+O\left({\epsilon }^2\right),\frac{d^2}{d{t}^2}=D_o^2+2\epsilon D_0{D}_1+O\left({\epsilon }^2\right),$

(4)

where $D_n=\frac{\partial }{\partial {\tau }_n};{\tau }_n={\epsilon }^{n}t\left(n=0,1\right).$

Considering that

$\alpha =\epsilon \stackrel{\sim }{\alpha },\beta =\epsilon \stackrel{\sim }{\beta },h=\epsilon \stackrel{\sim }{h},\alpha =\epsilon \stackrel{\sim }{\alpha },\gamma ={\epsilon }^2\stackrel{\sim }{\gamma },q=\epsilon \stackrel{\sim }{q},G_1=\epsilon {\stackrel{\sim }{G}}_1,G_2=\epsilon {\stackrel{\sim }{G}}_2$

(5)

Inserting equations (3)–(5) in (2) and equating the same power of coefficients $\epsilon$ , to obtain:

Order of ${\mathit{\epsilon }}^{\mathbf{0}}$

$\left(D_o^2+{\omega }_o^2\right)u_0=0.$

(6)

Order of $\mathit{\epsilon }$

$\left(D_o^2+1\right)u_1=\stackrel{\sim }{q}{\Omega }^2\hspace{0.17em}\mathit{cos}\left(\Omega T_0\right)+\stackrel{\sim }{h}{u}_0\hspace{0.17em}\mathit{cos}\left(\omega T_0\right)-\frac{1}{2}\stackrel{\sim }{q}{\Omega }^2\hspace{0.17em}\mathit{cos}\left(\Omega T_0\right)u_0^2+\stackrel{\sim }{\gamma }{u_0}^3+\stackrel{\sim }{\alpha }{D}_o{u}_0-{\stackrel{\sim }{G}}_1{D}_o{u}_0-\stackrel{\sim }{\beta }{u}_1^2{D}_o{u}_0-{\stackrel{\sim }{G}}_2{\left(D_o{u}_0\right)}^3-2D_o{D}_1{u}_0.$

(7)

Therefore, equations (6) and (7) can be solved subsequently. Consequently, we start with the homogenous solution of (6), which is given as

$u_0=A\left({\tau }_1\right)e^{i{\tau }_0}+c.c.$

(8)

Here, $c.c.$ stands for the complex conjugate of the preceding term, while $A$ is a complex function within ${\tau }_1$ . Substituting (8) into (7) with the removal of terms that yield an unbounded solution to get the desired asymptotic solution, the mentioned secular terms are symbolized as follows:

$3\stackrel{\sim }{\gamma }{A}^2\overline{A}-2i\frac{\partial {\rm A}}{\partial {\tau }_1}+i\left(\stackrel{\sim }{\alpha }-{\stackrel{\sim }{G}}_1\right)A-i\stackrel{\sim }{\beta }{A}^2\overline{A}-3i{\stackrel{\sim }{G}}_2{A}^2\overline{A}=0.$

(9)

Then, the second-order solutions may subsequently be expressed as:

$u_1=\frac{{\Omega }^2}{2\left(1-{\Omega }^2\right)}\stackrel{\sim }{q}{e}^{i\Omega \tau }+\frac{\stackrel{\sim }{h}}{2}\left\{\frac{A}{\left[1-{\left(1+\omega \right)}^2\right]}{e}^{i\left(1+\omega \right)\tau }+\frac{A}{\left[1-{\left(1-\omega \right)}^2\right]}{e}^{i\left(1-\omega \right)\tau }\right\}$

$-\frac{{\Omega }^2\stackrel{\sim }{q}}{4}\left\{\frac{A^2}{\left[1-{\left(2+\Omega \right)}^2\right]}{e}^{i\left(2+\Omega \right)\tau }+\frac{{\overline{A}}^2}{\left[1-{\left(\Omega -2\right)}^2\right]}{e}^{i\left(\Omega -2\right)\tau }+2\frac{A\overline{A}}{1-{\Omega }^2}{e}^{i\Omega \tau }\right\}-\frac{\stackrel{\sim }{\gamma }{A}^3}{8}{e}^{3i\tau }+\frac{\stackrel{\sim }{\beta }{A}^3}{8}i{e}^{3i\tau }-\frac{{\stackrel{\sim }{G}}_2{A}^3}{8}i{e}^{3i\tau }+c.c.$

(10)

The complex function $A$ can be determined from the condition (9) by using the starting conditions $A\left(0\right)=z_0$ where $z_0$ is known as initial quantity.

Now we study the worst operating mode of the system due to the resonance condition:

$\omega =2+\epsilon \sigma ,$

(11)

where $\sigma$ is the detuning parameter. Substituting (11) into (7), hence the secular terms will be:

$\frac{\stackrel{\sim }{h}\overline{A}}{2}{e}^{i\sigma {\tau }_1}-2i\frac{\partial {\rm A}}{\partial {\tau }_1}+3\stackrel{\sim }{\gamma }{A}^2\overline{A}+i\left(\stackrel{\sim }{\alpha }-{\stackrel{\sim }{G}}_1\right)A-i\stackrel{\sim }{\beta }{A}^2\overline{A}-3i{\stackrel{\sim }{G}}_2{A}^2\overline{A}=0.$

(12)

Converting $A$ to the polar form, we have

$A=\frac{a}{2}{e}^{i\psi },$

(13)

where $a$ and $\psi$ are, respectively, the amplitude and phase of the given system. Upon substituting (13) into (12) and matching the real and imaginary parts, we get

$\frac{da}{d\tau }=\frac{1}{4}ah\hspace{0.17em}\mathit{sin}\hspace{0.17em}\theta +\frac{a}{2}\left(\alpha -G_1\right)-\beta \frac{a^3}{8}-\frac{3a^3}{8}{G}_2,2a\frac{d\theta }{d\tau }=\frac{1}{4}ah\hspace{0.17em}\mathit{cos}\hspace{0.17em}\theta +\frac{a}{2}\sigma +\frac{3a^3}{8}\gamma ,$

(14)

where $\theta =\sigma {\tau }_1-2\psi$ .

To obtain a stationary solution for amplitude and phase, inserting $\dot{a}=\dot{\theta }=0$ into (14), squaring and adding both sides, ⁴⁵ yields:

$r_1{a}^4+r_2{a}^2+r_3=0,$

(15)

where

$r_1=\frac{1}{4h^2}\left({\beta }^2+9{\gamma }^2+9G_2^2+6\beta G_2\right),$

$r_2=\frac{1}{h^2}\left[2G_1\left(\beta +3G_2\right)-2\alpha \beta +6s\gamma -6\alpha G_2\right],r_3=\frac{1}{h^2}\left[-h^2+4\left(s^2+{\alpha }^2\right)-8\alpha G_1+4G_1^2\right].$

(16)

An important aspect of steady-state vibration is the examination of its stability in regions near the fixed points. To achieve this purpose, we take the below expression of $a$ and $\theta$ into account

$a=a_0+a_1,\theta ={\theta }_0+{\theta }_1,$

(17)

where $a_0$ and ${\theta }_0$ denote the unperturbed steady-state solutions, and $a_1$ and ${\theta }_1$ refer to the perturbed solutions. Inserting (17) into (14) to get

$\frac{d{a}_1}{d\tau }=\frac{h}{4}\left(a_1\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_0+a_0{\theta }_1\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_0\right)+\frac{a_1}{2}\left(\alpha -G_1\right)-\frac{3a_1{a}_0^2}{8}\left(\beta +3G_2\right),2a_0\frac{d{\theta }_1}{d\tau }=\frac{h}{4}\left(a_1\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_0-a_0{\theta }_1\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_0\right)+\frac{a_1}{2}\sigma +\frac{9a_1{a}_0^2}{8}\gamma .$

(18)

Since $a_1$ and ${\theta }_1$ are unknown functions of amplitude and phase, respectively, we express them as a linear combination of $k_q{e}^{\lambda \tau }\left(q=1,2\right)$ , where $\lambda$ is the eigenvalue associated with the unknown perturbation and $k_q$ is the constant. When the steady-state solutions $a_0$ and ${\theta }_0$ are stable, the characteristic equation derived from (18) has roots that possess a negative real part

${\lambda }^2+{\Gamma }_1\lambda +{\Gamma }_2=0,$

(19)

where ${\Gamma }_q$ can be get easily as

${\Gamma }_1=\frac{2}{a_{10}}\left(\frac{1}{8}h{\mathit{sin}\hspace{0.17em}\theta }_{10}{a}_{10}-\frac{\alpha a_{10}}{4}+\frac{3\beta a_{10}^3}{16}+\frac{a_{10}{G}_1}{4}+\frac{9}{16}{a}_{10}^3{G}_2\right),$

${\Gamma }_2=\frac{1}{16}h[4{\mathit{sin}\hspace{0.17em}\theta }_{10}{G}_1-2h-4\sigma {\mathit{cos}\hspace{0.17em}\theta }_{10}-4\alpha {\mathit{sin}\hspace{0.17em}\theta }_{10}-a_{10}^2(9\gamma {\mathit{cos}\hspace{0.17em}\theta }_{10}-3\beta {\mathit{sin}\hspace{0.17em}\theta }_{10}-9{\mathit{sin}\hspace{0.17em}\theta }_{10}{G}_2)].$

(20)

As per the Routh–Hurwitz criterion, ⁴⁶ the requirements for the stability of the steady-state solutions are outlined as follows

${\Gamma }_1>0,{\Gamma }_2>0.$

(21)

Results and discussions

This section demonstrates the system’s behavior under resonance conditions. It includes a comparison of several controllers to determine the most suitable one for the system and examines the influence of specific system parameters on its amplitude.

Time history

Figure 2 illustrates the system’s normal operation, where it operates without resonance or control (i.e. $G_1=G_2=0$ ). We observe that the amplitude has a value of $0.1125$ . Alternatively, Figure 3(a) displays the time response for the amplitude $a$ , while Figure 3(b) illustrates the system’s phase plane under resonance conditions without any control system applied $\left(G_1=G_2=0\right)$ , using the parameter values $q=0.001,\gamma =0.1,\Omega =25,h=0.5,\alpha =0.005,\omega =1.2,$ and $\beta =0.05.$

Figure 2.

(a) The behavior of $u\left(\tau \right)$ in normal operation (no resonance or control system applied), and (b) the system phase plane.

Figure 3.

(a) Time response of $u\left(\tau \right)$ under normal conditions, and (b) the system phase plane.

We observe that the amplitude has a value of $0.8609$ , that is, it increased by 7.65 times the normal case, which indicates a threat to the system, so it is important to reduce and limit these harmful vibrations by using a suitable control to maintain the system’s safety. Figure 4 shows the time history with the Negative derivative feedback (NDF) controller and NVF + NCVF. We find that due to the action of the NDF control and NVF + NCVF, the amplitude is decreased to $0.05531$ and $0.01929$ , respectively. In this system, the NVF + NCVF controller is employed due to its superior performance in vibration reduction compared to the NVF controller alone. A comparison of the time responses for the NVF and NVF + NCVF controllers is presented in Figure 5. When the NVF + NCVF controller is applied, and the amplitude is reduced to $0.01929$ , whereas the NVF controller alone reduces the amplitude to $0.01991$ . It is worth noting that, in general, the control effect on vibration and tension is very strong, while the impact on tension is very small. This is beneficial as the primary goal of control is to minimize dynamic vibrations, allowing us to take advantage of the resulting tension. As shown in Figure 6, the effective gain values $G_1\ge 0.1$ significantly reduce system vibrations. Notice that the amplitude value decreases as $G_1$ increases.

Figure 4.

A comparison between the time response of NDF and NVF + NCVF controller.

Figure 5.

A comparison between the time response of NVF and NVF + NCVF controller.

Figure 6.

The effect of the controller gain $G_1$ of the composite NVF + NCVF controller on the amplitude $u\left(\tau \right)$ .

Table 1 shows a comparison between different types of controllers. This comparison helps us to identify the optimal controller for specific applications. Due to the influence of the NDF controller, which is represented by the main system

$\ddot{u}\left(t\right)+(1-h\hspace{0.17em}\mathit{cos}\left(\omega t\right))u\left(t\right)-\left(\alpha -\beta u^3\left(t\right)\right)\dot{u}\left(t\right)-\gamma u^3\left(t\right)=q{\Omega }^2\hspace{0.17em}\mathit{cos}\left(\Omega t\right)$

$-\frac{q}{2}{\Omega }^2{u}^2\left(t\right)\mathit{cos}\left(\Omega t\right)+G_1\dot{v}\left(t\right),\ddot{v}\left(t\right)+{\eta }_2\dot{v}\left(t\right)+{\omega }_2^{2}v\left(t\right)=-G_2\dot{u}\left(t\right).$

(22)

Table 1.

Comparison between different types of controllers.

Controller type	Mathematical representation	Controller gains	Amplitude	Controller effectiveness %
NDF (blue curve)	$\ddot{v}\left(t\right)+{\eta }_2\dot{v}\left(t\right)+{\omega }_2^{2}v\left(t\right)=-G_2\dot{u}\left(t\right)$	$G_1=G_2=10$	$0.05531$	93.57
NVF (green curve)	$-G_1\dot{u}\left(t\right)$	$G_1=3$	0.01991	97.68
NVF + NCVF (red curve)	$-G_1\dot{u}\left(t\right)-G_2{\dot{u}}^3\left(t\right)$	$G_1=3,G_1=1$	0.01929	97.75

It was found that the amplitude decreased by 93.57% due to the NDF control. Moreover, it is observed that this controller has a minimal efficiency on amplitude compared to the other discussed ones. The NVF controller is represented by $-G_1\dot{u}\left(t\right)$ , in which the amplitude is decreased by 97.68%. Finally, the composite controller NVF + NCVF that is represented by $-G_1\dot{u}\left(t\right)-G_2{\dot{u}}^3\left(t\right)$ in which the amplitude is decreased by 97.75% which has a potent effect on mitigating vibrations, so we chose it in this study.

Frequency response

This section aims to examine the stability of steady-state oscillations.^47,48 As a result, the present part employs the nonlinear stability analysis approach to verify the stability of the system under consideration. This involves evaluating the impact of the acted parameters and conducting a stability simulation. It is observed that several crucial characteristics, such as parametric excitation force $h,$ the nonlinear stiffness coefficient $\gamma$ and finally, the feedback gains $G_1$ and $G_2$ are critical in the stability processes. The stability graphs for the modulation equations of the system were created using a particular action and several settings of the system (18). The changed amplitudes that fluctuate over time for different parametrical zones are illustrated, together with their features over phase plane pathways. The system’s potential fixed points are depicted. The relationship between the detuning parameter and the possible fixed points is examined in Figures 7–10. On the other hand, the continuous lines indicate the stable range, while the dotted curve indicates the unstable one. The amplitude versus the detuning parameter $\sigma$ is plotted in Figures 7 and 8 for the resonance case, with $G_1=2,3,5$ and $G_2=1,3,5$ , respectively. We observe that the amplitude decreased with the increase in the value of $G_1$ and $G_2.$ Figures 9 and 10, it is shown that system amplitude is inversely proportional to varying $h$ and $\gamma$ , respectively. The drawn curves indicate that the unstable fixed points correspond to $-10\le \sigma \le -1.9$ , while the stable fixed points correspond to $-1.9\le {\sigma }_1\le 10$ . Conversely, the stability of the fixed point is reversed at ${\sigma }_1=-1.9$ . As per the simulation above, the system exhibits a single critical fixed point for all parameter values, with no apparent qualitative change in behavior when the parameters are altered. Additionally, Figures 7–10 demonstrate the efficacy of the employed control.

Figure 7.

Frequency response curve at $G_1=2,3,5$ .

Figure 8.

Frequency response curve at $G_2=1,3,5$ .

Figure 9.

Frequency response curve at $h=0.1,0.9,1.9$ .

Figure 10.

Frequency response curve at $\gamma =0.1,0.5,0.9.$

Comparisons with Numerical Techniques

Here, we compare the analytical solution induced by the MSS and the NS of equation (1) using the RK-4 method. In Figure 11, the system’s modes amplitudes are graphed with time in both the analytical and numerical methods without control and with control. This comparison shows a high degree of agreement between the two solutions. To validate this comparison, we perform an error analysis by computing the absolute errors between AS and NS. Furthermore, we present this verification in Table 2, which highlights the accuracy and convergence of the MSS method.

Figure 11.

A comparison between the AS and the NS: (a) without control and (b) with control.

Table 2.

Absolute errors between AS and NS.

Without control				With control
t	AS	NS	error	AS	NS	error
101.8898	0.8337	0.8491	0.0154	0.0138	0.01	0.0038
105.2236	0.79121	0.8163	0.02509	0.0178	0.01	0.0078
108.5575	0.7693	0.7807	0.0114	0.0197	0.0099	0.0098
111.8913	0.7291	0.7526	0.0235	0.0196	0.0099	0.0097
115.2252	0.7719	0.7433	0.0286	0.0174	0.0099	0.0075
120.1116	0.7743	0.785	0.0107	0.0135	0.0098	0.0037
126.6267	0.8967	0.9184	0.0217	0.0085	0.0097	0.0012
130.9951	0.9931	0.9857	0.0074	0.0028	0.0095	0.0067
136.4743	0.9929	0.9936	0.0007	−0.0033	0.0094	0.0127
144.2448	0.8562	0.899	0.0428	−0.009	0.0092	0.0182
149.2757	0.8375	0.8229	0.0146	0.0055	0.0085	0.003
154.7351	0.7451	0.7678	0.0227	0.0109	0.0078	0.0031
160.1945	0.78931	0.7864	0.00291	0.0153	0.0072	0.0081
169.2158	0.9473	0.9552	0.0079	0.0183	0.0001	0.0182
174.2892	1.0138	1.004	0.0098	0.0198	0.0001	0.0197
178.4215	0.9795	0.9884	0.0089	0.0193	0.0001	0.0192
182.5539	0.9463	0.9397	0.0066	0.0164	0	0.0164
186.6863	0.8854	0.8767	0.0087	0.0114	0	0.0114
188.7525	0.8581	0.8448	0.0133	0.0051	0	0.0051
191.5868	0.8126	0.805	0.0076	−0.0007	0	0.0007
194.421	0.7673	0.7759	0.0086	0.00012	0	0.00012
197.2553	0.7671	0.7643	0.0028	0.00011	0	0.00011
200.0896	0.7864	0.7755	0.0109	−0.0001	0	0.0001

The chaotic response

Bifurcation curves graphically depict how a dynamical system’s behavior evolves with changes in a parameter. They highlight points where qualitative shifts occur, such as the emergence of new equilibrium states or periodic orbits. Combined with PM, bifurcation curves offer valuable insights into the dynamics and stability of systems, aiding in the prediction and understanding of complex phenomena like chaos and stability transitions.

In the mentioned system, the control gains influence the system’s behavior, leading to the emergence of various types of motion. The bifurcation graphs of $u$ , versus $\alpha$ were presented to demonstrate the kinds of motions.^49–51 In this study, we first examine the chaotic behavior of the system without control, as illustrated in Figure 12, which presents the bifurcation diagram of $u$ . Based on the figures below, the uncontrolled system exhibits two distinct complex behaviors across different ranges of $\alpha$ , with each range representing a specific degree of chaotic motion. In the first range at $\alpha \in \left[0,0.01\right]$ , the motion appears to follow a chaotic sequence. On the other hand, the second range is found at $\alpha >0.01$ , the system exhibits hyper-chaotic behavior. We validated the previous conclusion about the uncontrolled system by simulating its phase portraits and Poincare’ maps at different values of $\alpha$ . Specifically, for $\alpha =0.002$ , the system exhibits chaotic behavior, while $\alpha =0.02$ demonstrates a hyper-chaotic state, as shown in Figures 13 and 14, respectively. Notably, as $\alpha$ increases, the system’s complexity grows, transitioning from chaotic behavior to a hyper-chaotic state. It is important to note that the blue curves represent the phase portraits, while the red dots correspond to the PM.

Figure 12.

Bifurcation diagram of $\alpha$ without control.

Figure 13.

The phase portrait and PM at $\alpha =0.002$ .

Figure 14.

The phase portrait and PM at $\alpha =0.02$ .

Conversely, we analyzed the controlled system’s behavior through bifurcation analysis, as presented in the bifurcation diagram in Figure 15, and further supported by the phase portraits and PM in Figures 16 and 17. Our findings demonstrate that the control device significantly enhances the system’s stability. As shown in Figure 15, the bifurcation diagram indicates that the system maintains stable behavior within the range $0<\alpha \le 0.5$ . This stability is further confirmed in the PM in Figure 16, where a single red point denotes the system’s periodic motion at $\alpha =0.2$ . However, for $\alpha >0.5$ , the system moves from a stable state into a chaotic state. This is evident in Figure 17, where the PM displays red points distributed randomly, signifying chaotic behavior. These results underscore the critical role of the control device in regulating the system’s dynamics and mitigating chaotic responses.

Figure 15.

Bifurcation diagram of $\alpha$ with control.

Figure 16.

The phase portrait and PM at $\alpha =0.2$ .

Figure 17.

The phase portrait and PM at $\alpha =0.6$ .

Conclusions

A VMDO comprised with the NVF and NCVF controllers has been studied. Controllers are used to minimize any possible harmful vibrations. MSS are used to obtain the AS until the second-order approximation. These results are plotted and have been certified numerically utilizing the RK-4 method. The modulation equations are obtained; in addition, all resonance cases are determined. The influence of the control devices on the chaotic dynamics of the system is analyzed, which helps the system to behave steadily before transitioning to chaotic behavior. We demonstrated this through bifurcation graphs, phase portraits, and PM. Furthermore, a comparison between multi-control strategies was performed and it was discovered that, because of the composite NVF and NCVF controllers, the amplitude decreased by $97.75\%$ that is the most used and efficient one. Therefore, it is crucial, as the primary objective of the controller is to reduce dynamic vibrations. The system being studied is of considerable importance due to its engineering applications. The VMDO plays a vital role in many disciplines due to its rich dynamics and ability to model complex behaviors. This paper discusses the chaotic behavior and control of the oscillations of the van der Poel-Matthew-Duffing oscillator under parametric excitation and resonance. Despite the valuable insights obtained, much research remains to be done in the areas of advanced stability analysis, higher-order nonlinear effects, innovative vibration control strategies, experimental validation, and multi-degree-of-freedom systems. From mechanical engineering to biology to chaos theory, its application areas continue to expand as researchers explore new ways to analyze nonlinear systems.