A review of the frequency-amplitude formula for nonlinear oscillators and its advancements

Abstract

The frequency-amplitude relationship is pivotal in understanding oscillatory systems, dictating their dynamic behaviors and responses. This paper delves into the intricacies of the frequency amplitude formula, elucidating its foundational role in both linear and nonlinear oscillations. Through comprehensive analysis, we highlight the formula’s significance in predicting system behaviors, especially in environments with varying amplitudes. Our findings underscore the formula’s potential as a robust tool for enhanced system characterization, offering profound implications for diverse applications across scientific and engineering domains. This study delves deep into the formulation of the frequency-amplitude relationship, extending its application beyond the conventional cubic powers of the restoring force. We introduce three novel and equivalent formulations of the generalized frequency amplitude, breaking traditional boundaries by not confining the restoring force to just odd functions. The formula to determine the frequency of the singular oscillator has been set forth. Our approach, characterized by its simplicity, offers a robust tool for analyzing high nonlinearity vibrations, ushering in a richer, multidimensional perspective to vibration analysis.

Introduction

Indeed, the study of N/MEMS (Nano/Microelectromechanical Systems) oscillations and their periodicity feature is crucial due to their extensive use in microinstruments and various other domains.^1–3 These systems possess unique qualities such as small size, batch fabrication, low power consumption, and excellent reliability, making them indispensable for manufacturing a wide range of microdevices. ⁴ The dynamic oscillatory behavior of N/MEMS is complex and often influenced by numerous nonlinearities.^5,6 Nonlinear phenomena, such as mode coupling or amplitude-dependent frequency shifts, can lead to significant deviations from perfect periodicity in N/MEMS oscillators. Combining techniques like the variational iteration method (VIM) with the Laplace transform provides effective approaches to investigating the periodicity of oscillatory problems in N/MEMS.^7–10 This coupling allows for the transformation of the nonlinear system into a more manageable form, yielding an approximate analytic solution with reduced computational effort and simplifying the analysis. It facilitates the exploration of frequency-amplitude correlations and aids in understanding the fundamental dynamics of N/MEMS oscillations. The small size, low power consumption, and high sensitivity of N/MEMS oscillators make them valuable in various applications, including communication systems, sensors, and timing circuits. The periodicity feature of N/MEMS oscillators is essential for their proper functioning and for providing accurate timing and sensing capabilities required in many applications. Optimizing the performance and reliability of N/MEMS oscillators in real-world systems necessitates an understanding of and the ability to control their periodic behavior. Persistent research in this domain will enhance our comprehension of N/MEMS dynamics and facilitate innovative implementations across many domains.

Dynamic regulations, commonly encountered in physics and engineering, are typically described by a consistent nonlinear differential equation, the exact solution for which remains indeterminate. Conventional strategies, like perturbation approaches or numerical methodologies, are used to determine the appropriate frequency-amplitude connection and predict their dynamic reactions, facilitating a comprehensive understanding of their quantitative and qualitative behaviors.^11–14

Oscillatory systems are ubiquitous in nature and engineering, manifesting in phenomena ranging from the rhythmic beating of a heart to the vibrations of a bridge. While linear oscillations can often be described with relative simplicity, real-world systems frequently exhibit nonlinear behaviors that defy such straightforward analysis. Nonlinear oscillations, characterized by their dependence on amplitude and the potential for complex, unpredictable behaviors, have long intrigued and challenged scientists and engineers.^15–18

Traditional methods of analyzing nonlinear oscillations have often relied on perturbation techniques. These methods involve introducing small “perturbations” to a known, simple solution and iteratively refining the solution based on these perturbations. While perturbation methods have proven invaluable in many contexts, they come with inherent limitations, especially when dealing with systems that exhibit strong nonlinearity or when the perturbations themselves aren’t small.^19–22

He’s frequency formula, introduced by Ji-Huan He, is an integral part of the homotopy perturbation method (HPM) framework. This formula utilizes a combination of the homotopy perturbation method and an averaging technique to estimate the frequency of nonlinear oscillators.^23–26 Known for its simplicity and efficiency, He’s formula excels in handling weakly nonlinear oscillators without requiring extensive computational resources. It provides an analytical expression that establishes a relationship between the frequency of oscillation and the amplitude, particularly useful in systems with moderate nonlinearity. Widely employed in engineering and physics, especially in mechanical systems like beam oscillations and electrical circuits, He’s frequency formula has demonstrated its effectiveness across a broad spectrum of challenges in these fields.^23–26 Notably, it has been particularly valuable in obtaining closed-form analytical solutions for oscillators, especially those of the Duffing type.^23–26 The Duffing frequency, conceived by He, has emerged as a powerful mathematical tool for investigating the periodic solutions inherent in nonlinear oscillators. Recognized for its clarity, simplicity, and efficacy, He’s frequency formula has attracted considerable attention since its introduction in a seminal review study.^27–29 Its lucidity and empirical validation quickly captured the interest of the engineering community.^30–32 Subsequent research efforts have been dedicated to refining and enhancing the precision of He’s frequency formula.^32–34 Moreover, its applicability has been extended to encompass various oscillator types, including fractal oscillators and damping oscillators, underscoring its versatility and importance in nonlinear dynamics research.^35–38

Enter the non-perturbative approach. This method offers a fresh perspective, sidestepping the iterative refinement of perturbation techniques. Instead of relying on small deviations from a known solution, the non-perturbative approach seeks to directly capture the inherent complexities of the nonlinear system. By not being bound to the constraints of small perturbations, this approach can provide insights into a broader range of system behaviors, especially those that might be overlooked or inaccessible through traditional perturbative methods. The present El-Dib’s formula focuses on strong nonlinearities even for quadratic nonlinearity in oscillatory systems. ³⁹ El-Dib’s formula is designed to tackle systems with strong or singular nonlinearities, offering a robust solution where traditional methods might fail. This formula is particularly noted for its ability to deal with high-degree nonlinearities, providing accurate frequency estimations in systems where the nonlinearity significantly impacts the system’s behavior. It also offers El-Dib’s formula for the successive solutions ⁴⁰ while He’s formula, may not provide accurate results for this purpose. El-Dib’s formula can be used to evaluate the damping coefficient of the nonlinear oscillator. Further, this formula can introduce nonlinear quadratic forces into the frequency formula of the oscillating system. ⁴¹ It often involves more complex mathematical manipulations compared to He’s method. El-Dib’s method finds its utility in more complex systems, such as high-energy physics, advanced engineering applications, and systems where the assumptions of weak nonlinearity are invalid.

Several unresolved issues are still in search of an optimal solution. This study aims to address the pressing challenge of incorporating quadratic nonlinearity without resorting to perturbation techniques and proposes an adjustment to derive the periodic remedy. The original frequency method developed by He^33,42 is modified by the Hamiltonian-based frequency formulation.^43–46 Viewing this formulation through the lens of energy dynamics introduces the Hamilton principle, shedding light on the energy conservation inherent in the vibrational process. Deriving their closed-form solutions is the main goal of this study because such nonlinear vibrations are common in many physical and technical situations.

The Galerkin and Rayleigh–Ritz methods are fundamental approaches in applied mathematics and engineering, designed to derive approximate solutions to differential equations. The Galerkin method relies on the concept of weighted residuals, with basis functions serving as weights to ensure the orthogonality of the residual to the space spanned by these functions. ⁴⁷ This method is widely applied across various numerical analyses, including finite element methods, to transform continuous problems into discrete systems amenable to numerical solutions. Conversely, the Rayleigh–Ritz method is anchored in the variational principle, where approximate solutions are attained by minimizing (or maximizing) an energy function. This approach is notably effective in structural mechanics for determining natural frequencies and mode shapes of structures.^48,49 Both methods, while distinct in their foundational principles, are closely intertwined in their application, offering robust frameworks for tackling complex problems in structural mechanics, fluid dynamics, and a broad spectrum of physical sciences. ⁵⁰ By converting continuous problems into discrete systems, these methods facilitate the numerical solution of equations that are otherwise challenging to solve analytically, underscoring their pivotal role in modern engineering and applied mathematics.

The extended Galerkin method represents a sophisticated evolution of the traditional Galerkin approach, which is integral to numerical analysis for solving differential equations. In its classic form, the Galerkin method approximates a solution by expressing it through a series of basis functions and then minimizes the error by projecting this error onto the same basis functions. This process aims to reduce the error to a weighted average minimum, culminating in a series of algebraic equations that serve as an approximation to the original differential equation.^51–53 The extended version of the Galerkin method is particularly useful in fields such as fluid dynamics, structural analysis, and electromagnetics. These domains often present challenges such as complex boundary conditions, nonlinear material properties, or the need to address phenomena that span multiple scales, which can stretch the capabilities of traditional numerical methods. The extended Galerkin method addresses these challenges by offering a more refined approach that enhances the accuracy, efficiency, and robustness of simulations. This method allows for more sophisticated modeling of physical systems, providing researchers and engineers with a powerful tool for simulating complex scenarios with greater fidelity.^54–56

The non-perturbative approach to nonlinear oscillations is not just a mathematical novelty; it holds the promise of deeper insights into the fundamental nature of complex systems. Whether one is examining the chaotic motion of a double pendulum, the intricate patterns of a nonlinear electrical circuit, or the complex dynamics of biological systems, the non-perturbative approach offers a robust and comprehensive tool for understanding and predicting system behavior.

In this exploration, we will delve into the intricacies of nonlinear oscillations through the lens of the non-perturbative approach, uncovering the rich tapestry of behaviors that these systems can exhibit and the profound insights that this method can offer.

Frequency formula of the nonlinear oscillation

The frequency-amplitude relationship is a fundamental aspect of nonlinear oscillators. This relationship describes how the natural frequency of the oscillator varies with its amplitude, which is a distinctive feature of nonlinear systems. In linear systems, the natural frequency remains constant regardless of the amplitude of oscillation. The frequency-amplitude relationship can be derived using various methods, including energy methods, perturbation techniques, non-perturbative approaches, and numerical simulations. The frequency-amplitude relationship in nonlinear oscillators is a rich and complex topic that has seen significant developments over the years. Understanding this relationship is crucial for various applications, from engineering design to biological systems, and continues to be an active area of research.

One of the major scientific questions in the theoretical analysis of nonlinear oscillation is the investigation of bursting oscillations caused by frequency-domain multiscale effects. Given its numerous applications in science and engineering, the Duffing oscillator (or Duffing Equation) is one of the most significant and traditional nonlinear ordinary differential equations. The motion of a body subjected to a nonlinear spring power causes a Duffing oscillator.

One can look for the Duffing type in this fashion

$\ddot{u}+F\left(u\right)=0,$

(1)

where the over dot represents the derivative concerning time t and $F\left(u\right)$ is composed of all odd functions in u. We deal with the initial conditions as

$u\left(0\right)=A\hspace{0.17em}\text{and}\hspace{0.17em}\dot{u}\left(0\right)=0.$

(2)

Equation (1) can be reformulated as follows

$\ddot{u}+\left(\frac{F\left(u\right)}{u}\right)u=0,$

(3)

where the proportion $\frac{F\left(u\right)}{u}$ signifies the stiffness attributed to the odd nonlinearity. In its most basic depiction, He ¹⁴ reformulated equation (3) as

$\ddot{u}+{\omega }^2\left(A\right)u=0,$

(4)

where A represents the amplitude and ω indicates the frequency of the nonlinear oscillation. It’s pertinent to highlight that equation (4) serves as a linear approximation of the nonlinear oscillator depicted in equation (1). The integration of Galerkin method coupled with the mean square technique ⁵⁷ has proficiently determined the frequency ${\omega }^2\left(A\right)$ . A plethora of analytical methodologies to gauge the frequency are documented in existing literature.^58–63 Adhering to the conventional approach of He’s formula,^24,27 Younesiana and colleagues ⁴³ explored the extended Duffing vibration characterized by $F\left(u\right)={\displaystyle \sum }_{n=0}^N{a}_{2n+1}{u}^{2n+1}$ .

Utilizing the weighted residuals approach, the frequency was determined. Using the energy balance technique, the same researchers ⁴³ derived an alternative frequency formula. Drawing from the foundational principles of the Energy Balance Method, Mahmoud Bayat and colleagues ⁶⁴ formulated a different version of the frequency equation. Molla and Alam, ⁶⁵ following the procedures outlined in Refs. ⁴² identified the relationship between frequency and amplitude, suggesting a method based on the energy balance technique to procure a higher-order approximate solution. This approach resulted in complex interconnected nonlinear algebraic equations. After these advancements, He and his team ⁴⁴ applied the Hamiltonian-based approach to discern the frequency attributes of the nonlinear oscillator, particularly for the cubic Duffing equation. Building on the insights from Younesiana et al. ⁴³ He et al. ⁴⁵ and Mahmoud Bayat et al. ⁶⁴ Hongjin Ma ⁴⁶ derived two distinct formulations for the frequency-amplitude correlation. He, in collaboration with Liu, ⁶⁶ introduced a modified frequency formula for nonlinear oscillators, specifically for fractal vibrations in porous media

