Abstract

This study investigates the primary-resonance response and local stability of a reduced Duffing-type oscillator with fractional-order damping and an effective stiffness correction representing the dominant influence of a spatially periodic medium. The model includes cubic nonlinear stiffness, weak fractional damping of order $g$ , harmonic excitation, and a linear stiffness-modification parameter. The fractional derivative is interpreted in the Caputo sense and is represented under the long-time harmonic approximation for the steady-state analysis. Using the method of multiple scales, amplitude and phase modulation equations are derived near primary resonance. The resulting algebraic frequency-response equation is solved numerically to obtain the positive real steady-state amplitudes, while the eigenvalues of the Jacobian matrix of the modulation system are used to distinguish stable and unstable branches. The results show that positive cubic nonlinearity produces hardening-type bending of the resonance response. Variations in $g$ modify both the effective dissipative and effective detuning contributions through the factors $\text{sin}\left(\pi g/2\right)$ and $\text{cos}\left(\pi g/2\right)$ , respectively. The effects of the fractional order, damping coefficient, and excitation amplitude on the steady-state response are examined. The analytical frequency-response relation is also compared with direct fourth-order Runge-Kutta integration in the classical limiting case $g=1$ , where the fractional derivative reduces to the ordinary first derivative. This comparison is used only as a consistency check for the classical limit and is not a direct time-domain validation of the fractional model for $0<g<1$ . The results provide an analytical framework for studying resonance and local stability in nonlinear oscillators with memory-dependent damping and effective stiffness corrections associated with periodic media.

Keywords

Fractional Calculus, Duffing Oscillator, Fractional-Order Damping, Primary Resonance, Frequency Response, Stability Analysis, Multiple Scales, Periodic Medium

Share and Cite:

1. Introduction

The Duffing oscillator is one of the fundamental models used to describe nonlinear vibrations involving cubic stiffness effects, geometric nonlinearities, and nonlinear restoring forces. Its response may exhibit amplitude-dependent resonance frequencies, bending of frequency-response curves, jump phenomena, and the coexistence of multiple steady-state solutions. These features make the Duffing model a useful prototype for studying nonlinear resonance and stability in mechanical and structural systems [1] [2].

Classical Duffing models generally incorporate viscous damping through a first-order time derivative. Although this representation is suitable for many systems, it does not fully describe materials and structures whose dissipative response depends on their deformation history. Fractional calculus provides a convenient mathematical framework for modeling such memory-dependent and hereditary effects [3]-[14]. Fractional constitutive models have been used extensively in viscoelasticity, wave propagation, solid mechanics, and dynamic systems because they can represent frequency-dependent dissipation using relatively few parameters [5]-[14].

The application of fractional derivatives to nonlinear oscillators has shown that the fractional order may influence not only the amount of energy dissipation but also the phase relation and the effective dynamic stiffness. In harmonic analysis, the fractional operator introduces the complex factor ${\left(i\omega \right)}^g$ , whose real and imaginary parts contribute differently to the response. Consequently, changing the fractional order can modify the resonance amplitude, the location of the resonance peak, the shape of the frequency-response curve, and the stability of the resulting solution branches [15]-[19]. These effects are particularly relevant in systems involving viscoelastic components, fractional-order resonators, and metamaterial-inspired effective models [20] [21].

Spatially periodic media provide another important setting in which resonance and effective stiffness corrections arise. Periodic lattices and structured media may exhibit modified wave-propagation characteristics and effective dynamic parameters compared with homogeneous systems. In a reduced single-mode description, the influence of the underlying periodic medium may be represented by an effective correction to the linear stiffness. Such a reduction does not replace a full distributed wave model; rather, it provides a reduced-order oscillator that retains selected nonlinear resonance characteristics of the underlying periodic structure [22]-[24].

Several studies have investigated the primary resonance of fractional Duffing-type systems using averaging, perturbation, and harmonic-balance approaches [15] [16] [19]. Other investigations have addressed bifurcation, stability, and complex nonlinear behavior in fractional-order systems [17] [18] [25]-[30]. Numerical methods for fractional differential equations have also been developed to account explicitly for the nonlocal memory term [31] [32]. Nevertheless, a clear distinction is needed between an analytical steady-state treatment based on the harmonic representation of the fractional derivative and a direct time-domain numerical solution of a fractional differential equation. In particular, a classical Runge-Kutta method applied directly to the equation with $g=1$ validates only the classical limiting case; it does not constitute a numerical integration of the Caputo derivative for $0<g<1$ [31] [32].

The objective of the present work is to analyze the primary-resonance response and local stability of a reduced Duffing-type oscillator containing fractional-order damping and an effective stiffness correction. The main contribution is the systematic characterization of the effects of the fractional order on the analytical frequency-response curves and their stable and unstable branches. The analysis combines the method of multiple scales with numerical continuation of the steady-state amplitude equation [1] [2] [15] [16] [19]. The positive real roots of the response equation are computed over a prescribed detuning interval, and the local stability of each branch is evaluated from the Jacobian eigenvalues of the modulation system [19] [25] [26] [28]. The effects of the fractional order $g$ , damping coefficient $\mu$ , and excitation amplitude $F$ are then examined separately [3] [15] [16] [19] [20].

The remainder of the paper is organized as follows. Section 2 presents the reduced mathematical model and the assumptions used in the analysis. Section 3 introduces the harmonic representation of the fractional derivative and derives the modulation equations using the method of multiple scales. Section 4 presents the steady-state frequency-response equation, the numerical continuation procedure, and the stability criterion. Section 5 discusses the results corresponding to Figures 1-6. Section 6 summarizes the principal conclusions and limitations of the present formulation.

