A Note on the Solutions of the Van der Pol and Duffing Equations Using a Linearisation Method

We present a novel application of the successive linearisation method to the classical Van der Pol and Duffing oscillator equations. By recasting the governing equations as nonlinear eigenvalue problems we obtain accurate values of the frequency and amplitude. We demonstrate that the proposed method can be used to obtain the limit cycle and bifurcation diagrams of the governing equations. Comparison with exact and other results in the literature shows that the method is accurate and effective in finding solutions of nonlinear equations with oscillatory solutions, nonlinear eigenvalue problems, and other nonlinear problems with bifurcations.


Introduction
The Van der Pol equation is ẍ x ε x 2 − 1 ẋ 0, x 0 α, ẋ 0 0, 1.1 where x defines the displacement of the periodic solution, α is the amplitude of the oscillations, and ε > 0 is a classical problem in nonlinear dynamics that models systems with selfsustained oscillations.The equation is often studied for its rich set of oscillatory and chaotic solutions Atay 1 and as prototype for testing numerical schemes since one of its attributes is that for large ε, the equation is stiff with solutions exhibiting oscillations on a large time scale.
Finding solutions of the Van der Pol equation using the perturbation series expansion method generally fails since the solutions have multiple time scales.Nonetheless, solutions in the form of a Taylor series and numerical solutions are well documented in the literature, for example, Buonomo 2 presented a procedure for finding the periodic solutions as a power series in powers of ε.Recent attempts to solve the Van der Pol equation include using, among other methods, the decomposition method, Asadi Cordshooli and Vahidi 3 , the parameter expansion methods Darvishi and Kheybari 4 Kimiaeifar et al. 5 , and the homotopy analysis method, Li et al. 6 .
We also consider the so-called Duffing oscillator in space, which is given in dimensionless form as where the prime denotes differentiation with respect to the variable t and ε ≥ 0 is a given parameter.The bifurcation of this equation for some values of ε was examined using norm forms in Kahn and Zarmi 7 and the homotopy analysis method by Liao 8 .
In this paper, we present a novel application of the successive linearisation method to initial value problems and boundary value problems with oscillatory behaviour.We demonstrate the application of this technique by finding solutions of the classical Van der Pol and the Duffing oscillator equations.We calculate the amplitude and frequency of the limit cycle of the Van der Pol equation and compare our results with known exact solutions and other results from literature.In the case of the Duffing oscillator we obtain the bifurcation diagrams and compare with exact solutions.

SLM Solution of Van der Pol Equation
In this section we present the basic description of the successive linearisation method SLM that is used to obtain solutions of the nonlinear initial value problem IVP , 1.1 and the boundary value problem BVP , 1.2 .However, before the governing equation 1.1 is solved it is converted into a boundary value problem by using the definition of period functions.In its basic form, the SLM 10, 11 seeks to linearize and reduce the governing nonlinear equation to a system of linear differential equations.
Following 12, 13 we introduce the transformation z ωt, where ω is the frequency of the oscillations.In addition we set x t αu t so that the boundary conditions are independent of α.The frequency ω and amplitude α are considered to be unknown parameters which will be determined by the SLM.Under these transformations, 1.1 becomes where the primes denote differentiation with respect to z.Since ω and α are considered to be unknown parameters, additional boundary conditions are required to fully solve 2.1 .The extra boundary conditions are derived by noting that the Duffing-Van der Pol equation 1.1 has translational invariance, that is, if u z is a solution, then so is u z φ for an arbitrary constant φ.The limit cycle is automatically periodic with a period 2π.We therefore introduce the following conditions: We now seek a solution of the form: where

2.6
The linearized system 2.4 was solved using the Chebyshev collocation spectral method.We used the transformation z τ 1 /2 to map 0, 2π to −1, 1 .The solution space is discretized using the Chebyshev-Gauss-Lobatto collocation points Mathematical Problems in Engineering which are the extrema of the Nth order Chebyshev polynomial We approximated the derivatives of the unknown variables u s at the collocation points as the matrix vector product: where D D/π, with D being the Chebyshev derivative matrix and U is the vector function of u s at the collocation points τ j see, for example, 14, 15 .The entries of D can be computed in different ways.In this work we use the method proposed by Trefethen 15 in the cheb.mMATLAB m-file.If we denote the entries of the derivative matrix D by D jk , we can apply the spectral collocation method, with derivative matrices on the linearized equation 2.4 and the boundary conditions 2.5 to obtain the linear matrix system: where and a r,s−1 , r 1, 2, 3 are the diagonal matrices corresponding to a r,s−1 τ j .The vectors a r,s−1 , r 4, 5 are of size N 1 × 1 and they correspond to a r,s−1 τ j .The appropriate boundary conditions are After imposing the conditions 2.12 on 2.10 we obtain the following matrix system: The boundary conditions u s τ N 0 and u s τ 0 0 appearing in 2.12 have been imposed by deleting the first and last rows of A s−1 and R s−1 .Two additional rows are also added to accommodate the two derivative boundary conditions.Thus, starting from the initial approximations u 0 z , ω 0 , and α 0 , the subsequent solutions for u s z , ω s , α s s 1, 2, 3 . . .can be obtained by solving the matrix system 2.13 .

SLM Solution of the Duffing Equation
The exact solution of 1.2 is given see 7, 8 in terms of the relation between α and L as 2.14 which gives the exact solution where K ζ is the complete elliptic integral of the first kind.
To obtain an approximate solution of 1.2 , it is convenient to introduce the following transformations 8 : Substituting in 1.2 , we have Since α is unknown in 2.17 , we require an additional boundary condition to be able to solve 2.17 .The extra boundary condition can be obtained from 2.16 as u π 2 1.

