Poincaré–Lindstedt method

Example: the van der Pol oscillator

We solve the van der Pol oscillator only up to order 2. This method can be continued indefinitely in the same way, where the order-n term ${\displaystyle {ϵ}^n{x}_n}$ consists of a harmonic term ${\displaystyle a_n\cos (t)+b_n\cos (t)}$, plus some super-harmonic terms ${\displaystyle a_{n,2}\cos (2t)+b_{n,2}\cos (2t)+\cdots }$. The coefficients of the super-harmonic terms are solved directly, and the coefficients of the harmonic term are determined by expanding down to order-(n+1), and eliminating its secular term.

See chapter 10 of for a derivation up to order 3, and for a computer derivation up to order 164.

Consider the van der Pol oscillator with equation$${\displaystyle \ddot{x}+ϵ(x^2-1)\dot{x}+x=0}$$where ${\displaystyle ϵ}$ is a small positive number. Perform substitution to the second order:

${\displaystyle \tau =\omega t,}$ where ${\displaystyle \omega =1+ϵ{\omega }_1+{ϵ}^2{\omega }_2+O({ϵ}^3)}$

which yields the equation$${\displaystyle {\omega }^2\ddot{x}+\omega ϵ(x^2-1)\dot{x}+x=0}$$Now plug in ${\displaystyle x=x_0+ϵx_1+{ϵ}^2{x}_2+O({ϵ}^3)}$, and we have three equations, for the orders ${\displaystyle 1,ϵ,{ϵ}^2}$ respectively:$${\displaystyle \{\begin{array}{l}{\ddot{x}}_0+x_0=0\\ {\ddot{x}}_1+x_1+2{\omega }_1{\ddot{x}}_0+(x_0^2-1){\dot{x}}_0=0\\ {\ddot{x}}_2+x_2+({\omega }_1^2+2{\omega }_2){\ddot{x}}_0+2{\omega }_1{\ddot{x}}_1+2x_0{x}_1{\dot{x}}_0+{\omega }_1(x_0^2-1){\dot{x}}_0+{\dot{x}}_1(x_0^2-1)=0\end{array}}$$The first equation has general solution ${\displaystyle x_0=A\cos (\tau +\varphi )}$. Pick origin of time such that ${\displaystyle \varphi =0}$. Then plug it into the second equation to obtain (after some trigonometric identities)$${\displaystyle {\ddot{x}}_1+x_1+(A-A^3/4)\sin \tau -2{\omega }_{1}A\cos \tau -(A^3/4)\sin (3\tau )=0}$$To eliminate the secular term, we must set both ${\displaystyle \sin \tau ,\cos \tau }$ coefficients to zero, thus we have $${\displaystyle \{\begin{array}{l}A=A^3/4\\ 2{\omega }_{1}A=0\end{array}}$$yielding ${\displaystyle A=2,{\omega }_1=0}$. In particular, we found that when ${\displaystyle ϵ}$ increases from zero to a small positive constant, all circular orbits in phase space are destroyed, except the one at radius 2. Now solving ${\displaystyle {\ddot{x}}_1+x_1=2\sin (3\tau )}$ yields ${\displaystyle x_1=B\cos (\tau +\varphi )-\frac{1}{4}\sin (3\tau )}$. We can always absorb ${\displaystyle ϵB\cos (\tau +\varphi )}$ term into ${\displaystyle x_0}$, so we can WLOG have just ${\displaystyle x_1=-\frac{1}{4}\sin (3\tau )}$.

Now plug into the second equation to obtain$${\displaystyle {\ddot{x}}_2+x_2-(4{\omega }_2+1/4)\cos \tau -\frac{3}{4}\cos 3\tau -\frac{5}{4}\cos 5\tau =0}$$To eliminate the secular term, we set ${\displaystyle {\omega }_2=-\frac{1}{16}}$.

Thus we find that ${\displaystyle \omega =1-\frac{1}{16}{ϵ}^2+O({ϵ}^3)}$.

Example: Mathieu equation

This is an example of parametric resonance.

Consider the Mathieu equation ${\displaystyle \ddot{x}+(1+b{ϵ}^2+ϵ\cos (t))x=0}$, where ${\displaystyle b}$ is a constant, and ${\displaystyle ϵ}$ is small. The equation's solution would have two time-scales, one fast-varying on the order of ${\displaystyle t}$, and another slow-varying on the order of ${\displaystyle T={ϵ}^{2}t}$. So expand the solution as $${\displaystyle x(t)=x_0(t,T)+ϵx_1(t,T)+{ϵ}^2{x}_2(t,T)+O({ϵ}^3)}$$Now plug into the Mathieu equation and expand to obtain$${\displaystyle \{\begin{array}{l}{\mathrm{\partial }}_t^2{x}_0+x_0=0\\ {\mathrm{\partial }}_t^2{x}_1+x_1=-\cos (t)x_0\\ {\mathrm{\partial }}_t^2{x}_2+x_2=-b{x}_0-2{\mathrm{\partial }}_{tT}{x}_0-\cos (t)x_1\end{array}}$$As before, we have the solutions$${\displaystyle \{\begin{array}{l}{x}_0=A\cos (t)+B\sin (t)\\ x_1=-\frac{A}{2}+\frac{A}{6}\cos (2t)+\frac{B}{6}\sin (2t)\end{array}}$$The secular term coefficients in the third equation are $${\displaystyle \{\begin{array}{l}\frac{1}{12}\left(-12bA+5A-24B^{\prime }\right)\\ \frac{1}{12}\left(24A^{\prime }-12bB-B\right)\end{array}}$$Setting them to zero, we find the equations of motion:

${\displaystyle \frac{d}{dT}\left[\begin{array}{c}A\\ B\end{array}\right]=\left[\begin{array}{cc}0& \frac{1}{2}(\frac{1}{12}+b)\\ \frac{1}{2}(\frac{5}{12}-b)& 0\end{array}\right]\left[\begin{array}{c}A\\ B\end{array}\right]}$

Its determinant is ${\displaystyle \frac{1}{4}(b-5/12)(b+1/12)}$, and so when ${\displaystyle b\in (-1/12,5/12)}$, the origin is a saddle point, so the amplitude of oscillation ${\displaystyle \sqrt{A^2+B^2}}$ grows unboundedly.

In other words, when the angular frequency (in this case, ${\displaystyle 1}$) in the parameter is sufficiently close to the angular frequency (in this case, ${\displaystyle \sqrt{1+b{ϵ}^2}}$) of the original oscillator, the oscillation grows unboundedly, like a child swinging on a swing pumping all the way to the moon.

Shohat expansion

For the van der Pol oscillator, we have ${\displaystyle \omega \sim 1/ϵ}$ for large ${\displaystyle ϵ}$, so as ${\displaystyle ϵ}$ becomes large, the serial expansion of ${\displaystyle \omega }$ in terms of ${\displaystyle ϵ}$ diverges and we would need to keep more and more terms of it to keep ${\displaystyle \omega }$ bounded. This suggests to us a parametrization that is bounded:$${\displaystyle r:=\frac{ϵ}{1+ϵ}}$$Then, using serial expansions ${\displaystyle ϵ\omega =r+c_2{r}^2+c_3{r}^3+c_4{r}^4+\cdots }$ and ${\displaystyle x=x_0+r{x}_1+r^2{x}_2+\cdots }$, and using the same method of eliminating the secular terms, we find ${\displaystyle c_2=1,c_3=\frac{15}{16},c_4=\frac{13}{16}}$.

Because ${\displaystyle \underset{ϵ\to \mathrm{\infty }}{lim}r=1}$, the expansion ${\displaystyle ϵ\omega =r+c_2{r}^2+c_3{r}^3+c_4{r}^4+\cdots }$ allows us to take a finite number of terms for the series on the right, and it would converge to a finite value at ${\displaystyle ϵ\to \mathrm{\infty }}$ limit. Then we would have ${\displaystyle \omega \sim 1/ϵ}$, which is exactly the desired asymptotic behavior. This is the idea behind Shohat expansion.

The exact asymptotic constant is ${\displaystyle ϵ\omega \to \frac{2\pi }{3-2\ln 2}=3.8936\cdots }$, which as we can see is approached by ${\displaystyle 1+c_2+c_3+c_4=3.75}$.