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

Abstract

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.

1. Introduction

The Van der Pol equation is 𝑥̈+𝑥+𝜀𝑥²−1𝑥̇=0,𝑥(0)=𝛼,𝑥̇(0)=0, (1.1) where 𝑥 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 self-sustained 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 𝑥^′′+𝜀𝑥−𝑥³=0,𝑥(0)=0,𝑥(𝜋)=0, (1.2) where the prime denotes differentiation with respect to the variable 𝑡 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.

2. 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 𝑧=𝜔𝑡, where 𝜔 is the frequency of the oscillations. In addition we set 𝑥(𝑡)=𝛼𝑢(𝑡) 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 𝜔²𝑢^′′+𝑢−𝜀𝜔𝑢^′+𝜀𝜔𝛼²𝑢^′𝑢²=0,𝑢(0)=1,𝑢^′(0)=0, (2.1) where the primes denote differentiation with respect to 𝑧. 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 𝑢(𝑧) is a solution, then so is 𝑢(𝑧+𝜙) for an arbitrary constant 𝜙. The limit cycle is automatically periodic with a period 2𝜋. We therefore introduce the following conditions: 𝑢(2𝜋)=1,𝑢^′(2𝜋)=0. (2.2) We now seek a solution of the form: 𝑢(𝑧)=𝑢_𝑠(𝑧)+^𝑠−1_𝑚=0𝑢_𝑚(𝑧),𝜔=𝜔_𝑠+^𝑠−1_𝑚=0𝜔_𝑚,𝛼=𝛼_𝑠+^𝑠−1_𝑚=0𝛼_𝑚,𝑠=1,2,3,…, (2.3) where 𝑢_𝑠, 𝜔_𝑠, 𝛼_𝑠 are obtained iteratively by solving the linearized equations that result from substituting (2.3) in the governing equation (2.1) and using 𝑢₀, 𝜔₀, and 𝛼₀ as initial approximations. Substituting (2.3) in (2.1) and linearizing (neglecting nonlinear terms in 𝜔_𝑠, 𝛼_𝑠,𝑢_𝑠,  𝑢^′_𝑠 and 𝑢_𝑠^′′) yield 𝑎_1,𝑠−1(𝑧)𝑢_𝑠^′′+𝑎_2,𝑠−1(𝑧)𝑢^′_𝑠+𝑎_3,𝑠−1(𝑧)𝑢_𝑠+𝑎_4,𝑠−1(𝑧)𝜔_𝑠+𝑎_5,𝑠−1(𝑧)𝛼_𝑠=𝑅_𝑠−1(𝑧), (2.4) subject to the boundary conditions 𝑢_𝑠(0)=0,𝑢^′_𝑠(0)=0,𝑢_𝑠(2𝜋)=0,𝑢^′_𝑠(2𝜋)=0, (2.5) where 𝑎_1,𝑠−1=^𝑠−1_𝑚=0𝜔_𝑚²,𝑎_2,𝑠−1=−𝜀^𝑠−1_𝑚=0𝜔_𝑚+𝜀^𝑠−1_𝑚=0𝜔_𝑚^𝑠−1_𝑚=0𝑢_𝑚^2𝑠−1_𝑚=0𝛼_𝑚²,𝑎_3,𝑠−1=1+2𝜀^𝑠−1_𝑚=0𝜔_𝑚^𝑠−1_𝑚=0𝛼_𝑚^2𝑠−1_𝑚=0𝑢^′_𝑚^𝑠−1_𝑚=0𝑢_𝑚,𝑎_4,𝑠−1=2^𝑠−1_𝑚=0𝜔_𝑚^𝑠−1_𝑚=0𝑢_𝑚^′′−𝜀^𝑠−1_𝑚=0𝑢^′_𝑚+𝜀^𝑠−1_𝑚=0𝛼_𝑚^2𝑠−1_𝑚=0𝑢_𝑚^2𝑠−1_𝑚=0𝑢^′_𝑚,𝑎_5,𝑠−1=2𝜀^𝑠−1_𝑚=0𝜔_𝑚^𝑠−1_𝑚=0𝛼_𝑚^𝑠−1_𝑚=0𝑢^′_𝑚^𝑠−1_𝑚=0𝑢_𝑚²,𝑅_𝑠−1=−^𝑠−1_𝑚=0𝜔_𝑚^2𝑠−1_𝑚=0𝑢_𝑚^′′+𝜀^𝑠−1_𝑚=0𝜔_𝑚^𝑠−1_𝑚=0𝑢^′_𝑚−^𝑠−1_𝑚=0𝑢_𝑚−𝜀^𝑠−1_𝑚=0𝜔_𝑚^𝑠−1_𝑚=0𝑢_𝑚^2𝑠−1_𝑚=0𝑢^′_𝑚^𝑠−1_𝑚=0𝛼_𝑚². (2.6) The linearized system (2.4) was solved using the Chebyshev collocation spectral method. We used the transformation 𝑧=(𝜏+1)/2 to map [0,2𝜋] to [−1,1]. The solution space is discretized using the Chebyshev-Gauss-Lobatto collocation points 𝜏_𝑗=cos𝜋𝑗𝑁,𝑗=0,1,…𝑁, (2.7) which are the extrema of the 𝑁th order Chebyshev polynomial 𝑇_𝑁(𝜏)=cos𝑁cos⁻¹𝜏. (2.8) We approximated the derivatives of the unknown variables 𝑢_𝑠 at the collocation points as the matrix vector product: 𝑑𝑢_𝑠=𝑑𝑧^𝑁_𝑘=0𝐃_𝑗𝑘𝑢_𝑠𝜏_𝑘=𝐃𝐔_𝑠,𝑗=0,1,…,𝑁, (2.9) where 𝐃=𝐷/𝜋, with 𝐷 being the Chebyshev derivative matrix and 𝐔 is the vector function of 𝑢_𝑠 at the collocation points 𝜏_𝑗 (see, for example, [14, 15]). The entries of 𝐷 can be computed in different ways. In this work we use the method proposed by Trefethen [15] in the cheb.m MATLAB m-file. If we denote the entries of the derivative matrix 𝐷 by 𝐷_𝑗𝑘, 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: 𝐀_𝑠−1𝐔_𝑠+𝐚_4,𝑠−1𝜔_𝑠+𝐚_5,𝑠−1𝛼_𝑠=𝐑_𝑠−1, (2.10) where 𝐀_𝑠−1=𝐚_1,𝑠−1𝐃²+𝐚_2,𝑠−1𝐃+𝐚_3,𝑠−1, (2.11) and 𝐚_{𝑟,𝑠−1}, (𝑟=1,2,3) are the diagonal matrices corresponding to 𝑎_{𝑟,𝑠−1}(𝜏_𝑗). The vectors 𝐚_{𝑟,𝑠−1}, (𝑟=4,5) are of size (𝑁+1)×1 and they correspond to 𝑎_{𝑟,𝑠−1}(𝜏_𝑗). The appropriate boundary conditions are 𝑢_𝑠𝜏_𝑁=0,𝑢_𝑠𝜏₀=0,^𝑁_𝑘=0𝐷_0𝑘𝑢_𝑠𝜏_𝑘=0,^𝑁_𝑘=0𝐷_𝑁𝑘𝑢_𝑠𝜏_𝑘=0. (2.12)

