Numerical Investigation of the Fractional-Order Liénard and Duffing Equations Arising in Oscillating Circuit Theory

In this article, we present the Jacobi spectral colocation method to solve the fractional model of Liénard and Duffing equations with the Liouville–Caputo fractional derivative. These equations are the generalization of the spring–mass system equation and describe the oscillating circuit. The main reason for using this technique is high accuracy and low computational cost compared to some other methods. The main solution behaviors of these equations are due to fractional orders, which are explained graphically. The convergence analysis of the proposed method is also provided. A comparison is made between the exact and approximate solutions.


INTRODUCTION
The standard Liénard equation (LE) is a generalization of the damped pendulum equation or spring-mass system. Because this equation can be applied to describe the oscillating circuits, therefore, it is used in the development of radio and vacuum-tube technology. The LE was given by Liénard [1], and it is written as follows: where τ 1 (v) D ′ v is the damping force, τ 2 (v) is the restoring force, and τ 3 (t) is the external force. For different choices of the variable coefficients τ 1 (v) , τ 2 (v), and τ 3 (t), the LE is used in many phenomena. The Liénard Equation (1) becomes the van der Pol equation for τ 1 (v) = ε v 2 − 1 , τ 2 (v) = v, and τ 3 (t) = 0, which has many applications [2,3]. By usual way, we cannot find the exact solution for these equations [4]. Kong [5] studied the LE given as follows: where a, b, and c are real constants.
In particular, if we take c = 0 in the LE, then it reduces to the Duffing equation (DE). This special case of the LE is known as the DE and is given as follows; where a, d, and b are real constants. In recent years, fractional calculus has become an interesting and useful part of mathematical analysis and applied mathematics. The importance of fractional calculus arises because of its non-local nature. The real-life applications of fractional calculus are in fluid dynamics [6], signal processing [7], chemistry [8], viscoelasticity [9], and bioengineering [10]. For some other applications, see Srivastava et al. [11], Kilbas et al. [12], and Robinson [13]. Many physical problems are modeled by fractional-order LE (FLE) and DE. In addition, we are familiar with the fact that non-integer-order derivatives handle models accurately. So, for the accurate modeling of these equations, it is fundamentally needed to change integer-order equations to fractional-order equations.

Fractional-Order Liénard Equation
The FLE is given by (4) with the conditions: where ξ and η are real constants.

Fractional-Order DE
The fractional-order DE (FDE) is given by with the following conditions: where µ and σ are real constants. The innovator approach to solve the LE originates in the work by Kong [5], who provided an exact solution of these equations in some particular cases. For some particular choices of the involved real constants, Feng [14] obtained an exact solution of these equations, which were the generalization of Kong's [5] results. In 2008, Matinfar et al. [15] suggested a variational iteration method in order to obtain the approximate solutions of the LE. Subsequently, in 2011, a variational homotopy perturbation method was applied in order to solve LE (see Matinfar et al. [16]). Recently in 2017, a numerical method using homotopy analysis transform method (HATM) in order to solve fractional LE was proposed, and the uniqueness and existence of solutions were also given (see Kumar et al. [17]). Further, Singh [18,19] used Legendre polynomials and Chebyshev polynomials, respectively, to solve fractional models of these equations.
In this article, we propose an effective method for the FLE and DE. The proposed method is a spectral colocation method based on the applications of operational matrix of differentiation for the Jacobi polynomials. Spectral colocation method is used to solve many problems in differential calculus (see [20][21][22][23][24][25][26][27][28][29]). By using the spectral colocation method, these equations are converted into a system of non-linear algebraic equations whose solution gives approximate solution to these equations. The derived solution is discussed for different fractional orders. The obtained results are compared with the exact solution and presented in the form of numerical tables. Because fractional order derivatives are nonlocal in nature, and integer-order derivatives are a special case of fractional order derivative, it is important to study fractional order models. The proposed method is easy to implement because it is computer oriented. It is also a time-saving method. The integer as well as fractional order behavior of solution is shown in numerical section. The main solution behaviors of these equations are due to fractional orders, and using the proposed method, these behaviors of solution are explained clearly.

PRELIMINARIES
In this article, we have considered the non-integer-order differentiations in Liouville-Caputo (LC) sense, which are defined as follows: The LC non-integer derivative of order β is defined as follows [30,31]: In this article, we have used Jacobi polynomials as a basis for the approximation of unknown functions. The shifted Jacobi polynomial is given as follows [32][33][34]: (9) where e and f are parameters in Jacobi polynomials as given in Doha et al. [32].
The orthogonal property of Jacobi polynomials is as follows: where g (e,f ) (t) is weight function, and δ mn is kronecker delta function and given as A function f ∈ L 2 g [0, 1], with f ′′ (t) ≤ A, can be expanded as follows: where , and ., . denotes the usual inner product space.
Equation (12), for finite dimensional approximation, is written as follows: where C and p m (t) are (m + 1) × 1 matrices given by T be the shifted Jacobi vector and if v > 0, then where D (v) = q i, j is an (n + 1) × (n + 1) operational matrix of non-integer derivative of order v, and its entries are given by Proof. See Doha et al. [32], Ahmadian et al. [33], and Bhrawy et al. [34].

OUTLINE OF THE METHOD
Here, we will describe the algorithm for the construction of the solution for the fractional LE and DE using operational matrix and colocation method [27][28][29]. Let us take the following approximation: Then, by taking the derivative of order one on both sides of Equation (16), we get where D (1) is the operational matrix of differentiations for the Jacobi polynomials of order 1. Next, by taking the derivatives of orders α and β on both sides of Equation (16), we find that where D (α) and D (β) are the operational matrices of differentiations for the Jacobi polynomials of orders α and β, respectively.

