May 25, 2022

The Van Der Pol diﬀerential equation \[ x''(t) -\alpha \left ( 1-x^2(t) \right ) x'(t) +x(t) = 0 \] Was solved using perturbation with ﬁrst order approximation. Two diﬀerent solutions were obtained. The ﬁrst solution restricted the initial conditions to be \(x(0)^{2}+x^{\prime 2}=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.

The Van Der Pol equation is \[ \ddot {x}-\alpha \left ( 1-x^{2}\right ) \dot {x}+x=0 \] Initial conditions are \(x(0)=\varphi \) and \(\dot {x}(0)=\xi \).

If \(\alpha \lll 1\), then the equation becomes \(\ddot {x}_{0}+x_{0}=0\) which has an exact solution. We perturb the exact solution for \(\ddot {x}_{0}+x_{0}=0\) to approximate the solution of the nonlinear equation \(\ddot {x}-\alpha \left ( 1-x^{2}\right ) \dot {x}+x=0\) in terms of perturbation parameter \(\alpha \) Let the solution of the Van Der Pol equation be \(x(t)\) and the solution to the exact equation \(\ddot {x}_{0}+x_{0}=0\) be \(x_{0}(t)\), then \[ x(t)=x_{0}(t)+\alpha x_{1}(t)+\alpha ^{2}x_{2}(t)+\cdots \] First order approximation results in \[ x(t)=x_{0}(t)+\alpha x_{1}(t) \] Now we need to determine \(x_{0}\left ( t\right ) \) and \(x_{1}(t)\). Substituting (2) into (1) gives \[ \ddot {x}_{0}+\alpha \ddot {x}_{1}-\alpha \dot {x}_{0}+\alpha x_{0}^{2}\dot {x}_{0}+\alpha ^{3}x_{1}^{2}\dot {x}_{0}+2\alpha ^{2}x_{1}x_{0}\dot {x}_{0}-\alpha ^{2}\dot {x}_{1}+\alpha ^{2}x_{0}^{2}\dot {x}_{1}+\alpha ^{4}x_{1}^{2}\dot {x}_{1}+2\alpha ^{3}x_{1}x_{0}\dot {x}_{1}+x_{0}+\alpha x_{1}=0 \] Collecting terms with same power of \(\alpha \) gives \[ \alpha ^{0}\left ( \ddot {x}_{0}+x_{0}\right ) +\alpha \left ( \ddot {x}_{1}-\dot {x}_{0}+x_{0}^{2}\dot {x}_{0}+x_{1}\right ) +\alpha ^{2}\left ( x_{0}^{2}\dot {x}_{1}+2x_{1}x_{0}\dot {x}_{0}-\dot {x}_{1}\right ) +\alpha ^{3}\left ( x_{1}^{2}\dot {x}_{0}+2x_{1}x_{0}\dot {x}_{1}\right ) +\alpha ^{4}x_{1}^{2}\dot {x}_{1}=0 \] Setting terms which multiply by higher power \(\alpha \) to zero and since it is assumed that \(\alpha \) is very small therefore \[ \left ( \ddot {x}_{0}+x_{0}\right ) +\alpha \left ( \ddot {x}_{1}-\dot {x}_{0}+x_{0}^{2}\dot {x}_{0}+x_{1}\right ) =0 \] For the LHS to be zero implies that \[ \ddot {x}_{0}+x_{0}=0 \] And \[ \ddot {x}_{1}-\dot {x}_{0}+x_{0}^{2}\dot {x}_{0}+x_{1}=0 \] Equation (3) is solved for \(x_{0}\) and (4) is solved for \(x_{1}\) in order to ﬁnd the solution \(x(t)\). Solution to (3) is \[ x_{0}(t)=A_{0}\cos t+B_{0}\sin t \] Assuming initial conditions for \(x_{1}(t)\) are zero, then initial conditions for \(x_{0}(t)\) can be taken to be those given for \(x(t)\). Solving for \(A_{0},B_{0}\) gives \[ x_{0}(t)=\varphi \cos t+\xi \sin t \] Substituting the solution just found (and its derivative) into (4) gives \begin {align*} \ddot {x}_{1}+x_{1} & =-\varphi \sin t+\xi \cos t-\xi ^{3}\cos t\sin ^{2}t-2\xi ^{2}\varphi \cos ^{2}t\sin t+\xi ^{2}\varphi \sin ^{3}t\\ & -\xi \varphi ^{2}\cos ^{3}t+\allowbreak 2\xi \varphi ^{2}\cos t\sin ^{2}t+\varphi ^{3}\cos ^{2}t\sin t \end {align*}

Using \(\sin t\cos ^{2}t=\frac {1}{4}\left ( \sin t+\sin 3t\right ) \) and \(\cos t\sin ^{2}t=\frac {1}{4}(\cos t-\cos 3t)\) the above can be simpliﬁed to\begin {align*} \ddot {x}_{1}+x_{1} & =\left ( -\varphi -\frac {\xi ^{2}\varphi }{2}+\frac {\varphi ^{3}}{4}\right ) \sin t+\left ( \xi -\frac {\xi ^{3}}{4}+\frac {\allowbreak \xi \varphi ^{2}}{2}\right ) \cos t+\left ( \frac {\xi ^{3}}{4}-\frac {\allowbreak \xi \varphi ^{2}}{2}\right ) \cos 3t\\ & +\left ( -\frac {\xi ^{2}\varphi }{2}+\frac {\varphi ^{3}}{4}\right ) \sin 3t+\xi \varphi \left ( \xi \sin ^{3}t-\varphi \cos ^{3}t\right ) \end {align*}

Using \(\xi \sin ^{3}t-\varphi \cos ^{3}t=\frac {1}{4}\left ( -3\varphi \cos t-\varphi \cos 3t+3\xi \sin t-\xi \sin 3t\right ) \) the above can be simpliﬁed further to\begin {align} \ddot {x}_{1}+x_{1} & =\left ( -\varphi -\frac {\xi ^{2}\varphi }{2}+\frac {\varphi ^{3}}{4}+\frac {3\xi ^{2}\varphi }{4}\right ) \sin t+\left ( \xi -\frac {\xi ^{3}}{4}+\frac {\allowbreak \xi \varphi ^{2}}{2}-\frac {3}{4}\xi \varphi ^{2}\right ) \cos t\nonumber \\ & +\left ( \frac {\xi ^{3}}{4}-\frac {\allowbreak \xi \varphi ^{2}}{2}-\frac {\xi \varphi ^{2}}{4}\right ) \cos 3t+\left ( -\frac {\xi ^{2}\varphi }{2}+\frac {\varphi ^{3}}{4}-\frac {\xi ^{2}\varphi }{4}\right ) \sin 3t\tag {5} \end {align}

The above is the equation we will now solve for \(x_{1}\left ( t\right ) \), 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 ﬁrst two in the RHS shown above). If we however restrict the initial conditions such that \[ \left ( -\varphi -\frac {\xi ^{2}\varphi }{2}+\frac {\varphi ^{3}}{4}+\frac {3\xi ^{2}\varphi }{4}\right ) =0 \] And \[ \left ( \xi -\frac {\xi ^{3}}{4}+\frac {\allowbreak \xi \varphi ^{2}}{2}-\frac {3}{4}\xi \varphi ^{2}\right ) =0 \] Then those forcing function will vanish. The above two equation result in the same restriction, which can be found to be\[ 4=\xi ^{2}+\varphi ^{2}\] By restricting \(\xi ^{2}+\varphi ^{2}\) to be exactly \(4\), then the solution we obtain for \(x_{1}\left ( t\right ) \) from (5) will not blow up. Hence (4) can now be rewritten without the two forcing functions which caused resonance resulting in\[ \ddot {x}_{1}+x_{1}=\xi \left ( \frac {\xi ^{2}-3\varphi ^{2}}{4}\right ) \cos 3t+\varphi \left ( \frac {\varphi ^{2}-3\xi ^{2}}{4}\right ) \sin 3t \] But \(\varphi ^{2}=4-\xi ^{2}\) and \(\xi ^{2}=4-\varphi ^{2}\), hence the above becomes after some simpliﬁcation\begin {equation} \ddot {x}_{1}+x_{1}=\xi \left ( \xi ^{2}-3\right ) \cos 3t+\varphi \left ( \varphi ^{2}-3\right ) \sin 3t\tag {6} \end {equation} The homogeneous solution to the above is \[ x_{1,h}=c_{1}\cos t+c_{2}\sin t \] And we guess two particular solutions \(x_{1,p1}=a_{1}\cos 3t+a_{2}\sin 3t\) and \(x_{1,p2}=d_{1}\cos 3t+d_{2}\sin 3t\). By substituting each of these particular solutions into (6) one at a time, and comparing coeﬃcients, we can determine \(a_{1},a_{2},d_{1},d_{2}\). This results in \[ x_{1,p1}=\frac {\xi \left ( 3-\xi ^{2}\right ) }{8}\cos 3t \] Similarly, \[ x_{1,p2}=\frac {\varphi \left ( 3-\varphi ^{2}\right ) }{8}\sin 3t \] Hence the solution to (6) becomes\begin {align} x_{1}\left ( t\right ) & =x_{1,h}+x_{1,p1}+x_{1,p2}\nonumber \\ & =\left ( c_{1}\cos t+c_{2}\sin t\right ) +\frac {\xi \left ( 3-\xi ^{2}\right ) }{8}\cos 3t+\frac {\varphi \left ( 3-\varphi ^{2}\right ) }{8}\sin 3t\tag {7} \end {align}

