# Analytical and Numerical Study to a Forced Van der Pol Oscillator

Authors: Alvaro Salas; Lorenzo J. Martínez H; David L. Ocampo R
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

In this paper we show the way to apply some and analytical and numerical techniques in order to solve the forced Van der Pol oscillator. We illustrate the obtained results with examples. A comparison with Runge–Kutta numerical method is made in order to see the accuracy of the approximated analytical solution.


---

## Body

## 1. Introduction

We consider the following forced Van der Pol oscillator:(1)x¨−ε1−x2x˙+αx=FcosΩt,x0=x0,x′0=x˙0,In the case when α=1 and β=F=0 we obtain the well known Van der Pol oscillator [1].(2)x¨−ε1−x2x˙+x=0,x0=x0,x′0=x˙0.Oscillator (2) may be solved using perturbation methods as described.


## 2. The Lindstedt–Poincaré Method

Let us consider the i.v.p.(3)x¨−ε1−x2x˙+ω02x=0,x0=2,x′0=0.The Lindstedt–Poincaré method assumes the solution in the ansatz form(4)xt=y0ωt+εy1ωt+ε2y2ωt+ε3y3ωt+ε4y4ωt+⋯,where(5)ω2=ω02+εω1+ε2ω2+ε3ω3+ε4ω4+⋯.Using this method gives the following solution up toε4:(6)xt=2cosωt+34ω0sinωt−14ω0sin3ωtε+596ω02cos3ωt−596ω02cos5ωtε2+−7256ω03sinωt+21256ω03sin3ωt−35576ω03sin5ωt+7576ω03sin7ωtε3+−6301276480ω04cos3ωt+108527648ω04cos5ωt−2149110592ω04cos7ωt+6120480ω04cos9ωtε4,whereω=ω02−18ε2+231536ω02ε4.This solution is periodic with periodT=2π/ω.


## 3. He’s Homotopy Perturbation Method

The homotopy perturbation method has been shown to solve a large class of nonlinear differential problems effectively, easily, and accurately; generally one iteration is enough for engineering applications with acceptable accuracy, making the method accessible to nonmathematicians. However, in this method secularity terms appears. We perform iterations until such terms appear. We construct the following homotopy [2].(7)x¨+ω02x−pε1−x2x˙+αp1−px=0,where(8)xt=x0ωt+px1ωt+p2x2ωt+p3x3ωt+p4x4ωt+⋯,ω2=ω02+pω1+p2ω2+p3ω3+p4ω4+⋯.After some algebra one finds that(9)xt=4−c02sinτ+c0cosτ+εsin3τ−34−c02c02ε+34−c02ε−8c03ω0+24c0ω064ω02+εcos3τ−3c03ε+9c0ε+84−c02c02ω0−84−c02ω064ω02+54−c02c04−3c02+1ε2sin5τ192ω02−5c0c04−5c02+5ε2cos5τ192ω02.Here,τ=c1+ω02−ε28tα=0,ω1=0,ω2=−ε28.The constantsc0 and c1 are determined from nthe initial conditions.


## 4. The Krylov–Bogoliubov–Mitropolsky Method (KBM)

