 # Approximate solutions of one dimensional systems with fractional derivative

The fractional calculus is useful to model non-local phenomena. We construct a method to evaluate the fractional Caputo derivative by means of a simple explicit quadratic segmentary interpolation. This method yields to numerical resolution of ordinary fractional differential equations. Due to the non-locality of the fractional derivative, we may establish an equivalence between fractional oscillators and ordinary oscillators with a dissipative term.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The study of fractional derivatives for its application in classical and quantum physics has lately received a lot of attention [1, 2, 3]. Needless to say that one of the simplest and most studied of those systems is the one dimensional harmonic oscillator. Thus, it would be a good point of departure in the study of systems with fractional derivative, a task which has been carried out in . Damped oscillator with fractional derivative has been also the objective of some studies, see . Some extensions of the theory to other classical systems have been proposed, see for instance  and references therein, or in .

In many of these papers, it was noted an analogy between a fractional oscillator and a classical oscillator with a damping term. This could be an idea to be exploited in order to model quantum systems with dissipation, in which the second derivative of the wave function in the Schrödinger equation be replaced by a fractional derivative, see .

The present work has been inspired by the article by Narahari et. al. , where in addition to the study of the one dimensional harmonic oscillator with fractional derivative, they give a comparison with an equivalent dissipative oscillator described on the phase plane and analyze the stability of the solution.

Along the present manuscript, we show that it may be possible the determination of a time interval in which the solution of a fractional one dimensional oscillator may be approximated by the solution of a one dimensional ordinary equation with a dissipative term. The idea could be described by using a very simple example. Let us consider the Caputo derivative , defined in (1) below, and the fractional differential equation , with and initial conditions and . The solution is . Then, let us consider the equation , . The goal is the determination of a value of such that the solution of this equation with initial conditions and , i.e., the same initial conditions imposed to the fractional equation, be approximated by the solution of the fractional equation over a finite time interval. This is clear, since the solution of the equation on is

 z(t)=1p(−1+e−pt).

Therefore, on a time interval , with , we have . In this sense of having similar approximate solutions on a finite interval, we say that the fractional and the dissipative equations are equivalent. Here, we want to extend this idea.

Observe that in our notion of equivalence, we have discarded the asymptotic regime. This is essentially due to two reasons: i.) for large values of time, the fractional oscillator does not show oscillations; ii.) the behaviour of the oscillator from a strictly physical point of view, whether linear or non-linear but particularly the latter, has interest for finite times only. Its asymptotic behaviour is not measurable and has a mathematical interest only, and it is not the object of our study.

The present paper is organized as follows: On Section 2, we construct a method to obtain approximate solutions of fractional differential equations, with fractional Caputo derivative to be defined there, based on segmentary interpolation. This kind of interpolation has been used successfully to obtain approximate solutions of ordinary differential equations



. On Section 3, we apply this method to the fractional linear oscillator and to some other simple examples and make estimations on its precision. We compare results with those obtained replacing the fractional oscillator by the ordinary oscillator with a dissipative term. We present a similar analysis by replacing the equation of the oscillator by the van der Pol equation on Section 4. We close this paper with some concluding remarks.

## 2 Caputo fractional derivative and its evaluation by segmentary interpolation

Let be a real positive number and denote by the smaller integer bigger than . Let us define the Caputo fractional derivative, , of a times differentiable function of real variable, , as 

 Dαax(t)=1Γ(n−α)∫tax(n)(s)(t−s)α−n+1ds, (1)

where means the -th derivative of the function . Our objective, as mentioned at the header of the present section, is an evaluation of (1) using segmentary interpolation. Here, we consider that , so that the only choice for is , and this will be the case for some of our applications. Segmentary interpolation is a standard tool of wide use in the approximation of solutions of differential equations . Let us sketch the method here for completeness, using an approach that has been used in previous articles by our group [13, 14].

Let be a compact interval in the real axis . At regular intervals, we select nodes, , with , for all , so that . Let be a continuous function and use the notation and  , for all .

Then, a quadratic segmentary interpolator for the function , is a continuous function , with first continuous derivative, such that

1.- On each interval of the form , , we have that , where is a polynomial of order two, depending on the given interval.

2.- The function interpolates , in the sense that for any of the nodes , one has that

 Pk(tk−1)=xk−1,Pk(tk)=xk,k=1,2,…,n. (2)