Now applying initial conditions, which are \(x_{1}\left ( 0\right ) =0\) and \(\dot {x}_{1}\left ( 0\right ) =0\), to (7) and determining \(c_{1}\) and \(c_{2}\) results in\[ x_{1}\left ( t\right ) =\frac {\xi \left ( \xi ^{2}-3\right ) }{8}\cos t+\frac {3}{8}\varphi \left ( \varphi ^{2}-3\right ) \sin t+\frac {\xi \left ( 3-\xi ^{2}\right ) }{8}\cos 3t+\frac {\varphi \left ( 3-\varphi ^{2}\right ) }{8}\sin 3t \] Therefore we now have the ﬁnal perturbation solution, which is\begin {align*} x\left ( t\right ) & =x_{0}\left ( t\right ) +\alpha x_{1}\left ( t\right ) \\ & =\varphi \cos t+\xi \sin t+\alpha \left ( \frac {\xi \left ( \xi ^{2}-3\right ) }{8}\cos t+\frac {3}{8}\varphi \left ( \varphi ^{2}-3\right ) \sin t+\frac {\xi \left ( 3-\xi ^{2}\right ) }{8}\cos 3t+\frac {\varphi \left ( 3-\varphi ^{2}\right ) }{8}\sin 3t\right ) \end {align*}