${\omega }^2=\underset{0}{\overset{A}{{\displaystyle \int }}}{u}^{4}F\left(u\right)du/\underset{0}{\overset{A}{{\displaystyle \int }}}{u}^{5}du.$

(5)

He^23,67,68 introduced a straightforward framework for extracting the frequency from the nonlinear stiffness, as represented by

${\omega }^2\left(A\right)={\frac{dF\left(u\right)}{du}|}_{u=\overline{u}}.$

(6)

It’s essential to emphasize that the identified point $\overline{u}$ is determined based on the proposition of a tentative solution $U\left(t\right)$ that adheres to the stipulated initial conditions. It’s also worth noting that the criterion for the least mean square of the displacement is articulated as

${\overline{u}}^2=\frac{1}{2T}\underset{0}{\overset{T}{{\displaystyle \int }}}{U}^2\left(t\right)dt.$

(7)

Typically, the trial solution aligning with the initial conditions outlined in equation (1) is postulated as

$U\left(t\right)=A\hspace{0.17em}\cos \hspace{0.17em}\omega t.$

(8)

It is noted that if the initial conditions of the Duffing equation are such that u (0) $=0$ and $\dot{u}\left(0\right)=0$ , then the trial solution needs to be tailored to satisfy these conditions. One possible form of the trial solution in this case is $U\left(t\right)=A\left(1-\hspace{0.17em}\cos \hspace{0.17em}\omega t\right).$

The proposed solution (8) encompasses a singular unknown, specifically the frequency ω, which necessitates determination. Substituting (8) into (7) results in

${\overline{u}}^2=\frac{1}{4}{A}^2.$

(9)

By inserting (9) into (6), the frequency is ascertained as suggested by He.^33,67,68 Additionally, He ⁴² defined the frequency of the cubic oscillator as comparable to the function’s first derivative concerning the variable u evaluated at the designated point $\overline{u}=\frac{1}{2}A.$

${\omega }^2={\frac{F^{\prime }\left(u\right)}{u^{\prime }}|}_{u=\overline{u}=\frac{1}{2}A}=a_1+\frac{3}{4}{A}^2{a}_3,$

(10)

where the overline indicates the derivative concerning the variable u. This frequency is in complete harmony with outcomes obtained using perturbation techniques, as showcased in. ⁶⁹ Although He’s frequency-amplitude formulation is straightforward and apt for practical applications, especially for the cubic Duffing oscillator, it doesn’t perform as efficiently for higher-order Duffing oscillators, as will be further elucidated. To rectify this limitation and to explore a formulation better equipped to handle the presence of a robust component, El-Dib^70,71 introduced an optimized frequency formula, achieved by enhancing the stiffness $F\left(u\right)/u$ as described.

Recently, El-Dib^39,71 presented an improved formulation of the stiffness term for the nonlinear equation (3), as illustrated below

$\ddot{u}+\left(\frac{uF\left(u\right)}{u^2}\right)u=0.$

(11)

Then, independently integrate the numerator and denominator of the stiffens term related to the variable t from zero to T where T is the period. When compared to equation (4), the formula of the frequency will take the following shape

${\omega }^2=\frac{\underset{0}{\overset{T}{{\displaystyle \int }}}u\left(t\right)F\left(u\left(t\right)\right)\hspace{0.17em}dt}{\underset{0}{\overset{T}{{\displaystyle \int }}}{u}^2\left(t\right)dt}.$

(12)

The above formula can be established using the mean square error technology:

The following error was made while comparing equation (4) to equation (1)

$E\left({\omega }^2\right)=\left|F\left(u\right)-{\omega }^2\left(A\right)u\left(t\right)\right|.$

(13)

The mean square error is defined as

${\overline{E}}^2\left({\omega }^2\right)={\int }_0^T{\left(F\left(u\right)-{\omega }^2\left(A\right)u\left(t\right)\right)}^{2}dt;\hspace{0.17em}T=\frac{\pi }{2\omega }.$

(14)

The minimum value requires that

$\frac{d{\overline{E}}^2\left({\omega }^2\right)}{d\left({\omega }^2\right)}=\frac{d}{d\left({\omega }^2\right)}{\int }_0^T\left[{\omega }^4\left(A\right)u^2\left(t\right)-2{\omega }^2\left(A\right)u\left(t\right)F\left(u\right)+F{\left(u\right)}^2\right]dt=0.$

(15)

Since the function $F\left(u\right)$ is not related to the frequency ω the simplification leads to

${\omega }^2\left(A\right){\int }_0^T{u}^2\left(t\right)dt-{\int }_0^{T}u\left(t\right)F\left(u\right)dt=0.$

(16)

The result of solving equation (16) is formula (12). The frequency corresponding to the appropriate trial solution that matches the initial conditions can be obtained by performing the aforementioned integrals. The benefit of using this formulation is that it enables analysis of the frequency for successive approximations to the generalized Duffing oscillator. The formula presented above offers the most efficient and straightforward method for calculating frequency. The analysis procedure described above validates the frequency formula’s accuracy.

Enhanced He’s frequency formula

We shall explain how He’s frequency formula (6) can be included in El-Dib’s frequency formula (12) in the paragraphs that follow:

Let’s define the function $F\left(u\right)$ as an odd-power polynomial, which may be written as

$F\left(u\right)=a_{1}u+a_3{u}^3+a_5{u}^5+...a_{2n+1}{u}^{2n+1}={\displaystyle \sum }_{n=0}^N{a}_{2n+1}{u}^{2n+1},$

(17)

where a _j represents a constant coefficient and N refers to a general integer number.

Consequently, the generalized form of equation (1) is represented as

$\ddot{u}+{\displaystyle \sum }_{n=0}^N{a}_{2n+1}{u}^{2n+1}=0.$

(18)

The linear harmonic equation is discovered if n = 0 in equation (18), and its potential frequency calculated as

${\omega }^2={\frac{F\left(u\right)}{u}|}_{u=\overline{u}=\frac{1}{2}A}=a_1.$

(19)

In expansion (17) with n = 1, the cubic Duffing’s equation appears, where

$F\left(u\right)=a_{1}u+a_3{u}^3.$

(20)

The cubic Duffing’s frequency in this instance has been proved by He ³³ to correspond to the first derivative of the function $F\left(u\right)$ involving the variable u at the point of $\overline{u}=\frac{1}{2}A.$

${\omega }^2={\frac{F^{\prime }\left(u\right)}{u^{\prime }}|}_{u=\overline{u}=\frac{1}{2}A}=a_1+\frac{3}{4}{A}^2{a}_3,$

(21)

where the over dash denotes the derivative concerning the variable u.

Utilizing $n=2$ , the result of the cubic-quintic Duffing equation is

$F\left(u\right)=a_{1}u+a_3{u}^3+a_5{u}^5.$

(22)

To estimate the frequency ${\omega }^2$ , according to the first derivative concerned to the variable u results in

${\omega }^2={\frac{dF\left(u\right)}{du}|}_{u=\frac{1}{2}A}=a_1+\frac{3}{4}{A}^2{a}_3+\frac{5}{16}{A}^4{a}_5.$

(23)

This outcome perfectly matches He’s ⁶⁹ observations. It does not, however, precisely match the results obtained from the sophisticated Lindstedt–Poincaré solution. ⁶⁶ The frequency-amplitude connection is found to be as follows in Ref. ⁷²

${\omega }^2=a_1+\frac{3}{4}{A}^2{a}_3+\frac{5}{8}{A}^4{a}_5.$

(24)

Consequently, relying solely on the first derivative for $n>1$ does not ensure that the frequency aligns with that derived from a perturbation method. An examination of the aforementioned frequency indicates that it can be interpreted as

${\omega }^2={\frac{d}{du}\left(a_{1}u+a_3{u}^3\right)+\frac{1}{2}{a}_{5}u\frac{d^2}{d{u}^2}{u}^5|}_{u=\frac{1}{2}A}.$

(25)

Furthermore, $n=3$ the function $F\left(u\right)$ will incorporate the septic power, transforming it into

$F\left(u\right)=a_{1}u+a_3{u}^3+a_5{u}^5+a_7{u}^7.$

(26)

Once more, the frequency estimation, when derived through the derivative $F\left(u\right)$ concerning the variable u evaluated at the point $\frac{1}{2}A$ , is given as

${\omega }^2\left(A\right)={\frac{dF\left(u\right)}{du}|}_{u=\frac{1}{2}A}=a_1+\frac{3}{4}{A}^2{a}_3+\frac{5}{16}{A}^4{a}_5+\frac{7}{64}{a}_7{A}^6.$

(27)

This result is not in concordance with the findings derived from the Homotopy Perturbation Method (HPM) ⁶⁶ nor with those from the equivalent linearized method ⁷³

${\omega }^2\left(A\right)=a_1+\frac{3}{4}\hspace{0.17em}{A}^2{a}_3+\frac{5}{8}\hspace{0.17em}{A}^4{a}_5+\frac{35}{64}\hspace{0.17em}{a}_7{A}^6.$

(28)

As previously highlighted, the appropriate derivative methodology is

${\omega }^2\left(A\right)={\frac{d}{du}\left(a_{1}u+a_3{u}^3\right)+\frac{1}{2!}{a}_{5}u\frac{d^2}{d{u}^2}{u}^5+\frac{1}{3!}{a}_7{u}^2\frac{d^3}{d{u}^3}{u}^7|}_{u=\frac{1}{2}A}.$

(29)

For elevated values of n, the precise derivative formula to determine the frequency can be expressed as

$\begin{array}{l}{\omega }^2\left(A\right)=\underset{u\to \frac{1}{2}A\hspace{0.17em}}{\lim }\left(D_u\left(a_{1}u+a_3{u}^3\right)+\frac{1}{2!}{a}_{5}u{D}_u^2{u}^5+\frac{1}{3!}{a}_7{u}^2{D}_u^3{u}^7+\frac{1}{4!}{a}_9{u}^3{D}_u^4{u}^9+...+\frac{1}{r!}{a}_{2r+1}{u}^{r-1}{D}_u^r{u}^{2r+1}\right)\\ =\underset{u\to \frac{1}{2}A\hspace{0.17em}}{\lim }{\displaystyle \sum }_{n=0}^{\infty }\frac{1}{n!}{a}_{2n+1}{u}^{n-1}{D}_u^n{u}^{2n+1},\end{array}$

(30)

This is a development of He’s fundamental frequency formula, where $D_u^n$ represents nth order derivative concerning the variable u.

Frequency formula in terms of a power series form

When the infinite series (30), which results in, is expanded, the significance of this frequency may be appreciated

${\omega }^2\left(A\right)=a_1+\frac{3}{4}{a}_3{A}^2+\frac{5}{8}{a}_5{A}^4+\frac{35}{64}{a}_7{A}^6+\frac{63}{128}{a}_9{A}^8+\frac{231}{512}{a}_{11}{A}^{10}...$

(31)

The series given above can be formatted as a power series as

$\begin{array}{l}{\omega }^2\left(A\right)=a_1+\frac{3!}{1!\hspace{0.17em}2!}{a}_3{\left(\frac{A}{2}\right)}^2+\frac{5!}{2!\hspace{0.17em}3!}{a}_5{\left(\frac{A}{2}\right)}^4+\frac{7!}{3!\hspace{0.17em}4!}{a}_7{\left(\frac{A}{2}\right)}^6+\frac{9!}{4!\hspace{0.17em}5!}{a}_9{\left(\frac{A}{2}\right)}^8+\frac{11!}{5!\hspace{0.17em}6!}{a}_{11}{\left(\frac{A}{2}\right)}^{10}+...\\ ={\displaystyle \sum }_{n=0}^N{a}_{2n+1}\frac{\left(2n+1\right)\hspace{0.17em}!}{n!\times \left(n+1\right)\hspace{0.17em}!}{\left(\frac{A}{2}\right)}^{2n}.\end{array}$

(32)

This suggests that the interrelation between frequency and amplitude has been examined through the lens of a power series, consistent with the differential expression depicted in (30).

It’s important to note that the nonlinear equation (18) can be simplified to the form of a fundamental pendulum equation if the coefficient $a_{2n+1}$ is replaced with $\frac{{\left(-1\right)}^n}{\left(n+1\right)!}$ and the coefficient $b_{2n}$ is set to zero