The condition on the continuity of the derivative implies that

 P′k(tk)=P′k+1(tk),k=1,2,…,n−1. (3)

Thus, the construction of the segmentary interpolator relies in the construction of the interpolating polynomials . We propose the following form for the interpolating polynomials: For each of the intervals , let us define,

 Pk(t)=pk(t)+ak(t−tk−1)(t−tk), (4)

with

 pk(t)=t−tk−1hxk−t−tkhxk−1, (5)

and the complex coefficients are given by

 ak=n∑j=0ck,jxj. (6)

We still need to determine the values of the , which are

 cj,k=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩(−1)kh2η1,ifj=0,(−1)k+1h2(2η1+η2),ifj=1,(−1)k+jh2(ηj−1+2ηj+ηj+1),if1

where if and if .

Taking into account that is an approximation of , on each of the nodes the Caputo fractional derivative (1) is approximated by

 Dαax(tk)≈1Γ(1−α)k∑j=1∫tjtj−1P′j(s)(tk−s)αds,k=1,…,n. (8)

Then, on each of the intervals , we have that , with

 αk=2ak,βk=xk−xk−1h−ak(2tk−1+h), (9)

and, consequently, equation (8), takes the following form:

 Dαax(tk)≈1Γ(1−α)k∑j=1˜ck,jαj+˜dk,jβj. (10)

The new coefficients and are given by

 ˜ck,j=∫tjtj−1sds(tk−s)α,˜dk,j=∫tjtj−1ds(tk−s)α, (11)

which obviously depend solely of the partition.

It is customary to choose , which obviously does not restrict generality. Since the integrals in (11) are easily solvable and we know expressions for and , we can write the right hand side in (10) as

 Dαax(tk)≈1(−1+α)(−2+α)Γ(1−α)k∑j=1γk,jαj, (12)

where,

 γk,j=[h(−j+k)]−α[h(1−j+k)]−α{h(−1+j−k)[h(−j+k)]α ×[−2h(−2+j+k)+h(−3+2j)α−2(−2+α)tj−1] −h(j−k)[h(1−j+k)]α[−2h(−1+j+k)+h(−1+2j)α−2(−2+α)tj−1]}. (13)

This is a quite simple and workable receipt to obtain, once is given, the values of its Caputo fractional derivative at the nodes , so that we have an estimation of this fractional derivative.

It is interesting to remark that, due to the linear dependence on of the coefficients given in (6), then, the derivative in (12) can be explicitly determined from .

### 2.1 A type of differential equations with fractional derivative

Let be a differentiable real function of the real variable and . In addition, we assume that , where is one of the nodes , , and . Then, let us consider the following fractional differential equation:

 Dαax(t)=f(t,x(t)). (14)

The objective is to obtain an approximation for the solution of equation (14) under the condition . We already know how to obtain the identity (14) in the nodes . Take these nodes with the exception of . Then, (14) provides of an algebraic system of equations where the indeterminates are with and . This algebraic system may or may not be linear depending on the form of and is of order . A numerical solution of this system could be obtained by whatever method, which gives a segmentary solution , which is given once one has obtained the coefficients defined in (6).

In the particular case in which

contains an eigenvalue

and be linear with respect to this algebraic system is linear and homogeneous. The eigenvalue is determined in the usual way.

As the reader may easily understand, this method is more general than the usual way to obtain a solution knowing an initial value, since now could be any node. In particular, the restriction to the solution that replaces the initial value condition could be imposed at , and this represents a great advantage when compared with the shooting method worked out in [11, 15].

## 3 The fractional linear oscillator

A simple example of an equation of the type (14) is the linear oscillator with the fractional derivative, which is defined as

 Dα0x(t)=−ω2x(t). (15)

As in the standard harmonic oscillator the constant , where is the oscillator mass and a constant, being the order of derivation that in the present case we assume to be . Using the definition (1), taking into account that for some differentiable function (in our case or , where the dot means first derivative), we have that

 limα→0+1Γ(α+1)tαf(0)=f(0), (16)

and that is either 2 or 3, we may integrate by parts (15) using (1), which gives the following integral version of (2):

 x(t)=x(0)+˙x(0)−ω2Γ(−α)∫t0(t−s)−α−1x(s)ds. (17)

The general solution has the form

 x(t)=c1Eα,1(−ω2tα)+c2tEα,2(−ω2tα), (18)