2. Mathematical Model

We consider a reduced single-degree-of-freedom Duffing-type oscillator with fractional-order damping, cubic nonlinear stiffness, harmonic excitation, and an effective linear stiffness correction associated with the underlying spatially periodic medium. The reduced coordinate $q\left(t\right)$ represents the dominant modal displacement of the system. The governing equation is written as

$\ddot{q}\left(t\right)+q\left(t\right)+\epsilon \mu { }^C{D}_t^{g}q\left(t\right)+\epsilon \beta q^3\left(t\right)+\epsilon \frac{A_0}{2}q\left(t\right)=\epsilon F\cos \left(\Omega t\right),1<g\le 1$ (1)

where $q\left(t\right)$ is the generalized displacement, $\mu$ is the fractional damping coefficient, $\beta$ is the cubic nonlinear stiffness coefficient, $A_0$ is the effective stiffness-modification parameter, $F$ is the excitation amplitude, Ω is the excitation frequency, and $\epsilon \ll 1$ is a bookkeeping parameter that measures the weakness of damping, nonlinearity, stiffness correction, and external excitation [1]-[3] [24] [26].

The operator ${}^C{D}_t^g$ denotes the Caputo fractional derivative of order $g$ . For $0<g<1$ , it is defined by

${}^C{D}_t^{g}q\left(t\right)=\frac{1}{\Gamma \left(1-g\right)}{\displaystyle {\int }_0^t\frac{\dot{q}\left(\tau \right)}{{\left(t-\tau \right)}^g}\text{d}\tau }.$ (2)

where $\Gamma (\cdot )$ is the Gamma function [4] [9] [12] [13]. The Caputo formulation is convenient for mechanical applications because it permits the use of conventional initial conditions expressed in terms of integer-order derivatives [5]-[8] [14]. In the limiting case $g=1$ , the fractional derivative reduces to the ordinary first derivative [4] [12] [13],

$\begin{array}{c}{}^C{D}_t^{1}q\left(t\right)=\dot{q}\left(t\right).\end{array}$ (3)

Equation (1) is a reduced-order model rather than a complete distributed wave equation. The spatial periodicity of the original medium is represented here through the effective stiffness correction $A_0/2$ . This type of reduced description is suitable for investigating how periodic-medium effects modify the resonance characteristics of a dominant mode, whereas a full spatially distributed analysis would require a separate partial differential or lattice formulation [22] [24].

The term $\beta q^3$ describes the cubic restoring-force contribution. For $\beta >0$ , the oscillator has a hardening-type nonlinearity, which generally produces a resonance curve bending toward higher excitation frequencies [1] [2]. The fractional damping term introduces a memory-dependent response and, in the harmonic regime, modifies both the dissipative and reactive parts of the effective dynamic response [3] [6] [7] [14] [19].

To investigate primary resonance, the excitation frequency is assumed to be close to the linear natural frequency, which is normalized to unity. We therefore introduce the detuning relation

$\begin{array}{c}\Omega =1+\epsilon \sigma ,\end{array}$ (4)

where $\sigma$ is the detuning parameter. The sign convention in Equation (4) is used consistently throughout the analytical derivation and numerical continuation. The time scales are defined by

$T_0=t,T_1=\epsilon t.$ (5)

The displacement is expanded in the multiple-scales form

$q\left(t,\epsilon \right)=q_0\left(T_0,T_1\right)+\epsilon q_1\left(T_0,T_1\right)+O\left({\epsilon }^2\right),$ (6)

The derivative operators then become

$\frac{\text{d}}{\text{d}t}=D_0+\epsilon D_1+O\left({\epsilon }^2\right),$ (7)

and

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

where

$D_n=\frac{\partial }{\partial T_n}.$

At the leading order, $O\left(1\right)$ , Equation (1) gives

$D_0^2{q}_0+q_0=0.$ (9)

The general leading-order solution may be expressed in complex form as

$q_0\left(T_0,T_1\right)=A\left(T_1\right){\text{e}}^{i{T}_0}+\overline{A}\left(T_1\right){\text{e}}^{-i{T}_0},$ (10)

where $A\left(T_1\right)$ is a slowly varying complex amplitude and the overbar denotes complex conjugation. Equivalently, the solution can be written in amplitude-phase form as

$q_0\left(T_0,T_1\right)=a\left(T_1\right)\cos \left(T_0+\gamma \left(T_1\right)\right),$ (11)

where $a\left(T_1\right)$ and $\gamma \left(T_1\right)$ are the slowly varying amplitude and phase, respectively [1] [3] [12] [13] [25].

The perturbation assumptions underlying the analysis require that the amplitude remain within a moderate range such that the weak-perturbation ordering remains meaningful. Accordingly, the results obtained below describe the primary-resonance regime of the reduced model and should not be interpreted as a uniformly valid approximation for arbitrarily large amplitudes, strong damping, or strongly nonlinear parameter values [1] [2] [15] [16].

Figure 1 illustrates the reduced Duffing-type oscillator considered in this study, including the fractional damping element, the cubic nonlinear stiffness, the effective stiffness correction, and the harmonic excitation.

Figure 1. Schematic of the reduced Duffing-type oscillator with fractional-order damping, cubic nonlinear stiffness, effective stiffness correction, and harmonic excitation.

The schematic represents the reduced single-degree-of-freedom model used throughout the analytical and numerical analysis. The spatially periodic medium is not solved as a complete distributed wave system; its influence is represented through the effective stiffness correction $A_0/2$ .

2.1. Fractional Harmonic Representation