$\ddot{u}+\sin \hspace{0.17em}u=0.$

(33)

Using (32) as a starting point, the corresponding linearized form is

$\ddot{u}+\frac{2}{A}{J}_1\left(A\right)u=0,$

(34)

This suggests that the frequency articulating the simple pendulum’s oscillation is as follows

${\omega }^2\left(A\right)=\frac{2}{A}{J}_1\left(A\right).$

(35)

This aligns with the findings previously reported in Ref. ²⁴ The following format can be used to indicate a rough solution

$u\left(t\right)=A\hspace{0.17em}\cos \sqrt{\frac{2}{A}{J}_1\left(A\right)}\hspace{0.17em}t.$

(36)

where $J_1\left(A\right)$ is the first type of the order one Bessel functions in the amplitude A.

The integrative performance of the frequency formula

Utilizing the expansion denoted by (31) is beneficial for calculating the frequency for higher values of n in the generalized $F\left(u\right)$ . Drawing from the trial solution (8) and the subsequent formula

$\frac{2}{\pi }\underset{0}{\overset{\pi /2}{{\displaystyle \int }}}{\cos }^{2m}\hspace{0.17em}\xi d\xi =\frac{2m-1}{2m},$

(37)

where m is an integer number. The integrative setup can be sought using the power series (31) as

${\omega }^2\left(A\right)=\frac{2}{A^{2}T}{a}_1\underset{0}{\overset{T}{{\displaystyle \int }}}{U}^{2}dt+a_3\frac{2}{A^{2}T}\underset{0}{\overset{T}{{\displaystyle \int }}}{U}^{4}dt+a_5\frac{2}{A^{2}T}\underset{0}{\overset{T}{{\displaystyle \int }}}{U}^{6}dt+a_7\frac{2}{A^{2}T}\underset{0}{\overset{T}{{\displaystyle \int }}}{U}^{8}dt+...$

$\begin{array}{l}{\omega }^2\left(A\right)=\frac{2}{\hspace{0.17em}{A}^2}{\displaystyle \sum }_{n=0}^N{a}_{2n+1}\frac{1}{T}\underset{0}{\overset{T}{{\displaystyle \int }}}{U}^{2n+2}dt=\frac{2}{\hspace{0.17em}{A}^{2}T}\underset{0}{\overset{T}{{\displaystyle \int }}}{\displaystyle \sum }_{n=0}^N{a}_{2n+1}{U}^{2n+2}dt\\ =\frac{2}{\hspace{0.17em}{A}^{2}T}\underset{0}{\overset{T}{{\displaystyle \int }}}UF\left(U\right)dt=\frac{1}{2{\overline{u}}^{2}T}\underset{0}{\overset{T}{{\displaystyle \int }}}UF\left(U\right)dt;\hspace{0.17em}T=\frac{\pi }{2\omega },\end{array}$

(38)

where the phrase ${\overline{u}}^2$ refers to the displacement with the least mean square^37,38 which is estimated in (7). Insert (7) into (38) yields the frequency formula in the final form

${\omega }^2\left(A\right)=\frac{\underset{0}{\overset{T}{{\displaystyle \int }}}UF\left(U\right)dt}{{\int }_0^T{U}^{2}dt};\hspace{0.17em}T=\frac{\pi }{2\omega },$

(39)

For a generic trial solution, the aforementioned formula can be represented in the format of (12).

To elucidate the precision of the solution using formula (12), one might examine the Duffing oscillator characterized by cubic nonlinearity

$\ddot{u}+{\omega }_0^{2}u+Q{u}^3=0;\hspace{0.17em}u\left(0\right)=A,\dot{u}\left(0\right)=0,$

(40)

where ω₀² is the natural frequency while Q refers to the Duffing coefficient.

Deriving the analytical solution can be expressed as

$u\left(t\right)=A\hspace{0.17em}\cos \hspace{0.17em}\omega t,$

(41)

where ω is determined by formula (12). To visualize this, the aforementioned analytical solution can be depicted and juxtaposed with the numerical solution of equation (40), as depicted in Figure 1.

Figure 1.

The visualization of the approximate solution to equation (40) for ${\omega }_0=1,Q=0.1$ and $A=1$ . The associated relative error for the frequency $\omega$ is given by a value of $0.006$ .

This comparison shows that the analytical solution has an excellent agreement with the numerical solution.

The Helmholtz oscillator

The Helmholtz oscillator has garnered significant attention from the scientific community, as evidenced by the extensive body of research on the topic. ⁷³ Recognized as a foundational principle in both physics and engineering, the Helmholtz oscillator has been the focal point of myriad scholarly inquiries. ⁷⁴ Over the past decade, the nonlinear Helmholtz equation has been a particular area of emphasis. The practice of translating real-world challenges into mathematical models serves as a potent instrument for comprehension. It is posited that solutions derived from these models can offer insights into the dynamics of nonlinear oscillations across various scientific and engineering domains, including but not limited to, electronics, naval applications, finance, and other models. ⁷⁵ This manuscript introduces an analytical method tailored for crafting the frequency formula pertinent to the Helmholtz oscillator. Helmholtz-type nonlinearities, characterized by quadratic characteristics, owe their nomenclature to Helmholtz’s pioneering assertion that the eardrum functions akin to a linear and quadratic geometric component-based restoring force to an asymmetric oscillator. ⁷⁶ Contemporary research, however, has shifted its focus towards non-perturbative approaches in the pursuit of optimal approximate solutions. ³⁷

Established the frequency formula for the quadratically nonlinear oscillator

When the function $F\left(u\right)$ is comprised solely of odd powers, the oscillation is identified as the Duffing oscillator. Conversely, if the oscillation consists only of even powers, it is recognized as the Helmholtz type. A nonlinear oscillator, where the restoring force is confined to being an even function, can be articulated as

$\ddot{u}+{\omega }_0^{2}u+g\left(u\right)=0,$

(42)

where $g\left(u\right)$ denotes the nonlinearity having a function of the quadratic format such that

$g\left(u\right)=b_2{u}^2+b_4{u}^4+...+b_{2n}{u}^{2n}={\displaystyle \sum }_{n=1}^N{b}_{2n}{u}^{2n},$

(43)

where the coefficient b _i is a constant coefficient.

Let’s assume that the initial condition remains valid. If we posit that equation (42) possesses a precise direct frequency formula ${\sigma }^2\left(A\right)$ , then it can be linearly expressed as

$\ddot{u}+{\sigma }^2\left(A\right)u=0,$

(44)

where

${\sigma }^2\left(A\right)={\omega }_0^2+{\frac{g\left(u\right)}{u}|}_{u\to \overline{u}}.$

(45)

It’s important to highlight that when the function $g\left(u\right)$ adopts a quadratic form, the frequency ${\sigma }^2\left(A\right)$ can be directly determined by substituting the displacement u with $\frac{1}{2}A$ . However, this holds only accurate to a quadratic power. For even greater powers, there’s a necessity to formulate a more comprehensive equation. About the three methodologies previously mentioned for deriving the frequency formula (namely, the derivative method, the power series method, and the integrative method), the frequency ${\sigma }^2\left(A\right)$ can be deduced in analogous manners.

Given that an even function can be derived from the differentiation of an odd function concerning the variable u, the even function $g\left(u\right)$ can be represented as

$g\left(u\right)=\frac{d}{du}F\left(u\right).$

(46)

Differentiating the power series (17) concerned with the variable u gives

$\frac{d}{du}F\left(u\right)=\frac{d}{du}{\displaystyle \sum }_{n=0}^N{a}_{2n+1}{u}^{2n+1}={\displaystyle \sum }_{n=0}^N\left(2n+1\right)a_{2n+1}{u}^{2n}.$

(47)

Replacing $a_{2n+1}$ into (6) by $\frac{b_{2n}}{2n+1}$ leads to the power series (43). Employing this transformation into the expansion (32) yields

${\sigma }^2\left(A\right)={\omega }_0^2+{\displaystyle \sum }_{n=1}^{\infty }{b}_{2n}\frac{2n\hspace{0.17em}!}{n!\times \left(n+1\right)\hspace{0.17em}!}{\left(\frac{A}{2}\right)}^{2n}={\omega }_0^2+\frac{1}{4}{b}_2{A}^2+\frac{1}{8}{b}_4{A}^4+\frac{5}{64}{b}_6{A}^6+\frac{7}{128}{b}_8{A}^8+...\hspace{0.17em}.$

(48)

It should be observed that the corresponding differentiative form $\sigma \left(A\right)$ can be expressed as

${\sigma }^2\left(A\right)={\omega }_0^2+\underset{u\to \frac{1}{2}A\hspace{0.17em}}{\lim }{\displaystyle \sum }_{m=1}^M\frac{1}{\left(m+1\right)!}{b}_{2m}{u}^m{D}_u^m{u}^{2m},$

(49)

where $M$ represents a generic integer. Expanding the above formularization, the expansion given in (7) will arise.

To derive the integrative form of the parameter ${\sigma }^2\left(A\right)$ , one can leverage the understanding that the odd function $F\left(u\right)$ can be obtained as the integration of the even function $g\left(u\right)$ , as illustrated below

$F\left(u\right)={\int }_0^u\left({\omega }_0^2+g\left(x\right)\right)dx.$

(50)

Employing (50) with the integrative frequency formula given in (39) yields

${\sigma }^2\left(A\right)=\underset{0}{\overset{T}{{\displaystyle \int }}}u\left(t\right)\left(\underset{0}{\overset{u\left(t\right)}{{\displaystyle \int }}}\left({\omega }_0^2+g\left(x\right)\right)dx\right)dt/{\int }_0^T{u}^2\left(t\right)dt.$

(51)

To illustrate this objective, the solution for the quadratically nonlinear oscillator, as represented by (42), can be expressed as follows

$u\left(t\right)=A\hspace{0.17em}\cos \hspace{0.17em}\sigma \left(A\right)t.$

(52)

The frequency $\sigma \left(A\right)$ is determined either by (49) when adopting the differentiative method or by (51) when using the integrative approach. It’s noteworthy that these formulae are congruent. Subsequently, a numerical implementation of the analytical solution for the Helmholtz oscillator, as denoted by (52), is warranted. Given this context, we can examine the subsequent Helmholtz oscillator

$\ddot{u}+{\omega }_0^{2}u+b_2{u}^2+b_4{u}^4+b_6{u}^6+...=0.$

(53)

The frequency $\sigma \left(A\right)$ according to the differentiate formula (49) or due to the integrative expansion (51) is

$\sigma \left(A\right)=\sqrt{{\omega }_0^2+\frac{1}{4}{b}_2{A}^2+\frac{1}{8}{b}_4{A}^4+\frac{5}{64}{b}_6{A}^6}...\hspace{0.17em}.$

(54)

To ascertain the accuracy of the suggested approach, a numerical solution is employed as a reference point. A side-by-side comparison with the numerical solution of equation (53) is undertaken to assess the solution given by equation (52). The outcomes are showcased in Figure 2. The exact numerical solution of equation (53) is represented by a solid red line, whereas the blue-dashed line signifies the analytical solution (52). This juxtaposition underscores the accuracy of the analytical solution provided by (52).

Figure 2.

The visualization of the approximate solution to equation (53) for ${\omega }_0=1,b_2=b_4=b_6=0.01$ and $A=1$ . The associated relative error for the frequency $\omega$ is $-0.001$ .

Frequency formula for the Helmholtz–Duffing oscillator

An oscillator where the restoring force is not confined solely to odd functions or exclusively to even powers, but encompasses both odd and even powers, is termed the Helmholtz–Duffing oscillator, as illustrated below

$\ddot{u}+F\left(u\right)+g\left(u\right)=0.$

(55)

The functions $F\left(u\right)$ and $g\left(u\right)$ , as mentioned earlier, stand for the odd and even functions, respectively. The construction of the overall frequency corresponding to the oscillator above includes both of these functions. This oscillator’s linearized equation can be expressed as follows

$\ddot{u}+\mathrm{\Omega }\left(A\right)u=0.$

(56)

If the trial solution of this oscillator is denoted as

$U\left(t\right)=A\hspace{0.17em}\cos \hspace{0.17em}\mathrm{\Omega }t,$

(57)

then the frequency Ω is structured from both frequency ${\omega }^2\left(A\right)$ and the frequency ${\sigma }^2\left(A\right)$ as follows

${\mathrm{\Omega }}^2\left(A\right)={\omega }^2\left(A\right)+{\sigma }^2\left(A\right)-{\omega }_0^2.$

(58)

where ${\omega }^2\left(A\right)$ and ${\sigma }^2\left(A\right)$ are given by the formula (12) and (51), respectively. At this end, the approximate solution is given by