where is the so called Mittag-Leffler function

 Eα,β(z)=∞∑k=0zkΓ(αk+β). (19)

Thus, in order to obtain a particular solution, we have to impose some initial conditions. For instance, if we choose and , the solution to (17) with these initial conditions is given by

 x(t)=Eα,1(−ω2tα). (20)

Let us find a particular numerical solution to the fractional linear oscillator, using the method introduced in Section 2.1. We have to choose a particular value for and the simplest possibility is . This is developed in the forthcoming subsection.

### 3.1 Some numerical estimations.

First of all, it is not the objective here to give explicit expressions for the approximate solutions for the studied examples. It is not difficult to plot these solutions for different values of .

Let us start with equation (15) with on the interval , with and initial condition . As seen above, this equation has exact solution given by . The objective is now an estimation on the precision of the proposed method. As customary, this precision is measured by using the following parameter:

 en(α)=∫10[xexact(t)−xn(t)]2dt, (21)

Here, is the number of sub-intervals in which we partite , the number of nodes being . The dependence of this parameter on shows that the smaller is the value of , or equivalently the closer is to zero, the lower is the precision and, therefore, the slower is the convergence to the exact value. However, we do not observe significative variations on the precision when we increase the value of , i.e., as we make the sub-intervals smaller and smaller.

Table 1.- Values of the precision in terms of and .

This can be seen in Table 1, where we have chosen values of ranging from 5 to 40. The values of studied are 0.1, 0.5 and 0.9.

Let us study the precision of the method with another example different from the fractional oscillator, yet an equation of the form (14). Here, we have chosen,

 D1/20x(t)=sinx(t), (22)

on the interval , with the initial condition , which was already studied in , where the integration was performed by means of the iterative shot method. Contrarily to the previous example, here we do not know an exact solution. The way out is to define the precision as

 en=∫10[D1/20xn(t)−sinxn(t)]2dx, (23)

where is again the number of intervals and is the interpolating function for the studied case. After integration and using the boundary condition, we obtain . Along with (23), we introduce another parameter that measures the convergence and that we denote as . It represents the relative variation between the value of , obtained for a given value of , and the value given for the precedent value of as listed on Table 2.

Table 2 is just a sample of numerous numerical examples we have performed. This sample is significative, as it manifest an obvious convergence and shows that the result obtained for a small number of nodes is satisfactory.

Table 2: Values of , and for a given value of

Finally, we have performed another reliability test, which was the use of the value of obtained numerically as initial value and evaluate the value of . In all cases, we have recovered the value .

### 3.2 Damped oscillator with entire derivative.

As is well known, the damped oscillator with entire derivative is given by

 m¨y(t)+p˙y(t)+ky(y)=0, (24)

where , and are constants. Here, we assume that .

For , the limit for in (15) should give the solution for (24), which we denote here as . The solution is damped oscillatory on the transitory regime only [9, 10, 1]. Based on these notions, we propose the following Ansatz:

For each given , there exists such that the solution of (15) with is a good approximation of the solution of (24) with , in the transitory regime.

Using this Ansatz, let us obtain an approximate solution for (24) such that this and its corresponding solution for (15) fulfil the conditions and . This is:

 y(t)=exp(−pmt)(c−exp(−Δt)+c+exp(Δt)), (25)

where,

 Δ=√(p/m)2−4ω2,λ±=−p±Δ2,c±=±λ∓Δ, (26)

where was given in (15).

Then, the point is the determination of the value of being given the value of , or equivalently the determination of a function , in application to our Ansatz. This is an optimal control problem. We have to find the optimal solution, which minimizes the following functional:

 E(α):=1T∫T0[x(t)−y(t)]2dt, (27)

where is some time scale in which the amplitude of the oscillations are reduced by a factor of . On the interval , the transitory regime, we compare the solutions of the fractional derivative and of the damped oscillator . The functional measures the deviation between and . Then, go back to (20) and note that the function is not asymptotically oscillating as . This permits us to choose a value of , although not small, not very high either. Numerical experiments have shown that the choice is appropriate.

Let us give some numerical results. On Table 3, we give the dependence between values of , and for the values and .

Table 3: Comparison between the values of , and ,

for , and .

An explicit expression of the function may be obtained by the least-square method and this gives . This is depicted on Figure 2.

On Figure 1, we represent the usual behaviour in the transitory regime for and , when we choose and .