Fractional-Order LE
Grouping Equations (4) and (16)- (18), we get The residual for Equation (22) is given as follows: 5 . (23) Now, colocating Equation (23) at n−1 points given by t i = i n , i = 1, 2, . . . , n − 1, we find that Further, from Equations (5), (20), and (21), we can write where ξ and η are real constants. Using the colocation points in Equation (24), together with Equation (25), we get a system of non-linear algebraic equations with the same number of unknowns. The solution of this system leads the solution for FLE.

Fractional-Order DE
Grouping Equations (6), (16), (17), and (19), we get The residual for Equation (26) is given as follows: 3 . (27) Now, colocating Equation (27) at the n − 1 points given by t i = i n , i = 1, 2, . . . , n − 1, we get R n (t i ) = C T D (β) p n (t i ) + aC T D (1) p n (t i ) + dC T p n (t i ) Further, from Equations (7), (20), and (21), we can write where µ and σ are real constants. By using the colocation points in Equation (28), together with the Equation (29), we get a system of equations with the same number of unknowns. The solution of this system leads the approximate solution for the FDE. and v n (t) be the n th approximation obtained by using Jacobi polynomials, then

CONVERGENCE ANALYSIS
and the error vector in Equation (30) tends to zero as n → ∞.
Theorem 4.2. If E α, g D, n be the error vector for α order operational matrix integration, which is obtained using (n + 1) Jacobi polynomials. Then and the error vector in Equation (31) tends to zero as n → ∞.
Proof: See Kazem [38]. Let V n be the n-dimensional subspace generated by λ (e,f ) i 0≤i≤n for L 2 g [0, 1]. Let δ n is the minimum value of the functional on the space V n . We can write V n ⊂ V n+1 and δ n+1 ≥ δ n . Proof: See Ezz-Eldien [39].
Functional for FLE is given as follows: Using Equations (16)-(18), we get Residual for Equation (33), is given as Now, similar as in Equation (23), colocating Equation (37), at n − 1 points given by t i = i n , i = 1, 2, . . . , n − 1, we get Using the colocation points in Equation (37), together with Equation (25), we get a system of non-linear algebraic equations. The solution of this system leads the solution for FLE. Let this solution be denoted by δ * n (t). Using Theorems 4.1 and 4.2 and taking n → ∞, From Theorem 4.3 and Equation (39), we achieve that lim n→∞ δ * n (t) = δ(t).
Proof completed. Similar proof can be written for convergence of DE.

NUMERICAL SIMULATION OF RESULTS
In this section, we implement our proposed algorithm by testing it on some special cases of the LE and DE. We study the applicability and accuracy of our proposed computational method by applying it on the FLE and DE. The parameters in the LE and DE are chosen in such a way for which the exact solution is known. Case 1. For the particular choices of the parameters a = −1, b = 4 and c = 3 in Equation (4), the FLE is given as follows (see Singh [18] and Tohidi et al. [20]): where τ = 4 3a 2 3b 2 − 16ac and δ = −1 + √ 3b The exact solution for the FLE given by Equation (40), with conditions in Equation (41), is given by where τ and δ are as given in Equation (42). In Figures 1, 2, we have shown the approximate solution for different values of α for the FLE choosing different parameters in the Jacobi polynomials. In Figure 3, we have compared     approximate solution by our proposed method and solution obtained by the methods of Singh [18,19] for integer-order LE. Figures 1, 2 show that the period will be really affected by the non-integer-order values, and the solution varies continually from non-integer-order solution to integer-order solution and coincides with the integer-order solution at α = 2. The solution has some different behavior when the value of fractional order is 1.76, and this is because the main solution behavior of LE takes place when α is very close to 2. Figure 3 shows that solution has exact the behavior as the methods of Singh [18,19]. In Table 1, we have listed approximate and exact solutions for the integerorder equation. Table 1 shows a good accuracy of the achieved solution. In Table 2, we have listed approximate solution by our method and the methods of Singh [18,19]. Table 2 shows good agreement with these methods.
Case 2. For the particular choices of the parameters a = 0.5, b = 25, and c = 25 in Equation (6), the fractional DE is given as follows [see Singh [18,19] and Nourazar and Mirzabeigy [40]]:  The analytical solution using the differential transform method (DTM) for fractional DE given by Equation (44), with the initial conditions in Equation (45) Figure 6, we have compared approximate solution by our proposed method and solution obtained by the methods of Singh [18,19] for integer-order DE.   has some different behavior when the value of fractional order is 1.76, and this is because the main solution behavior of DE takes place when β is very close to 2. Figure 6 shows that solution has the exact behavior as the methods of Singh [18,19]. In Table 3, we have listed the approximate and exact solutions by the DTM method for the integer-order equation. Table 3 shows a good accuracy of the achieved solution.

CONCLUSIONS
In this article, we have presented numerical solution and simulation for fractional-order and integer-order LE and DE. The proposed algorithm is easy to implement because the construction of the operational matrix is sufficiently easy, which makes our method remarkably attractive for practical applications. In the numerical section, it is presented how the approximate solution varies continuously for different values of the fractional time derivatives and for the integer-order approximate solution is the same as the exact solution for the fractional LE and DE. Recently, many equations in science and engineering appear in the form of non-linear fractional differential equations, which makes it necessary to investigate the method of solution for such equations. The main advantage of the proposed method is that it works for such type of equations arising in science and engineering. In the future, we can use operational matrices of different orthogonal polynomials to achieve better accuracy.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/supplementary material.