Where \(x\left ( 0\right ) =\varphi \) and \(\dot {x}\left ( 0\right ) =\xi \) and the above solution is valid under the restriction that \[ x^{2}\left ( 0\right ) +\dot {x}^{2}\left ( 0\right ) =4 \] We notice that the above restriction implies that both \(x\left ( 0\right ) \) and \(\dot {x}\left ( 0\right ) \) 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 \(\alpha \). To illustrate this solution, We show the phase plot which results from the above solution, and also plot the solution \(x\left ( t\right ) \). Selecting \(x\left ( 0\right ) =1.5\) and \(\dot {x}\left ( 0\right ) =\sqrt {4-1.5^{2}}=1.3229\)

The same initial steps are repeated as in the ﬁrst solution, up the equation to solve for \(x_{1}\left ( t\right ) \)\begin {align} \ddot {x}_{1}+x_{1} & =\varphi \left ( \frac {\varphi ^{2}+\xi ^{2}}{4}-1\right ) \sin t+\xi \left ( 1-\frac {\xi ^{2}-\varphi ^{2}}{4}\right ) \cos t\tag {8}\\ & +\xi \left ( \frac {\xi ^{2}-3\varphi ^{2}}{4}\right ) \cos 3t+\varphi \left ( \frac {\varphi ^{2}-3\xi ^{2}}{4}\right ) \sin 3t\nonumber \end {align}

There exists four particular solutions relating to the above four forcing functions\begin {align*} f_{1}\left ( t\right ) & =\varphi \left ( \frac {\varphi ^{2}+\xi ^{2}}{4}-1\right ) \sin t\\ f_{2}\left ( t\right ) & =\xi \left ( 1-\frac {\xi ^{2}-\varphi ^{2}}{4}\right ) \cos t\\ f_{3}\left ( t\right ) & =\xi \left ( \frac {\xi ^{2}-3\varphi ^{2}}{4}\right ) \cos 3t\\ f_{4}\left ( t\right ) & =\varphi \left ( \frac {\varphi ^{2}-3\xi ^{2}}{4}\right ) \sin 3t \end {align*}

For each of the above, we guess a particular solution, and determine the solution parameters. We guess \(x_{1,p1}=t\left ( c_{1}\sin t+c_{2}\cos t\right ) \), \(x_{1,p2}=t\left ( c_{1}\sin t+c_{2}\cos t\right ) \), \(x_{1,p3}=c_{1}\sin 3t+c_{2}\cos 3t\), and \(x_{1,p4}=c_{1}\sin 3t+c_{2}\cos 3t\). By substituting each of the above into (8) one at a time and comparing coeﬃcients, we arrive at the following solutions \begin {align*} x_{1,p1} & =t\left ( \frac {1}{2}\varphi \left ( 1-\frac {\varphi ^{2}+\xi ^{2}}{4}\right ) \cos t\right ) \\ x_{1,p2} & =t\left ( \frac {1}{2}\xi \left ( 1-\frac {\xi ^{2}-\varphi ^{2}}{4}\right ) \sin t\right ) \\ x_{1,p3} & =\frac {1}{8}\xi \left ( \frac {\xi ^{2}-\varphi ^{2}}{4}-1\right ) \cos 3t\\ x_{1,p4} & =\frac {1}{8}\varphi \left ( \frac {3\xi ^{2}-\varphi ^{2}}{4}\right ) \sin 3t \end {align*}