It is interesting to remark that numerical experiments show that does not depend on any choice of initial values.

### 3.3 Non-linear oscillator

Following the discussion on the damped oscillator, we present a similar problem given by the following non-linear oscillator:

 Dα0x(t)=y(t),Dα0y(t)=−sinx(t), (28)

with . Let us choose the initial values given by and . Clearly, for small oscillations equation (28) becomes equation (15) with the replacement . Again, this is a non-linear problem having no analytic solution for non-integer. Then, we proceed by analogy with the damped oscillator. In this case, we consider the following system involving entire derivatives only:

 ˙z(t)=w(t),˙w(t)=−pw(t)−sinz(t),p>0. (29)

On the transitory regime, we compare the solutions of systems (28) and (29) under the conditions and . To do this, we need the previous determination of , which we assume that minimizes the following quadratic dispersion:

 E(α)=1T∫T0{[xn(t)−z(t)]2+[yn(t)−w(t)]2}dt. (30)

Obviously, this expression generalizes (27). Again, we adjust the value of by numerical experiments, which show that is, again, a convenient choice. On Table 4, we give some values for the dependence between , and after the choice and .

Table 4: Comparison between the values of , and ,

for the values and .

The above examples manifest an analogous behaviour between a fractional linear oscillator and a damped or even non-linear oscillator on some time interval. The solutions between the fractional and the integer equation are very similar on some time scale. This could be a rather general situation, so that in many practical cases and inside a time interval, we conjecture that a fractional operator might be replaced by the frictional additional term on the classical oscillator. The behaviour of the solutions is similar to that shown in Figure 1. Figure 1: The continuous line represents the solution, x(t) of the fractional equation (15), which is given by (25). The dashed line gives the solution, y(t), of the equation with ordinary derivative (25). Figure 2: Function p(α).

## 4 On the fractional van der Pol equation

The van der Pol equation is a second order non-linear ordinary differential equation [16, 18]. It has the following form:

 ¨x(t)−μ(1−x2(t))˙x(t)+x(t)=0, (31)

where is a constant. When , (31) is the equation of the ordinary harmonic oscillator. The van der Pol equation may be written in terms of a first order system as

 ˙x(t):=z(t),˙z(t)=−μ(x2(t)−1)z(t)−x(t),μ≥0. (32)

This equation has a unique limit cycle for , after the Liénard theorem . This suggested us the interest of considering the possible existence of a limit cycle for the fractional system analogous to (32) given by

 Dα0x(t)=z(t),Dα0z(t)=−μz(t)(x2(t)−1)−x(t),0<α≤1. (33)

One fractional van der Pol equation of the type

 Dα+10x(t)+μ(x2(t)−1)Dα0x(t)+x(t)=0, (34)

has been studied in , where a relation between the parameters and was given as a sufficient condition for the existence of a limit cycle, using the balance harmonic method .

We have studied the system (33) through numerical as well as analytic methods. We have performed a big amount of numerical experiments, which have shown the existence of a value of the parameter , here called , where the subindex stands for critical, which depends on and , such that

• For values of with , there is a fixed point , which remains stable at the limit , . Therefore, there is no stable limit cycle. In addition, there is no evidence of the existence of an unstable limit cycle.

• For values , the fixed point is unstable. We found a unique stable limit cycle. In this case a Hopf bifurcation emerges with as critical parameter.

• As shown in Figure 3, decreases with and .

The point here is to show the existence of the critical value for the parameter for a given value of , . This existence has been manifested by the numerical estimations above. Nevertheless, this existence may be also shown analytically. To this end, we use the following result :

Let us consider the following system, where represents the Caputo fractional derivative:

 Dα0x(t)=f(x,z),Dα0z(t)=g(x,z),0<α<1. (35)

A solution is in equilibrium if . It is asymptotically stable if the eigenvalues, , of the Jacobian matrix

 J(x,z):=(∂f/∂x∂f/∂z∂g/∂x∂g/∂z), (36)

when evaluated at the equilibrium point satisfies

 |arg(λ)|>απ2. (37)

A comparison between (33) and (35) gives the precise form of and for our particular case. This gives the precise form of (36) as

 J(x,z):=(01−1−2μz(t)x(t)−μ(x2(t)−1)). (38)

Taking into account that the fixed point is located at , we have that

 J(0,0)=(01−1μ), (39)

