Acessibilidade / Reportar erro

Fractional Dynamics: A Comprehensive Exploration of Non-integer Order Systems

Abstract

This article delves into the applications of fractional calculus, an extension of classical calculus that introduces non-integer derivative orders. The primary focus of this research is to present a methodology for simulating fractional differential equations and to explore the effects of fractional order in two well-known dynamic systems: the Van der Pol and Duffing systems. These systems are known for their nonlinear characteristics and, in certain cases, exhibit complex and rich dynamic behaviors. Initially, the Grunwald-Letnikov definition of fractional derivatives is introduced, followed by the general numerical solution for a fractional differential equation. The Van der Pol system is modified by the inclusion of a fractional time derivative of order q, reducing the integer order of the system to 1+q, while the Duffing system is modified in terms of viscous damping, by the add of a fractional damping, which is now related to the fractional variation in displacement. The dynamics of the systems are characterized using classical methods of nonlinear dynamics, such as time-response, Poincaré sections, bifurcation diagrams and fast Fourier transform, as well as more advanced approaches, such as the continuous wavelet transform (CWT) and the Hilbert-Huang transform (HHT).

Keywords:
Fractional calculus; Nonlinear dynamics; Time-frequency analysis


1. Introduction

Fractional calculus (FC) presents itself as a generalization of integration and differentiation for a fundamental operator of non-integer order, that is, it refers to differential calculus, but with derivatives and integrals of fractional order, being an extension of traditional calculus, represented by aDtq where t and a are the limits of the operation and q belongs to the set of real numbers [1[1] I. Podlubny, Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications (Elsevier, San Diego, 1998), 1 ed.]. The concept of differentiation and integration of non-integer order is not new, Fig. 1 presents the timeline of important discovers of fractional calculus. This subject has been proposed since the beginning of classical calculus, as Leibniz mentions in a letter to L’Hopital in 1695, where he questioned the interpretation of the meaning of the expression d12x.

Figure 1
Fractional calculus timeline.

Riemann proposed a different definition which involved a definite integral and was applicable to power series with non-integer exponents [2[2] K.S. Miller, An introduction to the fractional calculus and fractional differential equations (John Wiley & Sons, New York, 1993), 1 ed.]. Grunwald [3[3] A.K. Grunwald, Zangew Math und Phys. 12, 441 (1867).] and Letnikov [4[4] A.V. Letnikov, Sbornik: Mathematics 3, 1 (1868).], independently developed an approach to derivatives involving non-integer numbers. This approach is expressed in terms of a convergent series [5[5] V. Lakshmikantham and A.S. Vatsala, Nonlinear Analysis: Theory, Methods & Applications 69, 2677 (2008).]. Krug, working through Cauchy’s integral formula for ordinary derivatives, showed that Riemann’s definite integral had to be interpreted as having a finite lower bound, while Liouville’s definition, in which no distinguishable lower bound appeared, corresponded to a lower bound [6[6] K.B. Oldham, Advances in Engineering software 41, 9 (2010).], where the breakthrough led to numerous other definitions, some of which can be found in [7[7] E.C.D. Oliveira and J.A. Tenreiro, Mathematical Problems in Engineering 2014, 38459 (2014).].

Until the 1970s, a modest number of papers were published on the subject of FC, where a chronology documenting the period 1965-1974 can be seen in [6[6] K.B. Oldham, Advances in Engineering software 41, 9 (2010).], after which Miller and Ross covered more current issues until the 1990s [2[2] K.S. Miller, An introduction to the fractional calculus and fractional differential equations (John Wiley & Sons, New York, 1993), 1 ed.]. The study of the theory of fractional differential equations is relatively new; this is due to the intense computational work required to simulate this type of system due to the effect known as memory [1[1] I. Podlubny, Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications (Elsevier, San Diego, 1998), 1 ed., 8[8] I. Petráš, Fractional-order nonlinear systems: modeling, analysis and simulation (Springer Science & Business Media, Berlim, 2011), 1 ed.]. However, interest in this tool has not waned due to its applicability in various academic disciplines. Thus, differential equations involving fractional differential operators become important in the modeling of various physical and dynamic phenomena, since characterization and identification are one of the objectives in the study of dynamic systems.

