Poincaré–Lindstedt method explained

In perturbation theory, the Poincaré–Lindstedt method or Lindstedt–Poincaré method is a technique for uniformly approximating periodic solutions to ordinary differential equations, when regular perturbation approaches fail. The method removes secular terms—terms growing without bound—arising in the straightforward application of perturbation theory to weakly nonlinear problems with finite oscillatory solutions.^{[1] [2]}

The method is named after Henri Poincaré,^[3] and Anders Lindstedt.^[4]

The article gives several examples. The theory can be found in Chapter 10 of Nonlinear Differential Equations and Dynamical Systems by Verhulst.^[5]

Example: the Duffing equation

The undamped, unforced Duffing equation is given by

\ddot{x} + x + \varepsilon x ^{3} = 0

for t > 0, with 0 < ε ≪ 1.^[6] Consider initial conditions

x( 0 ) = 1 , • x( 0 ) = 0 .

A perturbation-series solution of the form x(t) = x₀(t) + ε x₁(t) + ... is sought. The first two terms of the series are

x(t) = \cos(t) + \varepsilon \left[ \tfrac{ 1 }{ 32 } \left( \cos( 3 t) - \cos(t) \right) - \tfrac{ 3 }{ 8 } t \sin(t) \right] + … .

This approximation grows without bound in time, which is inconsistent with the physical system that the equation models.^[7] The term responsible for this unbounded growth, called the secular term, is

t\sin(t) . The Poincaré–Lindstedt method allows for the creation of an approximation that is accurate for all time, as follows.

In addition to expressing the solution itself as an asymptotic series, form another series with which to scale time t:

\tau = \omega t, where \omega = \omega _{0} + \varepsilon \omega _{1} + … .

We have the leading order ω₀ = 1, because when

\epsilon = 0 , the equation has solution x = \cos (t) . Then the original problem becomes\omega ^{2} x''(\tau) + x(\tau) + \varepsilon x ^{3(\tau)} = 0

Now search for a solution of the form x(τ) = x₀(τ) + ε x₁(τ) + ... . The following solutions for the zeroth and first order problem in ε are obtained:

\begin{align} x _{0} & = \cos(\tau) \\ and x _{1} & = \tfrac{ 1 }{ 32 } \left(\cos( 3 \tau) - \cos(\tau)\right) + \left( \omega _{1} - \tfrac{ 3 }{ 8 } \right) \tau \sin(\tau). \end{align}

So the secular term can be removed through the choice: ω₁ = . Higher orders of accuracy can be obtained by continuing the perturbation analysis along this way. As of now, the approximation—correct up to first order in ε—is

x(t) ≈ \cos l (\left( 1 + \tfrac{ 3 }{ 8 } \varepsilon \right) t r ) + \tfrac{ 1 }{ 32 } \varepsilon \left[\cos l ( 3 \left( 1 + \tfrac{ 3 }{ 8 }\varepsilon \right) t r ) - \cos l (\left( 1 + \tfrac{ 3 }{ 8 }\varepsilon \right) t r )\right].

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

\epsilon ^{n} x _{n} consists of a harmonic term a _{n\cos(t)} + b _{n\cos(t)} , plus some super-harmonic terms a _{n,} \cos( 2 t) + b _{n,} \cos( 2 t) + … . 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 ^[8] for a computer derivation up to order 164.

Consider the van der Pol oscillator with equation$$\backslash ddot\; x\; +\; \backslash epsilon\; (x^2-1)\; \backslash dot\; x\; +\; x\; =\; 0$$where

\epsilon is a small positive number. Perform substitution to the second order:

\tau = \omega t, where \omega = 1 + \epsilon \omega _{1} + \epsilon ^{2} \omega _{2} + O(\epsilon ^{3)}

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

Now plug into the second equation to obtain$$\backslash ddot\; x\_2\; +\; x\_2\; -(4\backslash omega\_2\; +\; 1/4)\backslash cos\backslash tau\; -\; \backslash frac\; 34\; \backslash cos\; 3\backslash tau\; -\; \backslash frac\; 54\; \backslash cos\; 5\backslash tau\; =\; 0$$To eliminate the secular term, we set

\omega _{2} =- (1)/(16) .

Thus we find that

\omega = 1 - (1)/(16) \epsilon ^{2} + O(\epsilon ^{3)} .

Example: Mathieu equation

This is an example of parametric resonance.