which has the following eigenvalues:

 λ±=12(μ±√μ2−4). (40)

Obviously, if , then . On the other hand, if , one has that,

 arg(λ±)=arctan⎛⎝±√(2μ)2−1⎞⎠. (41)

From (37), the critical value, , of should obey the following relation:

 arctan⎛⎝±√(2μc)2−1⎞⎠=απ2, (42)

which gives

 μc=2√1+tan2(απ2). (43)

This is to say, if we fix and start from , as we increase , we go from a situation with a asymptotically stable fixed point to an unstable point. This happens when . The transition from the stability to the unstability drives to the emergency of a limit cycle. We may qualitatively interpret the limit cycle loss as follows: let us consider in (33), which may then be approximated by

 Dα0x(t)=z(t),Dα0z(t)=−x(t). (44)

This is the same than equation (28) with the paraxial approximation . Note that (44) does not show a limit cycle and, further, the trivial solution is an attractor. The second equation in (33) contains the term , which in the case of is not negligible. This fact outbalances the dissipation and this is precisely which makes it possible the existence of a limit cycle.

On Figure 4 and on the phase plane, the continuous and slashed curves represent the solution with entire and fractional derivative, respectively. Both trajectories are determined with same initial values and same parameter . In all cases, the fractional limit circle is enclosed by the trajectory of the limit cycle with entire derivative.

On Figure 3, we show the relation . In the region above the curve, there exists a stable cycle limit and, furthermore, the fixed point is unstable. Below the curve the limit cycle does not exist and the fixed point is asymptotically stable. There is an obvious difference with the results obtained in , which is due to the fact that the fractional equations (33) and (34) are not equivalent.

### 4.1 Equivalence between the fractional van der Pol equation and the same equation with entire derivative and dissipation.

On the previous section, we have compared the approximate solutions of a dissipative oscillator with entire derivative with those of the linear oscillator with fractional derivative. Now, we want to carry out a similar analysis with the fractional van der Pol equation and a van der Pol equation with entire derivative and a dissipative term of the form , . This system has the form,

 x′(t)=z(t),z′(t)=−z(t)(β+μ(x2(t)−1))−z(t),μ,β>0. (45)

Here, the fixed point is . To check its stability, we consider again equation (37), which in the present case gives at the fixed point the following expression

 J(0,0)=(01−1μ−β), (46)

which has the eigenvalues

 λ±=(μ−β)±√(μ−β)2−42. (47)

Therefore, the fixed point is stable if Re and unstable if Re, or equivalently, if and , respectively. Then, for each , there exists a , where the Hopf bifurcation of the fixed point appears and, consequently, the destruction of the limit cycle.

In any case, according to the Liénard theorem [17, 18], we may show that there exists a unique stable limit cycle if . In consequence, the van der Pol equations with entire derivative and dissipation and fractional have qualitatively the same properties.

We have checked numerically a qualitative equivalence, in the sense of having approximately the same solution, through a substantial number of numerical experiments, between equations (44) (with fractional derivative) and (45) (with entire derivative). Thus by trial and error, we have determined a value of giving the same cycle in both cases. For instance, if we give the values and , we obtained . In an analogous manner, we trial with values for which and obtained similar conclusions.

On Figure 5, we represent limit cycles for the fractional and dissipative entire van der Pol equations. Figure 4: Comparison between the entire (continuous curve) and the fractional (slashed curve) van der Pol solutions, with the same value of μ. Figure 5: Comparison between the damped (continuous curve) and the fractional (slashed curve) van der Pol solutions.

## 5 Concluding remarks

We have applied a quadratic spline method in order to obtain functions that approximate the result of applying a fractional derivative to a given function. This is suitable to obtain solutions to some differential equations with initial values or mixed conditions of potential interest in Physics or Engineering. We have tested our method with the fractional linear oscillator, where exact solutions are known and checked its degree of precision. Results of numerous numerical experiments show that for values of in the range , the higher is , the better is the precision of our method. Here, is the order of the fractional derivation, . However, there are not substantial differences when we increase the number of nodes on the interval under our consideration. Similar results have been obtained for non-linear oscillators.

Based on previous studies on the approximation of solutions of the fractional linear oscillator by solutions of a damped oscillator, we have used our method to confirm these results. We have shown that there exists a time interval for which the solutions of both equations are similar with a high degree of accuracy, under the condition that a relation is given between the coefficient of the dissipative term of the damped oscillator and the order of the fractional derivation, .