After imposing the conditions (2.12) on (2.10) we obtain the following matrix system: 𝐀_𝑠−1𝑎_4,𝑠−1𝜏₁⋮𝑎_4,𝑠−1𝜏_𝑁−1𝑎_5,𝑠−1𝜏₁⋮𝑎_5,𝑠−1𝜏_𝑁−1𝐷₀₁⋯𝐷_0𝑁−100𝐷_𝑁1⋯𝐷_𝑁𝑁−100𝑢_𝑠𝜏₁⋮𝑢_𝑠𝜏_𝑁−1𝜔_𝑠𝛼_𝑠=𝑅_𝑠−1𝜏₁⋮𝑅_𝑠−1𝜏_𝑁−100. (2.13) The boundary conditions 𝑢_𝑠(𝜏_𝑁)=0 and 𝑢_𝑠(𝜏₀)=0 appearing in (2.12) have been imposed by deleting the first and last rows of 𝐀_𝑠−1 and 𝑅_𝑠−1. Two additional rows are also added to accommodate the two derivative boundary conditions. Thus, starting from the initial approximations 𝑢₀(𝑧), 𝜔₀, and 𝛼₀, the subsequent solutions for 𝑢_𝑠(𝑧), 𝜔_𝑠, 𝛼_𝑠  (𝑠=1,2,3…) can be obtained by solving the matrix system (2.13).

2.1. SLM Solution of the Duffing Equation