For example, fractional calculus is extensively reviewed in its application to controllers, heat diffusion and digital circuits in [9[9] J.A. Tenreiro, M.F. Silva, R.S. Barbosa, I.S. Jesus, C.M. Reis, M.G. Marcos and A.F. Galhano, Mathematical problems in engineering 2010, 639801 (2010).], coupled nonlinear differential equations of arbitrary orders (integer and fractional orders) are implemented with electronic circuits by means of analogical integrators, adders, multipliers and the fractance circuits units with capacitors/resistors in [10[10] S.A. David, C. Fischer and C.D. Oliveira, in: COUPLED PROBLEMS 2019 At: SITGES (Barcelona, 2019).]. Its application to mass-spring-damper systems is discussed in [11[11] N. Sene and J.F.Gómez, Fractional mass-spring-damper system described by generalized fractional order derivatives 3, 39 (2019)., 12[12] S.S. Ray, S. Sahoo and S. Das, Advances in Mechanical Engineering 8, 1 (2016).]. Various formulations of the Van der Pol oscillator, with the implementation of fractional calculus, are analyzed in [13[13] R.S. Barbosa, J.A. Tenreiro, B.M. Vinagre and A.J. Calderon, Journal of Vibration and Control.13, 1291 (2007)., 14[14] J.A. Tenreiro, A.M. Lopes and V.E. Tarasov, Applications in Physics, Part A. 4, 1 (2019)., 15[15] I.H. Outram, Frequency Pulling of the van der Pol Oscillator Doctoral thesis, McMaster University, Hamilton (1967)., 16[16] M. Tsatsos, Theoretical and Numerical Study of the Van der Pol Equation. Doctoral Thesis, Aristotle University of Thessaloniki, Thessaloniki (2006)., 17[17] M. Borowiec, G. Litak and A. Syta, arxiv.org/abs/nlin/ 0601033 (2007).]. The dynamics of the Duffing system with fractional damping is presented in [18[18] M.V. Varanis, A.M. Tusset, J.M. Balthazar, G. Litak, C.D. Oliveira, R.T. Rocha, A. Nabarrete and V. Piccirillo, Journal of the Franklin Institute 357, 2067 (2020).]. Finally, [19[19] W.M. Ahmad and J.C. Sprott, Chaos, Solitons & Fractals 16, 339 (2003).] offers a more comprehensive study of the impact of fractional calculus on a wide variety of engineering problems.

The application of fractional calculus in dynamical systems fits the nonlinear dynamics problem class. When addressing nonlinear systems, whose identification and characterization still lack a universal and comprehensive theory, we encounter several limitations. For instance, temporal analysis presents significant challenges since nonlinear systems do not follow a direct cause-and-effect relationship, making it difficult to accurately predict the system’s behavior over time. Additionally, these systems are highly sensitive to initial conditions, where small inaccuracies or uncertainties can lead to significant deviations in results over time. A perspective on nonlinear dynamics can be gained from bifurcation theory, which provides an approximate analysis of mathematical models but may oversimplify the actual system’s behavior, failing to capture all aspects of the dynamics. Poincaré sections, which reduce a continuous dynamic system to a discrete map, can also simplify the analysis, but this may result in the loss of important information.

In this context, an alternative perspective can be gained through the use of time-frequency tools. The use of analysis techniques such as Continuous Wavelet Transform (CWT) and Hilbert-Huang Transform (HHT) emerges as a promising solution. The CWT, as described by [20[20] R. Benítez, V.J. Bolós and M.E. Ramírez, Computers & Mathematics with Applications 60, 634 (2010)., 21[21] J. Chen, Z. Li, J. Pan, G. Chen, Y. Zi, J. Yuan, B. Chen and Z. He, Mechanical systems and signal processing 70, 1 (2016).], is a time-frequency analysis technique that excels in detecting short-duration transients and analyzing signals with non-stationary dynamics. The HHT, in turn, introduces the concept of instantaneous frequency, serving as a complementary and alternative tool to the CWT, as discussed by [22[22] N.E. Huang, Hilbert-Huang transform and its applications (World Scientific, Singapura, 2014), 2 ed., 23[23] N.E. Huang, Interdiscip. Math. 5, 1 (2005).].

This paper focuses on presenting the basic theory of fractional differential equations involving fractional derivatives by the Grunwald-Letnikov (GL) definition, the application of fractional derivatives in nonlinear dynamic systems such as the Van der Pol system and the Duffing oscillator, in order to present a methodology for simulating fractional differential equations. The changes in system dynamics due to fractional order will also be observed using temporal tools such as phase space, Poincaré sections, bifurcations and in the frequency domain using fast Fourier transforms (FFT), continuous wavelet transforms (CWT) and Hilbert-Huang transforms (HHT) to explore non-linearity, chaos and periodicity in these dynamic systems.

2. Fractional Derivative and Solution of a Fractional ODE

Considering a continuous function f(t), its derivative can be defined as

(1) d d t f ( t ) = lim h 0 f ( t ) f ( t h ) h

This definition can be extended to any n-derivative to obtain the definition of Grunwald-Letnikov (GL), first generalize the derivative of order n for any natural integer value from an alternating series that follows the binomial coefficients [1[1] I. Podlubny, Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications (Elsevier, San Diego, 1998), 1 ed., 8[8] I. Petráš, Fractional-order nonlinear systems: modeling, analysis and simulation (Springer Science & Business Media, Berlim, 2011), 1 ed.],