The Krylov–Bogoliubov–Mitropolsky method (KBM) [2], [3] is a technique to give approximate analytical solution to the weakly nonlinear second-order equation(10)d2udt2+ω02u=εfu,dudt.Whenε=0 the solution of (10) may be expressed as(11)u=acosω0+θ.where a and θ are constants. For the case when ε>0 is small, Krylov and Bogoliubov (1947) assumed that the solution is still given by (11) but with time-varying a and θ, and subject to the condition(12)dudt=−aω0sinϕ,ϕ=ω0t+θ.In the general case, the solution is assumed in the ansatz form(13)u=acosψ+∑n=1Nεnuna,ψ+OεN+1.where each un is a periodic function of ψ with a period 2π, and a and ψ are assumed to vary with time according to(14)dadt=∑n=1NεnAna+OεN+1,dψdt=ω0+∑n=1Nεnψna+OεN+1.In order to uniquely determineAn and ψn, we require that no un contains cosψ. Let N=3. Then,(15)dudt=−aω0sinψ+ω0u1,ψ−aψ1sinψ+A1cosψε+A1u1,a+ω0u2,ψ+ψ1u1,ψ−aψ2sinψ+A2cosψε2+A2u1,a+A1u2,a+ω0u3,ψ+ψ2u1,ψ+ψ1u2,ψ−aψ3sinψ+A3cosψε3+⋯,(16)d2udt2=−aω02cosψ+ω02u1,ψψ−2aψ1ω0cosψ−2A1ω0sinψε+2A1ω0u1,aψ+2ψ1ω0u1,ψψ+ω02u2,ψψ+sinψ−aA1ψ˙1−2A1ψ1−2A2ω0+cosψA1A˙1−a2ψ2ω0+ψ12ε2+A˙1A1u1,a+A12u1,aa+2A1ψ1u1,aψ+2A1ω0u2,aψ+2A2ω0u1,aψ+A1ψ˙1u1,ψ+2ψ2ω0u1,ψψ+2ψ1ω0u2,ψψ+ψ12u1,ψψ+ω02u3,ψψ+sinψ−aA2ψ˙1−aA1ψ˙2−2A2ψ1−2A1ψ2−2A3ω0+cosψ−2aψ3ω0−2aψ1ψ2+A2A˙1+A1A˙2ε3+⋯.Here,(17)ψ˙1=ψ1′t,u2,ψ=∂u2∂ψ,u2,ψψ=∂2u2∂ψ2,u1,aψ=∂2u1∂ψ∂a,etc.Following is the solution obtained using KBM to accuracyOε3 [2]:(18)x=xt=acosψ−132a3sin3ψε−a2+8a31024cos3ψ+5a53072cos5ψε2,where a=at and ψ=ψt obey the odds(19)a′t=ε2at−ε8a3t,ψ′t=1−ε28−ε2at4256+ε2at232,so that(20)at=2eεt/2e8c1+eεt,ψt=116−e8c1εe8c1+eεt+εloge8c1+eεt−2ε2−8t+c2.The constantsc1 and c2 are determined from the initial conditions x0=x0 and x′0=x˙0. For an illustration, see Figure 1 and Figure 2.Figure 1
The Vander Pol oscillator forε  = 0.75.Figure 2
The Vander Pol oscillator forε  = 0.75.


### 4.1. A More Accurate Analytical Solution