$u\left(t\right)=A\hspace{0.17em}\cos \sqrt{{\omega }^2\left(A\right)+{\sigma }^2\left(A\right)-{\omega }_0^2}t.$

(59)

A straightforward methodology to derive the linearized version of an oscillator (55) can be articulated using El-Dib’s formula (12). This can be accomplished by reformatting equation (55) into the subsequent alternative expression ⁶¹

$\ddot{u}+\left(\frac{uF\left(u\right)}{u^2}\right)u+\frac{uG\left(u\right)}{u^2}=0.$

(60)

In this context, the function $G\left(u\right)=ug\left(u\right)$ is characterized as an odd nonlinear function. It’s important to note that the central term is employed to formulate the Duffing’s frequency, while the final segment is utilized to realize the non-homogeneous component. By incorporating the trial solution (8) into (60) and subsequently integrating both the numerator and the denominator of the nonlinear stiffness term, as well as the non-stiffness term concerning the variable t, we achieve integration over the interval from zero to T.

Accordingly, equation (60) will linearized in the form

$\ddot{u}+{\omega }^2\left(A\right)u=-{\varpi }^2\left(A\right),$

(61)

where

${\omega }^2\left(A\right)=\frac{{\int }_0^{T}uF\left(u\right)dt}{{\int }_0^T{u}^{2}dt}={\omega }_0^2+\frac{3}{4}{a}_3{A}^2+\frac{5}{8}{a}_5{A}^4+...={\omega }_0^2+{\displaystyle \sum }_{n=1}^N{a}_{2n+1}\frac{\left(2n+1\right)\hspace{0.17em}!}{n!\times \left(n+1\right)\hspace{0.17em}!}{\left(\frac{A}{2}\right)}^{2n},$

(62)

${\varpi }^2\left(A\right)=\frac{{\int }_0^{T}uG\left(u\right)dt}{{\int }_0^T{u}^{2}dt}=\frac{3}{4}{b}_2{A}^2+\frac{5}{8}{b}_4{A}^4+...={\displaystyle \sum }_{n=1}^N{b}_{2n}\frac{\left(2n+1\right)\hspace{0.17em}!}{n!\times \left(n+1\right)\hspace{0.17em}!}{\left(\frac{A}{2}\right)}^{2n}.$

(63)

which is the non-homogenous harmonic oscillator its solution can be expressed in the form

$u\left(t\right)=\left(A+\frac{{\varpi }^2}{{\omega }^2}\right)\cos \hspace{0.17em}\omega t-\frac{{\varpi }^2}{{\omega }^2}.$

(64)

Figure 3 showcases a comparison between the analytical solution of equation (55) and its numerical counterpart. Both the analytical solutions, denoted by (59) and (64), will be juxtaposed for evaluation. For this illustration, we will consider specific functions, accompanied by the initial conditions set out in (2)

$F\left(u\right)={\omega }_0^{2}u+Q{u}^3+P{u}^5,$

(65)

$g\left(u\right)=b_2{u}^2+b_4{u}^4+b_6{u}^6.$

(66)

Figure 3.

Comparison of the approximate solution $u\left(t\right)$ as given by (59) and the analytical solution given by (64) with numerical solution of equation (55) for a system of ${\omega }_0=1,Q=P=0.1$ , $b_2=b_4=b_6=0.01$ and $A=1$ .

The frequency corresponding to solution (59) is defined as

$\mathrm{\Omega }=\sqrt{{\omega }^2\left(A\right)+{\sigma }^2\left(A\right)-{\omega }_0^2}=\sqrt{{\omega }_0^2+\frac{3}{4}{A}^{2}Q+\frac{5}{8}{A}^{4}P+\frac{A^2{b}_2}{4}+\frac{A^4{b}_4}{8}+\frac{35A^6{b}_6}{64}}.$

(67)

The frequency associated with solution (64) and the value of the non-homogeneous component are, respectively, provided by

$\omega =\sqrt{{\omega }_0^2+\frac{3}{4}{A}^{2}Q+\frac{5}{8}{A}^{4}P},$

(68)

$\varpi =\sqrt{\frac{3}{4}{A}^2{b}_2+\frac{5A^2{b}_4}{8}+\frac{35A^6{b}_6}{64}}.$

(69)

The analytical solutions, as denoted by (59) and (64), are graphically represented alongside the numerical solution of equation (55) in Figure 3. The approximate solution, as outlined by (59), is depicted by a blue-dashed line, whereas the analytical solution (64) is shown with a green-dotted line. In contrast, the numerical solution corresponding to the Helmholtz–Duffing equation (55) is marked by a red solid line, specifically for a system characterized by parameters ${\omega }_0=1,Q=P=0.1$ , $b_2=b_4=b_6=0.01$ and $A=1$ . The computed relative error for the analytical solution (59) stands at 0.0407, while for the analytical solution (64), it is −0.0258. This juxtaposition suggests that the solution provided by (64) boasts greater accuracy compared to that of (59).

Application

The application focuses on a general vibratory system about N/MEMS, elucidating the theory outlined in the preceding formalism.

Consider the dynamics of microstructures, which can be depicted by a nonlinear ordinary differential equation. This equation embodies the typical form associated with a set of oscillators from N/MEMS,^76–79 commonly employed in the realms of nanoscience and nanotechnology

$\left({\alpha }_0+{\alpha }_{1}u+{\alpha }_2{u}^2+{\alpha }_3{u}^3+{\alpha }_4{u}^4\right)\ddot{u}+\left({\beta }_1+{\beta }_{2}u+{\beta }_3{u}^2+{\beta }_4{u}^3+{\beta }_5{u}^4+{\beta }_6{u}^5+{\beta }_7{u}^6\right)u+{\beta }_0=0$

(70)

Given the initial conditions: $u\left(0\right)=A,$ $\dot{u}\left(0\right)=0,$ where ${\alpha }_0,{\alpha }_1,...,{\alpha }_4$ and ${\beta }_0,{\beta }_1,...,{\beta }_7$ are constants derived from the transformation of a multivariable differential equation into an ordinary differential equation, utilizing the Galerkin method.

Equation (70) serves as an extension of the Helmholtz–Duffing oscillator, which is composed of a combination of both odd and even nonlinear terms. It can be expressed as

$\ddot{u}+F\left(u,\ddot{u}\right)+g\left(u,\ddot{u}\right)+\frac{{\beta }_0}{{\alpha }_0}=0,$

(71)

where $F\left(u,\ddot{u}\right)$ and $g\left(u,\ddot{u}\right)$ are odd and even functions, respectively, which are given by

$F\left(u,\ddot{u}\right)=\frac{1}{{\alpha }_0}\left(\left({\alpha }_2{u}^2+{\alpha }_4{u}^4\right)\ddot{u}+{\beta }_{1}u+{\beta }_3{u}^3+{\beta }_5{u}^5+{\beta }_7{u}^7\right),$

(72)

$g\left(u,\ddot{u}\right)=\frac{1}{{\alpha }_0}\left(\left({\alpha }_{1}u+{\alpha }_3{u}^3\right)\ddot{u}+{\beta }_2{u}^2+{\beta }_4{u}^4+{\beta }_6{u}^6\right).$

(73)

Equation (71) can be reformulated as ⁴⁸

$\ddot{u}+\left(\frac{uF\left(u,\ddot{u}\right)}{u^2}\right)u+\left(\frac{u^{2}g\left(u,\ddot{u}\right)}{u^2}\right)+\frac{{\beta }_0}{{\alpha }_0}=0.$

(74)

By integrating both the numerator and the denominator independently over the interval (0, T), we get

$\ddot{u}+\left(\frac{{\int }_0^{T}uF\left(u,\ddot{u}\right)dt}{{\int }_0^T{u}^{2}dt}\right)u+\left(\frac{{\int }_0^T{u}^{2}g\left(u,\ddot{u}\right)}{{\int }_0^T{u}^{2}dt}\right)+\frac{{\beta }_0}{{\alpha }_0}=0.$

(75)

By utilizing the trial solution given in (8) and proceeding with the integration, equation (75) can be reformulated as

$\ddot{u}+{\omega }^{2}u+{\varpi }^2+\frac{{\beta }_0}{{\alpha }_0}=0,$

(76)

where

${\omega }^2=\frac{{\int }_0^{T}UF\left(U,\ddot{U}\right)dt}{{\int }_0^T{U}^{2}dt}=\frac{1}{{\alpha }_0}\left(-\left(\frac{3A^2{\alpha }_2}{4}+\frac{5A^4{\alpha }_4}{8}\right){\omega }^2+{\beta }_1+\frac{3A^2{\beta }_3}{4}+\frac{5A^4{\beta }_5}{8}+\frac{35A^6{\beta }_7}{64}\right),$

(77)

${\varpi }^2=\frac{{\int }_0^T{u}^{2}g\left(u,\ddot{u}\right)}{{\int }_0^T{u}^{2}dt}=\frac{1}{{\alpha }_0}\left(-\left(\frac{3A^2{\alpha }_1}{4}+\frac{5A^4{\alpha }_3}{8}\right){\omega }^2+\frac{3A^2{\beta }_2}{4}+\frac{5A^4{\beta }_4}{8}+\frac{35A^6{\beta }_6}{64}\right).$

(78)

Solving equation (77) for ${\omega }^2$ gets

${\omega }^2=\frac{\left({\beta }_1+\frac{3A^2{\beta }_3}{4}+\frac{5A^4{\beta }_5}{8}+\frac{35A^6{\beta }_7}{64}\right)}{\left({\alpha }_0+\frac{3}{4}{A}^2{\alpha }_2+\frac{5}{8}{A}^4{\alpha }_4\right)}.$

(79)

The solution of the non-homogenous second-order harmonic equation (76) is given by

$u\left(t\right)=\left(A+\frac{1}{{\omega }^2}\left({\varpi }^2+\frac{{\beta }_0}{{\alpha }_0}\right)\right)\cos \hspace{0.17em}\omega t-\frac{1}{{\omega }^2}\left({\varpi }^2+\frac{{\beta }_0}{{\alpha }_0}\right).$

(80)

The solution is visually represented in Figure 4, juxtaposed with the numerical solution for equation (70). A scrutiny of Figure 4 reveals that the analytical solution aligns exceptionally well with the numerical counterpart, showing a relative error estimated to be 0.027.

Figure 4.

Comparison of the approximate solution $u\left(t\right)$ as given (80) with a numerical solution of equation (70) for a system of ${\alpha }_0=1,{\alpha }_1={\alpha }_2={\alpha }_3={\alpha }_4=.1,{\beta }_1=3,{\beta }_0=1,{\beta }_2={\beta }_3={\beta }_4={\beta }_5={\beta }_6={\beta }_7=.1,$ and $A=1$ .

The Van der Pol oscillator, augmented with extended damping nonlinearity

Van der Pol is celebrated for the eponymous equation that delineates a rudimentary self-oscillating triode circuit. ⁸⁰ This equation has been the cornerstone for a myriad of electro-technical apparatuses, where the triode is replaced with entities like an arc lamp, a glow discharge tube, or an electronic vacuum tube. When augmented with a driving term, this equation stands as a benchmark system manifesting chaotic dynamics. Moreover, it’s associated with a prominent class of oscillations termed relaxation oscillations. Here, “relaxation” refers to the discharging activity of a capacitor. Van der Pol played a pivotal role in championing this nomenclature, which is deeply rooted in the nascent stages of nonlinear oscillation theory. ⁸¹

Van der Pol has demonstrated a dissipative equation that exhibits prolonged spontaneous oscillations without the use of force. The equation describing the amplitude of an oscillating current produced by a triode can be streamlined as follows

$\ddot{u}+\mu \left(1-u^2\right)\dot{u}+u=0,$

where the damping coefficient μ is in the original technical problem required to be smaller than 1.

In the ensuing section, our focus will be on the Van der Pol oscillator, with a particular emphasis on damping forces that incorporate higher-order nonlinearity. The goal is to distill this oscillator down to a fundamental harmonic oscillator. Drawing from the knowledge accumulated in our previous deliberations, we aim to refine the quadratic order of nonlinear forces.

Suppose we have a Van der Pol oscillator given in the form

$\ddot{u}+\mu g\left(u\right)\dot{u}+{\omega }_0^{2}u=0,$

(81)

where μ represents the coefficient associated with the linear damping force and

$g\left(u\right)=1+b_2{u}^2+b_4{u}^4+...$

(82)