\begin {align*} x_{1}\left ( t\right ) & =x_{1,h}+x_{1,p4}+x_{2,p4}+x_{3,p4}+x_{4,p4}\\ & =\left ( c_{1}\cos t+c_{2}\sin t\right ) +t\left ( \frac {1}{2}\varphi \left ( 1-\frac {\varphi ^{2}+\xi ^{2}}{4}\right ) \cos t\right ) +t\left ( \frac {1}{2}\xi \left ( 1-\frac {\xi ^{2}-\varphi ^{2}}{4}\right ) \sin t\right ) \\ & +\frac {1}{8}\xi \left ( \frac {\xi ^{2}-\varphi ^{2}}{4}-1\right ) \cos 3t+\frac {1}{8}\varphi \left ( \frac {3\xi ^{2}-\varphi ^{2}}{4}\right ) \sin 3t \end {align*}

Now applying initial conditions which are \(x_{1}\left ( 0\right ) =0\) and \(\dot {x}_{1}\left ( 0\right ) =0\) and determining \(c_{1}\) and \(c_{2}\) results in\begin {align*} x_{1}\left ( t\right ) & =\left ( -\frac {1}{32}\xi ^{3}+\frac {1}{32}\xi \varphi ^{2}+\frac {1}{8}\xi \right ) \cos t-\frac {1}{32}\varphi \left ( 5\xi ^{2}-7\varphi ^{2}+16\right ) \sin t\\ & +t\left ( \frac {1}{2}\varphi \left ( 1-\frac {\varphi ^{2}+\xi ^{2}}{4}\right ) \cos t\right ) +t\left ( \frac {1}{2}\xi \left ( 1-\frac {\xi ^{2}-\varphi ^{2}}{4}\right ) \sin t\right ) \\ & +\frac {1}{8}\xi \left ( \frac {\xi ^{2}-\varphi ^{2}}{4}-1\right ) \cos 3t+\frac {1}{8}\varphi \left ( \frac {3\xi ^{2}-\varphi ^{2}}{4}\right ) \sin 3t \end {align*}

We can now write the ﬁnal solution for \(x\left ( t\right ) \) as\begin {align*} x\left ( t\right ) & =x_{0}\left ( t\right ) +\alpha x_{1}\left ( t\right ) \\ & =\varphi \cos t+\xi \sin t+\\ & \alpha \left ( -\frac {1}{32}\xi ^{3}+\frac {1}{32}\xi \varphi ^{2}+\frac {1}{8}\xi \right ) \cos t-\alpha \frac {1}{32}\varphi \left ( 5\xi ^{2}-7\varphi ^{2}+16\right ) \sin t\\ & +\alpha t\left ( \frac {1}{2}\varphi \left ( 1-\frac {\varphi ^{2}+\xi ^{2}}{4}\right ) \cos t\right ) +\alpha t\left ( \frac {1}{2}\xi \left ( 1-\frac {\xi ^{2}-\varphi ^{2}}{4}\right ) \sin t\right ) \\ & +\alpha \frac {1}{8}\xi \left ( \frac {\xi ^{2}-\varphi ^{2}}{4}-1\right ) \cos 3t+\alpha \frac {1}{8}\varphi \left ( \frac {3\xi ^{2}-\varphi ^{2}}{4}\right ) \sin 3t \end {align*}

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

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

The eﬀect becomes more clear if we plot the solution \(x\left ( t\right ) \) itself, and compare it to the solution earlier with no resonance. Below, We show \(x(t)\) from both solutions. We clearly see the eﬀect 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\left ( 0\right ) \) and \(\dot {x}\left ( 0\right ) \). This is the phase plot for \(x\left ( 0\right ) =5\) and \(\dot {x}\left ( 0\right ) =2\). This would not be possible using the ﬁrst solution, due to the restriction on \(x\left ( 0\right ) ^{2}\) \(+\) \(\dot {x}\left ( 0\right ) ^{2}\) having to be equal to \(4\)

Van Der Pol was solved using perturbation for ﬁrst order approximation. It was solved using two methods. In the ﬁrst 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\left ( 0\right ) \) and \(\dot {x}\left ( 0\right ) \) can both be nonzero.

In the ﬁrst method, it was found that the we must restrict the initial conditions to be such that \(x\left ( 0\right ) ^{2}\) \(+\dot {x}\left ( 0\right ) ^{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 suﬀer from ﬂuttering 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 ﬂexibility of choosing any initial conditions.

It was also observed that increasing the value of the perturbation parameter \(\alpha \) gradually resulted in an
inaccurate solution as would be expected, as the solutions derived here all assumed a very
small^{1} \(\alpha \)
was used in generating the solutions and plots.

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

^{1}\(\alpha =0.01\) value was used