The Caputo derivative in Equation (2) is a nonlocal operator and therefore cannot be replaced, for general transient motion, by an ordinary local derivative. In the present steady-state multiple-scales analysis, the fractional derivative is evaluated using its long-time harmonic representation. For a harmonic component with unit leading frequency [4]-[7] [9]-[14],

$\begin{array}{c}{}^C{D}_t^g{\text{e}}^{i{T}_0}~{\left(i\right)}^g{\text{e}}^{i{T}_0},\end{array}$ (12)

where the principal complex power is [6] [9] [11] [14]-[16] [19]

${\left(i\right)}^g=\exp \left(i\frac{\pi g}{2}\right)=\cos \left(\frac{\pi g}{2}\right)+i\sin \left(\frac{\pi g}{2}\right),$ (13)

For a harmonic component with frequency $\omega$ , the corresponding relation is

${\left(i\omega \right)}^g={\omega }^g\left[\cos \left(\frac{\pi g}{2}\right)+i\sin \left(\frac{\pi g}{2}\right)\right],$ (14)

Since the leading-order frequency is normalized to unity and the detuning is introduced as a small perturbation, Equation (12) is used in deriving the first-order steady-state modulation equations. The real and imaginary parts of Equation (13) have different physical roles. The term involving

$\text{cos}\left(\frac{\pi g}{2}\right)$

contributes to the effective reactive or detuning component, whereas the term involving

$\text{sin}\left(\frac{\pi g}{2}\right)$

contributes to the effective dissipative component. Thus, changing $g$ does not merely scale the damping magnitude; it changes the phase and frequency-dependent balance between energy dissipation and effective stiffness [6] [11] [14]-[16] [19].

This approximation is appropriate for the analytical steady-state frequency-response calculation considered here. It does not represent a direct numerical discretization of the Caputo memory integral. A direct time-domain computation for $0<g<1$ would require a suitable fractional numerical method, such as a convolution-based, predictor-corrector, or other memory-preserving scheme [12] [31] [32].

2.2. Limiting Classical Case

When $g=1$ , Equation (13) becomes

${\left(i\right)}^1=i,$

because

$\cos \left(\frac{\pi }{2}\right)=0,\sin \left(\frac{\pi }{2}\right)=1.$

Consequently,

${}^C{D}_t^{1}q=\dot{q},$

and Equation (1) reduces to the classical weakly damped Duffing equation [4] [9] [12]-[14]

$\ddot{q}\left(t\right)+q\left(t\right)+\epsilon \beta q^3\left(t\right)+\epsilon \frac{A_0}{2}q\left(t\right)=\epsilon F\cos \left(\Omega t\right),$ (15)

This limiting case is used later only for the consistency check shown in Figure 6, where the resulting ordinary differential equation is integrated using the classical fourth-order Runge-Kutta method. The RK4 calculation is therefore not a direct numerical validation of the fractional model for $0<g<1$ [12] [31] [32].

2.3. Model Parameters

The parameters used in the numerical figures are selected as follows:

$\mu =0.2,\beta =0.4,A_0=0.5,F=1.5,$ (16)

unless another value is specified in the corresponding figure description. The fractional-order study uses

$g\in \left\{1.0,0.7,0.5\right\}.$ (17)

The reference response in Figure 2 is computed for one fixed fractional order, $g=0.7$ , while Figure 3 compares the three values in Equation (17). The continuation calculations use

$-2\le \sigma \le 2,0<a\le 8.$ (18)

These numerical bounds define the computational domain used to locate the positive real steady-state amplitudes. They do not represent mathematical restrictions on the model itself; rather, they specify the parameter and amplitude ranges examined in the present study.

3. Multiple-Scales Analysis and Modulation Equations

We apply the method of multiple scales to Equation (1) in order to derive the slow evolution equations governing the amplitude and phase near primary resonance. The independent time scales are

$T_0=t,T_1=\epsilon t,$

and the displacement is expanded as

$q\left(t,\epsilon \right)=q_0\left(T_0,T_1\right)+\epsilon q_1\left(T_0,T_1\right)+O\left({\epsilon }^2\right),$ (19)

The time-derivative operators are

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

and

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

where $D_j=\partial /\partial T_j$ .

Using the detuning relation

$\Omega =1+\epsilon \sigma ,$

the excitation term becomes

$\cos \left(\Omega t\right)=\cos \left(T_0+\sigma T_1\right),$ (21)

Substitution of Equations (19)-(21) into Equation (1) gives a sequence of problems in increasing powers of $\epsilon$ [1] [2] [15] [16] [19] [23].

3.1. Leading-Order Problem

At $O\left(1\right)$ , we obtain

$\begin{array}{c}{D}_0^2{q}_0+q_0=0.\end{array}$ (22)

The general real solution is written as

$q_0\left(T_0,T_1\right)=a\left(T_1\right)\cos \theta ,\theta =T_0+\gamma T_1$ (23)

where $a\left(T_1\right)$ and $\gamma \left(T_1\right)$ denote the slowly varying amplitude and phase [1] [2] [15] [16].

Equivalently, the leading-order solution can be expressed as

$q_0\left(T_0,T_1\right)=A\left(T_1\right){\text{e}}^{i{T}_0}+\overline{A}\left(T_1\right){\text{e}}^{-i{T}_0},$ (24)

The amplitude-phase representation in Equation (23) is used below to obtain the real modulation equations [1] [2] [15] [16].

3.2. First-Order Problem

At $O\left(\epsilon \right)$ , Equation (1) yields

$D_0^2{q}_1+q_1=-2D_0{D}_1{q}_0-\mu {}^C{D}_{T_0}^g{q}_0-\beta q_0^3-\frac{A_0}{2}{q}_0+F\cos \left(T_0+\sigma T_1\right).$ (25)