The exact solution of (1.2) is given (see [7, 8]) in terms of the relation between 𝛼 and 𝐿 as 𝐿=2^𝛼₀𝑑𝑡𝛼²−𝑧²−𝛼⁴−𝑧⁴,/2 (2.14) which gives the exact solution 𝜀^1/2=2𝜋1−𝛼²𝐾𝛼/22−²𝛼², (2.15) where 𝐾(𝜁) 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]: 𝜋𝛼=𝑣2,𝑣(𝑡)=𝛼𝑢(𝑡). (2.16) Substituting in (1.2), we have 𝑢^′′+𝜀𝑢−𝛼²𝑢³=0,𝑢(0)=𝑢(𝜋)=0. (2.17) 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 𝑢𝜋2=1. (2.18) Equations (2.17) and (2.18) form an eigenvalue problem with 𝛼 being the unknown eigenvalue. We now look for solutions 𝑢(𝑡) and 𝛼 in the form 𝑢(𝑡)=𝑢_𝑠(𝑡)+^𝑠−1_𝑚=0𝑢_𝑚(𝑡),𝛼=𝛼_𝑠+^𝑠−1_𝑚=0𝛼_𝑚,𝑠=1,2,3,…, (2.19) where 𝑢_𝑠 and 𝛼_𝑠 are unknown functions. Substituting (2.19) in (2.17) and neglecting nonlinear terms give 𝑢_𝑠^′′+𝑏_1,𝑠−1𝑢_𝑠+𝑏_2,𝑠−1𝛼_𝑠=𝑃_𝑠−1, (2.20) subject to the boundary conditions 𝑢_𝑠(0)=0,𝑢_𝑠𝜋2=0,𝑢_𝑠(𝜋)=0, (2.21) where 𝑏_1,𝑠−1=𝜀−3𝜀^𝑠−1_𝑚=0𝑢_𝑚^2𝑠−1_𝑚=0𝛼_𝑚²,𝑏_2,𝑠−1=−2𝜀^𝑠−1_𝑚=0𝛼_𝑚^𝑠−1_𝑚=0𝑢_𝑚³,𝑃_𝑠−1=−^𝑠−1_𝑚=0𝑢_𝑚^′′−𝜀^𝑠−1_𝑚=0𝑢_𝑚+𝜀^𝑠−1_𝑚=0𝑢_𝑚^3𝑠−1_𝑚=0𝛼_𝑚². (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 𝑧=𝜋(𝜏+1)/2 to map [0,𝜋] to [−1,1] gives 𝐁_𝑠−1𝐔_𝑠+𝐛_2,𝑠−1𝛼_𝑠=𝐏_𝑠−1, (2.23) with the boundary conditions 𝑢_𝑠𝜏_𝑁=0,𝑢_𝑠𝜏_𝑁/2+1=0,𝑢_𝑠𝜏₀=0, (2.24) where 𝐁_𝑠−1=𝐃²+𝐛_1,𝑠−1, (2.25) and 𝐛_{𝑟,𝑠−1}, (𝑟=1,2) are the diagonal matrices corresponding to 𝑏_{𝑟,𝑠−1}(𝜏_𝑗). After imposing the conditions (2.24) on (2.23) we obtain the following matrix system: 𝐁_𝑠−1𝑏_2,𝑠−1𝜏₁⋮𝑏_2,𝑠−1𝜏_𝑁/2−20𝑏_2,𝑠−1𝜏_𝑁/2⋮𝑏_2,𝑠−1𝜏_𝑁−10⋯010⋯00𝑢_𝑠𝜏₁⋮𝑢_𝑠𝜏_𝑁/2−2𝑢_𝑠𝜏_𝑁/2−1𝑢_𝑠𝜏_𝑁/2⋮𝑢_𝑠𝜏_𝑁−1𝛼_𝑠=𝑃_𝑠−1𝜏₁⋮𝑃_𝑠−1𝜏_𝑁/2−20𝑃_𝑠−1𝜏_𝑁/2⋮𝑃_𝑠−1𝜏_𝑁−10. (2.26) Starting from suitable initial approximations, 𝑢₀ and 𝛼₀, (2.26) can be used to find approximate 𝑢 and 𝛼.

3. 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 𝑢₀=cos(𝑧) while for the Duffing equation (1.2) 𝑢₀=sin(𝑡) was used. In both cases the initial values of all eigenvalues were set to be 1. The number of collocations points used was 𝑁=200 and 𝑁=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 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 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 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.

Table 2 

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

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.

Figure 1 

Limit cycle when 𝜀=0.5. Solid line: 6th order SLM solution. Circles: numerical solution.

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).

Table 3 

Comparison of the SLM results and the HAM and Homotopy-Padé [8] approximate results for 𝛼 at different orders 𝑝 for 𝛼 for 𝜀=10 and 𝜀=25.

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.

Figure 2 

Comparison of the exact solution (dots) with the 4th order SLM solution (line).

4. 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.