It’s important to highlight that the initial conditions, as denoted by (2), remain applicable. The primary objective is to linearize this oscillator. It’s observed that the damping force’s nonlinearity manifests in a quadratic form. The fundamental principle revolves around transforming squared powers into odd powers. This transformation can be achieved in one of two methods. Firstly, it’s recognized that odd powers can be derived by integrating the quadratic powers related to the variable u. Subsequently, one can utilize formula (51) to derive a constant damping term. The alternative approach involves expressing these powers as a rational fraction. This is done by multiplying the numerator with the variable u and then dividing by the same variable. Consequently, both the numerator and denominator evolve into odd functions. In such scenarios, the El-Dib integral formula, denoted by (12), becomes applicable. This formula yields a constant that symbolizes a function in terms of amplitude A. The result is a straightforward harmonic equation.

Accordingly, the first scope of the nonlinear oscillator (81) can be rewritten in the form

$\ddot{u}+\mu {\sigma }^2\left(A\right)\dot{u}+{\omega }_0^{2}u=0,$

(83)

where the constant function ${\sigma }^2\left(A\right)$ ³⁹ is given by

${\sigma }^2\left(A\right)=\frac{{\int }_0^{T}U\left({\int }_0^{U}g\left(x\right)dx\right)dt}{{\int }_0^T{U}^{2}dt}=1+\frac{1}{4}{b}_2{A}^2+\frac{1}{8}{b}_4{A}^4+\frac{5}{64}{b}_6{A}^6+...$

(84)

The damping linear oscillator (83) has the following solution

$u\left(t\right)=A{e}^{-\frac{1}{2}\mu {\sigma }^2\left(A\right)t}\left(\cos \left(\sqrt{{\omega }_0^2-\frac{1}{4}{\mu }^2{\sigma }^4\left(A\right)}\right)t+\frac{\mu {\sigma }^2\left(A\right)}{2\sqrt{{\omega }_0^2-\frac{1}{4}{\mu }^2{\sigma }^4\left(A\right)}}\sin \left(\sqrt{{\omega }_0^2-\frac{1}{4}{\mu }^2{\sigma }^4\left(A\right)}\right)t\right),$

(85)

Given the second scope, equation (81) can be reformulated as

$\ddot{u}+\mu \left(\frac{F\left(u\right)}{u}\right)\dot{u}+{\omega }_0^{2}u=0,$

(86)

where $F\left(u\right)=ug\left(u\right)$ is used. It is noted that equation (86) can be enhanced ³⁰ to becomes

$\ddot{u}+\mu \left(\frac{uF\left(u\right)}{u^2}\right)\dot{u}+{\omega }_0^{2}u=0.$

(87)

Employing the trial solution (8) and using the El-Dib formula (12) to the above equation which is reduced to

$\ddot{u}+\mu {\omega }^2\left(A\right)\dot{u}+{\omega }_0^{2}u=0,$

(88)

where

${\omega }^2\left(A\right)=\frac{\underset{0}{\overset{T}{{\displaystyle \int }}}U\left(t\right)F\left(U\left(t\right)\right)\hspace{0.17em}dt}{\underset{0}{\overset{T}{{\displaystyle \int }}}{U}^2\left(t\right)dt}=1+\frac{3}{4}{A}^2{b}_2+\frac{5}{8}{A}^4{b}_4+\frac{35}{64}{A}^6{b}_6+...$

(89)

As mentioned before, the solution of equation (88) can be written as

$u\left(t\right)=A{e}^{-\frac{1}{2}\mu {\omega }^2\left(A\right)t}\left(\cos \left(\sqrt{{\omega }_0^2-\frac{1}{4}{\mu }^2{\omega }^4\left(A\right)}\right)t+\frac{\mu {\omega }^2\left(A\right)}{2\sqrt{{\omega }_0^2-\frac{1}{4}{\mu }^2{\omega }^4\left(A\right)}}\sin \left(\sqrt{{\omega }_0^2-\frac{1}{4}{\mu }^2{\omega }^4\left(A\right)}\right)t\right),$

(90)

To shed light on the analytical solutions of equation (81), we undertook a numerical comparison of solutions (85) and (90) against the numerical solution of equation (81), as previously delineated. By consolidating the numerical representations of these three solutions into one unified format, we aimed to simplify the comparative analysis and pinpoint the most precise solutions. The graphical representation reveals a commendable alignment of both solutions with their numerical counterpart. Figure 5 Specifically, the relative error for solution (85) stands at 0.02601, while the error associated with solution (90) is measured at 0.06954.

Figure 5.

Comparison of the approximate solution $u\left(t\right)$ as given by (85) and the analytical solution given by (90) with the numerical solution of equation (81) for a system of ${\omega }_0=1.5,\mu =0.1$ , $b_2=b_4=b_6=0.1$ and $A=1$ .

with the numerical solution of equation (81) for a system of ${\omega }_0=1.5,\mu =0.1$ , $b_2=b_4=b_6=0.1$ and $A=1$ .

Fast and simple approach to establish the equivalent linearization oscillation

In the present section, our focus is on broadening the application of the El-Dib formula (12) to address oscillators exhibiting pronounced nonlinearity. The fundamental premise of El-Dib’s approach operates as follows:

Suppose we possess a function represented by a differential equation. This function exhibits nonlinearity of the variable u and includes both its second and first derivatives concerning the independent variable t, presented in the format

$F\left(u,\dot{u},\ddot{u}\right)=0.$

(91)

The target is to discern how to derive an equivalent linear differential equation from the nonlinear version, characterized as follows

$a\ddot{u}+b\dot{u}+cu=0,$

(92)

where the undetermined constant coefficients a, b, and c need to be ascertained. It’s widely recognized that - the aforementioned equation can be approached from the subsequent perspective

$\frac{dF}{d\ddot{u}}\ddot{u}+\frac{dF}{d\dot{u}}\dot{u}+\frac{dF}{du}u=0.$

(93)

where the coefficients $\frac{dF}{d\ddot{u}}$ , $\frac{dF}{d\dot{u}}$ and $\frac{dF}{du}$ are nonlinear functions in the variable u(t) and its traditional derivatives. It is better to recall equation (93) in the form

$f_2\left(u\right)\hspace{0.17em}\ddot{u}+f_1\left(u\right)\dot{u}+f_0\left(u\right)u=0,$

(94)

where

$f_2\left(u\right)=\frac{dF}{d\ddot{u}},f_1\left(u\right)=\frac{dF}{d\dot{u}}\hspace{0.17em}\text{and}\hspace{0.17em}{f}_0\left(u\right)=\frac{dF}{du}.$

(95)

The comparison of equations (92) and (94), leads to express the following errors

$E_2=\left|\left(f_2\left(u\right)-a\right)\ddot{u}\left(t\right)\right|,E_1=\left|\left(f_1\left(u\right)-b\right)\dot{u}\left(t\right)\right|\hspace{0.17em}\text{and}\hspace{0.17em}{E}_0=\left|\left(f_0\left(u\right)-c\right)u\left(t\right)\right|.$

(96)

The total mean square errors are convenient to use here, which are defined as

$\begin{array}{l}{\overline{E}}^2={\overline{E}}_2^2+{\overline{E}}_1^2+{\overline{E}}_0^2\\ =\underset{T\to \infty }{\lim }\frac{1}{2T}{\int }_0^T{\left(f_2\left(u\right)-a\right)}^2{\ddot{u}}^2\left(t\right)dt+\underset{T\to \infty }{\lim }\frac{1}{2T}{\int }_0^T{\left(f_1\left(u\right)-b\right)}^2{\dot{u}}^2\left(t\right)dt\\ +\underset{T\to \infty }{\lim }\frac{1}{2T}{\int }_0^T{\left(f_0\left(u\right)-c\right)}^2{u}^2\left(t\right)dt;\hspace{0.17em}T=\frac{\pi }{2\omega }.\end{array}$

(97)

The minimum value requires that

$\frac{d{\overline{E}}^2}{da}=\underset{T\to \infty }{\lim }\frac{1}{2T}\frac{d}{da}{\int }_0^T\left({\left(f_2\left(u\right)-a\right)}^2{\ddot{u}}^2\left(t\right)+{\left(f_1\left(u\right)-b\right)}^2{\dot{u}}^2\left(t\right)+{\left(f_0\left(u\right)-c\right)}^2{u}^2\left(t\right)\right)dt=0,$

(98)

$\frac{d{\overline{E}}^2}{db}=\underset{T\to \infty }{\lim }\frac{1}{2T}\frac{d}{db}{\int }_0^T\left({\left(f_2\left(u\right)-a\right)}^2{\ddot{u}}^2\left(t\right)+{\left(f_1\left(u\right)-b\right)}^2{\dot{u}}^2\left(t\right)+{\left(f_0\left(u\right)-c\right)}^2{u}^2\left(t\right)\right)dt=0,$

(99)

$\frac{d{\overline{E}}^2}{dc}=\underset{T\to \infty }{\lim }\frac{1}{2T}\frac{d}{dc}{\int }_0^T\left({\left(f_2\left(u\right)-a\right)}^2{\ddot{u}}^2\left(t\right)+{\left(f_1\left(u\right)-b\right)}^2{\dot{u}}^2\left(t\right)+{\left(f_0\left(u\right)-c\right)}^2{u}^2\left(t\right)\right)dt=0.$

(100)

Upon simplification, we arrive at the subsequent equations

$a{\int }_0^T{\ddot{u}}^{2}dt={\int }_0^T{\ddot{u}}^2{f}_2\left(u\right)dt,$

(101)

$b{\int }_0^T{\dot{u}}^{2}dt={\int }_0^T{\dot{u}}^2{f}_1\left(u\right)dt,$

(102)

$c{\int }_0^T{u}^{2}dt={\int }_0^T{u}^2{f}_0\left(u\right)dt.$

(103)

To derive the equivalent linear second-order differential equation using the conventional derivative, one must introduce the trial solution U that adheres to the initial conditions specified for the problem at hand. By aptly selecting the trial solution, the corresponding approximate solution can be achieved. Moreover, the aforementioned formulas can be employed to generate successive approximate solutions, as detailed in Ref. ⁸²

Let’s assume we possess a series of trial solutions for the nonlinear oscillation: U₀(t), U₁ (t), U₂ (t), …, and so forth. These are denoted as the zero-order, first-order, second-order, and subsequent orders of approximate solutions. For each instance, there are distinct coefficients corresponding to the order of the approximate solution. Consequently, the linear equation (92) transforms into

$a_i\ddot{u}+b_i\dot{u}+c_{i}u=0;\hspace{0.17em}i=0,1,2,...$

(104)

where $a_i,b_i$ and $c_i$ are

$a_i={\int }_0^T{\ddot{U}}_i^2{f}_2\left(U_i\right)dt/{\int }_0^T{\ddot{U}}_i^{2}dt,$

(105)

$b_i={\int }_0^T{\dot{U}}_i^2{f}_1\left(U_i\right)dt/{\int }_0^T{\dot{U}}_i^{2}dt,$

(106)

$c_i={\int }_0^T{U}_i^2{f}_0\left(U_i\right)dt/{\int }_0^T{U}_i^{2}dt.$

(107)

Assuming the partial derivatives provided in equation (95) are expressed in terms of the trial solution $U\left(t\right)$ , it can be represented as

$f_2\left(U\right)={\frac{dF}{d\ddot{u}}|}_{u=U},f_1\left(U\right)={\frac{dF}{d\dot{u}}|}_{u=U},f_0\left(U\right)={\frac{dF}{du}|}_{u=U}.$

(108)

Employing (108) into equation (94) which can be rearranged as

$\left(\frac{{\ddot{U}}^2{f}_2\left(U\right)}{{\ddot{U}}^2}\right)\hspace{0.17em}\ddot{u}+\left(\frac{{\dot{U}}^2{f}_1\left(U\right)\hspace{0.17em}}{{\dot{U}}^2}\right)\dot{u}+\left(\frac{U^2{f}_0\left(U\right)\hspace{0.17em}}{U^2}\right)u=0.$

(109)

By integrating both the numerator and the denominator for each coefficient of the aforementioned second-order equation across the domain $\left(0,T\right)$ , we obtain

$\left(\frac{{\int }_0^T{\ddot{U}}^2{f}_2\left(U\right)dt}{{\int }_0^T{\ddot{U}}^{2}dt}\right)\ddot{u}+\left(\frac{{\int }_0^T{\dot{U}}^2{f}_1\left(U\right)dt}{{\int }_0^T{\dot{U}}^{2}dt}\right)\dot{u}+\left(\frac{{\int }_0^T{U}^2{f}_0\left(U\right)dt}{{\int }_0^T{U}^{2}dt}\right)u=0.$

(110)

This can be rewritten in the familiar situation of the damping linear second-order equation as follows

$\ddot{u}+{\mu }_{eq}\dot{u}+{\omega }_{eq}^{2}u=0,$

(111)