For the steady-state harmonic analysis, the fractional derivative is represented using [1] [2] [15] [16] [19]

${}^C{D}_{T_0}^g{\text{e}}^{i{T}_0}~{\left(i\right)}^g{\text{e}}^{i{T}_0},$

where

${\left(i\right)}^g=\cos \left(\frac{\pi g}{2}\right)+i\sin \left(\frac{\pi g}{2}\right),$ (26)

For compactness, define [4]-[7] [9]-[14]

$c_g=\cos \left(\frac{\pi g}{2}\right),s_g=\sin \left(\frac{\pi g}{2}\right).$ (27)

The first harmonic of the cubic term is obtained from

$q_0^3=a^3{\cos }^3\theta =\frac{3a^3}{4}\cos \theta +\frac{a^3}{4}\cos \left(3\theta \right).$ (28)

Only the component proportional to $\cos \theta$ is resonant with the homogeneous solution of Equation (22). The third-harmonic contribution does not generate secular growth at this order and is omitted from the solvability condition [1] [2] [15] [16].

Using the harmonic representation of the fractional derivative, the resonant part of Equation (25) can be collected in terms of $\sin \theta$ and $\cos \theta$ . Eliminating the resonant terms gives the slow-flow equations [1] [2] [15] [16] [19]

$\frac{\text{d}a}{\text{d}{T}_1}=-\frac{\mu }{2}a{s}_g+\frac{F}{2}\sin \gamma ,$ (29)

and

$\frac{\text{d}\gamma }{\text{d}{T}_1}=\sigma +\frac{A_0}{2}-\frac{3\beta a^2}{8}-\frac{\mu }{2}{c}_g+\frac{F}{2a}\cos \gamma ,$ (30)

Here, $\gamma$ is defined consistently with the phase relations

$\sin \gamma =\frac{\mu a{s}_g}{F},$ (31)

and

$\cos \gamma =\frac{2a}{F}\left[\frac{3\beta a^2}{8}+\frac{\mu }{2}{c}_g-\sigma -\frac{A_0}{2}\right].$ (32)

Equations (29) and (30) constitute the autonomous modulation system used for the numerical continuation and stability calculations [1] [2] [15] [16] [19].

Important phase convention

The phase $\gamma$ in Equations (29)-(32) is defined so that the steady-state relations take the form given in Equations (31) and (32). An alternative definition of the phase, such as $\tilde{\gamma }=\gamma +\pi /2$ , would interchange the sine and cosine terms. Such formulations are mathematically equivalent if used consistently, but mixing phase conventions changes the apparent signs in the modulation equations and leads to an incorrect Jacobian. The present manuscript uses only the convention in Equations (29)-(32).

3.3. Steady-State Conditions

At steady state,

$\frac{\text{d}a}{\text{d}{T}_1}=0,\frac{\text{d}\gamma }{\text{d}{T}_1}=0.$ (33)

Equation (29) then gives

$\frac{F}{2}\sin \gamma =\frac{\mu }{2}a{s}_g,$

or

$\sin \gamma =\frac{\mu a{s}_g}{F}.$ (34)

Equation (30) gives

$\frac{F}{2a}\cos \gamma =-\sigma -\frac{A_0}{2}+\frac{3\beta a^2}{8}+\frac{\mu }{2}{c}_g,$

which leads to

$\cos \gamma =\frac{2a}{F}\left[\frac{3\beta a^2}{8}+\frac{\mu }{2}{c}_g-\sigma -\frac{A_0}{2}\right].$ (35)

4. Steady-State Frequency-Response Equation

Squaring Equations (34) and (35) and using

${\sin }^2\gamma +{\cos }^2\gamma =1,$

we obtain

${\left(\frac{\mu a{s}_g}{F}\right)}^2+{\left\{\frac{2a}{F}\left[\frac{3\beta a^2}{8}+\frac{\mu }{2}{c}_g-\sigma -\frac{A_0}{2}\right]\right\}}^2=1.$ (36)

Multiplying Equation (36) by $F^2$ , the algebraic amplitude equation becomes

${\mu }^2{a}^2{s}_g^2+4a^2{\left[\frac{3\beta a^2}{8}+\frac{\mu }{2}{c}_g-\sigma -\frac{A_0}{2}\right]}^2-F^2=0.$ (37)

Using the definitions in Equation (27), Equation (37) may be written as

$\Phi \left(a;\sigma ,g\right)=a^2\left[{\mu }^2{s}_g^2+{\left(\frac{3\beta a^2}{4}+\mu c_g-2\sigma -A_0\right)}^2\right]-F^2=0.$ (38)

Equation (38) is the steady-state frequency-response equation used in the numerical computations. For each prescribed value of the detuning parameter $\sigma$ , the positive real roots of this equation are calculated over the amplitude interval

$0\le a\le 8.$

The detuning interval is

$-2\le \sigma \le 2,$

and it is sampled at $N_{\sigma }=301$ equally spaced values, giving

${\Delta }_{\sigma }=\frac{2-\left(-2\right)}{301-1}=\mathrm{0.0133333.}$

The amplitude interval is discretized using $N_a=4001$ grid points, corresponding to

${\Delta }_a=\frac{8-0}{4001-1}=\mathrm{0.002.}$

Sign-changing intervals of the residual function are identified and refined using a bracketed fzero procedure. Near-zero residual values are also retained to reduce the possibility of missing tangential roots. Roots separated by less than 10⁻⁶ in amplitude are treated as duplicate roots. The resulting roots define the steady-state frequency-response branches used in Figures 2-5 [16].