(2) f ( n ) ( t ) = lim h 0 1 h n j = 0 n ( 1 ) j ( j n ) f ( t j h )

where,

(3) ( j n ) = n ( n 1 ) . . . ( n j + 1 ) j ! = n ! j ! ( n j ) !

this equation, in turn, can be extended to the positive real numbers from the Euler gamma function, given by,

(4) Γ ( n ) = 0 t n 1 e t d t

which is a generalization of the factorial as follows,

(5) Γ ( n ) = ( n 1 ) !

thus obtaining a method for determining the binomial coefficients for any real value greater than zero for the order of the derivative q,

(6) ( j q ) = q ! j ! ( q j ) ! = Γ ( q + 1 ) Γ ( j + 1 ) Γ ( q j + 1 )

This defines the fractional order derivative of GL, where a is the memory length, t is the time and h is the time step,

(7) D t q a = lim h 0 1 h q j = 0 [ t a h ] ( 1 ) j ( j q ) f ( t j h )

For the numerical calculation of fractional order derivatives, using the relation below, where applying the short memory principle [1[1] I. Podlubny, Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications (Elsevier, San Diego, 1998), 1 ed., 8[8] I. Petráš, Fractional-order nonlinear systems: modeling, analysis and simulation (Springer Science & Business Media, Berlim, 2011), 1 ed.], which states that for a large time t, the behavior of the function f(t) near the lower terminal t = a can be disregarded under certain assumptions, such as taking into account only the recent past, which is the interval [tLm,t] obtaining the expression below,

(8) D t k q ( k L m / h ) f ( t ) h q j = 0 k ( 1 ) j ( j q ) f ( t k j )

Where Lm is the memory length, for tk=kh, and the equation below presents a method for calculating the binomial coefficients,

(9) c 0 ( q ) = 1 , c j ( q ) = ( 1 1 + q j ) c j 1 ( q )

Given this, the general numerical solution of a first-order fractional differential equation can be derived [8[8] I. Petráš, Fractional-order nonlinear systems: modeling, analysis and simulation (Springer Science & Business Media, Berlim, 2011), 1 ed.], thus obtaining equation (10):

(10) y ( t k ) = f ( y ( t k ) , t k ) h q j = v k c j ( q ) y ( t k j )

3. Formulation of Fractional Dynamic Systems

This article explores two well-known dynamical systems in the field of non-linear oscillations: the Duffing system and the Van der Pol oscillator. The focus is on their fractional models and the various applications of fractional derivatives. In the Van der Pol system, the total order of the system is reduced to 1 + q, while in the Duffing system, fractional damping is applied.

3.1. Fractional Van der Pol

The Van der Pol model is a self-excited oscillator model with limit cycles. Regardless of what the initial conditions of the system are, the oscillator always converges to an isolated cycle, its standard form is shown by the equation below.

(11) x ¨ + α ( x 2 1 ) x ˙ + x = δ c o s ( ω t )

Analyzing the equation, is possible see the similarity to a harmonic oscillator with a non-linear damping term α(x21)x˙, where α is the non-linear parameter. This nonlinear term has positive damping for |x|>1 and negative damping for |x|<1. From this it’s possible see that, for cases where the amplitude is very small, the damping increases and for cases where the amplitude is very large, the damping decreases, causing the system to remain in a cycle [24[24] S.H. Strogatz,Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering (CRC Press, Florida, 2000), 1 ed.]. In addition, there are the parameters of the external force, where δ is the amplitude of the external force and ω is its excitation frequency.

The fractional model of the Van der Pol oscillator was initially proposed by [25[25] R.S. Barbosa, J.A. Tenreiro, I.M. Ferreira and J.K. Tar, in: In Anais Second IEEE International Conference on Computational Cybernetics (Vienna, 2004).], where it was later revisited, and the addition of an external force was addressed by [13[13] R.S. Barbosa, J.A. Tenreiro, B.M. Vinagre and A.J. Calderon, Journal of Vibration and Control.13, 1291 (2007).]. The modification is the introduction of a fractional time derivative of order q in the system, which due to the representation of the state space is also applied to the next derivative, where the following equation is obtained,

(12) x 1 + q + α ( x 2 1 ) x q + x = 0 0 < q < 1

applying the state space formulation obtain,

