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

Authors: Sandile S. Motsa; Precious Sibanda
Journal: Mathematical Problems in Engineering
 (2012)

---

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


---

## Body

## 1. Introduction

The Van der Pol equation is(1.1)x¨+x+ε(x2-1)x˙=0,x(0)=α,x˙(0)=0,
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 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(1.2)x′′+ε(x-x3)=0,x(0)=0,x(π)=0,
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.


## 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 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
(2.1)ω2u′′+u-εωu'+εωα2u′u2=0,u(0)=1,u′(0)=0,
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:
(2.2)u(2π)=1,u′(2π)=0.
We now seek a solution of the form:
(2.3)u(z)=us(z)+∑m=0s-1‍um(z),ω=ωs+∑m=0s-1‍ωm,α=αs+∑m=0s-1‍αm,s=1,2,3,…,
where us, ωs, αs are obtained iteratively by solving the linearized equations that result from substituting (2.3) in the governing equation (2.1) and using u0, ω0, and α0 as initial approximations. Substituting (2.3) in (2.1) and linearizing (neglecting nonlinear terms in ωs, αs,us,  us′ and us′′) yield
(2.4)a1,s-1(z)us′′+a2,s-1(z)us′+a3,s-1(z)us+a4,s-1(z)ωs+a5,s-1(z)αs=Rs-1(z),
subject to the boundary conditions
(2.5)us(0)=0,us′(0)=0,us(2π)=0,us′(2π)=0,
where
(2.6)a1,s-1=(∑m=0s-1‍ωm)2,a2,s-1=-ε∑m=0s-1‍ωm+ε(∑m=0s-1‍ωm)(∑m=0s-1‍um)2(∑m=0s-1‍αm)2,a3,s-1=1+2ε(∑m=0s-1‍ωm)(∑m=0s-1‍αm)2(∑m=0s-1‍um′)(∑m=0s-1‍um),a4,s-1=2∑m=0s-1‍ωm∑m=0s-1‍um′′-ε∑m=0s-1‍um′+ε(∑m=0s-1‍αm)2(∑m=0s-1‍um)2(∑m=0s-1‍um′),a5,s-1=2ε(∑m=0s-1‍ωm)(∑m=0s-1‍αm)(∑m=0s-1‍um′)(∑m=0s-1‍um)2,Rs-1=-(∑m=0s-1‍ωm)2∑m=0s-1‍um′′+ε(∑m=0s-1‍ωm)(∑m=0s-1‍um′)-∑m=0s-1‍um-ε(∑m=0s-1‍ωm)(∑m=0s-1‍um)2(∑m=0s-1‍um′)(∑m=0s-1‍αm)2.
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
(2.7)τj=cos(πjN),j=0,1,…N,
which are the extrema of the Nth order Chebyshev polynomial
(2.8)TN(τ)=cos(Ncos-1τ).
We approximated the derivatives of the unknown variables us at the collocation points as the matrix vector product:
(2.9)dusdz=∑k=0N‍Djkus(τk)=DUs,j=0,1,…,N,
where D=D/π, with D being the Chebyshev derivative matrix and U is the vector function of us 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.m MATLAB m-file. If we denote the entries of the derivative matrix D by Djk, 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:
(2.10)As-1Us+a4,s-1ωs+a5,s-1αs=Rs-1,
where
(2.11)As-1=a1,s-1D2+a2,s-1D+a3,s-1,
and ar,s-1, (r=1,2,3) are the diagonal matrices corresponding to ar,s-1(τj). The vectors ar,s-1, (r=4,5) are of size (N+1)×1 and they correspond to ar,s-1(τj). The appropriate boundary conditions are
(2.12)us(τN)=0,us(τ0)=0,∑k=0N‍D0kus(τk)=0,∑k=0N‍DNkus(τk)=0.After imposing the conditions (2.12) on (2.10) we obtain the following matrix system:
(2.13)[As-1a4,s-1(τ1)⋮a4,s-1(τN-1)a5,s-1(τ1)⋮a5,s-1(τN-1)D01⋯D0N-100DN1⋯DNN-100][us(τ1)⋮us(τN-1)ωsαs]=[Rs-1(τ1)⋮Rs-1(τN-1)00].
The boundary conditions us(τN)=0 and us(τ0)=0 appearing in (2.12) have been imposed by deleting the first and last rows of As-1 and Rs-1. Two additional rows are also added to accommodate the two derivative boundary conditions. Thus, starting from the initial approximations u0(z), ω0, and α0, the subsequent solutions for us(z), ωs, αs  (s=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 L as
(2.14)L=2∫0α‍dtα2-z2-(α4-z4)/2,
which gives the exact solution
(2.15)ε1/2=2π1-α2/2K(α22-α2),
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]:
(2.16)α=v(π2),v(t)=αu(t).
Substituting in (1.2), we have
(2.17)u′′+ε(u-α2u3)=0,u(0)=u(π)=0.
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.18)u(π2)=1.
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
(2.19)u(t)=us(t)+∑m=0s-1‍um(t),α=αs+∑m=0s-1‍αm,s=1,2,3,…,
where us and αs are unknown functions. Substituting (2.19) in (2.17) and neglecting nonlinear terms give
(2.20)us′′+b1,s-1us+b2,s-1αs=Ps-1,
subject to the boundary conditions
(2.21)us(0)=0,us(π2)=0,us(π)=0,
where
(2.22)b1,s-1=ε-3ε(∑m=0s-1‍um)2(∑m=0s-1‍αm)2,b2,s-1=-2ε(∑m=0s-1‍αm)(∑m=0s-1‍um)3,Ps-1=-∑m=0s-1‍um′′-ε∑m=0s-1‍um+ε(∑m=0s-1‍um)3(∑m=0s-1‍αm)2.
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
(2.23)Bs-1Us+b2,s-1αs=Ps-1,
with the boundary conditions
(2.24)us(τN)=0,us(τN/2+1)=0,us(τ0)=0,
where
(2.25)Bs-1=D2+b1,s-1,
and br,s-1, (r=1,2) are the diagonal matrices corresponding to br,s-1(τj). After imposing the conditions (2.24) on (2.23) we obtain the following matrix system:
(2.26)[Bs-1b2,s-1(τ1)⋮b2,s-1(τN/2-2)0b2,s-1(τN/2)⋮b2,s-1(τN-1)0⋯010⋯00][us(τ1)⋮us(τN/2-2)us(τN/2-1)us(τN/2)⋮us(τN-1)αs]=[Ps-1(τ1)⋮Ps-1(τN/2-2)0Ps-1(τN/2)⋮Ps-1(τN-1)0].
Starting from suitable initial approximations, u0 and α0, (2.26) can be used to find approximate u and α.