The positive real roots of Equation (38) represent the possible steady-state amplitudes for a prescribed detuning parameter. Depending on the parameter values, the equation may have one or several positive roots. Multiple roots correspond to coexisting steady-state responses and are associated with the nonlinear bending and multistability of the Duffing response. The roots are evaluated independently at each detuning value and are ordered according to their amplitude for branch visualization. The resulting branches therefore represent a grid-based algebraic root continuation procedure rather than a pseudo-arclength continuation algorithm [19] [23].

For each positive real root $a$ , the corresponding phase $\gamma$ is reconstructed from Equations (34) and (35). The Jacobian matrix of the modulation system is then evaluated at the corresponding steady-state point. A branch is classified as locally asymptotically stable when

$\text{Re}\left({\lambda }_i\right)<0\text{for}\text{all}i,$

where ${\lambda }_i$ are the eigenvalues of the Jacobian matrix. If at least one eigenvalue satisfies

$\text{Re}\left({\lambda }_i\right)>0,$

the branch is classified as unstable.

4.1. Classical Limiting Case

For $g=1$ ,

$s_g=1,c_g=0,$

and Equation (38) reduces to

$a^2\left[{\mu }^2+{\left(\frac{3\beta a^2}{4}-2\sigma -A_0\right)}^2\right]-F^2=0.$ (39)

This is the frequency-response equation associated with the classical viscously damped Duffing equation in the present normalization [1] [2] [4] [12] [13]. The analytical curve obtained from Equation (39) is the curve compared with the RK4 solution in Figure 6.

4.2. Numerical Root Continuation

For each value of $\sigma$ , Equation (38) is evaluated over the amplitude interval

$0<a\le 8.$

The amplitude interval is discretized using a sufficiently fine grid. Intervals in which $\Phi \left(a;\sigma ,g\right)$ changes sign are identified, and each bracket is refined using MATLAB’s fzero algorithm. Near-zero grid values are also retained to reduce the possibility of missing roots that occur close to grid points [12] [31] [32].

The procedure is repeated over the detuning interval

$-2\le \sigma \le 2.$

The resulting pairs $\left(\sigma ,a\right)$ form the steady-state frequency-response branches. The root-continuation procedure is algebraic and is based on Equation (38); it does not perform direct time integration of the fractional differential equation [12] [31] [32].

4.3. Stability of the Modulation Solutions

The local stability of each steady-state solution is determined from the Jacobian matrix of the slow-flow system in Equations (29) and (30). Let

$f\left(a,\gamma \right)=-\frac{\mu }{2}a{s}_g+\frac{F}{2}\sin \gamma ,$

and

$h\left(a,\gamma \right)=\sigma +\frac{A_0}{2}-\frac{3\beta a^2}{8}-\frac{\mu }{2}{c}_g+\frac{F}{2a}\cos \gamma .$

Then the Jacobian is [1] [2] [15] [16] [19] [25] [26] [28]

$J=\left(\begin{array}{cc}-\frac{\mu }{2}{s}_g& \frac{F}{2}\cos \gamma \\ -\frac{3\beta a}{4}-\frac{F}{2a^2}\cos \gamma & -\frac{F}{2a}\sin \gamma \end{array}\right).$ (40)

The solution is classified as locally asymptotically stable if both eigenvalues of $J$ have negative real parts. For this two-dimensional continuous system, the equivalent criterion is [1] [2] [19] [25] [26] [28]

$\text{tr}\left(J\right)<0,\det \left(J\right)>0.$ (41)

The stable and unstable branches are plotted using filled and open markers, respectively.

4.4. Consistency of the Analytical and Computational Formulations

The numerical calculations are performed using the same amplitude equation, phase relations, and slow-flow Jacobian derived in Sections 3 and 4. Equation (38) is used to calculate the positive steady-state amplitudes, Equations (34) and (35) are used to reconstruct the corresponding phases, and the Jacobian of the modulation system is used to classify the local stability of each solution. Figures 2-5 are obtained by algebraic continuation of the steady-state equation, whereas Figure 6 uses direct RK4 integration only in the classical limiting case $g=1$ . Maintaining the same phase convention and slow-flow formulation throughout the analysis prevents inconsistencies in the stability classification [16] [19].

5. Numerical Results and Discussion

This section presents the steady-state frequency-response results obtained from the algebraic amplitude equation derived in Section 4. For each prescribed value of the detuning parameter $\sigma$ , all positive real roots of Equation (38) are calculated over the amplitude interval $0<a\le 8$ . The local stability of each root is then determined from the eigenvalues of the Jacobian matrix in Equation (40). Filled markers denote locally asymptotically stable solutions, whereas open markers denote unstable solutions [1] [2] [15] [16] [19] [23] [25] [26] [28].

The calculations in Figures 2-5 are based on the long-time harmonic representation of the fractional derivative. Therefore, these figures describe the steady-state response predicted by the modulation equations and should not be interpreted as direct time-domain simulations of the Caputo model for $0<g<1$ . Direct numerical integration of the fractional memory equation would require a dedicated fractional time-stepping method [6] [11]-[16] [19] [31] [32].

Unless otherwise stated, the reference parameters are

$\mu =0.2,\beta =0.4,A_0=0.5,F=\mathrm{1.5.}$ (42)

The detuning parameter is varied over

$-2\le \sigma \le 2,$

and the positive amplitude roots are searched over

$0<a\le 8.$

These intervals are numerical search domains and do not constitute restrictions on the underlying mathematical model.

5.1. Reference Frequency Response

Figure 2 presents the reference steady-state frequency response for the fixed fractional order $g=0.7$ , with $\mu =0.2$ , $\beta =0.4$ , $A_0=0.5$ , and $F=1.5$ . The plotted points are obtained from the positive real roots of Equation (38). Filled markers denote locally asymptotically stable solutions, whereas open markers denote unstable solutions.

