How to find the period of periodic solutions of the van der Pol equation?

Question

The equation $$y''+1.115(y^2-1)y'+y=0$$ has solutions that tend towards periodic solutions and I am asked to enter the period of the periodic solutions. How can I find the period without any boundary conditions? And what is the period?

Answer 1

For a numerical solution, initial numerical experiments give a period between $6.7$ and $6.8$ and that the point $(y,y')=(2,0)$ is close to the periodic solution.

Now use Newton's method for the function $f(t)=y'(t)$ to find the zero crossing of the derivative close to $t=6.75$. The value found is an approximation for the period $T$.

Apply this twice to find the sequence $(a_0,a_1,a_2)=\bigl(y(0),y(T),y(2T)\bigr)$ and apply the Aitken's delta-squared iteration formula $$ a=a_0-\frac{(a_1-a_0)^2}{a_2-2a_1+a_0} $$ to get the next better estimate for the initial point $(y,y')=(a,0)$ on the periodic solution.

An implementation of this in Python, aiming at an error level of $10^{-9}$ in $a$, gives the computational steps

  Newton T=6.750000000, y'(T)=-0.028197, y''(T)=-1.9144
  Newton T=6.735271592, y'(T)=0.0007177, y''(T)=-2.0126
  Newton T=6.735628185, y'(T)=4.3368e-07, y''(T)=-2.0102
  Newton T=6.735628401, y'(T)=1.962e-12, y''(T)=-2.0102
Aitken: y(0)=2.000000000, y(T)-y(0)=0.010215, y(2T)-2y(T)+y(0)=-0.010215, a=2.010215386
  Newton T=6.735628401, y'(T)=0.032431, y''(T)=-2.1199
  Newton T=6.750926957, y'(T)=0.00082388, y''(T)=-2.013
  Newton T=6.751336236, y'(T)=5.7135e-07, y''(T)=-2.0102
  Newton T=6.751336521, y'(T)=5.1818e-13, y''(T)=-2.0102
Aitken: y(0)=2.010215386, y(T)-y(0)=3.2853e-06, y(2T)-2y(T)+y(0)=-3.2853e-06, a=2.010218672
  Newton T=6.751336521, y'(T)=1.0205e-05, y''(T)=-2.0103
  Newton T=6.751341597, y'(T)=8.7319e-11, y''(T)=-2.0102
Aitken: y(0)=2.010218672, y(T)-y(0)=1.0227e-09, y(2T)-2y(T)+y(0)=-1.0227e-09, a=2.010218673
  Newton T=6.751341597, y'(T)=3.1768e-09, y''(T)=-2.0102
  Newton T=6.751341599, y'(T)=5.0462e-13, y''(T)=-2.0102
Aitken: y(0)=2.010218673, y(T)-y(0)=2.9843e-13, y(2T)-2y(T)+y(0)=-2.5757e-13, a=2.010218673

---

For large values of the parameter $\mu$ (the $1.115$ in the equation above) this numerical method fails as the zero crossings are phase change points in the fast-slow system, so that there is no good tangent for the Newton method. It is better to use the second coordinate $v=\frac{y'}{\mu}+\frac{y}3(y^2-3)$ so that the resulting first order system adapted to the fast-slow dynamic is \begin{align} y'&=μ\left(v-\frac{y}3(y^2-3)\right)\\ v'&=-\frac{y}{μ} \end{align}

Answer 2

Answer to the first question "How can I find the period without any boundary conditions? " :

You don't need boundary condition since the limit cycle doesn't depend of them. You can chose any initial conditions.

Answer to the second question " What is the period? " (of the limit cycle) : $$y''+\mu (y^2-1)y'+y=0$$ For small values of $\mu$ the equation is approximately $y''+y\simeq 0$ which solution is $y=C\:\sin\left(2\pi\frac{t}{T}+\varphi \right)$ where the period $T=2\pi$.

For a first approximate, you can take $\quad T\simeq 2\pi\quad$ until $\mu=1.115$ is not large. See the empirical graph below.

Semi-empirical formula from M.Cartwright : $$T\simeq \left(3-2\ln(2) \right)\mu+3\frac{2.2338}{\mu^{1/3}} $$ This formula, which isn't convenient for small $\mu$, cannot be used in the present case.

An updated equation is represented on the figure below.

In case of not large $\mu$, expending in series of power of $\mu$ leads to : $$\frac{dt}{dy}\simeq \frac{1}{\sqrt{4-y^2}}+\frac{\mu}{4}y+\frac{\mu^2}{96}(5y^2-2)\sqrt{4-y^2} +...$$ $$T\simeq 2\int_{-2}^2 \left( \frac{1}{\sqrt{4-y^2}}+\frac{\mu}{4}y+\frac{\mu^2}{96}(5y^2-2)\sqrt{4-y^2} \right)dy$$ $$T\simeq 2\pi\left(1+\frac{\mu^2}{16}\right)\qquad\text{not large }\mu .$$ This analytic approximate formula gives $\quad T\simeq 6.77$

Direct numerical solving of the ODE gives $\quad T\simeq 6.75$

More simply, the rough approximate $\quad T\simeq 2\pi\simeq 6.28\quad$ is not too bad in case of $\mu=1.115$.