Let(21)u¨−ε1−u2u˙+ω02u=0,u0=u0,u′0=u˙0.In view of the KBM, we assume the ansatz form(22)u=acosψ+∑n=13εnuna,ψ+Oε4.Hence,(23)dadt=∑n=13εnAna+Oε4,dψdt=ω0+∑n=13εnψna+Oε4.Making use of (15) and (16), we will have(24)u¨−ε1−u2u˙+ω02u=H1ε+H2ε2+H3ε3+Oε4,where(25)H1=144ω02u1,ψψ+4u1ω02+14sinψ−ω0a3+4ω0a−8A1ω0−14a3ω0sin3ψ−2aψ1ω0cosψ,(26)H2=142a2ω0u1,ψ+8A1ω0u1,aψ+8ψ1ω0u1,ψψ−4ω0u1,ψ+4ω02u2,ψψ+4u2ω02+12a2ω0cos2ψu1,ψ+14sinψa3−ψ1−4aA1ψ˙1+4aψ1−8A1ψ1−8A2ω0−14a3ψ1sin3ψ+14cosψ3a2A1−8aψ2ω0−4aψ12−4A1+4A1A˙1+14a2A1cos3ψ−a2u1ω0sin2ψ,(27)H3=cosψ2au1ω0u1,ψ+3a2A24−2aψ3ω0−2aψ1ψ2−A2+A2A˙1+A1A˙2+cos2ψ12a2A1u1,a+12a2ω0u2,ψ+12a2ψ1u1,ψ+aA1u1+12a2A1u1,a+12a2ω0u2,ψ+12a2ψ1u1,ψ−A1u1,a+A1A˙1u1,a+A12u1,aa+2A1ψ1u1,aψ+2A2ω0u1,aψ+2A1ω0u2,aψ+A1ψ˙1u1,ψ+2ψ2ω0u1,ψψ+2ψ1ω0u2,ψψ+ψ12u1,ψψ−ω0u2,ψ−ψ1u1,ψ+ω02u3,ψψ+sinψ−14a3ψ2−aA2ψ˙1−aA1ψ˙2−au12ω0+aψ2−2A2ψ1−2A1ψ2−2A3ω0−14a3ψ2sin3ψ+14a2A2cos3ψ+sin2ψa2−u1ψ1−a2u2ω0+aA1u1+u3ω02.Equating to zero the coefficients ofcosψ and sinψ in (25) gives(28)A1=184a−a3,ψ1=0.From the conditionH1=0, we get the second order linear ode(29)144ω02u1,ψψ+4u1ω02−14a3ω0sin3ψ=0.The general solution to the ode (29) reads(30)u1a,ψ=32c1ω0232ω02cosψ+32c2ω02−2a3ω032ω02sinψ−a332ω0sin3ψ.Define(31)c1=0,c2=116ω0a3.Then,(32)u1=u1a,ψ=−132ω0a3sin3ψ.Replacing the expressions (31) and (32) into (26), we obtain(33)H2=128ω04u2a,ψ+128ω04u00,2a,ψ128ω02−5a5ω02128ω02cos5ψ−256A2ω03128ω02sinψ−132a3a2+2cos3ψ+316a3−14a−5128a5−2aω0ψ2cosψ,Equating to zero the coefficients ofsinψ and sinψ in (33), we get(34)A2=0,ψ2=−1256ω05a4−24a2+32.From the conditionH2=0 we obtain the following linear second order ode(35)ω02u2a,ψ+u20,2a,ψ−5a5ω02128ω02cos5ψ−132a3a2+2cos3ψ=0.The general solution to the ode (35) is given by(36)u2=u2a,ψ=−5a5ω023072ω04cos5ψ+c4sinψ−1256a3ω02a2+2cos3ψ−1128ω02a5+2a3−128c3ω02cosψ.Define(37)c3=1128a3ω02a2+2,c4=0.Then,(38)u2=u2a,ψ=−5a5ω023072ω04cos5ψ−1256a3ω02a2+2cos3ψ.Finally, we must haveH3=0. Proceeding in a similar manner as before, we get(39)A3=−18192ω0225a7−184a5+320a3,ψ3=0.u3=28a7sin7ψ−671a2+280a5sinψ+53a2+56a5sin5ψ−929a4−168a2+64a3sin3ψ294912ω03.We thus have that(40)u=acosψ−a3sin3ψ32ω0ε−5a5cos5ψ3072ω02+a2+8a3cos3ψ1024ω02ε2+7a773728ω03sin7ψ+53a2+56a5294912ω03sin5ψ−29a4−168a2+64a332768ω03sin3ψε3+Oε4.Also, from (23),(41)a′t=12εat−5ε3128ω02+ε8at3+23ε31024ω02at5−25ε38192ω02at7,Ode (41) is hard to solve in closed form. We use the approximation(42)aε2−5ε3128ω02+ε8a3+23ε31024ω02a5−25ε38192ω02a7≈ra+sa3fora≤M,where,r=ε2+ε325M2−184M465536ω02,s=−ε8−ε3175M4−1472M2+256065536ω02.Solving the ode(43)dadt=ra+sa3,one gets(44)a=at=±a0rertr+a02s1−a02se2rt/a02s+r.The expression forψ=ψt is obtained from the ode ψ′t=ω0−ε2/8ω0+ε2/32ω0at2−ε2/256ω0at4:(45)ψt=1512s2ω0r−64s2t−r+8slog1−a02sre2rt−1ε2−ε2ra02s+rr−a02se2rt−1+sa02ε2+512sω0ω0t+ψ0.Example 1.
Let(46)u0=u0=−0.104,u′0=u˙0=−0.25,ε=0.75,ω0=1.
The i.v.p to be solved reads(47)u¨−0.751−u2u˙+u=0,u0=−0.104,u′0=−0.25.
From (43),(48)r=0.366348,s=−0.0903511.
We have(49)at=0.251949e0.366348t1+0.0156555e0.732697t,ψt=0.239256−0.09035110.0351563−46.2598t+2−0.074330.00564694e0.732697t−1+0.366348+−0.522453t−0.56250.356461−log0.00564694e0.732697t−1+0.366348−1.00417+0.366348.
See Figure3.Figure 3
Comparison between analytical and numerical solution.


## 4.1. A More Accurate Analytical Solution