Figure 2. Reference steady-state frequency response for $g=0.7$ . Filled and open markers denote stable and unstable solutions, respectively.

For $\beta >0$ , the resonance curve displays the characteristic hardening-type behavior of the Duffing oscillator. The nonlinear frequency response bends toward higher detuning values as the amplitude increases. Depending on the detuning, more than one positive steady-state amplitude may exist. The simultaneous presence of stable and unstable branches indicates the possibility of coexisting responses and jump transitions under a slowly varying excitation frequency [1] [2] [15] [16] [19] [23].

It is important to distinguish between the algebraic existence of a root and its stability. Equation (38) determines the possible steady-state amplitudes, whereas the Jacobian criterion in Equation (41) determines the local stability of the corresponding modulation solution. Consequently, a complete frequency-response diagram requires both root calculation and stability classification [1] [2] [19] [25] [26] [28].

5.2. Effect of the Fractional Order

Figure 3 compares the steady-state frequency responses for $g=1.0$ , 0.7, and 0.5, while $\mu =0.2$ , $\beta =0.4$ , $A_0=0.5$ , and $F=1.5$ are kept fixed. The stable and unstable branches are identified from the eigenvalues of the Jacobian matrix of the modulation system.

Figure 3. Effect of the fractional order $g$ on the steady-state frequency response. Filled and open markers denote stable and unstable branches, respectively.

The effect of $g$ enters the response equation through the two quantities

$s_g=\text{sin}\left(\frac{\pi g}{2}\right),c_g=\text{cos}\left(\frac{\pi g}{2}\right).$

The first quantity modifies the effective dissipative contribution, while the second modifies the effective detuning contribution [6] [11] [14]-[16] [19]. Consequently, a change in $g$ may alter both the response amplitude and the horizontal position of the resonance curve. The effect should therefore not be interpreted solely as a monotonic change in damping strength [6] [14]-[16] [19].

At $g=1$ ,

$s_g=1,c_g=0,$

and the fractional term reduces to the classical viscous damping contribution. At $g=0.7$ and $g=0.5$ , both $s_g$ and $c_g$ contribute to the response. Thus, the fractional cases differ from the classical case through a combined modification of dissipation, phase, and effective frequency shift [4] [6] [9]-[14].

The changes in the stable and unstable branches demonstrate that the fractional order can influence not only the peak response but also the topology of the steady-state solution set. In particular, the extent of regions containing multiple positive roots may change as $g$ varies. These features are relevant to vibration control in systems with memory-dependent damping and fractional-order resonators [14] [17]-[21] [25]-[29].

As shown in Figure 3, varying $g$ changes the resonance amplitude and may also shift the response curve along the detuning axis. These changes result from the simultaneous modification of the effective dissipative and detuning contributions. Therefore, the fractional order acts as a response-shaping parameter rather than as a damping coefficient alone [6] [14]-[16] [19].

5.3. Effect of the Fractional Damping Coefficient

Figure 4 examines the influence of the fractional damping coefficient $\mu$ on the steady-state response. The values of $g$ , $\beta$ , $A_0$ , and $F$ are held fixed, while $\mu$ is varied over the values specified in the numerical code.

Figure 4. Effect of the fractional damping coefficient $\mu$ on the steady-state frequency response for $\mu =0.05$ , 0.10, 0.15, and 0.20. The fixed parameters are $g=0.7$ , $\beta =0.4$ , $A_0=0.5$ , and $F=1.5$ .

The parameter $\mu$ multiplies both fractional contributions in Equation (38). Specifically, it appears in ${\mu }^2{s}_g^2$ and in $\mu c_g$ .

Therefore, changing $\mu$ modifies both the effective dissipative component and the effective reactive or detuning component of the steady-state response. For sufficiently large damping coefficients, the response peak is generally reduced and the resonance curve becomes less sharply pronounced. The precise horizontal displacement of the curve depends on the value of $g$ , because the reactive contribution is proportional to $c_g$ [3] [6] [7] [14] [19].

The stability classification in Figure 4 is obtained from the Jacobian of the same modulation system used to calculate the amplitude roots. This consistency is important: the stable and unstable branches cannot be inferred reliably from the amplitude curve alone and must be determined from the linearized slow-flow equations [1] [2] [15] [16] [19] [25] [26] [28].

5.4. Effect of the Excitation Amplitude

Figure 5 illustrates the effect of the excitation amplitude $F$ on the steady-state response. In this calculation, $\mu$ , $\beta$ , $A_0$ , and $g$ are fixed, while $F$ is varied over the selected values.

Figure 5. Effect of the excitation amplitude $F$ on the steady-state frequency response for $F=0.5$ , 1.0, 1.5, and 2.0. The fixed parameters are $\mu =0.2$ , $\beta =0.4$ , $A_0=0.5$ , and $g=0.7$ .

The forcing amplitude appears on the right-hand side of the amplitude equation through $F^2$ . An increase in $F$ generally raises the response amplitude and can enlarge the portion of the resonance curve affected by the cubic nonlinearity. For $\beta >0$ , this may make the hardening-type bending more pronounced and may alter the range of detuning values over which multiple steady-state solutions coexist [1] [2] [15] [16] [19].

The effect of $F$ is not a fractional effect by itself; rather, it reveals how external forcing interacts with the fractional damping and nonlinear stiffness. The resulting response is governed by the combined balance represented in Equation (38) [1] [2] [15] [16] [19].

5.5. Classical-Limit Consistency Check

