Non-perturbative analysis of weakly nonlinear oscillators: Stability,chaos,and bifurcation investigations

Abstract

The nonlinear oscillations frequently occur in real-world challenges, such as mechanical structures and biological cycles, and can give rise to complex behaviour like bifurcations and chaos that are vital in understanding and predicting dynamical systems. Advances in analytical and numerical methods have enabled novel applications in fields including engineering, medicine, and materials science. We analyse oscillators with pronounced nonlinear features, incorporating both damping and restoring forces, using a blend of theoretical and computational approaches. Two examples are presented from diverse scientific and technical areas. The innovative methodology described significantly reduces computation time and resources when compared to conventional perturbation methods widely used in this field. Based on He’s frequency formula, the proposed non-perturbative approach transforms weakly nonlinear oscillator of ordinary differential equation into linear one, allowing us to establish a new frequency corresponding to the linearized one. The theoretical outcomes are validated via numerical simulations using Mathematica Software, with results revealing excellent regularity between the two ordinary differential equations. An inclusive investigation of system stability can be conducted with this approach, expanding the capabilities beyond those available with previous techniques. Consequently, the non-perturbative approach offers a more practical and reliable framework via numerical solutions of weakly nonlinear oscillators. Moreover, it stands out as a flexible tool for use in applied research and engineering due to its flexibility to a variety of nonlinear situations. We also examine the impact of different parameters on stability, with results approving the approach’s simplicity, efficiency, and reliability. Adjusting bifurcation parameters alters the structure of bifurcation diagrams and their associated Poincaré maps, and we further illustrate these effects by mapping the Lyapunov exponent curves.

Keywords

nonlinear systems he’s formulation formula non-perturbative approach stability analysis differential equations chaotic motion 65P40 70k20 70Q05 74H45

1. Introduction

Nonlinear oscillation is an increasingly noticeable and attractive topic in mathematical and engineering disciplines. In the prototypes of physics, biology, and chemistry, nonlinear ODEs in progressive form constitute a natural characteristic. They were employed to address problems of dynamical systems, encompassing networking mathematics and optimization, in addition to engineering challenges related to materials, energy, and electrical structure.^1,2 An assembly of analytical and numerical methods is formulated and employed to address these issues, with extensive research on their examination in the literature. Therefore, employing the coupled homotopy-variational approach to the nonlinear DO revealed novel frequency-amplitude connections.³ Numerous technical applications employed actual springs linked in succession of flexible components designed to impersonator springs. The system’s linear or nonlinear performance is determined by the operational range of its components. The relevant literature elucidated the methodology of ascertaining the equivalent spring of linearly connected springs in series. A system including a mass supported by a linear spring and a nonlinear spring was reported.⁴ The current methodology markedly diverges from prior studies by introducing a novel framework that amalgamates adaptive modelling, real-time data, and iterative validation into a cohesive strategy, thereby enhancing accuracy, flexibility, and providing profound insights beyond traditional methods.

The Aboodh integral transform utilizing HPM was employed to derive an approximate analytical solution of classified issue.⁵ A nonlinear oscillator was augmented with two substantial linear components.⁶ An approximate response is calculated through a power series procedure. This concise statement functions as a framework for various applications expected in determining the period of nonlinear oscillators, as HPM has demonstrated its significance as a mathematical tool in understanding nonlinear oscillators.⁷ The Hamiltonian method was an analytical approximation technique employed to analyse extremely nonlinear dynamical systems.⁸ This system encompasses free oscillatory motion of a conservative autonomous oscillator exhibiting inertial and static fifth-order nonlinear effects, as well as motion of a non-deformable rod oscillating along a circular surface without slipping. A unique approach integrating HPM through the use of a variational formula was devised to get an accurate solution derived by analysis of a nonlinear ODE characterized by inertia and static nonlinearity.⁹ The produced findings are contrasted with alternative analytical and exact results to demonstrate the unexpected level of exactitude and generality of the approximation analytical technique. The obtained first-order approximation is shown to be virtually equivalent and possesses exact responses. A revised approach of global error minimization was proposed to address nonlinear third-degree Jerk ODE by concentrating on the determination of unknown parameters up to the third order.¹⁰ The nonlinear fractional Drinfeld–Sokolov system was investigated using the conformable fractional derivative.¹¹ The fractional traveling wave transformation and Galilean transformation were utilized to change the nonlinear fractional Drinfeld–Sokolov system into a planar dynamical system, followed by an investigation of its bifurcation and chaos, illustrated through 3D and 2D phase portraits. The fractional Yajima–Oikawa equation utilizing a conformable fractional derivative was investigated.¹² The fractional Yajima–Oikawa equation serves as a significant model in physics, utilized to characterize ion sound waves influenced by Langmuir action waves. The nonlinear fractional Schrödinger–Hirota equation is a significant model in physics, utilized to elucidate the propagation of optical solitons in optical fibers. The variational concept of the nonlinear fractional Schrödinger–Hirota equation was effectively formulated utilizing the semi-inverse method.¹³ The fractional Kaup–Newell system was examined, serving as a significant model in optical fibers.¹⁴ New soliton solutions were successfully generated utilizing two advanced mathematical methodologies and a novel fractional auxiliary equation method. The intricate and dynamic nature of epidemic modeling, especially with the inclusion of stochastic factors, presents a significant challenge in developing accurate and effective numerical methods for solving ODEs. The endeavor by presenting a novel finite difference approach for both linear and nonlinear stochastic and deterministic ODEs was documented.¹⁵ A numerical method for addressing stochastic time-dependent partial differential equations was proposed.¹⁶ This offers the benefit of addressing issues with affirmative resolutions. The scheme establishes criteria of acquiring positive answers, which the current Euler–Maruyama approach is incapable of achieving. A stochastic finite difference method, which was an explicit scheme, was proposed.¹⁷ The method can be utilized to discretize temporal variables in the examined stochastic parabolic equations. Nonlinear stochastic modeling is crucial in fields like psychology, finance, physical sciences, engineering, econometrics, and biological sciences. Dynamical consistency, positivity, and boundedness are essential characteristics of stochastic modelling. A stochastic model of coronavirus was analysed using transition probability approaches and parametric perturbation methods.¹⁸ Nonlinear stochastic modelling is significant in various scientific disciplines, including environmental science, materials science, engineering, chemistry, physics, and biomedical engineering, among others. The computational dynamics of the stochastic dengue model utilizing actual model data was examined.¹⁹