2.18
Equations 2.17 and 2.18 form an eigenvalue problem with α being the unknown eigenvalue.We now look for solutions u t and α in the form where u s and α s are unknown functions.Substituting 2.19 in 2.17 and neglecting nonlinear terms give subject to the boundary conditions where

2.22
The linearized system 2.20 was solved using the Chebyshev collocation spectral method described in the previous section.Using the spectral method and the transformation z π τ 1 /2 to map 0, π to −1, 1 gives with the boundary conditions where and b r,s−1 , r 1, 2 are the diagonal matrices corresponding to b r,s−1 τ j .After imposing the conditions 2.24 on 2.23 we obtain the following matrix system: . . .

2.26
Starting from suitable initial approximations, u 0 and α 0 , 2.26 can be used to find approximate u and α.

Results and Discussion
In this section we present solutions of the Van der Pol equation 1.1 and Duffing oscillatory equation 1.2 obtained using the SLM solution procedure outlined in the previous sections.The initial approximations were obtained using physical considerations and the given boundary conditions.For 1.1 , the initial approximation used is u 0 cos z while for the Duffing equation 1.2 u 0 sin t was used.In both cases the initial values of all eigenvalues were set to be 1.The number of collocations points used was N 200 and N 100 for the Van der Pol and Duffing equations, respectively.To determine the accuracy and performance of our method, the results were compared against previous results in literature and exact solutions for selected parameter values.In particular, for the Van der Pol equation, the present results were compared to the exact solution and the homotopy analysis method results reported in 9 and against the power series solution results of 2 when ε is small.For the Duffing equation we compared with the exact solution 2.15 .
Table 1 shows a comparison of the amplitude and frequency of the Van der Pol equation obtained using the successive linearisation method compared to the exact results and recent solutions obtained using the homotopy analysis method HAM .The parameter ε 1 has been kept fixed.We note firstly that only six iterations of the linearisation method are sufficient to match the exact result to ten decimal places compared to at least fifty iterations in the case of the homotopy analysis method, Chen and Liu 9 .This shows that the rate of convergence of the present SLM approach is higher than that of the HAM.
Table 2 shows the effect of increasing the parameter ε on the accuracy and efficiency of the successive linearisation method.The solutions are compared with the power series solutions of Buonomo 2 and the homotopy analysis results of Chen and Liu 9 .For Mathematical Problems in Engineering  values of ε ≤ 0.6 only two iterations of the SLM are sufficient to give accuracy of up to 6 decimal places.However, as ε increases, more iterations are required to achieve the same level of accuracy.We also note that the SLM requires fewer iterations to converge to the exact solutions compared to the HAM. Figure 1 shows the limit cycle generated using the SLM when ε 0.5 compared with numerical results generated using the MATLAB ode45 solver.It can be seen from the graph that there is good agreement between the two results.
In Table 3 we give the SLM results of the unknown parameter α at different orders when ε 10 and ε 15.The results are compared with the HAM and homotopy-Padé results of 8 .The SLM converges much faster than both the HAM and the homotopy-Padé technique.When ε 10, convergence to 5 decimal places of accuracy is achieved after only 4 iterations in the SLM compared to 20 iterations and 6, 6 in the HAM and homotopy-Padé techniques, respectively.We also note that convergence to 8 digits of accuracy is achieved after only 5 iterations when ε 10.When ε 25, convergence to 5 decimal places is achieved after only 3 iterations in the SLM compared to 10, 10 in the homotopy-Padé.The HAM solution will not have converged to a 5-digit accurate solution after 30 iterations.This is an advantage of SLM method over other nonperturbation techniques when solving BVPs of type 1.2 .Figure 2 gives the bifurcation diagram in the α, ε plane and a further comparison between the exact solution 2.15 and the fourth-order SLM result.It can be see from the graph that there is good agreement between the two results.The two examples considered in this paper show that the SLM is sufficiently robust to handle solutions of nonlinear problems with oscillatory solutions.In Motsa and Sibanda 16 , it has been shown that, with minor modifications, the SLM can be easily extended to chaotic and hyperchaotic systems of initial value problems.

Conclusion
We have presented a new application of the successive linearisation method to find solutions of oscillatory systems.The application of the method has been demonstrated by solving the Van der Pol and Duffing oscillatory equations.Using suitable transformations, the governing equations were transformed to eigenvalue boundary value problems and solved iteratively using the SLM.Comparison was made between exact analytical solutions of the governing equations and against results in the literature.Our results show that the SLM rapidly converges to the exact solutions and can be used as an efficient method for finding solutions of boundary value problems, requiring only a few iterations to give accurate solutions.The method can also be easily extended to other nonlinear oscillating systems and nonlinear problems with bifurcations.

Figure 2 :
Figure 2: Comparison of the exact solution dots with the 4th order SLM solution line .

Table 1 :
Comparison of the results of the frequency and amplitude by the proposed SLM approach with the exact results and the homotopy analysis method results of Chen and Liu 9 when ε 1.

Table 2 :
Comparison of the frequency by the SLM approach with Buonomo 2 and the HAM, Chen and Liu 9 .

Table 3 :
Comparison of the SLM results and the HAM and Homotopy-Padé 8 approximate results for α at different orders p for α for ε 10 and ε 25.