A similar study compares a fractional and a damped van der Pol equations, written as a system of two equations on phase space, with similar results. In addition, we have considered the behaviour of limit cycles and fixed points in terms of and a parameter characteristic of the van der Pol equation. Using analytic as well as numerical arguments, we show the existence of a critical value for the parameter , , such that if the origin of phase space is stable and if is unestable. This limit value depends on and we give the exact relation.

## Acknowledgements

This research has been financed by the Projects No. ING 19/i 402 and ING 80020180100064 of the Universidad Nacional de Rosario, the Spanish MINECO (Project No. MTM2014-57129) and the Junta de Castilla y León (Project Nos. BU229P18 and VA137G18).

## References

•  A. Kilbas, H. Srivastava, J. Trujillo. Theory and Applications of Fractional differential equations, Elsevier, Amsterdam (2006).
•  R. Herrmann, Fractional Calculus: An Introduction for Physicists. World Scientific, (2011).
•  V.V. Uchaikin, Fractional derivatives for physicists and engineers. Volume I Background and theory volume, Volume II Applications. Berlin Heidelberg: Springer Science & Business Media (2013).
•  F. Mainardi, Fractional relaxation-oscillation and fractional diffusion-wave phenomena, Chaos, Solitons and Fractals, 7 (9), 1461-1477 (1996).
•  F. Olivar-Romero, O. Rosas-Ortiz, Fractional driven damped oscillator, J. Phys. Conf. Ser., 839, 012010 (2017).
•  F. Olivar-Romero, O. Rosas-Ortiz, Transition from the Wave Equation to Either the Heat or the Transport Equations through Fractional Differential Expressions, Symmetry, 10, 524 (2018).
•  F. Olivar-Romero, O. Rosas-Ortiz, Factorization of the quantum fractional oscillator, J. Phys. Conf. Ser., 698, 012025 (2016).
•  Can Evren Yarman, Approximating fractional derivative of Faddeeva function, Gaussian function, and Dawson’s integral, Math. Met. Appl. Scie., DOI: 10.1002/mma.5679 (2019).
•  B.N. Narahari Archar, J.W. Hanneken, T. Enck, T. Clarke, Dynamics of the fractional oscillator, Physica A, 297, 361-367 (2001).
•  K. Diethelm, The Analysis of Fractional Differential Equations. An Application Oriented Exposition Using Differential Operators of Caputo Type, Springer Verlag, Berlin (2010).
•  K. Diethelm, W. McLean, Volterra integral equations and fractional calculus: do neighboring solutions intersect?, Journal of Integral Equations and Applications, 24 (1), 25-37 (2012).
•  E. Hairer, S.P. Nørsett and G. Wanner, Solving Ordinary Differential Equations I, Springer, Berlin and New York, 1993.
•  M. Gadella, L.P. Lara, A study of periodic potentials based in quadratic splines, Int. J. Mod. Phys. C, 29, 1850067 (2018).
•  A. Ferrari, L.P. Lara, E Santillan-Marcus, Convergence analysis and parity conservation of a new form of quadratic splines, arXiv:1906.10559v1 (2019).
•  H. Demir, Y. Balturk, On numerical solutions of fractional order boundary value problem with shooting method, ITM Web of Conferences, 13, 01032 (2017).
•  B. van der Pol, On relaxation equations, The London, Edinburgh and Dublin Phil. Mag. & J. of Sci., 2 (7), 978-992 (1927).
•  S. Strogatz, Nonlinear Dynamics and Chaos, CRC Press, Taylor and Francis, Boca Raton, London and New York, 2015.
•  M. Farkas, Periodic Motions, Springer, New York and Berlin (1994).
•  Z. Guo, A.Y.T. Leung, H.X. Yang, Osillatory region and asymptotic solution of fractional van der Pol oscillator via residue harmonic balance technique, Apply Mathematical Modelling, 35, 3918-3925 (2011).
•  P. Deuflhard, Newton Methods for Nonlinear Problems, Springer, Berlin, 2006.
•  E. Ahmed, A. El-Sayed, A.A. El-Saka, On some Routh-Hurwitz conditions for fractional order differential equations and their applications in Lorenz, Rössler, Chua and Chen systems, Phys. Lett. A, 358, 1-4 (2006).