where ${\mu }_{eq}$ denotes the equivalent damping coefficient and ${\omega }_{eq}^2$ signifies the equivalent conservative frequency. These coefficients are defined as

${\mu }_{eq}=\left(\frac{{\int }_0^T\left({\dot{U}}^2{f}_1\right)dt}{{\int }_0^T{\dot{U}}^{2}dt}\right)/\left(\frac{{\int }_0^T\left({\ddot{U}}^2{f}_2\right)dt}{{\int }_0^T{\ddot{U}}^{2}dt}\right),$

(112)

${\omega }_{eq}^2=\left(\frac{{\int }_0^T\left(U^2{f}_0\right)dt}{{\int }_0^T{U}^{2}dt}\right)/\left(\frac{{\int }_0^T\left({\ddot{U}}^2{f}_2\right)dt}{{\int }_0^T{\ddot{U}}^{2}dt}\right).$

(113)

The linear damping equation, represented by (111), possesses a solution in the format

$u\left(t\right)=e^{-\frac{1}{2}{\mu }_{eq}t}\left(A\hspace{0.17em}\cos \hspace{0.17em}\omega t+B\hspace{0.17em}\sin \hspace{0.17em}\omega t\right),$

(114)

where A and B are arbitrary constants determined by applying the appropriate initial conditions. The frequency $\omega$ is given by

${\omega }^2={\omega }_{eq}^2-\frac{1}{4}{\mu }_{eq}^2.$

(115)

Application

One can apply the formula presented in equation (110) to the example outlined in equation (70). It’s evident that this equation is formulated as

$F_1\left(u\right)\ddot{u}+F_0\left(u\right)u+g\left(u,\ddot{u}\right)+{\beta }_0=0,$

(116)

where $F_1\left(u\right),F_0\left(u\right),$ and $g\left(u,\ddot{u}\right)$ are even functions defined below

$F_1\left(u\right)={\alpha }_0+{\alpha }_2{u}^2+{\alpha }_4{u}^4,$

(117)

$F_0\left(u\right)={\beta }_1+{\beta }_3{u}^2+{\beta }_5{u}^4+{\beta }_7{u}^6,$

(118)

$g\left(u,\ddot{u}\right)=\left({\alpha }_{1}u+{\alpha }_3{u}^3\right)\ddot{u}+\left({\beta }_2{u}^2+{\beta }_4{u}^4+{\beta }_6{u}^6\right).$

(119)

Utilizing El-Dib’s enhancement, ⁴⁰ equation (116) can be restructured as

$\left(\frac{{\ddot{u}}^2{F}_1\left(u\right)}{{\ddot{u}}^2}\right)\ddot{u}+\left(\frac{u^2{F}_0\left(u\right)}{u^2}\right)u+\left(\frac{u^{2}g\left(u,\ddot{u}\right)}{u^2}\right)+{\beta }_0=0.$

(120)

Applying the integral form yields

$\left(\frac{{\int }_0^T{\ddot{u}}^2{F}_1\left(u\right)dt}{{\int }_0^T{\ddot{u}}^{2}dt}\right)\ddot{u}+\left(\frac{{\int }_0^T{u}^2{F}_0\left(u\right)dt}{{\int }_0^T{u}^{2}dt}\right)u+\left(\frac{{\int }_0^T{u}^{2}g\left(u,\ddot{u}\right)dt}{{\int }_0^T{u}^{2}dt}\right)+{\beta }_0=0.$

(121)

By utilizing the trial solution given in (8) and completing the integrations, the linearized version corresponding to the original equation (116) is derived as

$a_{eq}\ddot{u}+b_{eq}u+c_{eq}+{\beta }_0=0,$

(122)

where

$a_{eq}=\frac{{\int }_0^T{\ddot{U}}^2{F}_1\left(U\right)dt}{{\int }_0^T{\ddot{U}}^{2}dt}={\alpha }_0+\frac{3A^2{\alpha }_2}{4}+\frac{5A^4{\alpha }_4}{8},$

(123)

$b_{eq}=\frac{{\int }_0^T{U}^2{F}_0\left(U\right)dt}{{\int }_0^T{U}^{2}dt}={\beta }_1+\frac{3A^2{\beta }_3}{4}+\frac{5A^4{\beta }_5}{8}+\frac{35A^6{\beta }_7}{64},$

(124)

$c_{eq}=\frac{{\int }_0^T{U}^{2}g\left(U,\ddot{U}\right)dt}{{\int }_0^T{U}^{2}dt}=-\left(\frac{3}{4}{A}^2{\alpha }_1+\frac{5}{8}{A}^4{\alpha }_3\right){\omega }^2+\frac{3A^2{\beta }_2}{4}+\frac{5A^4{\beta }_4}{8}+\frac{35A^6{\beta }_6}{64}.$

(125)

By dividing each side of equation (122) by the non-zero coefficient $a_{eq}$ , the equation simplifies to a non-homogeneous simple harmonic form as

$\ddot{u}+{\omega }_{eq}^{2}u=-\frac{{\beta }_0+c_{eq}}{a_{eq}},$

(126)

where ${\omega }_{eq}^2=b_{eq}/a_{eq}$ represents the frequency associated with the subsequent solution

$u\left(t\right)=\left(A+\frac{{\beta }_0+c_{eq}}{b_{eq}}\right)\cos \hspace{0.17em}{\omega }_{eq}t-\frac{{\beta }_0+c_{eq}}{b_{eq}}.$

(127)

From (124) and (123), it’s evident that the frequency ${\omega }_{eq}^2$ aligns with what was previously derived in (79). Similarly, the coefficient $c_{eq}$ matches the value ${\varpi }^2$ presented in (78). However, a comparison between the linearized equation (76) and the linearized equation (126) reveals discrepancies in the inhomogeneity components. These variations result in only minor effects.

Implementation of the Toda oscillator

The Toda oscillator is a well-known nonlinear system that arises from the Toda lattice, a one-dimensional chain of particles with an exponential interaction potential. The Toda lattice and its associated oscillator have been studied extensively due to their integrable nature and the presence of soliton solutions. The Toda oscillator and lattice have applications in various fields, including condensed matter physics, nonlinear optics, and even biology. The integrable nature of the system makes it a valuable tool for studying more complex, non-integrable systems.

This example is dedicated to providing a concise overview of the Toda oscillator as referenced in Refs^66,83,84

$\ddot{u}+e^u-1=0;\hspace{0.17em}u\left(0\right)=A\hspace{0.17em}\&\hspace{0.17em}\dot{u}\left(0\right)=0.$

(128)

where $\ddot{u}$ is the acceleration, $\dot{u}\left(0\right)$ is the initial velocity at t = 0. For the Toda oscillator, the frequency-amplitude relationship can be derived using various methods, such as perturbation techniques or direct numerical simulations. Due to the exponential nonlinearity, the relationship can exhibit complex behaviors, especially at higher amplitudes.

It’s advantageous to bifurcate the nonlinearity into odd and even nonlinear forces. Substituting the exponential function in equation (80) with its counterparts, the hyperbolic functions $\cosh \hspace{0.17em}u$ and $\sinh \hspace{0.17em}u$ , results in

$\ddot{u}+\sinh \hspace{0.17em}u=1-\hspace{0.17em}\cosh \hspace{0.17em}u,$

(129)

Selecting

$F\left(u\right)=\sinh \hspace{0.17em}u\hspace{0.17em}\text{and}\hspace{0.17em}g\left(u\right)=\cosh u-1.$

(130)

Equation (129), when expressed in terms of the Taylor expansion, can be represented in the form of a power series. This utilizes expansion (32) for the odd function and expansion (48) for the even function.

The Toda oscillator will have the form if the coefficient $a_{2n+1}$ is exchange by $1/\left(2n+1\right)!$ and the coefficient $b_{2n}$ is reciprocity by $1/2n!$

$\ddot{u}\left(t\right)+\left({\displaystyle \sum }_{n=0}^{\infty }\frac{1}{n!\times \left(n+1\right)\hspace{0.17em}!}{\left(\frac{A}{2}\right)}^{2n}\right)u\left(t\right)+{\displaystyle \sum }_{n=0}^{\infty }\frac{1}{n!\times \left(n+1\right)\hspace{0.17em}!}{\left(\frac{A}{2}\right)}^{2n}-1=0.$

(131)

When the Bessel function is updated, equation (131) undergoes a modification

$\ddot{u}\left(t\right)+\frac{2}{A}{I}_1\left(A\right)u\left(t\right)+\frac{2}{A}{I}_1\left(A\right)-1=0.$

(132)

This is the equivalent linearized form of the Toda oscillator given by equation (128), which has the following solution

$u\left(t\right)=\left(A-\frac{A-2I_1\left(A\right)}{2I_1\left(A\right)}\right)\cos \sqrt{\frac{2}{A}{I}_1\left(A\right)}\hspace{0.17em}t+\frac{A-2I_1\left(A\right)}{2I_1\left(A\right)}.$

(133)

where $I_1\left(A\right)$ is the first kind of modified Bessel function concerning A.

El-Dib’s formula (12) can be used to derive the previously indicated solution. Equation (129) can be rewritten as

$\ddot{u}+\left(\frac{u\hspace{0.17em}\sinh \hspace{0.17em}u}{u^2}\right)u=\frac{u^2\left(1-\hspace{0.17em}\cosh \hspace{0.17em}u\right)}{u^2}.$

(134)

Utilizing the trial solution given by equation (8), integrate the denominators of the two fractions across the time interval (0, T), as well as the numerators. The resulting integral representation of the first kind of Bessel functions can be expressed as

$\frac{1}{\hspace{0.17em}\pi }\underset{0}{\overset{\pi }{{\displaystyle \int }}}\cos \hspace{0.17em}t\hspace{0.17em}\sinh \left(x\hspace{0.17em}\cos \hspace{0.17em}t\right)dt=I_1\left(x\right).$

(135)

Equation (134) transformed to the linearized form

$\ddot{u}+{\omega }^2\left(A\right)u={\varpi }^2\left(A\right),$

(136)

where the expressions ${\omega }^2\left(A\right){\varpi }^2\left(A\right)$ for frequency and inhomogeneity, respectively, manifest in the configuration

${\omega }^2\left(A\right)=\frac{\underset{0}{\overset{T}{{\displaystyle \int }}}U\hspace{0.17em}\sinh \left(U\right)dt}{\underset{0}{\overset{T}{{\displaystyle \int }}}{U}^{2}dt}=\frac{2}{A}{I}_1\left(A\right),$

(137)

${\varpi }^2\left(A\right)=\frac{\underset{0}{\overset{T}{{\displaystyle \int }}}U\left(t\right)(\underset{0}{\overset{U\left(t\right)}{{\displaystyle \int }}}\cosh \hspace{0.17em}x-1)dx)dt}{\underset{0}{\overset{T}{{\displaystyle \int }}}{U}^{2}dt}=\frac{\underset{0}{\overset{T}{{\displaystyle \int }}}U\hspace{0.17em}\left(\sinh \left(U\right)-U\right)dt}{\underset{0}{\overset{T}{{\displaystyle \int }}}{U}^{2}dt}=2I_0\left(A\right)-\frac{2}{A}{I}_1\left(A\right)-1.$

(138)

where the functions $I_0\left(A\right)$ and $I_1\left(A\right)$ are the modified Bessel function in A. Consequently, the Toda oscillator’s solution is represented by equation (133).

The Toda oscillator is a rich and complex nonlinear system that offers a plethora of interesting dynamics and behaviors. Its study can provide valuable insights into the general nature of nonlinear oscillatory systems and their associated phenomena. Figure 6 shows the numerical representation of Toda’s frequency formula. The comparison with the numerical solution demonstrates that the analytic periodic solution (133), represented by the blue-dot-curve, and the numerical solution, represented by the red-dash-curve, have a good agreement.

Figure 6.

Contrasting between the numerical solution of equation (128) (red-dash curve) and the analytical solution (133) (blue-dot curve), for a system having $A=0.1$ .

Implementing the Tangent (or hyperbolic tangent) oscillator

The Tangent (hyperbolic tangent) oscillator is a specific type of nonlinear oscillator that involves the tangent function in its formulation. To implement or analyze the Tangent oscillator, one would typically start by defining its differential equation and then proceed to solve it using various methods. For a small amplitude of oscillation, you might be able to use perturbation methods to find an approximate solution. For larger amplitudes or more complex scenarios, numerical methods method might be more appropriate. Here’s a general analytical non-perturbative approach to implementing the Tangent oscillator:

The oscillator characterized by the tangent function is presented as^85,86

$\ddot{u}+{\omega }_0^2\hspace{0.17em}\tan \hspace{0.17em}u=0;\hspace{0.17em}u\left(0\right)=A,\dot{u}\left(0\right)=0,$