Figure 6 compares the analytical frequency-response relation with direct fourth-order Runge-Kutta integration in the classical limiting case $g=1$ . In this limit, ${}^C{D}_t^{1}q=\dot{q}$ , the Caputo derivative reduces to the ordinary first derivative, and the governing equation becomes a classical viscously damped Duffing equation [12] [13].

$\ddot{q}+q+\epsilon \mu \dot{q}+\epsilon \beta q^3+\epsilon \frac{A_0}{2}q=\epsilon F\cos \left(\Omega t\right),$ (43)

Therefore, the RK4 calculation is applied only to this ordinary differential equation and not directly to the fractional model for $0<g<1$ .

Equation (43) is integrated using the classical fourth-order Runge-Kutta method with zero initial conditions, a time step of $\text{Δ}t=0.01$ , and a final integration time of $t_{\text{end}}=800$ . The steady-state amplitude is determined from the final 20% of the computed time history.

Figure 6. Consistency check of the analytical frequency-response relation against direct fourth-order Runge-Kutta integration of the classical ordinary differential equation for $g=1$ . The RK4 calculation uses zero initial conditions, $\text{Δ}t=0.01$ , and $t_{\text{end}}=800$ .

The comparison in Figure 6 is therefore a consistency check for the analytical response relation in the classical limit $g=1$ , rather than a direct time-domain validation of the fractional model for $0<g<1$ . Direct fractional validation would require a time-integration method that preserves the history-dependent contribution of the Caputo operator, such as a convolution-based or predictor-corrector method [32].

The numerical amplitude is estimated from the final 20% of the computed time history using the peak-to-peak definition,

$a_{\text{RK}4}=\frac{\max \left(q_{\text{tail}}\right)-\min \left(q_{\text{tail}}\right)}{2},$ (44)

A time-step refinement test was performed using $\text{Δ}t=0.02$ , 0.01, 0.005, and 0.0025, while keeping $t_{\text{end}}=800$ and all model parameters fixed. Taking $\text{Δ}t=0.0025$ as the reference, the relative amplitude error at $\text{Δ}t=0.01$ was $1.53\times {10}^{-6}\text{\%}$ , this confirms that $\text{Δ}t=0.01$ is adequate for the present RK4 calculation.

The analytical amplitude is obtained from the positive real roots of Equation (39). Since the time integration is initialized at $q\left(0\right)=\dot{q}\left(0\right)=0$ , the trajectory may converge to one dynamically selected stable branch when multiple algebraic solutions coexist. Therefore, the RK4 points are interpreted as a consistency check for the classical limiting case $g=1$ , rather than as a simultaneous validation of all algebraic branches.

5.6. Discussion of the Numerical Procedure

The response curves in Figures 2-5 are generated by solving the algebraic equation pointwise in $\sigma$ . This approach is adequate for locating the positive roots of the steady-state equation over the selected parameter range. However, the procedure is more precisely described as a grid-based root continuation rather than a full pseudo-arclength continuation method [19] [23] [25] [26] [28].

A true pseudo-arclength method would track a branch through folds by treating both amplitude and detuning as unknowns and introducing an additional arclength constraint. Such a method is often useful when folds or turning points must be followed without relying on a predefined detuning grid. In the present work, the algebraic root search is retained because the response equation is scalar in the amplitude for each prescribed $\sigma$ . The method is therefore sufficient for the plotted parameter ranges, but this computational limitation should be acknowledged [19] [23] [25] [26] [28].

The root-search resolution may also affect the detection of saddle-node points. A sign-change search can miss a tangential root where the residual touches zero without changing sign. The implementation therefore includes near-zero grid-point detection; nevertheless, the precise location of a fold should preferably be obtained by solving

$\Phi \left(a;\sigma ,g\right)=0,\frac{\partial \Phi }{\partial a}=0,$ (45)

simultaneously. This refinement is recommended if the manuscript makes quantitative claims about saddle-node bifurcation boundaries [19] [23] [25] [26] [28].

5.7. Summary of Numerical Findings

The numerical results may be summarized as follows:

1) The positive cubic coefficient $\beta >0$ produces hardening-type bending of the steady-state resonance curve [1] [2].

2) Multiple positive roots of the amplitude equation may occur over portions of the detuning interval, indicating coexisting steady-state responses [1] [2] [15] [16] [19].

3) The fractional order $g$ modifies both the dissipative and reactive contributions through $\text{sin}\left(\pi g/2\right)$ and $\text{cos}\left(\pi g/2\right)$ , respectively [6] [14]-[16] [19].

4) The damping coefficient $\mu$ affects the response amplitude and the effective frequency shift because it multiplies both fractional contributions [3] [6] [7] [14] [19].

5) The excitation amplitude $F$ controls the overall response level and influences the prominence of the nonlinear bending [1] [2] [15] [16] [19].

6) Stability labels are obtained from the Jacobian eigenvalues of the modulation system and are not assigned solely from the visual shape of the response curve [1] [2] [19] [25] [26] [28].

7) The RK4 comparison in Figure 6 is restricted to the classical limit $g=1$ and is not a direct numerical validation of the fractional model for $0<g<1$ [12] [31] [32].

6. Conclusions, Limitations, and Future Work

6.1. Conclusions

This work analyzed the primary-resonance response and local stability of a reduced Duffing-type oscillator with fractional-order damping, cubic nonlinear stiffness, harmonic excitation, and an effective stiffness correction associated with a spatially periodic medium. The analytical treatment was based on the method of multiple scales and the long-time harmonic representation of the fractional derivative [1] [2] [6] [15] [16] [19] [24].

The main conclusions are as follows:

1) The method of multiple scales provides a systematic approximation for deriving the slow amplitude-phase dynamics near primary resonance [1] [2] [15] [16].