Let(21)u¨−ε1−u2u˙+ω02u=0,u0=u0,u′0=u˙0.In view of the KBM, we assume the ansatz form(22)u=acosψ+∑n=13εnuna,ψ+Oε4.Hence,(23)dadt=∑n=13εnAna+Oε4,dψdt=ω0+∑n=13εnψna+Oε4.Making use of (15) and (16), we will have(24)u¨−ε1−u2u˙+ω02u=H1ε+H2ε2+H3ε3+Oε4,where(25)H1=144ω02u1,ψψ+4u1ω02+14sinψ−ω0a3+4ω0a−8A1ω0−14a3ω0sin3ψ−2aψ1ω0cosψ,(26)H2=142a2ω0u1,ψ+8A1ω0u1,aψ+8ψ1ω0u1,ψψ−4ω0u1,ψ+4ω02u2,ψψ+4u2ω02+12a2ω0cos2ψu1,ψ+14sinψa3−ψ1−4aA1ψ˙1+4aψ1−8A1ψ1−8A2ω0−14a3ψ1sin3ψ+14cosψ3a2A1−8aψ2ω0−4aψ12−4A1+4A1A˙1+14a2A1cos3ψ−a2u1ω0sin2ψ,(27)H3=cosψ2au1ω0u1,ψ+3a2A24−2aψ3ω0−2aψ1ψ2−A2+A2A˙1+A1A˙2+cos2ψ12a2A1u1,a+12a2ω0u2,ψ+12a2ψ1u1,ψ+aA1u1+12a2A1u1,a+12a2ω0u2,ψ+12a2ψ1u1,ψ−A1u1,a+A1A˙1u1,a+A12u1,aa+2A1ψ1u1,aψ+2A2ω0u1,aψ+2A1ω0u2,aψ+A1ψ˙1u1,ψ+2ψ2ω0u1,ψψ+2ψ1ω0u2,ψψ+ψ12u1,ψψ−ω0u2,ψ−ψ1u1,ψ+ω02u3,ψψ+sinψ−14a3ψ2−aA2ψ˙1−aA1ψ˙2−au12ω0+aψ2−2A2ψ1−2A1ψ2−2A3ω0−14a3ψ2sin3ψ+14a2A2cos3ψ+sin2ψa2−u1ψ1−a2u2ω0+aA1u1+u3ω02.Equating to zero the coefficients ofcosψ and sinψ in (25) gives(28)A1=184a−a3,ψ1=0.From the conditionH1=0, we get the second order linear ode(29)144ω02u1,ψψ+4u1ω02−14a3ω0sin3ψ=0.The general solution to the ode (29) reads(30)u1a,ψ=32c1ω0232ω02cosψ+32c2ω02−2a3ω032ω02sinψ−a332ω0sin3ψ.Define(31)c1=0,c2=116ω0a3.Then,(32)u1=u1a,ψ=−132ω0a3sin3ψ.Replacing the expressions (31) and (32) into (26), we obtain(33)H2=128ω04u2a,ψ+128ω04u00,2a,ψ128ω02−5a5ω02128ω02cos5ψ−256A2ω03128ω02sinψ−132a3a2+2cos3ψ+316a3−14a−5128a5−2aω0ψ2cosψ,Equating to zero the coefficients ofsinψ and sinψ in (33), we get(34)A2=0,ψ2=−1256ω05a4−24a2+32.From the conditionH2=0 we obtain the following linear second order ode(35)ω02u2a,ψ+u20,2a,ψ−5a5ω02128ω02cos5ψ−132a3a2+2cos3ψ=0.The general solution to the ode (35) is given by(36)u2=u2a,ψ=−5a5ω023072ω04cos5ψ+c4sinψ−1256a3ω02a2+2cos3ψ−1128ω02a5+2a3−128c3ω02cosψ.Define(37)c3=1128a3ω02a2+2,c4=0.Then,(38)u2=u2a,ψ=−5a5ω023072ω04cos5ψ−1256a3ω02a2+2cos3ψ.Finally, we must haveH3=0. Proceeding in a similar manner as before, we get(39)A3=−18192ω0225a7−184a5+320a3,ψ3=0.u3=28a7sin7ψ−671a2+280a5sinψ+53a2+56a5sin5ψ−929a4−168a2+64a3sin3ψ294912ω03.We thus have that(40)u=acosψ−a3sin3ψ32ω0ε−5a5cos5ψ3072ω02+a2+8a3cos3ψ1024ω02ε2+7a773728ω03sin7ψ+53a2+56a5294912ω03sin5ψ−29a4−168a2+64a332768ω03sin3ψε3+Oε4.Also, from (23),(41)a′t=12εat−5ε3128ω02+ε8at3+23ε31024ω02at5−25ε38192ω02at7,Ode (41) is hard to solve in closed form. We use the approximation(42)aε2−5ε3128ω02+ε8a3+23ε31024ω02a5−25ε38192ω02a7≈ra+sa3fora≤M,where,r=ε2+ε325M2−184M465536ω02,s=−ε8−ε3175M4−1472M2+256065536ω02.Solving the ode(43)dadt=ra+sa3,one gets(44)a=at=±a0rertr+a02s1−a02se2rt/a02s+r.The expression forψ=ψt is obtained from the ode ψ′t=ω0−ε2/8ω0+ε2/32ω0at2−ε2/256ω0at4:(45)ψt=1512s2ω0r−64s2t−r+8slog1−a02sre2rt−1ε2−ε2ra02s+rr−a02se2rt−1+sa02ε2+512sω0ω0t+ψ0.Example 1.
Let(46)u0=u0=−0.104,u′0=u˙0=−0.25,ε=0.75,ω0=1.
The i.v.p to be solved reads(47)u¨−0.751−u2u˙+u=0,u0=−0.104,u′0=−0.25.
From (43),(48)r=0.366348,s=−0.0903511.
We have(49)at=0.251949e0.366348t1+0.0156555e0.732697t,ψt=0.239256−0.09035110.0351563−46.2598t+2−0.074330.00564694e0.732697t−1+0.366348+−0.522453t−0.56250.356461−log0.00564694e0.732697t−1+0.366348−1.00417+0.366348.
See Figure3.Figure 3
Comparison between analytical and numerical solution.