Recent advancements were made in asymptotic approaches of both weakly and highly nonlinear oscillators.²⁰ Three assessments of nonlinear oscillators were performed utilizing fundamental techniques: max-min methodology, the HFF, and HPM.²¹ The weighted average is employed to enhance the accuracy of its frequency-domain representation in a mathematical model. It was proposed that severely nonlinear oscillators could be completely and effortlessly controlled.²² The findings showed that the strategy produced responses with a high degree of accuracy. Both an experimental micro-electro-mechanical system and a packaging system demonstrated the use of the frequency-amplitude correlation in nonlinear vibration systems, as documented in prior research.²³ A direct frequency prediction method is introduced of nonlinear oscillators with impulsive ICs. Recent contributions were presented in the domains of Dynamical Systems in alignment with NPA were examined by Moatimid et al.^24–27 and El-Dib.^28–31 The following details must emphasize that the innovative method produces a distinct linear ODE that corresponds to the existing nonlinear one, along with the following advantages:

1. An innovative approach is employed to derive a new linear ODE that is equivalent to the original nonlinear one.

2. In light of the new method, there is a perfect connection between these two ODEs.

3. The shortcomings of earlier techniques are successfully addressed by this strategy.

4. Unlike other perturbation approaches, NPA allows, simply, experts to assess the problem’s stability examination.

5. One of the methods presents itself as a simple, effective, and compelling tool.

6. The NPA can be extended to include various combinations of interconnected dynamical systems regarded as meaningful, efficient, and convincing.

7. Taylor expansion assists in calculating these restoring forces in all perturbation methods, including the core technique known as the MSTM. However, NPA effectively away from this limitation.

8. The NPA employs a distinctive strategy in handling restoring forces, setting it apart from conventional perturbation methods and classifying it outside the perturbation technique framework.

The NPA possesses three primary limitations. The following is a summary of these restrictions:

i. It pertains to a weakly oscillatory second-order nonlinear ODE.

ii. The initial conditions remain unaltered.

iii. The initial amplitude must be less than one to get greater precision.

The subsequent sections of the paper will delineate the methodology of our inquiry in greater depth. § 2 examines two unique issues as highly nonlinear ODEs, in which the solutions are obtained, and the chaotic behaviors of the examined models are discussed through the BDs, PMs, and LEs. Moreover, the achieved results and debates on stability analysis are also calculated. § 3 presents the ultimate principal key outcomes.

2. Applications

This section aims to examine two specific instances of highly nonlinear systems utilizing NPA, as previously referenced by Moatimid et al.^24–27 and El-Dib.^28–31

2.1 example 1