2) The fractional-order contribution enters the modulation equations through the spectral factors

$\text{sin}\left(\frac{\pi g}{2}\right)\text{and}\text{cos}\left(\frac{\pi g}{2}\right).$

The former modifies the effective dissipative contribution, whereas the latter modifies the effective detuning or reactive contribution [6] [11] [14]-[16] [19].

3) The positive cubic coefficient $\beta >0$ produces hardening-type behavior, reflected by the bending of the steady-state response toward higher detuning values. This is consistent with the classical behavior of hardening Duffing oscillators [1] [2].

4) The steady-state response is determined by the positive real roots of the algebraic amplitude equation. Multiple roots may occur for the same detuning value, indicating the possible coexistence of different steady-state responses [1] [2] [15] [16] [19] [23].

5) The stability classification based on the Jacobian eigenvalues distinguishes stable and unstable branches of the modulation system. Therefore, the response topology cannot be described adequately by plotting amplitude roots alone [1] [2] [19] [25] [26] [28].

6) Varying the fractional order $g$ changes both the amplitude and the position of the resonance response. The effect of $g$ should not be interpreted as a simple monotonic increase or decrease in damping strength because the fractional contribution modifies dissipation and effective detuning simultaneously [6] [14]-[16] [19].

7) Increasing the fractional damping coefficient $\mu$ changes both terms associated with fractional damping in the amplitude equation. Its influence therefore includes both amplitude reduction and a possible displacement of the resonance curve [3] [6] [7] [14] [19].

8) Increasing the forcing amplitude $F$ raises the overall response level and enhances the influence of the cubic nonlinear stiffness [1] [2] [15] [16] [19].

9) The comparison with direct fourth-order Runge-Kutta integration is restricted to the classical limiting case $g=1$ . In this limit, the Caputo derivative becomes the ordinary first derivative, and the fractional model reduces to a classical viscously damped Duffing equation [1] [2] [4] [12] [13].

10) The numerical results for $0<g<1$ are obtained from the steady-state harmonic approximation and algebraic continuation. They are not obtained from direct time-domain integration of the Caputo memory equation [12] [31] [32].

6.2. Scope and Limitations

The present formulation has several limitations that should be stated explicitly.

First, the spatially periodic medium is represented by an effective stiffness correction $A_0/2$ within a single-degree-of-freedom reduction. Thus, the model does not resolve the complete spatial field, dispersion relation, mode coupling, or energy transport of the original periodic medium. A full wave or lattice model would be required to study these effects directly [22] [24].

Second, the fractional derivative is treated analytically through its long-time harmonic representation,

${}^C{D}_t^g{\text{e}}^{it}~{\left(i\right)}^g{\text{e}}^{it}.$

This approximation is appropriate for the steady-state harmonic response considered here, but it does not describe the complete transient memory evolution. In particular, the numerical continuation results for $0<g<1$ should not be interpreted as direct numerical solutions of the Caputo initial-value problem [4] [6] [11]-[14] [32].

Third, the RK4 calculation in Figure 6 is performed only for $g=1$ . Since classical RK4 is designed for ordinary differential equations, it cannot by itself provide a direct numerical solution of the fractional equation for $0<g<1$ . A direct fractional validation would require a method that retains the history dependence of the Caputo operator, such as a convolution-based or predictor-corrector method [12] [31] [32].

Fourth, the root calculation is performed over finite numerical intervals in $\sigma$ and $a$ . Roots outside these ranges are not represented in the figures. In addition, a grid-based sign-change method may require refinement near tangential roots and fold points. Quantitative determination of saddle-node boundaries should therefore be performed by solving the simultaneous conditions [19] [23] [25] [26] [28]

$\Phi \left(a;\sigma ,g\right)=0,\frac{\partial \Phi }{\partial a}=0.$

Finally, the analysis is restricted to primary resonance and moderate response amplitudes. Strongly nonlinear regimes, secondary resonances, chaotic responses, and transient phenomena require additional analysis and may not be captured by the first-order modulation equations [1] [2] [14]-[18] [29].

6.3. Future Work

Future work will focus first on direct time-domain validation of the fractional cases $0<g<1$ . The Caputo memory term should be discretized using a suitable fractional numerical scheme, and the resulting periodic amplitudes should be compared with the steady-state response predicted by the harmonic approximation [12] [31] [32].

A second direction is to replace the reduced single-degree-of-freedom model with a distributed periodic lattice or continuum formulation in order to examine dispersion, mode coupling, localization, and band-gap effects explicitly [22] [24]. More accurate continuation methods, such as pseudo-arclength continuation, may also be used to trace response branches through folds and determine saddle-node points with higher accuracy [19] [23] [25] [26] [28].

Finally, experimental validation using a viscoelastic oscillator, a fractional-order resonator, or a metamaterial-inspired mechanical element would help assess the physical interpretation of the effective fractional order and its role in frequency-dependent damping and vibration control. Extensions to secondary resonances, higher-order perturbation corrections, and alternative fractional constitutive formulations may be considered after these primary numerical and experimental validations [3] [6] [7] [14] [20] [21].

The present results establish a consistent analytical framework for examining the primary-resonance response and local stability of a Duffing-type oscillator with fractional-order damping. The principal contribution is the identification of the dual role of the fractional order in modifying both effective dissipation and effective detuning. The classical RK4 comparison confirms consistency only in the limiting case $g=1$ , while the fractional-order response for $0<g<1$ remains based on the long-time harmonic approximation. This distinction clarifies the scope of the results and provides a well-defined basis for future direct fractional numerical and experimental validation [1] [2] [6] [12] [14]-[16] [19] [31] [32].

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References