## 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 L as
(2.14)L=2∫0α‍dtα2-z2-(α4-z4)/2,
which gives the exact solution
(2.15)ε1/2=2π1-α2/2K(α22-α2),
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]:
(2.16)α=v(π2),v(t)=αu(t).
Substituting in (1.2), we have
(2.17)u′′+ε(u-α2u3)=0,u(0)=u(π)=0.
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.18)u(π2)=1.
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
(2.19)u(t)=us(t)+∑m=0s-1‍um(t),α=αs+∑m=0s-1‍αm,s=1,2,3,…,
where us and αs are unknown functions. Substituting (2.19) in (2.17) and neglecting nonlinear terms give
(2.20)us′′+b1,s-1us+b2,s-1αs=Ps-1,
subject to the boundary conditions
(2.21)us(0)=0,us(π2)=0,us(π)=0,
where
(2.22)b1,s-1=ε-3ε(∑m=0s-1‍um)2(∑m=0s-1‍αm)2,b2,s-1=-2ε(∑m=0s-1‍αm)(∑m=0s-1‍um)3,Ps-1=-∑m=0s-1‍um′′-ε∑m=0s-1‍um+ε(∑m=0s-1‍um)3(∑m=0s-1‍αm)2.
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
(2.23)Bs-1Us+b2,s-1αs=Ps-1,
with the boundary conditions
(2.24)us(τN)=0,us(τN/2+1)=0,us(τ0)=0,
where
(2.25)Bs-1=D2+b1,s-1,
and br,s-1, (r=1,2) are the diagonal matrices corresponding to br,s-1(τj). After imposing the conditions (2.24) on (2.23) we obtain the following matrix system:
(2.26)[Bs-1b2,s-1(τ1)⋮b2,s-1(τN/2-2)0b2,s-1(τN/2)⋮b2,s-1(τN-1)0⋯010⋯00][us(τ1)⋮us(τN/2-2)us(τN/2-1)us(τN/2)⋮us(τN-1)αs]=[Ps-1(τ1)⋮Ps-1(τN/2-2)0Ps-1(τN/2)⋮Ps-1(τN-1)0].
Starting from suitable initial approximations, u0 and α0, (2.26) can be used to find approximate u 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 u0=cos(z) while for the Duffing equation (1.2) u0=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).Table1 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.
Chen and Liu [9]
SLM results
Order
ω
α
Order
ω
α
10
0.94295340115
2.00861561331
2
0.93782084070
2.01344273793
20
0.94295586899
2.00861988872
4
0.94295585006
2.00861985316
50
0.94295584744
2.00861986087
6
0.94295584744
2.00861986087
100
0.94295584744
2.00861986087
8
0.94295584744
2.00861986087
Exact
0.94295584744
2.00861986087
Exact
0.94295584744
2.00861986087Table2 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].
(SLM)i
ε
0.3
0.6
1.2
1.5
1.8
2
0.9944195093
0.9781496961
0.8970288929
0.7933396502
0.7274495579
4
0.9944198444
0.9782168788
0.9211217519
0.8871147446
0.7394634260
6
0.9944198444
0.9782168788
0.9211220869
0.8854079093
0.8498705704
8
0.9944198444
0.9782168788
0.9211220869
0.8854079098
0.8482341694
10
0.9944198444
0.9782168788
0.9211220869
0.8854079098
0.8482341695
Buonomo [2]
0.9944198444
0.9782168788
0.9211220869
0.8854079098
0.8482341695
(HAM [9])
10
0.9944198397
0.8855524710
0.8485246885
20
0.9944198444
0.8854094218
0.8482487089
50
0.9944198444
0.8854079098
0.8482341635
100
0.9944198444
0.8854079098
0.8482341695Figure1 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 Table3 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 p for α for ε=10 and ε=25.
HAM results [8]
Homotopy-Padé [8]
SLM results
p
ε
=
10
ε
=
25
[
p, p]
ε
=
10
ε
=
10
p
ε
=
10
ε
=
25
5
0.99833
1.01046
[2, 2]
0.99914
1.01167
1
1.01834916
1.00009888
10
0.99588
1.00313
[4, 4]
0.99651
1.00113
2
1.00142002
0.99994805
15
0.99624
1.00117
[6, 6]
0.99644
1.00012
3
0.99684045
0.99994004
20
0.99644
1.00049
[8, 8]
0.99644
0.99996
4
0.99643870
0.99993998
25
0.99644
1.00017
[10, 10]
0.99644
0.99994
5
0.99643534
0.99993998
30
0.99644
1.00000
[12, 12]
0.99644
0.99994
6
0.99643534
0.99993998Figure2 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.



---

*Source: 693453-2012-09-27.xml*