The first example can be formulated as a smooth, frictionless joint A moves along a circular path of radius r. Attached to it is the arm AB, whose end B pulls a nonlinear spring connected to a mass P as shown in Figure 1. The path of point P is a fixed horizontal straight line that intersects circular path perpendicularly at point D. The main dynamic equation of block P of mass m, restricted to move along a straight line, can be derived from this setup. Assuming that the spring is subjected to an external force $F_{y} = m a (l - l_{0})$ , and the oscillating as block experiences a resistive force denoted by $F_{R e s} = m c y^{'} + m b y^{3}$ . Here the quantities $a, b, c$ are physical parameters that maintain the force’s units of the two previously mentioned forces.

Figure 1.

Displays a constrained pendulum.

From the indicated Figure 1, it can be easily stated that:

l - l_{0} = - y + L - \sqrt{L^{2} - 2 r^{2} (1 - \cos φ)};

(1)

let

φ = t

(2)

The total force that affects the moving block is

F_{T o t} = F_{y} - F_{R e s}

(3)

where

F_{y} = m a (l - l_{0})

, and

F_{R e s} = m c y^{'} + m b y^{3}

The equation that controls the motion is known as:

m y^{″} = F_{T o t}

(4)

Inserting Eqs. (1)-(3) in Eq. (4), it follows that the controlling dynamical equation of this system may be written as follows:

y^{″} + c y^{'} + a y + b y^{3} = a L - a \sqrt{L^{2} - 2 r^{2} (1 - \cos t)}

(5)

Using Taylor expansion, the righ hand side of ODE (5) can be approximately simplied under the condition $r ≪ L$ and written as follows:

y^{″} + c y^{'} + a y + b y^{3} = \frac{a r^{2}}{L} (1 - \cos t)

(6)

Within perturbation methods, quadratic terms can be expressed as various odd-order terms through operations such as multiplication, division, differentiation, and integration. The application of the NPA is primarily based on distinguishing odd functions generated by damping effects from those associated with stiffness nonlinearities. In this framework, quadratic components are introduced to represent both damping and stiffness forces. The first principal integral is employed to evaluate the equivalent damping, whereas the second is used to determine the equivalent frequency. The contribution of even-order terms may be incorporated into either the main integral or the alternative one, depending on the formulation adopted. The NPA methodology has been extensively developed and examined by El-Dib in a number of publications.^29,30 Therefore, ODE as given in Eq. (5) can be expressed in the following form:

\ddot{y} + G (\dot{y}) + H (y) = 0

(7)

where

H (y) = a (1 - \frac{r^{2}}{L} (1 - \cos t)) y + b y^{3}, G (\dot{y}) = c \dot{y}

(8)

Figure 2 shows the curves given in Eq. 6, 7 in the part (a), and the curve of absolute error between them in part (b).

Figure 2.

(a) Illustrates a comparison for $y (t)$ with respect to the nonlinear ODE (6) and its comparable nonlinear ODE (7). (b): Shows the error curve for $y (t)$ with respect to the nonlinear ODE (6) and its corresponding nonlinear ODE (7) a maximum value $0.0021$ .

The guessing solution of Eq. (5) that satisfies the given ICs can be expressed as follows after the preceding NPA, see Moatimid et al.^24–27:

\tilde{u} = A \cos Ω t

(9)

wherein

Ω

represents the total frequency. It is established latter.

The equivalent frequencies are evaluated using MS as follows³¹:

ω_{e q v}^{2} = \int_{0}^{2 π / Ω} u H (u) d t / \int_{0}^{2 π / Ω} u^{2} d t = a + \frac{3}{4} b A^{2} + \frac{a r^{2}}{L (1 - 4 Ω^{2})}

(10)

On the other hand, see Moatimid et al.,^24–27 the equivalent damping term is given by

σ_{e q v} = \int_{0}^{2 π / Ω} \dot{u} G (\dot{u}) d t / \int_{0}^{2 π / Ω} {\dot{u}}^{2} d t = c

(11)

Therefore, the comparable linear ODE might be formulated as detailed below:

\ddot{u} + σ_{e q v} \dot{u} + ω_{e q v}^{2} u = 0

(12)

Figure 3 Illustrates a comparison for $(t)$ , $u (t)$ in part (a), and the absolute error between them in part (b).

Figure 3.

(a) Illustrates a comparison for $y (t)$ with respect to the nonlinear ODE (6) and its comparable linear ODE (12). (b): Illustrates the error curve for $y (t)$ with respect to the nonlinear ODE (6) and its corresponding linear ODE (12) with a maximum value $0.041$ .

To remove the damping term $c \dot{u}$ from linear ODE (12), one assumes that $u (t) = e^{- c t / 2} v (t)$ . Formerly, the corresponding standard normal form reads as:

\ddot{v} + Ω^{2} v = 0

where the total frequency

Ω

is defined as follows:

Ω^{2} = ω_{e q v}^{2} - \frac{1}{4} σ_{e q v}^{2} = \frac{(1 + 4 a + 3 A^{2} b - c^{2}) + \sqrt{{(1 - 4 a - 3 A^{2} b + c^{2})}^{2} - 16 a r^{2} / L^{2}}}{8} > 0

(13)

For an additional ease, the MS’s via the command NDSolve can be used to draw nonlinear ODE solutions as provided in Eq. (5) and their linear equivalent linear ODE as given in Eq. (10). For this motivation, suppose that the data are in the following given chosen sample:

a = 0.5, b = 0.4, c = 0.45, L = 1, r = 0.01 . and A = 0.2

Now, we are going to perform the achieved results in view of their graphical representations. The analysis of the first problem is presented throughout Figures 4–9. To ensure clarity, each figure contains no more than two parts. Let us break them down one by one:

Figure 4.

Shows the time-dependent solution $u$ at: (a) $a (= 1.0, 1.1, 1.2, 1.3)$ , (b) $b (= 2, 3, 4, 5)$ .

Figure 5.

Shows the time-dependent solution $u$ at $c (= 0.1, 0.3, 0.5, 0.6)$ .

Figure 6.

Presents the corresponding curves in the plane $u u^{'}$ at the same values of the parameters in Fig. 4.

Figure 7.

Presents the corresponding curves in the plane $u u^{'}$ at the same values of the parameters in Fig. 5.

Figure 8.

Demonstrates the stability areas according to condition (13) at the same values of the considered parameters in Fig. 4.

Figure 9.

Demonstrates the stability areas according to condition (13) at the same values of the considered parameters in Fig. 5.

The solution’s response over time (displacement) $u$ of the first problem is shown in Figure 4(a) and (b) and 4. Figure 4(a) and (b) displays the responses $u$ over time for two different parameter cases $a, b$ , while Figure 5 presents results for a third parameter $c$ setting. It is noted that as an increase of $a$ (from $1.0$ to $1.3$ ), the amplitude of the oscillations stabilizes more quickly, and the system appears more stable. As $b$ increases (from 3 to 6) the oscillations reduce in amplitude, leading to further stabilization. Additionally, with increasing $c$ (from $0.1$ to $0.6$ ), the oscillations diminish more rapidly, suggesting increasing stability. Higher $c$ values are associated with more rapid settling of oscillations, indicating increased resistance to displacement changes and a quicker stabilization.

The related phase plane diagrams of the visualized curves are provided in Figure 6(a) and (b) and Figure 7. Figure 6(a) and (b) show the phase diagrams (velocity $u^{'}$ versus displacement $u$ ) for the first two parameter cases, and Figure 7 corresponds to the third case. With increasing $a$ , the trajectories form more stable closed loops, indicating stable periodic motion. Smaller, more consistent loops as an increase represent more efficient cycling through states, without significant external influence. With increasing $b$ , the trajectories shrink in size, indicating the system’s progression towards a stable state. Higher $c$ values result in tighter loops, reflecting damped behaviour and stable equilibrium points. The tightening of loops suggests a fading oscillatory response, while larger $c$ values lead to phase spirals converging more tightly toward the centre, visual evidence of a highly damped and stabilized system.

A closer look at the stability regions according to condition (13) is given in Figure 8(a) and (b) and Figure 9. Figure 8(a) and (b) display the stability regions as defined by the parameters $A$ and $Ω^{2}$ for the first two parameter sets, while Figure 9 illustrates these regions for the third case. Stable and unstable regions are clearly delineated. As $a$ increases, the boundary for the transition to stability shifts, potentially increasing the stable region. The parameter $b$ also shifts the boundary between stability and instability. Increasing $b$ further expands the stable zone. Another parameter, $c$ , has a significant effect on the stable/unstable boundary. As $c$ increases, stability is further promoted, with substantial shifts toward a stable regime. Larger $c$ values reflect greater mechanisms for stabilization, leading to a more substantial stable region and better control over oscillatory behaviour.

Overall, these figures collectively illustrate how parameters $a, b,$ and $c$ influence the system’s stability and dynamic response. An increase in each parameter triggers stabilization mechanisms, as demonstrated by time histories, phase planes, and stability region boundaries. In summary, higher values generally result in more favourable conditions of stability, reducing undesirable oscillations and expanding operational safety margins.

2.1.1. Sensitivity analysis of bifurcations and chaotic dynamics to parameters

In this section, we analyze plotted figures that show how a system’s behavior changes as we adjust its settings. These graphs highlight the critical points where the system undergoes a major shift, revealing new stable states or periodic patterns. We also use PMs to gauge the system’s stability and potential for chaos. Together, these procedures reveal a deeper understanding of how complex systems evolve.²⁴ The “butterfly effect' is a famous way to describe chaos: minuscule changes at the start can completely alter the future, making it impossible to predict what will happen in the long run.²⁵ While the path a chaotic system takes seems random, it is actually following strict, deterministic rules, playing out within a complex, fractal-like shape known as a strange attractor.

We can measure this inherent unpredictability using LE. In simple terms, the LE calculates the average speed at which two nearly identical starting points will drift apart over time. A positive LE is indicative of chaos, indicating a constant divergence of trajectories. A negative LE, on the other hand, indicates stability and convergence.²⁵ Because of this, LE acts as a fundamental numeric signature for chaos, helping us to distinguish it from simple irregularity. Adjusting parameter $Ω$ changes the system’s behaviour. The following diagrams map out these different motion types and measure the resulting chaos.

Figures 10 and 11 are calculated at $a = 2$ and $a = 1$ , respectively. Both figures represent the maximal value of $y (t)$ as a function of $L$ over the interval $[2, 5]$ . In Figure 10, the system starts with higher values of $y (t)$ , decreasing smoothly as $L$ increases, and then flattens at about $y \approx 0.5$ for $L > 3.4$ . The BD is continuous and shows no complex features such as period-doubling or chaos. In Figure 11, the maximal $y (t)$ values start lower and decrease very rapidly between $L = 2$ and $L \approx 2.9$ , then quickly settle to a constant value $y \approx 0.525$ . The transition to the steady state occurs over a narrower range of $L$ compared to Figure 9.

Figure 10.

Shows BD as $L$ varies $(2 to 5)$ , $a = 2$ .

Figure 11.

Displays BD as $L$ varies $(2 to 5)$ , $a = 1$ .

Based on these graphs, one can interpret this discussion physically as follows: Increasing $a$ from $1$ to $2$ appears to delay and smooth out the onset of steady-state behaviour and allows for higher maximal response amplitudes, suggesting a more robust response for larger $a$ , while smaller $a$ leads the system rapidly and directly into steady-state with little transient oscillatory behaviour.

Both Figure 12 ( $a = 2, L = 2.7$ ) and Figure 13 ( $a = 2, L = 3.7$ ) show the Poincaré section in the $(y - y^{'})$ plane for their respective parameter sets. In Figure 12, the points form a nearly closed curve with a few sparse outliers, suggesting quasi-periodic motion with a clear dominant cycle and some minor irregularities. The amplitude and spread are larger, indicating higher energy or larger oscillations.

Figure 12.

Displays BM for $L = 2.7$ .

Figure 13.

Shows PM for $L = 3.7$ .

In Figure 13, the closed curve is much tighter and smaller, with fewer and less pronounced, energetic, and potentially slightly less regular oscillation pronounced outliers, demonstrating a more regular periodic motion with a smaller amplitude and reduced dynamical range, possibly indicating stronger damping or weaker forcing. Physically, one can state that the system with $a = 2$ (a larger) displays most, while with $a = 1$ , the motion is more strictly periodic and constrained. This agrees with the BDs: a larger a supports richer and larger oscillatory behaviour.

Figure 14 at $a = 2$ and Figure 15 at $a = 1$ show the two largest LEs over the interval $L \in [2, 5]$ . In Figure 14, both exponents are negative throughout, indicating stable (non-chaotic) behaviour, but there is a clear variation with L: the exponents separate, then level off (where the BD also levels off), and then diverge again. This suggests varying rates of convergence/divergence in different $L$ ranges, possibly corresponding to transient changes in the system’s response.

Figure 14.

Shows the LEs when $a = 2$ .

In Figure 15, both LEs are always negative and decrease only slightly. The higher exponent $λ_{1}$ approaches $- 0.5$ for large $L$ , while the lower $λ_{2}$ is more negative. There is much less structure-this matches the rapid transition to steady state seen in the bifurcation and PMs. Therefore, the system is always stable (no chaos for either set), but the exponents’ behaviour for $a = 2$ shows more sensitivity to $L$ , which is consistent with the more extended non-steady-state region in the other plots. For $a = 1$ , the LEs show more uniform damping, explaining the rapid settling seen elsewhere.

Figure 15.

Displays the LEs when $a = 1$ .

2.2. example 2

The second example involves a mass $m$ moving without friction along a circular path of radius $R$ which rotates at a constant angular velocity $Ω$ about its physical diameter, as illustrated in Figure 16. The forces acting on the particle include the gravitational force $m g$ , the centrifugal force $m Ω^{2} R \sin θ$ , and the reaction force $N$ . By taking moments about the center of the circle and equating their total to the rate of change of the particle’s angular momentum around $O$ , we obtain:

Figure 16.

Displays a particle moving on a circular wire that rotates smoothly.

As previously shown,^32,33 the kinetic component of the energy in the particle is T:

T = 0.5 {(R \dot{θ})}^{2}

(14)

The gravitational potential becomes: $V g = - m g R \cos θ$ .

The potential energy from fictitious centrifugal force: $V c e n t r i f u g a l = - 0.5 m Ω^{2} R^{2} si n^{2} θ$

So, the total potential energy will be:

V = - m g R \cos θ - 0.5 m Ω^{2} R^{2} si n^{2} θ

(15)

Formerly, the Lagrangian may be written as:

L = T - V = 0.5 m {(R \dot{θ})}^{2} + m g R \cos θ + 0.5 m Ω^{2} R^{2} si n^{2} θ

(16)

Applying Euler–Lagrange equation, one gets

\frac{d}{d t} (\frac{\partial L}{\partial \dot{θ}}) - (\frac{\partial L}{\partial θ}) = 0

(17)

Inserting Eq (16). in Eq. (17), the equation describing the motion is given by:

m R^{2} \ddot{θ} = m Ω^{2} R^{2} \sin θ \cos θ - m g R \sin θ

(18)

Dividing both sides over $m R^{2}$ , one gets:

\ddot{θ} - Ω^{2} \sin θ \cos θ + \frac{g}{R} \sin θ = 0

(19)

The ICs is represented as:

θ (0) = C, and \dot{θ} (0) = 0

(20)

Eq (19). can be stated as follows:

\ddot{θ} + f (θ) = 0

(21)

where

f (θ) = - Ω^{2} \sin θ \cos θ + \frac{g}{R} \sin θ

The initial solution may be represented as follows:

ζ = C \cos Π t

(22)

where

Π

is the total frequency, which will be originated later.

With similar ICs as given in Eqs. (20).

Because of the NPA’s benefit over all other conventional perturbation techniques, using Taylor expansion’s drawback is not necessary in this case.

As previously shown by Moatimid et al.^24–27 and El-Dib,^28–31 the equivalent frequency could be evaluated as follows:

Π^{2} = \frac{1}{C R} (2 g J_{1} (C) - R Π^{2} J_{1} (2 C))

(23)

where

J_{1} (C)

, and

J_{1} (2 C)

denote Bessel functions for the arguments

C

and

2 C

, in respective order.

Eq (23). can be rewritten as follows:

Π^{2} = \frac{2 g J_{1} (C)}{R (C + J_{1} (2 C))}

(24)

The stability criterion requires:

\frac{2 g J_{1} (C)}{R (C + J_{1} (2 C))} > 0

(25)

Accordingly, the equivalent linear equation may be represented as:

\ddot{ζ} + Π^{2} ζ = 0

(26)

For more convenience, MS via command NDsolve to match the nonlinear ODE as given in Eq. (19) based on the linear comparative one as reflected in Eq. (26). Choose a sample system with the data:

Ω = 1.0, g = 10.0, R = 5.0, and C = 0.2

A significant coincidence arises between two the plane curves of the oscillators, one described by a nonlinear ODE, as specified in Eq. (19), and the other following a linear form, as indicated in Eq. (26). This coincidence occurs when their trajectories in phase space or over time align evidently, despite their essentially distinct mathematical frameworks. The phenomenon appears under specific conditions where the impact of nonlinearity is negligible or sufficiently controlled by designated parameters, resulting in motion that closely resembles the linear case. Under specific conditions, the nonlinear system may display harmonic or nearly harmonic oscillations with frequency and amplitude attributes akin to those of its linear equivalent. This behaviour indicates a notable underlying symmetry within the nonlinear dynamics. The MS calculates the Absolute Error between the two solutions, yielding a value of 0.0440441 Figure 17.

Figure 17.

Illustrates a strong correspondence between the two solutions of Eqs (19). and (26).

A plot of the absolute error, $| θ - ζ |$ , over the time interval $t = 0$ to $100$ (see Figure 18). It displays an oscillatory error with an amplitude that grows linearly over time, reaching a maximum value of approximately $0.040$ by the end of the simulation.

Figure 18.

Shows absolute error $| θ - ζ |$ .

Here, we discuss the obtained results for the considered dynamical model in Example 2. Let us perform and interpret Figures 19–21 in the context referring to the movement of a mass $m$ in the absence of friction on a rotating circle.

Figure 19.

Presents the behaviour of the solution $u$ and the related phase plane plots at: (a) $R (= 5.0, 5.05, 5.1, 5.15)$ , and (b) $C (= 0.2, 0.3, 0.4, 0.5)$ .

In view of the plotted curves in Figure 19, one can perform the panels (a) and (b), which express the time histories of $ζ$ , as $R$ increases, the oscillation frequency increases (i.e.,the cycles become more numerous within the same time window), but the amplitude stays nearly the same. The increase in oscillation frequency with $R$ suggests that the system’s restoring force or effective “springiness” becomes stronger with larger $R$ . Increasing the radius of the ring lengthens the path along which the mass can move. The amplitude consistency implies energy input/output or damping is not strongly affected by $R$ in this regime. One can physically interpret that the rotation with angular velocity $Ω$ creates centrifugal and Coriolis forces. As $R$ increases, the effect of centrifugal force (which grows with radius) changes the “potential landscape' in which the mass moves constriction the restoring force and thus speeding up the oscillations. The amplitude remains similar, in which the energy remains roughly constant. So, the maximum displacement does not dramatically change. Furthermore, as $C$ increases, both amplitude and frequency of oscillation increase, meaning the mass swings farther from the centre. This parameter likely increases the energy in the system, or it could represent a parameter related to the system’s initial conditions, enhancing the oscillations. Higher $C$ leads to less “restraint,” allowing larger excursions in $ζ$ .

In Figure 20, panels (a) and (b) show how the phase portrait for this system plots ( $ζ^{'}$ vs. $ζ$ ) of the represented curves in the panels Figure 19 a, Figure 19 b, respectively. All trajectories are closed loops (indicative of undamped oscillations or simple harmonic motion), as in (a). The loops are concentric (slightly shifting outwards with $R$ ), representing increased frequency but similar amplitude. This means that the system’s phase space behaviour is typical of a simple oscillator, with $R$ primarily affecting the frequency, not the amplitude. Regarding the phase plots for $C$ variation in panels (b), we can observe that the plotted loops grow larger as $C$ increases, reflecting higher amplitude oscillations for larger $C$ , as before, $C$ amplifies the system response, so energy or coupling represented by $C$ directly translates to larger oscillations.

Figure 20.

Presents the behaviour of the phase plane plots at: (a) $R (= 5.0, 5.05, 5.1, 5.15)$ , and (b) $C (= 0.2, 0.3, 0.4, 0.5)$ .

Analysis of the stability curves axes $Π^{2}$ versus $C$ , of several $R$ values is presented in Figure 21. The plot clearly separates parameter space into “Stable” (upper, shaded) and “Unstable” (lower, shaded) regions. Each curve (red, blue, green, black) corresponds to a constant $R$ . Based on this description, one can interpret these curves as: For each $R$ , there is a critical curve: Above the curve, the system is stable (oscillations do not grow); below the curve, the system is unstable (perturbations grow without bound). As $R$ increases, the curve shifts downward and right-meaning:

Figure 21.

Illustrates the stability regions considering condition (25) at $R (= 5.0, 5.05, 5.1, 5.15)$ .

For larger $R$ values, a larger $C$ and/or smaller $Π^{2}$ is required before the system becomes unstable. For a given $C$ , a lower $Π^{2}$ (i.e.,a lower frequency or energy parameter) is needed to destabilize the system as $R$ increases. In other words, the system is more stable (against perturbations) for smaller $R$ , and more susceptible to instability for larger $R$ or higher $C$ .

Generally, to keep the system stable, you need to manage the size of $R$ , the excitation/couplingC, and keep $Π^{2}$ above its critical value relative to $C$ and $R$ . Instabilities arise when system parameters cross the outlined thresholds, potentially leading to runaway oscillations.

2.2.1. Chaos and bifurcation analysis with parameter sensitivity

Here, we investigate how the system responds to changes in its parameters by analyzing a series of plots. These figures illustrate key transitions, such as the emergence of new steady states or periodic behavior. PMs are used to assess stability and detect signs of chaos, offering deeper insights into the evolution of complex systems.

Figure 22 Figure 22 present the BDs illustrating the maximum angular displacement, $\max [θ (t)]$ as a function of the driving frequency parameter $Ω$ . The system parameters and ICs were varied as follows: for Figure 22, $c = 0.1$ , and $Ω$ varies in the range $[1, 3]$ ; for the Figure 23, $c = 0.15$ , and $Ω$ varies in the range $[0.5, 1.5]$ at $g = 10$ and $R = 5 .$

Figure 22.

Shows BM as $Ω$ varies $(1 to 3)$ .

Figure 23.

Shows BM as $Ω$ varies $(0.5 to 1.5)$ .

Based on the plotted graphs in Figure 22, one concludes that when $Ω < 1.4$ , the system exhibits stable, small-amplitude oscillations, with $\max [θ (t)]$ remaining close to the IC ( $c = 0.1$ ). Around $Ω \approx 1.4$ , a clear bifurcation occurs, manifested by a sharp increase in the maximum amplitude. This transition marks the threshold for parametric resonance, where driving frequency matches a natural instability of the system. Furthermore, beyond this threshold, $\max [θ (t)]$ continuously increases with increasing $Ω$ , indicating amplified regular oscillatory motion. The smooth, single-valued nature of plot in this regime suggests that the system remains periodic and stable without evidence of chaos or period-doubling in the scanned parameter range.

In contrast, Figure 23 is plotted, where a higher initial angle ( $c = 0.15$ ) and lower $Ω$ range $[0.5, 1.5]$ were chosen. It is shown that the maximum angular displacement remains nearly constant near the initial value of all examined $Ω$ . There is no indication of a bifurcation or the onset of large initial angular velocity motion in this regime. The slight fluctuations and isolated spike observed are likely due to numerical sensitivity or isolated resonant effects, but the overall response is characterized by robust stability and absence of resonance.

The emergence of bifurcation and parametric resonance is highly sensitive to both the initial displacement and the range of

Ω

. Figure 23 captures the typical transition from stable, low-amplitude oscillations to large, regular parametric excitation as the driving frequency increases; meanwhile, Figure 23 demonstrates that, away from resonance, the system returns to near its original state. These results underscore the importance of the driving parameter

Ω

in governing the nonlinear dynamics and highlighting the system’s ability to undergo abrupt transitions, a hallmark of bifurcation phenomena in nonlinear oscillatory systems. In the following table 1 shows a comparison between Figures 22 and 23.

Table 1.

Shows a comparison of two Figures 22 and 23.

Fig.	$c$	Range of $Ω$	Behavior	Remarks
22	$0.10$	$[1, 3]$	Abrupt bifurcation; smooth amplitude increase after threshold	Onset of parametric resonance near $Ω \sim 1.4$
23	$0.15$	$[0.5, 1.5]$	Tiny, stable oscillations; no bifurcation	System remains near the initial condition

Figures 24 and 25 display PMs for two sets of parameters: $c = 0.1, Ω = 2$ and $c = 0.15, Ω = 1$ , respectively, when $g = 10$ and $R = 5$ . The Poincaré section is constructed by sampling the phase space $({θ, θ}^{'})$ at intervals of the driving period. When $Ω = 2$ , as seen in Figure 23, the PM forms a smooth, closed curve, indicating the presence of regular, bounded motion that may be periodic or quasi-periodic. The shape of the orbit is distinct from a simple ellipse or circle, suggesting significant nonlinear effects and involvement of more than one frequency in the system dynamics. Furthermore, the distribution of points shows no indication of chaos (no scattered cloud), but marks a transition away from strictly linear response.

Figure 24.

Displays PM for $Ω = 2$ .

Figure 25.

Shows PM for $Ω = 1$ .

At $Ω = 1$ , as drawn in Figure 25, one concludes that the PM forms a nearly perfect circle in phase space. This indicates a simple, regular periodic motion characteristic of a linear or weakly nonlinear oscillator. The system oscillates harmonically with minimal influence from nonlinear effects, and the motion is highly predictable.

The interpretation of these results can be stated as follows; increasing the driving frequency from

Ω = 1

Ω = 2

bring out nonlinearities and richer phase-space structures, but the system remains in a regime of regular (non-chaotic) motion for both cases. The contrast between the maps highlights how system parameters influence the transition from linear, harmonic oscillations to more complex periodic or quasi-periodic behaviour. No chaotic attractor is observed in these regimes, which is consistent with the BDs described previously for these parameter values. The Table 2 below explore a comparison between Figures 22 and 23.

Table 2.

Shows a comparison of two Figures 23 and 24.

Fig.	$c$	$Ω$	Behavior	System behavior
24	$0.10$	2	Smooth closed curve	Regular, large-amplitude, possibly quasi-periodic
25	$0.15$	1	Circle	Simple periodic (harmonic-like)

Figures 26 and 27 show the computed Lyapunov spectra for the rotating system as a function of the driving frequency Ω, with parameters $g = 10.0, R = 5$ , and initial angular velocity $c = 0.1$ (Ω from 1.5 to 4.0, Figure 25) and $c = 0.15$ (Ω from 0.5 to 1.5, Figure 26).

Figure 26.

Shows LEs when $g = 10.0, R = 5, c = 0.1$ .

Figure 27.

Shows LEs at $g = 10.0, R = 5, c = 0.15$ .

Based on Figure 26, we can write; the curves indicate the two largest LEs, $λ_{1}$ (blue) and $λ_{2}$ (red). For all values of Ω, $λ_{1}$ is slightly positive (ranging roughly from 0.004 to 0.009) while $λ_{2}$ is approximately the same magnitude but negative (−0.004 to −0.007). A positive largest LE $λ_{1} > 0$ is a classic indicator of chaotic dynamics: infinitesimal perturbations in ICs grow exponentially with time. The relatively low values of $λ_{1}$ suggest weak chaos, in which the system is sensitive but not in a strongly turbulent regime. The negative second exponent is typical of two-dimensional dissipative systems, where the phase-space volume contracts despite local stretching.

According to Figure 27, we can see that, for this parameter regime, $λ_{1}$ is again consistently positive (about 0.006 to 0.009) and $λ_{2}$ is negative (−0.003 to −0.007). The spectrum remains in a regime indicative of chaos or at least complex dynamics, with sensitivity to ICs. Similar structure and amplitude of exponents in both figures suggest that the system tends to display weak chaos throughout both frequency intervals, independent of the initial angle parameter c, though the precise rate of divergence of nearby trajectories varies with Ω.

3. Concluding remarks

Nonlinear oscillations commonly manifest in real-world problems, including mechanical structures and biological cycles, leading to intricate behaviors such as bifurcations and chaos, which are essential in comprehending and forecasting dynamical systems. Progress in analytical and numerical techniques has facilitated innovative applications in engineering, medicine, and materials research. We examined oscillators exhibiting significant nonlinear characteristics, integrating both damping and restoring forces, through a combination of theoretical and computational methodologies. Two instances were provided from various scientific and technical domains. The novel methodology presented markedly decreases computing time and resource utilization relative to traditional perturbation methods commonly employed in this domain. Utilizing HFF, the suggested NPA converts a weakly nonlinear oscillator described by an ODE into a linear form, hence enabling the derivation of a new frequency associated with the linearized ODE. The theoretical conclusions were corroborated through NSs conducted via MS, demonstrating a strong correlation between the two ODEs. This approach enables a comprehensive examination of system stability, enhancing capabilities beyond those offered by prior methods. Accordingly, the NPA provides a more pragmatic and dependable framework through NSs of weakly nonlinear oscillators. Furthermore, it is distinguished as a versatile instrument for applied research and engineering owing to its adaptability to many nonlinear scenarios. We also analyze the influence of various parameters on stability, with results validating the approach’s simplicity, efficiency, and reliability. Modifying bifurcation parameters changes the configuration of BDs and their corresponding PMs, which we further demonstrated by mapping the LE curves. The principal key outcomes can be encapsulated as follows:

1. Two nonlinear oscillators governed by weak ODEs were analysed.

2. The stability study was conducted for each oscillator.

3. Dangerous operating zones characterized by unpredictable and unstable nonlinear oscillations are detected through chaos analysis.

4. The critical parameter thresholds that induce abrupt changes in system dynamics are determined using bifurcation analysis.

5. The LEs quantitatively assess stability and susceptibility to ICs.

Investigating linked weakly nonlinear oscillators is crucial as numerous physical and engineering systems comprise interacting components whose movement’s exhibit nonlinearity, particularly at moderate amplitudes. Weak nonlinearity engenders phenomena including amplitude-dependent frequencies, energy transfer across modes, synchronization, internal resonance, bifurcations, and the emergence of chaotic behavior, which linear theory alone cannot predict. Coupling enables oscillators to transfer energy and information, rendering these models crucial for comprehending vibrations in mechanical structures, electrical circuits, MEMS and NEMS devices, aeronautical systems, rotating machinery, biological rhythms, laser dynamics, and power grids. Analyzing coupled weakly nonlinear oscillators from an engineering standpoint facilitates the design of stable and efficient systems, the control of undesirable vibrations and resonance, the enhancement of signal processing and communication technologies, the improvement of energy harvesting devices, and the prediction of failure mechanisms in complex interconnected systems.

Footnotes

Acknowledgments

The Researchers would like to thank the Deanship of Graduate Studies and Scientific Research at Qassim University (www.qu.edu.sa) for financial support (QU-APC-2026).

ORCID iDs

Tarek S. Amer

Abdallah Galal

Author contributions

Ahmad Almutlg: Conceptualization, Methodology, Formal Analysis, and Validation.

Galal Moatimid: Conceptualization, Resources, Methodology, Formal analysis, Validation, Visualization, Writing- Original draft preparation, Reviewing.

Tarek S. Amer: Investigation, Validation, Methodology, Data duration, Conceptualization, Reviewing.

Abdallah Galal: Conceptualization, Formal analysis, Validation, Data duration, Software, Reviewing and Editing.

Funding

The Researchers would like to thank the Deanship of Graduate Studies and Scientific Research at Qassim University (www.qu.edu.sa) forfinancial support (QU-APC-2026).

Declaration of conflicting interests

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

Data Availability Statement

This study was not supported by any data.*

Appendix

References

Ӧcalan

Ӧzkan

. Oscillations of dynamic equations on time scales with advanced arguments. International Journal of Dynamical Systems and Differential Equations 2016; 6(4): 275–284. https://doi.org/10.1504/ijdsde.2016.081813

Biçer

Tunc

. Tunç, On the asymptotic stability behaviours of solutions of non-linear differential equations with multiple variable advanced arguments. Journal of Applied Nonlinear Dynamics 2019; 8(2): 239–249. https://doi.org/10.5890/jand.2019.06.007

Ismail

Abul-Ez

Zayed

, et al. Analytical accurate solutions of nonlinear oscillator systems via coupled homotopy-variational approach. Alexandria Engineering Journal 2022; 61(7): 5051–5058. https://doi.org/10.1016/j.aej.2021.09.021

Kopmaz

. Free vibrations of a mass grounded by linear and non-linear springs in series. Journal of Sound and Vibration 2006; 289(4): 689–710.

Manimegalai

Sagar Zephania

Bera

, et al. Study of strongly nonlinear oscillators using the Aboodh transform and the homotopy perturbation method. The European Physical Journal Plus 2019; 134: 462. https://doi.org/10.1140/epjp/i2019-12824-6

Abu-As’ad

Asad

. Power series approach to nonlinear oscillators. Journal of Low Frequency Noise,Vibration and Active Control 2023; 43(1): 220–238. https://doi.org/10.1177/14613484231188756

Shou

D-H

. The homotopy perturbation method for nonlinear oscillators. Computers and Mathematics with Applications 2009; 58: 2456–2459. https://doi.org/10.1016/j.camwa.2009.03.034

Hermann

Saravi

Khah

. Analytical study of nonlinear oscillatory systems using the Hamiltonian approach technique. Journal of Theoretical and Applied Physics 2014; 8: 133. https://doi.org/10.1007/s40094-014-0133-9

Ismail

. An analytical coupled homotopy-variational approach for solving strongly nonlinear differential equation. Journal of the Egyptian Mathematical Society 2017; 25: 434–437. https://doi.org/10.1016/j.joems.2017.07.006

10.

Ismail

Abu-Zinadah

. Analytic approximations to non-linear third order jerk equations via modified global error minimization method. Journal of King Saud University - Science 2021; 33: 101219. https://doi.org/10.1016/j.jksus.2020.10.016

11.

Kang-Le

Bifurcation

. variational principle and new soliton solutions for the nonlinear fractional drinfeld–sokolov system. World Scientific Publishing Co, 2026. https://doi.org/10.1142/S0218348X26500623

12.

Kang-Le

. Dynamical behaviors and diverse soliton solutions for the fractional Yajima–Oikawa equation. Fractals 2026; 34(5): 2650020.

13.

Liu

Wang

K-L

. Variational principle, Hamiltonian, sensitivity analysis and new soliton solutions for the nonlinear fractional Schrödinger–Hirota equation. Fractals 2026; 34(03): 2650002. https://doi.org/10.1142/s0218348x26500027

14.

Kang-Le

. Investigation of soliton structures and dynamical behaviors for fractional Kaup–Newell system. World Scientific Publishing Co, 2026. https://doi.org/10.1142/S0218348X26500817

15.

Arif

Abodayeh

Nawaz

. Precision in disease dynamics: Finite difference solutions for stochastic epidemics with treatment cure and partial immunity. Partial Differential Equations in Applied Mathematics 2024; 9: 100660. https://doi.org/10.1016/j.padiff.2024.100660

16.

Arif

Abodayeh

Nawaz

. A reliable computational scheme for stochastic reaction–diffusion nonlinear chemical model. Axioms 2023; 12(5): 460. https://doi.org/10.3390/axioms12050460

17.

Baazeem

Nawaz

Arif

, et al. Modelling infectious disease dynamics: A robust computational approach for stochastic SIRS with partial immunity and an incidence rate. Mathematics 2023; 11(23): 4794. https://doi.org/10.3390/math11234794

18.

Shatanawi

Raza

Arif

, et al. An effective numerical method for the solution of a stochastic coronavirus (2019-nCovid) pandemic model. Computers, Materials & Continua 2021; 66(2): 1121–1137. https://doi.org/10.32604/cmc.2020.012070

19.

Shatanawi

Raza

Arif

, et al. Essential features preserving dynamics of stochastic dengue model. Computer Modeling in Engineering & Sciences 2021; 126(1): 201–215. https://doi.org/10.32604/cmes.2021.012111

20.

. Some asymptotic methods for strongly nonlinear equations. International Journal of Modern Physics B 2006; 20: 1141–1199. https://doi.org/10.1142/s0217979206033796

21.

J-H

. The simpler, the better: analytical methods for nonlinear oscillators and fractional oscillators. Journal of Low Frequency Noise, Vibration and Active Control 2019; 38: 1252–1260. https://doi.org/10.1177/1461348419844145

22.

Qie

Hou

J-H

. The fastest insight into the large amplitude vibration of a string. Reports in Mechanical Engineering 2020; 2: 1–5.

23.

Almutlg

Moatimid

Amer

, et al. Analytical modelling of a nonlinear 2DOF Mathieu–Duffing oscillator. Journal of Low Frequency Noise, Vibration and Active Control 2026; 10: 1. https://doi.org/10.1177/14613484261460109

24.

Moatimid

Amer

Galal

. Studying highly nonlinear oscillators using the non-perturbative methodology. Scientific Reports 2023; 13: 20288. https://doi.org/10.1038/s41598-023-47519-5

25.

Moatimid

Amer

Galal

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

26.

Moatimid

Amer

Galal

AAA

. Non-Perturbative framework in analyzing weakly nonlinear oscillators and their chaotic dynamics. Journal of machines 2026; 14: 267. https://doi.org/10.3390/machines14030267

27.

Moatimid

Amer

Mohamed

MAA

, et al. Insights of Nonlinear Vibration of Oscillators Linked with Two-Degrees-of-Freedom. Journal of Computational Applied Mechanics 2025; 56(2): 424–456. https://doi.org/10.22059/JCAMECH.2025.390125.1357

28.

El-Dib

. Insightful and comprehensive formularization of frequency-amplitude formula for strong or singular nonlinear oscillators. Journal of Low Frequency Noise, Vibration and Active Control 2023; 42(1): 89–109. https://doi.org/10.1177/14613484221118177

29.

El-Dib

. The masking technique for forced nonlinear oscillator stability behavior analysis using the non-perturbative approach. Journal of Low Frequency Noise, Vibration and Active Control 2024; 43(4): 1481–1497. https://doi.org/10.1177/14613484241253178

30.

El‐Dib

. A comprehensive study of stability analysis for nonlinear Mathieu equation without a non-perturbative technique. ZAMM‐Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik 2024; 104(9): e202400047. https://doi.org/10.1002/zamm.202400047

31.

El-Dib

. Insights into transferal to fractal space modeling: delayed forced Helmholtz–Duffing oscillator with the non-perturbative approach. Communications in Theoretical Physics 2024; 77(1): 015002. https://doi.org/10.1088/1572-9494/ad7ceb

32.

Strogatz

. Nonlinear Dynamics and Chaos: With applications to Physics, Biology. Chemistry, and Engineering (Studies in nonlinearity). CRC Press, 2000.

33.

Amer

Moatimid

Zakria

, et al. Vibrational and stability analysis of planar double pendulum dynamics near resonance. Nonlinear Dynamics 2024; 112: 21667–21699. https://doi.org/10.1007/s11071-024-10169-x