home

PDF (letter size)

Solving the Van Der Pol nonlinear differential equation using first order approximation perturbation method

Nasser M. Abbasi

Sometime in 2009   Compiled on January 30, 2024 at 4:42am

Abstract

The Van Der Pol differential equation x(t)α(1x2(t))x(t)+x(t)=0 Was solved using perturbation with first order approximation. Two different solutions were obtained. The first solution restricted the initial conditions to be x(0)2+x2=4 which resulted in forcing function that caused resonance to be eliminated. This gave a stable solution but with initial conditions restricted to be near the origin of the phase plane space.

The second solution allowed arbitrary initial conditions any where in the phase plane but the the resulting forcing function caused resonance resulting in a solution which became unstable after some time.

Phase plane plots are used to compare the two solutions.

1 First solution: Restriction on initial condition. No resonance

The Van Der Pol equation is x¨α(1x2)x˙+x=0 Initial conditions are x(0)=φ and x˙(0)=ξ.

If α1, then the equation becomes x¨0+x0=0 which has an exact solution. We perturb the exact solution for x¨0+x0=0 to approximate the solution of the nonlinear equation x¨α(1x2)x˙+x=0 in terms of perturbation parameter α Let the solution of the Van Der Pol equation be x(t) and the solution to the exact equation x¨0+x0=0 be x0(t), then x(t)=x0(t)+αx1(t)+α2x2(t)+ First order approximation results in x(t)=x0(t)+αx1(t) Now we need to determine x0(t) and x1(t). Substituting (2) into (1) gives x¨0+αx¨1αx˙0+αx02x˙0+α3x12x˙0+2α2x1x0x˙0α2x˙1+α2x02x˙1+α4x12x˙1+2α3x1x0x˙1+x0+αx1=0 Collecting terms with same power of α gives α0(x¨0+x0)+α(x¨1x˙0+x02x˙0+x1)+α2(x02x˙1+2x1x0x˙0x˙1)+α3(x12x˙0+2x1x0x˙1)+α4x12x˙1=0 Setting terms which multiply by higher power α to zero and since it is assumed that α is very small therefore (x¨0+x0)+α(x¨1x˙0+x02x˙0+x1)=0 For the LHS to be zero implies that x¨0+x0=0 And x¨1x˙0+x02x˙0+x1=0 Equation (3) is solved for x0 and (4) is solved for x1 in order to find the solution x(t). Solution to (3) is x0(t)=A0cost+B0sint Assuming initial conditions for x1(t) are zero, then initial conditions for x0(t) can be taken to be those given for x(t). Solving for A0,B0 gives x0(t)=φcost+ξsint Substituting the solution just found (and its derivative) into (4) gives x¨1+x1=φsint+ξcostξ3costsin2t2ξ2φcos2tsint+ξ2φsin3tξφ2cos3t+2ξφ2costsin2t+φ3cos2tsint

Using sintcos2t=14(sint+sin3t) and costsin2t=14(costcos3t) the above can be simplified tox¨1+x1=(φξ2φ2+φ34)sint+(ξξ34+ξφ22)cost+(ξ34ξφ22)cos3t+(ξ2φ2+φ34)sin3t+ξφ(ξsin3tφcos3t)

Using ξsin3tφcos3t=14(3φcostφcos3t+3ξsintξsin3t) the above can be simplified further tox¨1+x1=(φξ2φ2+φ34+3ξ2φ4)sint+(ξξ34+ξφ2234ξφ2)cost(5)+(ξ34ξφ22ξφ24)cos3t+(ξ2φ2+φ34ξ2φ4)sin3t

The above is the equation we will now solve for x1(t), which has four forcing functions, hence four particular solutions. two of these particular solutions will cause a complete solution blows up as time increases (the first two in the RHS shown above). If we however restrict the initial conditions such that (φξ2φ2+φ34+3ξ2φ4)=0 And (ξξ34+ξφ2234ξφ2)=0 Then those forcing function will vanish. The above two equation result in the same restriction, which can be found to be4=ξ2+φ2 By restricting ξ2+φ2 to be exactly 4, then the solution we obtain for x1(t) from (5) will not blow up. Hence (4) can now be rewritten without the two forcing functions which caused resonance resulting inx¨1+x1=ξ(ξ23φ24)cos3t+φ(φ23ξ24)sin3t But φ2=4ξ2 and ξ2=4φ2, hence the above becomes after some simplification(6)x¨1+x1=ξ(ξ23)cos3t+φ(φ23)sin3t The homogeneous solution to the above is x1,h=c1cost+c2sint And we guess two particular solutions x1,p1=a1cos3t+a2sin3t and x1,p2=d1cos3t+d2sin3t. By substituting each of these particular solutions into (6) one at a time, and comparing coefficients, we can determine a1,a2,d1,d2. This results in x1,p1=ξ(3ξ2)8cos3t Similarly, x1,p2=φ(3φ2)8sin3t Hence the solution to (6) becomesx1(t)=x1,h+x1,p1+x1,p2(7)=(c1cost+c2sint)+ξ(3ξ2)8cos3t+φ(3φ2)8sin3t