(139)

where $\ddot{u}$ is the second derivative of u concerning time (acceleration) and ${\omega }_0^2$ is a constant related to the natural frequency of the oscillator.

It makes sense to assume that the weighted residuals methodology and the frequency formula of He, which are two different approaches that can be used to determine the frequency formula for the nonlinear oscillator, would yield consistent results. The frequency formula of He is specifically used in ⁸⁶ for the odd function and is described as follows

${\omega }^2={\frac{\tan \hspace{0.17em}u}{u}|}_{u=\frac{\sqrt{3}}{2}A}.$

(140)

As was previously established, the frequency extracted from the aforementioned peculiar non-polynomial function is not accurate. It is well-known that if the function isn’t expanded, the implementation of Galerkin’s integral formula is impractical. The methods presented in this study can be used to derive a more accurate frequency formula. When the function is replaced with the hyperbolic tangent, the oscillator takes on the subsequent form, as detailed in reference ⁸⁷

$\ddot{u}+{\omega }_0^2\hspace{0.17em}\tanh \hspace{0.17em}u=0;\hspace{0.17em}u\left(0\right)=A,\dot{u}\left(0\right)=0,$

(141)

Equation (141) can be reformulated as

$\ddot{u}+{\omega }_0^2\left(\frac{\sinh \hspace{0.17em}u}{\cosh \hspace{0.17em}u}\right)=0.$

(142)

Considering the fraction in the preceding equation, which features an odd numerator and an even denominator, one can expect this equation to be transformed into the subsequent linear form

$\ddot{u}+\left(\frac{{\omega }^2\left(A\right)}{{\varpi }^2\left(A\right)}\right)u=0.$

(143)

To evaluate both the numerator and the denominator of equation (143), one can restructure equation (142) in the following manner

$\ddot{u}+\left({\omega }_0^2\left(\frac{u\hspace{0.17em}\sinh \hspace{0.17em}u}{u^2}\right)u/\frac{u^2\hspace{0.17em}\cosh \hspace{0.17em}u}{u^2}\right)=0.$

(144)

Applying the integrative approach for the numerator as well as the denominator should be treated separately, as follows

$\ddot{u}+\left(\frac{{\omega }_0^2\underset{0}{\overset{T}{{\displaystyle \int }}}u\hspace{0.17em}\sinh \hspace{0.17em}u\hspace{0.17em}dt/\underset{0}{\overset{T}{{\displaystyle \int }}}{u}^{2}dt}{\underset{0}{\overset{T}{{\displaystyle \int }}}{u}^2\hspace{0.17em}\cosh \hspace{0.17em}u\hspace{0.17em}dt/\underset{0}{\overset{T}{{\displaystyle \int }}}{u}^{2}dt}\right)u=0.$

(145)

In terms of the trial solution (8) and El-Dib’s Formula (12), equation (143) will be derived where

${\omega }^2\left(A\right)={\omega }_0^2\underset{0}{\overset{T}{{\displaystyle \int }}}U\hspace{0.17em}\sinh \hspace{0.17em}U\hspace{0.17em}dt/\underset{0}{\overset{T}{{\displaystyle \int }}}{U}^{2}dt={\omega }_0^2{I}_1\left(A\right),$

(146)

${\varpi }^2\left(A\right)=\underset{0}{\overset{T}{{\displaystyle \int }}}{U}^2\hspace{0.17em}\cosh \hspace{0.17em}U\hspace{0.17em}dt/\underset{0}{\overset{T}{{\displaystyle \int }}}{U}^{2}dt=I_1\left(A\right)+A{I}_2\left(A\right).$

(147)

Accordingly, equation (145) becomes

$\ddot{u}+\left(\frac{{\omega }_0^2{I}_1\left(A\right)}{I_1\left(A\right)+A{I}_2\left(A\right)}\right)u=0.$

(148)

Its solution seeking in the form

$u\left(t\right)=A\hspace{0.17em}\cos \left(\sqrt{\frac{{\omega }_0^2{I}_1\left(A\right)}{I_1\left(A\right)+A{I}_2\left(A\right)}}\right)t.$

(149)

Should the function $\tanh \hspace{0.17em}u$ be replaced with $\tan \hspace{0.17em}u$ , the resultant solution, as previously indicated, takes on the ensuing configuration

$u\left(t\right)=A\hspace{0.17em}\cos \left(\sqrt{\frac{{\omega }_0^2{J}_1\left(A\right)}{J_1\left(A\right)-A{J}_2\left(A\right)}}\right)t,$

(150)

Once we have a solution, we can analyze it to understand the behavior of the Tangent oscillator. This might involve plotting u(t) versus time. Fig.(7) illustrates a comparison between the analytical solutions outlined in equation (149) and the precise numerical solution corresponding to equation (141). This visualization spans various amplitude values to facilitate a comprehensive comparison. Significantly, the solution delineated by equation (149) aligns impressively with the numerical solution across the entire amplitude spectrum.

Figure 7.

Juxtaposition the numerical solution of equation (141) (represented by the red dashed curve) with the solution given by equation (149) (depicted by the blue dotted curve), for ${\omega }_0=1,A=0.1,0.4,0.7,1,1.5.$

Relative restoring forces oscillators

Oscillators with relative restoring forces are a subset of nonlinear oscillators where the restoring force is not just a function of the displacement but also depends on a fractional composed of two parts of a numerator function and another denominator function in variables of the system. The behavior of oscillators with relative restoring forces can be quite complex due to the interplay between the various terms. They can exhibit a range of behaviors, from simple periodic motion to chaotic dynamics.

This section explores nonlinear oscillators characterized by an even denominator present within the odd numerator of the restoring force. The literature has presented various restoring forces characterized by an even denominator. For an oscillator with a specific relative restoring force, it can be expressed as

$F\left(u\right)=\frac{f\left(u\right)}{g\left(u\right)}=\frac{{\displaystyle \sum }_{n=0}^N{a}_{2n+1}{u}^{2n+1}}{{\displaystyle \sum }_{m=1}^M{b}_{2m}{u}^{2m}},$

(151)

where the values of $N$ and $M$ are integers. It is convenient to utilize the adapted He’s formula (30) to produce the frequency from (151). This approach aids in transforming the nonlinear structure into a more manageable linear format. Here, $f\left(u\right)$ is an odd function in polynomial form, and $g\left(u\right)$ is an even function. If this function is raised to an even power, one can refine it by converting the polynomial to an odd form, achieved by multiplying both the numerator and denominator by the variable u. At this juncture, formula (30) ought to be applied independently to the denominator and the numerator.

Utilizing the expanded He’s frequency formula (30), the frequency can be determined in the following manner

${\omega }^2=\frac{F\left(u\right)}{u}=\frac{f\left(u\right)}{ug\left(u\right)}=\underset{u\to \overline{u}}{\lim }\frac{{\displaystyle \sum }_{n=0}^N\frac{1}{n!}{a}_{2n+1}{u}^{n-1}{D}_u^n{u}^{2n+1}}{{\displaystyle \sum }_{m=1}^M\frac{1}{m!}{b}_{2m}{u}^{m-1}{D}_u^m{u}^{2m+1}},$

(152)

where ${\overline{u}}^2$ represents the mean square of the displacement u, which is determined based on the selection of an appropriate trial solution U(t)

${\overline{u}}^2=\frac{1}{2T}{\int }_0^T{U}^{2}dt.$

(153)

Application to vibration mode of a tapered beam

The vibration of tapered beams is a topic of interest in structural dynamics and engineering. Tapered beams, unlike uniform beams, have a non-constant cross-sectional area along their length. This variation in cross-sectional area introduces complexities in the vibration behavior of the beam. The following mathematical expression describes the principal vibrational formula of a tapered beam as mentioned in Ref ⁸⁸

$\ddot{u}+u+\alpha u^3+\beta \left(u^2\ddot{u}+u{\dot{u}}^2\right)=0,$

(154)

where $\ddot{u}$ is the acceleration, $\alpha$ is a coefficient related to the linear stiffness of the beam, and $\beta$ is a coefficient associated with the nonlinear stiffness due to the tapering of the beam. Assuming the ensuing preliminary circumstances

$u\left(0\right)=A,\dot{u}\left(0\right)=0.$

(155)

To linearize the above nonlinear equation, one can employ the modified He’s formula (30). By doing so, the nonlinear terms can be approximated, and the equation can be transformed into a linear form. This linearized equation can then be solved using standard techniques to obtain the vibration mode shapes and frequencies of the tapered beam. Upon applying the modified He’s formula (30), the vibration characteristics of the tapered beam can be determined. ⁸⁹ These characteristics, such as the natural frequencies and mode shapes, are crucial for design and analysis purposes. ⁹⁰

It is possible to rewrite this oscillator as

$\ddot{u}+f\left(u\right)=0,$

(156)

where the configuration $f\left(u\right)$ is constructed as

$f\left(u\right)=\frac{u+\alpha u^3+\beta u{\dot{u}}^2}{1+\beta u^2}.$

(157)

Since the odd numerator is present, the frequency can be roughly expressed as follows

${\omega }^2=\underset{\begin{array}{l}x\to \frac{1}{2}A,\\ \dot{x}\to \frac{1}{2}A\omega \end{array}}{\lim }\frac{f\left(u\right)}{u}=\underset{\begin{array}{l}u\to \frac{1}{2}A,\\ \dot{u}\to \frac{1}{2}A\omega \end{array}}{\lim }\frac{\frac{d}{du}\left(u+\alpha u^3+\beta u{\dot{u}}^2\right)}{\frac{d}{du}\left(u+\beta u^3\right)}=\frac{1+\frac{3}{4}{A}^2\alpha +\frac{1}{4}{A}^2\beta {\omega }^2}{1+\frac{3}{4}{A}^2\beta }.$

(158)

As a result, the frequency is expressed in the format

${\omega }^2=\frac{1+\frac{3}{4}{A}^2\alpha }{1+\frac{1}{2}{A}^2\beta }.$

(159)

We arrive at the same frequency formula by using Galerkin’s approach and the trial solution (8), which is presented as

$\frac{2}{\hspace{0.17em}{A}^{2}T}\underset{0}{\overset{T}{{\displaystyle \int }}}\left(\ddot{u}\left(1+\beta u^2\right)+u+\alpha u^3+\beta u{\dot{u}}^2\right)\cos \hspace{0.17em}\omega tdt=1+\frac{3}{4}{A}^2\alpha -{\omega }^2\left(1+\frac{1}{2}{A}^2\beta \right)=0.$

(160)

The rough answer to equation (154) is as follows

$u\left(t\right)=A\hspace{0.17em}\cos \hspace{0.17em}\omega t.$

(161)

It’s worth noting that the results obtained using He’s formula should be validated against experimental data or numerical simulations to ensure accuracy. ⁹¹ In conclusion, the application of He’s formula (30) to the vibration analysis of tapered beams provides a simplified approach to understanding their dynamic behavior. This method offers a quick and efficient way to analyze such complex structures without delving deep into intricate numerical simulations.

The free vibrations with type fifth-order nonlinearities ⁹²

The fifth-order nonlinearities introduce additional complexities to the vibration behavior of the oscillator. Such nonlinearities can arise in various physical systems, especially when the amplitude of vibration becomes large.

Mathematical representation for the differential equation representing the free vibrations of a conservative oscillator with fifth-order nonlinearities is expressed by

$\ddot{x}+{\omega }_0^{2}x+a_1{x}^3+a_2\left(x^2\ddot{x}+x{\dot{x}}^2\right)+b_1{x}^5+b_2\left(x^4\ddot{x}+2x^3{\dot{x}}^2\right)=0,$

(162)

where $\ddot{x}$ is the acceleration of the oscillator (162), ${\omega }_0^2$ is a coefficient representing the linear stiffness of the system, $a_1,a_2$ represents the cubic nonlinearity, $b_1$ and $b_2$ are coefficient associated with the fifth-order nonlinearity. The starting circumstances as defined in (155). The presence of these nonlinearities indicates that the system exhibits nonlinear behavior, especially at large amplitudes. These terms can lead to phenomena such as softening or hardening behavior, bifurcations, and other complex dynamical responses.

To analyze the system, one could employ techniques like the modified He’s formula (30), which has been discussed earlier. This would involve approximating the fifth-order term to transform the equation into a more tractable form. Once linearized, standard techniques can be used to determine the system’s response. Solving the linearized equation will provide insights into the vibration characteristics of the oscillator.

Equation (162) can be rearranged as follows

$\ddot{x}\left(1+a_2{x}^2+b_2{x}^4\right)+\left({\omega }_0^{2}x+a_1{x}^3+b_1{x}^5\right)+\left(a_{2}x+2b_2{x}^3\right){\dot{x}}^2=0.$