(13) { L D t q x ( t ) = y ( t ) L D t 1 y ( t ) = α ( x 2 ( t ) 1 ) y ( t ) x ( t ) + δ c o s ( ω t ) .

The discretization of this system of equations into first-order approximations using the Forward Euler integration scheme for integer integrals and equation (10) for fractional integrals is given by equations (1415).

(14) x ( t k ) = y ( t k 1 ) h q j = 1 N 1 c j ( q ) x ( t k j )

(15) y ( t k ) = y ( t k 1 ) + [ α ( x 2 ( t k 1 ) 1 ) y ( t k 1 ) x ( t k 1 ) + δ c o s ( ω ( t k 1 ) ) ] h

3.2. Fractional duffing

The Duffing oscillator is a second-order non-linear differential equation commonly used to model non-linear (cubic) stiffness that deviates from Hooke’s law. It is typically employed to observe the system’s real behavior in natural settings. The model is expressed by equation (14).

(16) x ¨ + β x ˙ x + x 3 = δ c o s ( ω t )

In this system, the restoring force is described by a linear term x, similar to a standard linear oscillator, but with the addition of the term x3, which represents the non-linearity of the system, where the two make up a non-linear restoring force. In this case β is a constant that represents the total damping in the system and like the Van der Pol system presented, δ is the amplitude of the external force and ω the excitation frequency.

The modification presented in the equation is that the viscous damping related to the velocity (the integer temporal rate of change of the displacement) is now related to the fractional rate of change of the displacement as presented by [17[17] M. Borowiec, G. Litak and A. Syta, arxiv.org/abs/nlin/ 0601033 (2007).], thus obtaining equation (17), which can be written in the form of the following system of equation (18).

(17) d 2 x d t 2 + β d q x d t q x + x 3 = δ c o s ( ω t )

(18) { L D t 1 x ( t ) = y ( t ) L D t q x ( t ) = w ( t ) L D t 1 y ( t ) = x ( t ) β w ( t ) x 3 ( t ) + δ c o s ( ω t ) .

The discretization of this system of equations, in first-order approximations using the Forward Euler integration scheme for the integer integrals and equation (10) for the fractional equations, can be seen by the equations (1921).

(19) x ( t k ) = x ( t k 1 ) + y ( t k 1 ) h

(20) x ( t k ) = w ( t k 1 ) h q j = 1 N 1 c j ( q ) x ( t k j )

(21) y ( t k ) = y ( t k 1 ) + [ x ( t k 1 ) β w ( t k 1 ) x 3 ( t k 1 ) + δ c o s ( ω ( t k 1 ) ) ] h

It is worth noting that for the system to have fractional damping, w is solved by isolating equation (20) and replacing x with the definition given in equation (19), thus obtaining the solution of w given by equation (22).

(22) w ( t k 1 ) = 1 h q [ x ( t k 1 ) + y ( t k 1 ) h + j = 1 N 1 c j q x ( t k j ) ]

4. Numerical Simulations

For the simulations, the Python programming language with its main scientific libraries NumPy and matplotlib were used, with the add of the Numba library for optimizing the codes, the SciPy library for signal processing with FFTs and CWTs and the EMD library for using the HHT [26[26] A.J. Quinn, V. Lopes-dos-Santos, D. Dupret, A.C. Nobre and M.W. Woolrich, Journal of Open Source Software 6, 2977 (2021).]. It is worth noting that in the numerical simulations, a substitution was made in the parameters, replacing x with x1 and y with x2.

The simulations are presented in the following order: temporal response, phase space with Poincaré sections, bifurcation, FFT, CWT and HHT, with the first responses being purely temporal and the last three being frequency responses. For the frequency responses, it was used the signals presented in the time response with a greater interval, where for the Van der Pol system it was took the time interval around 2900 t 2950 s with a time step of dt = 2πω250 and for the Duffing system 2750 t 2950 s with a time step of dt = 2πω200.

To obtain the bifurcation graphs, the time interval 2100 t 3000 s, was used, where the first 150 peaks of the displacement time series were taken and plotted against a control parameter, so that for cases where the system has a periodic response, it will mark the same point several times on the x-axis (without taking into account the numerical error), and for cases where non periodic/chaotic responses are obtained, this is characterized by several points on the same parameter on the x-axis. In addition, in order to observe the signals over time, the signals were plotted using their amplitude as a colormap, thus helping to observe the transient region and which parameters lead the system to different responses in conjunction with the bifurcation.

To investigate the system in frequency domain, fast Fourier transforms were used to observe the different frequencies and their harmonics. To observe the frequency response over time, it was used continuous Wavelet transforms, which is a different approach in the analysis of non-stationary signals, because they present a response in the time and frequency domain, thus allowing us to observe important phenomena that are not visible in the purely time or purely frequency domain, having good temporal precision and are widespread in signal analysis tools, as seen in [27[27] S. Mallat, A wavelet tour of signal processing (Elsevier, San Diego, 1999), 1 ed., 28[28] M.V. Varanis, C.D. Oliveira, M.A. Ribeiro, W.B. Lenz, A.M. Tusse and J.M. Balthazar, Nonlinear Vibrations Excited by Limited Power Sources 116, 175 (2022)., 29[29] M.V. Varanis, S. Hemmati, M.C. Filipus, F.L.D. Abreu, J.M. Balthazar and C. Nataraj, Advances in Nonlinear Dynamics Volume III, (Springer,Cham, Switzerland, 2024), v. 2., 30[30] F.L.D. Abreu, M.V. Varanis, C.D. Oliveira, J.M. Balthazar and A.M. Tusset, in: In Anais XLIII Ibero-Latin American Congress on Computational Methods in Engineering (Foz do Iguaçu, 2022).].

Finally, as an alternative and complementary approach to CWT, the Hilber-Huang transform will be used. Unlike the two previous tools, the decomposition does not take place in the frequency domain but in the time domain using the sifting algorithm. In this algorithm, the signal is initially decomposed into a set of intrinsic mode functions which are more basic oscillatory functions, which are then transformed into the frequency domain. One of the interesting points of this method is that as the decomposition is not tied to any fixed basis such as the FFT and CWT, in addition to not presenting convolutions, the characterization of the signals in the frequency domain is given without loss of information, as this is not affected by the uncertainty principle. This allows the frequency of the system to be described without the need for harmonics, but from an oscillating frequency over time, which becomes interesting when dealing with non-linear systems, thus creating an instantaneous frequency [31[31] N.E. Huang, Z. Shen, S.R. Long, M.C. Wu, H.H. Shih, Q. Zheng, N.C. Yen, C.C. Tung and H.H. Liu, Proceedings of the Royal Society of London. Series A: mathematical, physical and engineering sciences 454, 903 (1998)., 32[32] N.E. Huang, Z. Wu, S.R. Long, K.C. Arnold, X. Chen and K. Blank, Advances in adaptive data analysis 1, 177 (2009)., 33[33] M.V. Varanis, A.L. Silva, J.M. Balthazar and R. Pederiva, Brazilian Journal of Physics 51, 859 (2021).].

5. Results and Discussions

5.1. Fractional Van der Pol

For the simulations of the fractional Van der Pol system, the damping and external force frequency parameters α = 5 and ω = 2.46 [rad/s], were used, respectively. In this case, the fractional order q of the system will be set to 0.85, where the analysis of the system will be based on the variation of the parameter δ which corresponds to the amplitude of the external force, that will assume the values of: 2.0, 1.5 and 5.5, these being the same parameters as in the analysis of [25[25] R.S. Barbosa, J.A. Tenreiro, I.M. Ferreira and J.K. Tar, in: In Anais Second IEEE International Conference on Computational Cybernetics (Vienna, 2004).].

Figure 2 shows the temporal response of the system. It gives an initial idea of the periodicity of each signal related to a change in the amplitude of the external force, where its possible to note that for δ = 2 and 5.5 cases, the system respond in a periodic regime, where for the δ = 1.5 case, the time response appears to be aperiodic.

Figure 2
Time response (a) δ = 2 (b) δ = 1.5 (c) δ = 5.5.

Figure 3 shows the phase space together with the Poincaré sections for the system. It can be seen that when δ = 1.5, the system exhibits aperiodic behavior, which is evidenced both by the various lines crossing each other in the phase space and by the various scattered points in the Poincaré sections. For values of δ = 2 and 5.5, the behavior is stable, where the phase space shows lines that do not cross, together with the Poincaré sections which show points that remain close, where periodic responses of periods 4 and 1 are observed respectively.

Figure 3
Phase space and Poincaré sections (a) δ = 2 (b) δ = 1.5 (c) δ = 5.5.

By analyzing the bifurcation diagram, Fig. 4, can see that by setting the fractional order parameter to 0.85 and varying the amplitude of the external force, obtain a behavior composed initially of aperiodic solutions in the interval 0 δ 3, where in the middle of these solutions there are periodic windows, and after δ > 3, the responses converge to a periodic solution. To observe the behavior of the system in the time domain, the solutions of x1 were plotted again in the interval 0t3000 s, where the response amplitudes were plotted with a color map, where it is possible to observe the behavior of changing the parameter δ presented in the bifurcation diagram in the time domain.

Figure 4
Varying the parameter δ (a) bifurcation with respect to x1 (b) time response.

A possible explanation of the observed effect can be found in [34[34] M. Cercek, M. Stanojevic and T. Gyergyek, Nuclear Energy in Central Europe 1, 562 (1996).], which presents the phenomenon known as frequency pulling. This phenomenon occurs in non-linear systems when an oscillator has its natural frequency adjusted or “pulled” to coincide with the frequency of a stronger oscillator or an external signal. However, unlike the case observed by the author where the external frequency is altered, here it is the amplitude of the external signal that varies significantly to change the system’s dynamics. What happens is that the external force is not enough to control the system, as shown in the bifurcation of the force. Larger forces lead to a periodic result, and intermediate forces lead to a competition between the natural response and the externalresponse.

In the FFT analysis (Fig. 5) there are small visually noticeable differences from the case of δ = 2 to δ = 1.5, where when analyzing the peaks with the highest amplitude, they have very close values, but when using other analyses, it is possible to identify the change in system dynamics much more clearly. When δ = 2, a periodic region of period 4 is observed, as seen in the time response, the Poincaré diagram and the bifurcation diagram through the region with 4 lines, then when δ = 1.5 this is located in a region that is under the action of both the natural response of the system and the external force, which does not have the amplitude to command the response on its own, thus generating zones of multiple orbits, which are described in the FFT by the small peaks, the problem is that analyzing only the peaks with the greatest amplitude could lead to an incoherent analysis of the response of these signals since the two cases are very close numerically. Finally, when analyzing the parameter δ = 5.5, the main frequency that is presented by the FFT is the frequency of the external force of 0.39 Hz (or 2.46 Rad/s), together with its superharmonics generated due to the non-linearity of the system, it is worth noting that the natural response in this case has been completely suppressed, in addition as presented by the Poincaré section, there is a case of period 1.

Figure 5
Fast Fourier Transform (a) δ = 2 (b) δ = 1.5 (c) δ = 5.5.

Next, the analysis in the Time-Frequency domain is obtained from the CWT, shown in Fig. 6. The observation made by the FFT for the case δ = 5.5 is confirmed, showing a main frequency and harmonics that are invariant in time, but the spectral content of the cases where δ assumes the values of 2 and 1.5 shows characteristics of a non-stationary dynamics. The main characteristic being the oscillating frequency along with its harmonics over time, where it is possible to observe a greater variation in the characteristics of the signals than observed by the FFT.

Figure 6
Continuos Wavelet Transform (a) δ = 2 (b) δ = 1.5 (c) δ = 5.5.

As an alternative to the analysis in the Time-Frequency domain, the HHT shown in Fig. 7 confirms the observations made by the CWT, but through its description in instantaneous frequency it is possible to observe a certain distinct dynamic characteristic for each case, where here the description in frequency shows a big difference from the case of δ = 2 and 1.5. In the case of δ = 2, three separate oscillating frequencies can be identified, but in the case of δ = 1.5 only one oscillating frequency with several peaks can be observed. Furthermore, in this analysis for δ = 5.5 there is only one oscillating frequency to describe the frequency now commanded by the external force and its harmonics.

Figure 7
Hilbert–Huang Transform (a) δ = 2 (b) δ = 1.5 (c) δ = 5.5.

5.2. Duffing with fractional damping

For the simulations of the Duffing system with fractional damping, β = 0.15 and ω = 1.0 [rad/s] were used as the damping and external force frequency parameters, respectively. In this case, the system will be analyzed by varying the q parameter, which corresponds to the order of the fractional derivative related to the damping parameter, which in turn will take on the values of 0.6, 0.8 and 1, as proposed in [35[35] A. Syta, G. Litak, S. Lenci and M. Scheffler, Chaos: An Interdisciplinary Journal of Nonlinear Science 24, 013107 (2014).].

Figure 8 shows the temporal response of the system. For q = 1 the time domain signal presents a non periodic response, characteristic of chaotic systems. However, with the application of fractional damping, as seen in the cases of q = 0.6 and 0.8 the response becomes periodic. This demonstrates that the introduction of fractional damping alters the system’s behavior from non-periodic to periodic.

Figure 8
Time response (a) q = 0.6 (b) q = 0.8 (c) q = 1.

Figure 9 shows the phase space with Poincaré sections. It can be seen that as the fractional term of the system approaches 1, the similarity to the chaotic behavior characteristic of the Duffing model increases, showing an attractor. When q reaches a value of 1, the system faithfully replicates the chaotic behavior predicted by the model, including the presence of numerous periods. Consistent with the time-domain response, the introduction of fractional damping leads the system to exhibit periodic orbits of periods 1 and 2.

Figure 9
Phase space and Poincaré sections (a) q = 0.6 (b) q = 0.8 (c) q = 1.

By analyzing the bifurcation diagram shown in Fig. 10, it can be seen that for values of 0.1 q 0.8, the system is composed of periodic solutions, followed by an aperiodic regime corresponding to the region 0.8 q 1.4, and finally another periodic region where 1.4 q 2. This is related to how the fractional term influences the system. For cases where the fractional term is equal to an integer term, the system tends to behave like a chaotic duffing oscillator, but when it moves away from this integer value the system cannot efficiently reproduce the chaotic behavior of the system.

Figure 10
Varying the order of the fractional derivative (a) bifurcation with respect to x1 (b) time response.

Figure 11 shows the FFT of the system for the 3 cases mentioned, it can be seen that for the system with q values of 0.6 the system has few frequencies due to the distance that the fractional term has in relation to the integer value, making the system have the characteristics of a periodic system, in this case presenting only a few harmonics, which are necessary to compose the Duffing response, which is non-linear. For the case where q = 0.8 the appearance of more frequencies, until at q = 1, there is a grouping of multiple frequencies, showing the chaos presented in the system.

Figure 11
Fast Fourier Transform (a) q = 0.6 (b) q = 0.8 (c) q = 1.

Figure 12 shows the response of the CWT. For q = 0.6 and 0.8 the frequencies remain constant over time, with the frequency in 0.159 Hz presenting the external forcing and weaker superharmonics that apears in both cases to compose the signal, finally in the case of q = 1 it is possible to notice a spread of frequencies, characteristic of chaotic systems.

Figure 12
Continuos Wavelet Transform (a) q = 0.6 (b) q = 0.8 (c) q = 1.

The HHT confirms the observations made earlier, but with an alternative view as shown in Fig. 13. For the case with q = 0.6 a periodic behavior is apparent due to the fractional term, where only 1 frequency oscillates around the external excitation frequency of 0.159 Hz. For q = 0.8, 3 oscillating frequencies appears, where the frequency that oscillates around the external excitation now has three peaks, there is a lower frequency that oscillates with a small variation below and finally an almost imperceptible frequency around zero, due to the non-zero average of the signal. Finally, there is the frequency spread at q = 1, showing the chaos present in the system.

Figure 13
Hilbert–Huang Transform (a) q = 0.6 (b) q = 0.8 (c) q = 1.

6. Conclusion

In summary, a simple methodology for solving fractional differential equations has been presented, starting with the definition of GL, along with the presentation of two nonlinear systems widely discussed in the nonlinear dynamics literature, where the influence of fractional order on their dynamics is explored, altering the total order of the system and by adding an fractional damping. In addition, these systems were simulated using the equations briefly introduced and characterized using classical methods of nonlinear dynamics, both in the time domain using phase space, Poincaré sections and bifurcation diagrams and in the frequency domain using FFT and emergent methods such as CWTand HHT.

With regard to the systems presented, it is worth noting the difference in the mathematical description of the two cases analyzed. The fractional Van der Pol system has a total order of 1 + q, due to its derivative not being taken into account separately, where the second derivative is of integer order. The Duffing with fractional damping has a total order of 2, because only the damping has a fractional order, so it is necessary to create an extra equation to take into account both the first and second derivatives of integer order, as well as the derivative of fractional order, necessary to simulate the system.

In terms of characterizing the dynamics, these methods proved to be excellent qualitative indicators for characterizing both the Van der Pol and Duffing systems. It is worth noting that the Poincaré diagram was able to characterize the orbits presented well, and in terms of frequency, the tools presented were also able to extract good information. As a highlight, in the analysis of the Van der Pol system, the characterization carried out by the HHT was able to better distinguish the responses presented and may be more suitable for extracting signal characteristics.

Finally, observing the “frequency pulling” phenomenon in the Van der Pol oscillator offers insights into the distinct dynamic behaviors observed across different regimes, particularly when considering variations in the amplitude of the external signal. In the case of the Duffing system, the introduction of fractional damping demonstrated effective chaos control as the order of the damping was modified. The methods presented in this study have proven effective and suggest that they can serve as valuable tools for the qualitative characterization of similar dynamical systems.

Acknowledgements

The authors gratefully acknowledge the financial support from CNPq (National Council for Scientific and Technological Development), as well as the support of the National Laboratory for Scientific Computing (LNCC/MCTI, Brazil), through the Ambassador Program (UFGD), MNSE subproject, for providing the HPC resources of the SDumont supercomputer, which contributed to the research results reported in this paper. URL:http://sdumont.lncc.br.

References

  • [1]
    I. Podlubny, Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications (Elsevier, San Diego, 1998), 1 ed.
  • [2]
    K.S. Miller, An introduction to the fractional calculus and fractional differential equations (John Wiley & Sons, New York, 1993), 1 ed.
  • [3]
    A.K. Grunwald, Zangew Math und Phys. 12, 441 (1867).
  • [4]
    A.V. Letnikov, Sbornik: Mathematics 3, 1 (1868).
  • [5]
    V. Lakshmikantham and A.S. Vatsala, Nonlinear Analysis: Theory, Methods & Applications 69, 2677 (2008).
  • [6]
    K.B. Oldham, Advances in Engineering software 41, 9 (2010).
  • [7]
    E.C.D. Oliveira and J.A. Tenreiro, Mathematical Problems in Engineering 2014, 38459 (2014).
  • [8]
    I. Petráš, Fractional-order nonlinear systems: modeling, analysis and simulation (Springer Science & Business Media, Berlim, 2011), 1 ed.
  • [9]
    J.A. Tenreiro, M.F. Silva, R.S. Barbosa, I.S. Jesus, C.M. Reis, M.G. Marcos and A.F. Galhano, Mathematical problems in engineering 2010, 639801 (2010).
  • [10]
    S.A. David, C. Fischer and C.D. Oliveira, in: COUPLED PROBLEMS 2019 At: SITGES (Barcelona, 2019).
  • [11]
    N. Sene and J.F.Gómez, Fractional mass-spring-damper system described by generalized fractional order derivatives 3, 39 (2019).
  • [12]
    S.S. Ray, S. Sahoo and S. Das, Advances in Mechanical Engineering 8, 1 (2016).
  • [13]
    R.S. Barbosa, J.A. Tenreiro, B.M. Vinagre and A.J. Calderon, Journal of Vibration and Control.13, 1291 (2007).
  • [14]
    J.A. Tenreiro, A.M. Lopes and V.E. Tarasov, Applications in Physics, Part A. 4, 1 (2019).
  • [15]
    I.H. Outram, Frequency Pulling of the van der Pol Oscillator Doctoral thesis, McMaster University, Hamilton (1967).
  • [16]
    M. Tsatsos, Theoretical and Numerical Study of the Van der Pol Equation Doctoral Thesis, Aristotle University of Thessaloniki, Thessaloniki (2006).
  • [17]
    M. Borowiec, G. Litak and A. Syta, arxiv.org/abs/nlin/ 0601033 (2007).
  • [18]
    M.V. Varanis, A.M. Tusset, J.M. Balthazar, G. Litak, C.D. Oliveira, R.T. Rocha, A. Nabarrete and V. Piccirillo, Journal of the Franklin Institute 357, 2067 (2020).
  • [19]
    W.M. Ahmad and J.C. Sprott, Chaos, Solitons & Fractals 16, 339 (2003).
  • [20]
    R. Benítez, V.J. Bolós and M.E. Ramírez, Computers & Mathematics with Applications 60, 634 (2010).
  • [21]
    J. Chen, Z. Li, J. Pan, G. Chen, Y. Zi, J. Yuan, B. Chen and Z. He, Mechanical systems and signal processing 70, 1 (2016).
  • [22]
    N.E. Huang, Hilbert-Huang transform and its applications (World Scientific, Singapura, 2014), 2 ed.
  • [23]
    N.E. Huang, Interdiscip. Math. 5, 1 (2005).
  • [24]
    S.H. Strogatz,Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering (CRC Press, Florida, 2000), 1 ed.
  • [25]
    R.S. Barbosa, J.A. Tenreiro, I.M. Ferreira and J.K. Tar, in: In Anais Second IEEE International Conference on Computational Cybernetics (Vienna, 2004).
  • [26]
    A.J. Quinn, V. Lopes-dos-Santos, D. Dupret, A.C. Nobre and M.W. Woolrich, Journal of Open Source Software 6, 2977 (2021).
  • [27]
    S. Mallat, A wavelet tour of signal processing (Elsevier, San Diego, 1999), 1 ed.
  • [28]
    M.V. Varanis, C.D. Oliveira, M.A. Ribeiro, W.B. Lenz, A.M. Tusse and J.M. Balthazar, Nonlinear Vibrations Excited by Limited Power Sources 116, 175 (2022).
  • [29]
    M.V. Varanis, S. Hemmati, M.C. Filipus, F.L.D. Abreu, J.M. Balthazar and C. Nataraj, Advances in Nonlinear Dynamics Volume III, (Springer,Cham, Switzerland, 2024), v. 2.
  • [30]
    F.L.D. Abreu, M.V. Varanis, C.D. Oliveira, J.M. Balthazar and A.M. Tusset, in: In Anais XLIII Ibero-Latin American Congress on Computational Methods in Engineering (Foz do Iguaçu, 2022).
  • [31]
    N.E. Huang, Z. Shen, S.R. Long, M.C. Wu, H.H. Shih, Q. Zheng, N.C. Yen, C.C. Tung and H.H. Liu, Proceedings of the Royal Society of London. Series A: mathematical, physical and engineering sciences 454, 903 (1998).
  • [32]
    N.E. Huang, Z. Wu, S.R. Long, K.C. Arnold, X. Chen and K. Blank, Advances in adaptive data analysis 1, 177 (2009).
  • [33]
    M.V. Varanis, A.L. Silva, J.M. Balthazar and R. Pederiva, Brazilian Journal of Physics 51, 859 (2021).
  • [34]
    M. Cercek, M. Stanojevic and T. Gyergyek, Nuclear Energy in Central Europe 1, 562 (1996).
  • [35]
    A. Syta, G. Litak, S. Lenci and M. Scheffler, Chaos: An Interdisciplinary Journal of Nonlinear Science 24, 013107 (2014).

Publication Dates

  • Publication in this collection
    15 Nov 2024
  • Date of issue
    2024

History

  • Received
    14 May 2024
  • Reviewed
    06 Sept 2024
  • Accepted
    21 Sept 2024
Sociedade Brasileira de Física Caixa Postal 66328, 05389-970 São Paulo SP - Brazil - São Paulo - SP - Brazil
E-mail: marcio@sbfisica.org.br