\ddot x + ( 1 + b\epsilon ^{2} + \epsilon \cos(t))x = 0 , where b is a constant, and \epsilon is small. The equation's solution would have two time-scales, one fast-varying on the order of t , and another slow-varying on the order of T = \epsilon ^{2} t . So expand the solution as $$x(t)\; =\; x\_0(t,\; T)\; +\; \backslash epsilon\; x\_1(t,\; T)\; +\; \backslash epsilon^2\; x\_2(t,\; T)\; +\; O(\backslash epsilon^3)$$Now plug into the Mathieu equation and expand to obtain$$\backslash begin\backslash partial\_t^2\; x\_0\; +\; x\_0\; =\; 0\; \backslash \backslash \backslash partial\_t^2\; x\_1\; +\; x\_1\; =\; -\backslash cos(t)x\_0\; \backslash \backslash \backslash partial\_t^2\; x\_2\; +\; x\_2\; =\; -bx\_0\; -\; 2\backslash partial\_x\_0\; -\; \backslash cos(t)x\_1\backslash end$$As before, we have the solutions$$\backslash beginx\_0\; =\; A\backslash cos(t)\; +\; B\; \backslash sin(t)\; \backslash \backslash x\_1\; =\; -\backslash frac\; A2\; +\; \backslash frac\; A6\; \backslash cos(2t)\; +\; \backslash frac\; B6\; \backslash sin(2t)\backslash end$$The secular term coefficients in the third equation are $$\backslash begin\backslash frac\; \backslash left(-12\; b\; A+5\; A-24\; B\text{'}\backslash right)\; \backslash \backslash \backslash frac\; \backslash left(24\; A\text{'}-12\; b\; B-B\backslash right)\backslash end$$Setting them to zero, we find the equations of motion:(d)/(dT) \begin{bmatrix} A \ B \end{bmatrix} = \begin{bmatrix} 0 & (12)/(( (1)/(12)) + b) \\ (12)/(( (5)/(12)) - b) & 0 \\ \end{bmatrix} \begin{bmatrix} A\\B \end{bmatrix}

Its determinant is

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

In other words, when the angular frequency (in this case,

1 ) in the parameter is sufficiently close to the angular frequency (in this case, \sqrt{ 1+ b\epsilon ^{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

^[9]

For the van der Pol oscillator, we have

\omega \sim 1 /\epsilon for large \epsilon , so as \epsilon becomes large, the serial expansion of \omega in terms of \epsilon diverges and we would need to keep more and more terms of it to keep \omega bounded. This suggests to us a parametrization that is bounded:$$r\; :=\; \backslash frac$$Then, using serial expansions \epsilon \omega = r + c _{2} r ^{2} + c _{3} r ^{3} + c _{4} r ^{4} + … and x = x _{0} + r x _{1} + r ^{2} x _{2} + … , and using the same method of eliminating the secular terms, we find c _{2} = 1 , c _{3} = (15)/(16) , c _{4} = (13)/(16) .

Because

\lim _{\epsilon} r = 1 , the expansion \epsilon \omega = r + c _{2} r ^{2} + c _{3} r ^{3} + c _{4} r ^{4} + … allows us to take a finite number of terms for the series on the right, and it would converge to a finite value at \epsilon \to inf ty limit. Then we would have \omega \sim 1 /\epsilon , which is exactly the desired asymptotic behavior. This is the idea behind Shohat expansion.

The exact asymptotic constant is

\epsilon\omega \to (2 \pi)/(3-2ln2) = 3 . 8936 … , which as we can see is approached by 1 + c _{2} + c _{3} + c _{4} = 3 . 75 .

References and notes

1. , pp. 181–186.
2. Book: Strogatz, Steven . Nonlinear dynamics and chaos : with applications to physics, biology, chemistry, and engineering . 2019 . 978-0-367-09206-1 . 2nd . Boca Raton . Exercise 7.6.19, 7.6.21 . 1112373147.
3. , §123–§128.
4. A. Lindstedt, Abh. K. Akad. Wiss. St. Petersburg 31, No. 4 (1882)
5. Book: Verhulst, Ferdinand . Nonlinear Differential Equations and Dynamical Systems . 1996 . Springer Berlin Heidelberg . 978-3-540-60934-6 . Universitext . Berlin, Heidelberg . 10.1007/978-3-642-61453-8.
6. J. David Logan. Applied Mathematics, Second Edition, John Wiley & Sons, 1997. .
7. The Duffing equation has an invariant energy \scriptstyle E = \tfrac 12 • x ^{2} + \tfrac 12 x ^{2} + \tfrac 14 \varepsilon x ^{4}  = constant, as can be seen by multiplying the Duffing equation with \scriptstyle • x and integrating with respect to time t. For the example considered, from its initial conditions, is found: E =  +  ε.
8. Andersen . C. M. . Geer . James F. . June 1982 . Power Series Expansions for the Frequency and Period of the Limit Cycle of the Van Der Pol Equation . SIAM Journal on Applied Mathematics . en . 42 . 3 . 678–693 . 10.1137/0142047 . 0036-1399.
9. Book: Bellman, Richard . Perturbation techniques in mathematics, engineering & physics . 2003 . Dover Publications . 0-486-43258-0 . Mineola, N.Y. . 2.5. The Shohat Expansion . 51942387.