(163)

This leads to the following formulation of the regaining force $f\left(x,\dot{x}\right)$ as follows

$f\left(x,\dot{x}\right)=\frac{\left({\omega }_0^{2}x+a_1{x}^3+b_1{x}^5\right)+\left(a_{2}x+2b_2{x}^3\right){\dot{x}}^2}{1+a_2{x}^2+b_2{x}^4}.$

(164)

Hence, the frequency is determined to be

$\begin{array}{l}{\omega }^2=\underset{\begin{array}{l}x\to \frac{1}{2}A,\\ \dot{x}\to \frac{1}{2}A\omega \end{array}}{\lim }\frac{f\left(x,\dot{x}\right)}{x}=\underset{\begin{array}{l}u\to \frac{1}{2}A,\\ \dot{u}\to \frac{1}{2}A\omega \end{array}}{\lim }\left(\frac{D_x\left({\omega }_0^{2}x+a_1{x}^3\right)+\frac{1}{2!}u{D}_x^2{b}_1{x}^5+{\dot{x}}^2{D}_x\left(a_{2}x+2b_2{x}^3\right)}{D_x\left(x+a_2{x}^3\right)+\frac{1}{2!}u{D}_x^2{b}_2{x}^5}\right)\\ =\frac{{\omega }_0^2+\frac{3}{4}{A}^2{a}_1+\frac{5}{8}{A}^4{b}_1+\frac{1}{4}{A}^2{\omega }^2\left(a_2+\frac{3}{2}{b}_2{A}^2\right)}{1+\frac{3}{4}{A}^2{a}_2+\frac{5}{8}{A}^4{b}_2}.\end{array}$

(165)

Upon further simplification, the frequency is deduced as

${\omega }^2=\frac{{\omega }_0^2+\frac{3}{4}{A}^2{a}_1+\frac{5}{8}{A}^4{b}_1}{\left(1+\frac{1}{2}{A}^2{a}_2+\frac{2}{8}{A}^4{b}_2\right)}.$

(166)

By opting for the suggested solution (8) as previously indicated, the frequency derived from the Galerkin formula takes the shape

$\frac{2}{\hspace{0.17em}{A}^{2}T}\underset{0}{\overset{T}{{\displaystyle \int }}}\left(\ddot{x}\left(1+a_2{x}^2+b_2{x}^4\right)+\left({\omega }_0^{2}x+a_1{x}^3+b_1{x}^5\right)+\left(a_{2}x+2b_2{x}^3\right){\dot{x}}^2\right)\cos \hspace{0.17em}\omega tdt=0,$

(167)

which yields

${\omega }^2=\frac{{\omega }_0^2+\frac{3}{4}{A}^2{a}_1+\frac{5}{8}{A}^4{b}_1}{\left(1+\frac{1}{2}{A}^2{a}_2+\frac{3}{8}{A}^4{b}_2\right)}.$

(168)

However, it’s essential to remember that the linearized solution is an approximation that becomes significant

$x\left(t\right)=A\hspace{0.17em}\cos \hspace{0.17em}\omega t.$

(169)

In conclusion, while fifth-order nonlinearities introduce complexities, they also offer a richer dynamical behavior. Understanding these systems is crucial in various engineering applications where such nonlinearities might be present.

The simplest way to find nonlinear singular oscillator frequency

The frequency of nonlinear singular oscillators is of importance in the field of nonlinear dynamics. ⁹³ Singular oscillator equations of motion contain terms that, under certain conditions, become unbounded, or infinite. ⁹³ Determining the frequency of such oscillators is challenging due to their singularities. ⁹⁴ Since singularities might lead to infinite reactions, several traditional methods are not suitable. The frequency may vary in tandem with the amplitude, especially in highly nonlinear systems.

Finding the nonlinear singular oscillator’s frequency is the main goal of this section. We shall look at the single oscillator’s frequency-amplitude connection. Connecting the frequency and amplitude of nonlinear oscillators is one way to use this technology. It may be necessary to modify the conventional method for singular oscillators in order to accommodate the singularities.

For the special attention to the oscillator

Although there are various approaches to finding the frequency of nonlinear singular oscillators, the best approach for a given system often depends on its particular features, such as its singularity and nonlinearity. Examine the system below that exhibits strong singularity as

$\ddot{u}+\frac{1}{u}=0;\hspace{0.17em}u\left(0\right)=A,\dot{u}\left(0\right)=0.$

(170)

Applying He’s frequency formula directly yields the following frequency calculation

${\omega }^2={f^{\prime }|}_{u=\frac{1}{2}A}=-\hspace{0.17em}{\frac{1}{u^2}|}_{u=\frac{1}{2}A}=-\frac{4}{A^2}.$

(171)

This result makes ω imaginary and shows that there is no oscillation. The distinctive development of the restorative force $f\left(u\right)=\frac{1}{u}$ can be rephrased as follows rather than being applied directly

$f\left(u\right)=\frac{u}{u^2}.$

(172)

One interesting word for the restorative energy shown in (172) is the enhancement of the original function. The frequency is determined directly from the original potential function; the enhancement technique is solely needed to produce an odd numerator and an even denominator. Taking (172) as a ballpark estimate for the frequency, we arrive at

${\omega }^2={\frac{f\left(u\right)}{u}|}_{u=\frac{1}{2}A}.$

(173)

Combining (172) and (173) yields a valid formula with two odd polynomials in the numerator and denominator. As such, He’s differentiative technique will be applied separately to the numerator and denominator. The frequency is then determined to be

${\omega }^2={\frac{D_{u}u}{D_u{u}^3}|}_{u=\frac{1}{2}A}=\frac{4}{3A^2}.$

(174)

Equation (170) can be rewritten to adapt the integrative frequency formula to the particular requirements and characteristics of the nonlinear system being studied

$\ddot{u}+\left(\frac{u^2}{u^4}\right)u=0.$

(175)

With the trial solution given by (8), the frequency can be ascertained

${\omega }^2=\underset{0}{\overset{T}{{\displaystyle \int }}}{u}^2\hspace{0.17em}dt/\underset{0}{\overset{T}{{\displaystyle \int }}}{u}^4\hspace{0.17em}dt=\frac{1}{\frac{3}{4}{A}^2}.$

(176)

Interestingly, the outcomes of this approach agree with the results of using He’s formula. This congruence demonstrates the effectiveness of the proposed method in solving special issues in nonlinear oscillatory systems. In addition to confirming the methodology, the consistency of the results confirms El-Dib’s frequency formula (12) as a trustworthy analytical tool. This type of connection is significant because it demonstrates how multiple analytical methods can converge to provide comparable outcomes. This boosts confidence in the theoretical understanding and practical implementation of these strategies in the handling of difficult dynamic problems.

Remark: for the anti-cubic nonlinearity

$\ddot{u}+\frac{1}{u^3}=0;\hspace{0.17em}u\left(0\right)=A,\dot{u}\left(0\right)=0.$

(177)

By making the restorative force stronger $f\left(u\right)=\frac{u^3}{u^6}$ , the frequency is explored as

${\omega }^2=\frac{{\int }_0^T{u}^{2}dt}{{\int }_0^T{u}^{6}dt}=\frac{8}{5A^4}.$

(178)

In the general singular case

$\ddot{u}+\frac{1}{u^{2n+1}}=0;\hspace{0.17em}u\left(0\right)=A,\dot{u}\left(0\right)=0.$

(179)

The stronger restorative force is $\left(u\right)=\frac{u^3}{u^{2n+4}}$ ; consequently, one may estimate the frequency as

${\omega }^2=\frac{{\int }_0^T{u}^{2}dt}{{\int }_0^T{u}^{2n+4}dt}.$

(180)

Examine the singular oscillator composed of anti-cubic and anti-quintic equations that follow

$\ddot{u}+\frac{1}{u}+\frac{1}{u^3}+\frac{1}{u^5}=0;\hspace{0.17em}u\left(0\right)=A,\dot{u}\left(0\right)=0.$

(181)

It is convenient to rewrite equation (181) in the form

$\ddot{u}+\left(\frac{1+u^2+u^4}{u^5}\right)=0.$

(182)

The improved restorative force has the following format

$\ddot{u}+\left(\frac{\left(1+u^2+u^4\right)u^2}{u^8}\right)u=0.$

(183)

Thus, the frequency is executed as

${\omega }^2=\frac{{\int }_0^T\left(1+u^2+u^4\right)u^{2}dt}{{\int }_0^T{u}^{8}dt}=\frac{1+\frac{3}{4}{A}^2+\frac{5}{8}{A}^4}{\frac{35}{64}{A}^6}.$

(184)

Hence, the following is a representation of the approximate solution

$u\left(t\right)=A\hspace{0.17em}\cos \sqrt{\frac{64+48A^2+40A^4}{35A^6}}t.$

(185)

Examine this Duffing using an Anti-Cubic Oscillator

$\ddot{u}+u+u^3+\frac{1}{u^3}=0;\hspace{0.17em}u\left(0\right)=A,\dot{u}\left(0\right)=0.$

(186)

There are two distinct parts to the restorative force

$f\left(u\right)=f_1\left(u\right)+f_2\left(u\right),$

(187)

where $f_2\left(u\right)$ represents the anti-cubic function and is selected to symbolize the promise of the original Duffing $f_1\left(u\right)$ . Thus, a simpler version of equation (186) can be expressed

$\ddot{u}+\left(\frac{1+u^4+u^6}{u^3}\right)=0.$

(188)

Enhancing the restoring force as

$\ddot{u}+\left(\frac{\left(1+u^4+u^6\right)u^2}{u^6}\right)u=0.$

(189)

The frequency can therefore be calculated as

${\omega }^2=\frac{{\int }_0^T\left(1+u^4+u^6\right)u^{2}dt}{{\int }_0^T{u}^{6}dt}=\frac{8}{5A^4}+\frac{6}{5A^2}+\frac{7A^2}{8}.$

(190)

The following form of the solution is identified

$u=A\hspace{0.17em}\cos \left(\sqrt{\frac{8}{5A^4}+\frac{6}{5A^2}+\frac{7A^2}{8}}\right)t.$

(191)

Conclusion

This study presents refined methodologies adept at formulating periodic solutions for the generalized Duffing and Helmholtz oscillators, as well as their related variants. The nonlinear facets of the oscillator are categorized into distinct odd and even components. Both these components play pivotal roles in shaping the frequency. Notably, the even components necessitate integration concerning the system variable for frequency determination. Collectively, they can be synthesized to craft a non-homogeneous linear harmonic equation. The oscillator’s frequency is discerned from the secular function, which embodies the odd component. In contrast, the non-secular function, representing the even component, aids in elucidating the non-homogeneity of the harmonic oscillator. This strategy culminates in three unique nonlinear frequency styles and their corresponding inhomogeneity components. The frequency, in its differential guise, is showcased as an extension of He’s formula. Simultaneously, a frequency rendition in the vein of a power series is elucidated. Moreover, a frequency interpretation in an integrative mold is also delineated. This technique finds applicability in the Toda, Tangent, and Van der Pol oscillators. The proposition’s scope is further expanded to encapsulate the singular oscillator. In the realms of engineering and science, it’s commonplace to encounter even functions, quadratic forms, or singularities. This research augments the guidance on the utilization of frequency formulas for all-encompassing nonlinear oscillators. It’s imperative to underscore that the El-Dib frequency integrative formula, encapsulated in equation (12), epitomizes the most streamlined articulation of the amplitude-frequency relationship. As the saying goes, “In nonlinear vibration theory, elegance lies in simplicity.” This manuscript proffers a potentially streamlined methodology for the expeditious evaluation of a nonlinear oscillator’s frequency attributes. The outcome of this singular approach is laudably precise, underscoring its commendable accuracy.

Future directions in this research area will focus on several key aspects. Firstly, there will be a continued emphasis on refining and enhancing the present approach, particularly in terms of the advanced trial solution. This will involve developing more accurate approximations and incorporating additional factors to improve the precision of the analytical solutions obtained. Furthermore, there will be a concerted effort to address the discussion of the solution in the presence of the periodic force. This entails exploring the effects of the periodic force on the stability and dynamics of the system in greater detail. Additionally, there will be a focus on understanding how various parameters, such as the amplitude and frequency of the periodic force, influence the behavior of the system, and its response to external stimuli. Overall, future research endeavors will seek to deepen our understanding of nonlinear oscillatory systems and their responses to external forces. By refining analytical techniques and exploring the effects of different factors on system dynamics, researchers aim to provide valuable insights into the behavior of these complex systems and pave the way for advancements in various fields of science and engineering.