Now applying initial conditions, which are x1(0)=0 and x˙1(0)=0, to (7) and determining c1 and c2 results inx1(t)=ξ(ξ23)8cost+38φ(φ23)sint+ξ(3ξ2)8cos3t+φ(3φ2)8sin3t Therefore we now have the final perturbation solution, which isx(t)=x0(t)+αx1(t)=φcost+ξsint+α(ξ(ξ23)8cost+38φ(φ23)sint+ξ(3ξ2)8cos3t+φ(3φ2)8sin3t)

Where x(0)=φ and x˙(0)=ξ and the above solution is valid under the restriction that x2(0)+x˙2(0)=4 We notice that the above restriction implies that both x(0) and x˙(0) have to be less than 4 to avoid getting a quantity which is complex. This implies that his solution is valid near the center of the phase plane only. In addition, it is valid only for small α. To illustrate this solution, We show the phase plot which results from the above solution, and also plot the solution x(t). Selecting x(0)=1.5 and x˙(0)=41.52=1.3229

2 Second solution. No restriction on initial conditions. Resonance present in the solution

The same initial steps are repeated as in the first solution, up the equation to solve for x1(t)(8)x¨1+x1=φ(φ2+ξ241)sint+ξ(1ξ2φ24)cost+ξ(ξ23φ24)cos3t+φ(φ23ξ24)sin3t

There exists four particular solutions relating to the above four forcing functionsf1(t)=φ(φ2+ξ241)sintf2(t)=ξ(1ξ2φ24)costf3(t)=ξ(ξ23φ24)cos3tf4(t)=φ(φ23ξ24)sin3t

For each of the above, we guess a particular solution, and determine the solution parameters. We guess x1,p1=t(c1sint+c2cost),  x1,p2=t(c1sint+c2cost), x1,p3=c1sin3t+c2cos3t, and x1,p4=c1sin3t+c2cos3t. By substituting each of the above into (8) one at a time and comparing coefficients, we arrive at the following solutions x1,p1=t(12φ(1φ2+ξ24)cost)x1,p2=t(12ξ(1ξ2φ24)sint)x1,p3=18ξ(ξ2φ241)cos3tx1,p4=18φ(3ξ2φ24)sin3t

x1(t)=x1,h+x1,p4+x2,p4+x3,p4+x4,p4=(c1cost+c2sint)+t(12φ(1φ2+ξ24)cost)+t(12ξ(1ξ2φ24)sint)+18ξ(ξ2φ241)cos3t+18φ(3ξ2φ24)sin3t

Now applying initial conditions which are x1(0)=0 and x˙1(0)=0 and determining c1 and c2 results inx1(t)=(132ξ3+132ξφ2+18ξ)cost132φ(5ξ27φ2+16)sint+t(12φ(1φ2+ξ24)cost)+t(12ξ(1ξ2φ24)sint)+18ξ(ξ2φ241)cos3t+18φ(3ξ2φ24)sin3t

We can now write the final solution for x(t) asx(t)=x0(t)+αx1(t)=φcost+ξsint+α(132ξ3+132ξφ2+18ξ)costα132φ(5ξ27φ2+16)sint+αt(12φ(1φ2+ξ24)cost)+αt(12ξ(1ξ2φ24)sint)+α18ξ(ξ2φ241)cos3t+α18φ(3ξ2φ24)sin3t

To illustrate this solution, we show the phase plot which results from the above solution, and also plot the solution x(t). We select x(0)=1.5 and x˙(0)=1.3229. We note that the same initial conditions are used as in the earlier solution to compare both solutions.

We see clearly the effect of including resonance particular solutions. The limit cycle boundary is increasing with time.

The effect becomes more clear if we plot the solution x(t) itself, and compare it to the solution earlier with no resonance. Below, We show x(t) from both solutions. We clearly see the effect of including the resonance terms. This is the solution with no resonance terms

This is the solution with resonance terms

In the second solution (with resonance), there is no restriction on initial conditions. Hence we can for example look at the solution starting from any x(0) and x˙(0). This is the phase plot for x(0)=5 and x˙(0)=2. This would not be possible using the first solution, due to the restriction on x(0)2 + x˙(0)2 having to be equal to 4

3 Conclusion

Van Der Pol was solved using perturbation for first order approximation. It was solved using two methods. In the first method, the solution was restricted to not include forcing functions which leads to resonance.

In the second method, no such restriction was made. In both methods, this problem was solved for general initial conditions (i.e. both x(0) and x˙(0) can both be nonzero.

In the first method, it was found that the we must restrict the initial conditions to be such that x(0)2 +x˙(0)2=4. This is the condition which resulted in the resonance terms vanishing from the solution. This leads to a solution which generated a limit cycle which did not blow up as the solution is run for longer time. In other words, once the solution enters a limit cycle, it remains in the limit cycle.

In the second method, no restriction on the initial conditions was made. This allowed the solution to start from any state. However, the limit cycle would grow with time, and the solution will suffer from fluttering due to the presence of the resonance terms. However, even though the solution was not stable in the long term, this second approach allowed one to examine the solution for a shorter time but with the flexibility of choosing any initial conditions.

It was also observed that increasing the value of the perturbation parameter α gradually resulted in an inaccurate solution as would be expected, as the solutions derived here all assumed a very small1 α was used in generating the solutions and plots.

For future research, it would be interesting to consider the effect of higher order approximation on the solutions and compare with accurate numerical solutions.

1α=0.01 value was used