## 5. Analytical Solution to a Forced Van der Pol Oscillator

Let us consider the i v p.(50)x¨−ε1−x2x˙+ω02x=FcosΩt.given that(51)x0=x0,x′0=x˙0.We will assume the ansatz(52)xt=ut+CcosΩt+DsinΩt.Here, the functionu=ut is the solution to the i v p(53)u¨−ε1−u2u˙+ω02u=0,u0=u0=x0−C,u′0=u˙0=x˙0−DΩ.We now will determine suitable values for the constantsC and D. We have(54)x¨−ε1−x2x˙+ω02x−FcosΩt=ω02ut−εsin2ΩtC2Ωut−CDu˙0−D2Ωut+12εcos2ΩtC2u˙0+4CDΩut−D2u˙0+12εu˙0C2+D2−14cosΩt⋅−4εut2DΩ−8Cεu˙0ut−C2DεΩ−4Cω02+4CΩ2−D3εΩ+4DεΩ+4F−14εsin3ΩtC3Ω−3CD2Ω+14εcos3Ωt3C2DΩ−3βCD2−D3Ω−14sinΩt⋅4CΩεut2−8Dεu˙0ut+C3εΩ−3βC2Dε+CD2εΩ−4Dω02+4DΩ2−4CεΩ.We choose the values ofC and D so that(55)−C2DεΩ−4Cω02+4CΩ2−D3εΩ+4DεΩ+4F=0,C3εΩ−3βC2Dε+CD2εΩ−4Dω02+4DΩ2−4CεΩ=0.From (55), it follows that(56)ε2F2Ω2C3+8ε2FΩ2Ω2−ω02C2+16Ω2−ω022ε2Ω2+Ω4−2Ω2ω02+ω04C+16FΩ2−ω023=0.We choose the least in magnitude real rootC to (56) and the least in magnitude real root D to (57).Example 2.
Let0≤t≤T=130, ε=0.1, ω0=1, F=1,Ω=5.2, x0=0 and x˙0=−0.183. Consider the i.v.p.(57)x¨t−0.11−xt2x˙t+xt=cos5.2t∧x0=0∧x′0=−0.183.
Our calculations give(58)a0=0.185483,ψ0=4.50402,C=−0.0383872.D=−0.000766281.r=0.049979,s=−0.0124918.at=−0.186053e0.049979t1+0.00865188e0.099958t.ψt=0.00625262log0.0495503+0.000428704e0.099958t+0.99875t+=1−1.38308e0.099958t−159.859+4.52869.
The [4] expression for u=ut is [5] obtained from (40). In Figure 4, we compare[6] the approximate [7] analytical solution (dashed curve) with the Runge–Kutta numerical solution. The error equals 0.0756086.Figure 4
Comparison between the approximate analytical solution (dashed curve) and the Runge–Kutta numerical.



---

*Source: 9736427-2022-03-25.xml*