Example 1: This problem is discussed in the book of Lin & Segel, pp. 233. A projectile is sent vertically up from a planet without atmosphere. The motion is described by the position x^✶(t^✶) above the planet’s surface, where t✶ is the time and the initial conditions are \[ x^{\ast}(0) = 0, \quad \frac{{\text d}x^{\ast}}{{\text d}t^{\ast}}\left( 0 \right) = v . \tag{1.1} \] The projectile will be affected by a force given by Newton’s law of gravitation, that is, \[ F \left( x^{\ast} \right) = -G\,\frac{M\,m}{(R + x^{\ast})^2} , \] where G = 6.6743 × 10⁻¹¹ m³/(kg s²) is the gravitational constant, M is the planet’s mass, R the planet’s radius, and m the projectile’s mass. Similar to Earth, the gravity force on the planet’s surface may be written as F(0) = −mg, so that g = GM/R². Thus, it follows that \[ m\,\frac{{\text d}^2 x^{\ast}}{{\text d} t^{{\ast}2}} = - \frac{R^2 gm}{(R + x^{\ast})^2} . \tag{1.2} \] The mathematical model thus consists of the non-linear differential equation (1.2) subject to the initial conditions stated in (1.1).

We are going to study a situation where the initial velocity v is much smaller than the planet’s escape velocity. If v is larger than the escape velocity, the projective will leave the planet permanently. For Earth the escape velocity is about 11.2km/s. However, here the assumption implies that x^✶(t^✶) ⪡ R for the whole trip of the projectile. Under this assumption, the equation simplifies to \[ \frac{{\text d}^2 x^{\ast}}{{\text d} t^{{\ast}2}} = - \frac{R^2 g}{(R + x^{\ast})^2} = - \frac{g}{(1 + x^{\ast}/R)^2} \approx -g . \tag{1.3} \] This equation can be solved with the given initial conditions: \[ x^{\ast} (t^{\ast}) \approx - \frac{1}{2}\, gt^{{\ast}2} + vt^{\ast} . \]

The approximate maximum height follows from the formula above by observing that the time to maximum height is approximately given by \[ \frac{{\text d} x^{\ast}}{{\text d} t^{\ast}} \approx -g t^{\ast} + v = 0, \] so \[ t_{\max} \approx \frac{v}{g} . \] Thus, \[ x_{\max} \approx - \frac{1}{2} \,g \left( \frac{v}{g} \right)^2 + v \left( \frac{v}{g} \right) = \frac{1}{2}\,\frac{v^2}{g} . \] Reasonable scaling of the equation of motion is \[ X = \frac{v^2}{g} , \quad T = \frac{v}{g} . \] This leads to the differential equation \[ \frac{{\text d}^2 \left( \frac{v^2}{g}\,x \right)}{{\text d} \left( \frac{v}{g}\,t \right)^2} = - \frac{R^2 g}{\left( R + \frac{v^2}{g}\,x \right)^2} , \] with the initial conditions \[ \frac{v^2}{g}\, x(0) = 0 , \quad \frac{{\text d} \left( \frac{v^2}{g}\,x \right)}{{\text d} \left( \frac{v}{g}\,t \right)} \left( 0 \right) = v . \] Upon cleaning the expression, we get the scaled problem: \[ \ddot{x} = - \frac{1}{(1 + \varepsilon x)^2} , \tag{1.4} \] subject to \[ x(0) = 0, \quad \dot{x}(0) = =1 , \quad \varepsilon = \frac{v^2}{R} . \tag{1.5} \] It turns out that it is not possible to express the solution of this initial value problem, x = x(t, ε), by means of elementary functions (this is not quite obvious!). Note that since the parameter ε is approximately equal to 2x_max/R, it is indeed small under the assumption we made above.

Solution by the perturbation method:    We are looking for solution of the IVP (1.4), (1.5) in the form of infinite series \[ x (t) = x_0 (t) + x_1 (t)\,\varepsilon + x_2 (t) \,\varepsilon^2 + \cdots . \] Substitution of this series into Eq.(1.4) yields

\begin{align*} \ddot{x} &= \ddot{x}_0 + \ddot{x}_1 \varepsilon + \ddot{x}_2 \varepsilon^2 + \cdots = - \left( 1 + \varepsilon x \right)^{-2} \\ &= -\left[ 1 + (-2) \varepsilon x + \frac{(-2)(-3)}{2} \left( \varepsilon x \right)^2 + \cdots \right] \\ &= -1 + 2 \varepsilon \left( x_0 + \varepsilon x_1 + \cdots \right) - 3 \varepsilon^2 x_0^2 + \cdots \\ &= -1 + \varepsilon 2 x_0 + \varepsilon^2 \left( 2 x_1 - 3 x_0^2 \right) + \cdots . \end{align*} By collecting the coefficients in front of each power of ε, we find the system \begin{align*} \ddot{x}_0 &= -1 , \\ \ddot{x}_1 &= 2 x_0 . \\ \ddot{x}_2 &= 2 x_1 - 3 x_0^2 , \\ \ddot{x}_3 &= 2 x_2 + 2 x_0 x_1 - 2 x_0 \left( 2 x_1 + x_0^2 \right) - 2 \left( 2 x_1 - 3 x_0^2 \right) x_0 , \end{align*} and so on. We solve these differential equations under the following initial conditions: \begin{align*} x_0 &= 0, \quad \dot{x}_0 = 1, \\ x_1 &= 0, \quad \dot{x}_1 = 0, \\ x_2 &= 0, \quad \dot{x}_2 = 0. \end{align*} Thus, x₀ takes care of the initial conditions (1.5), which are consequently satisfied no matter where we stop the series expansion. The solution for x₀ follows immediately from the equation \[ x_0 (t) = t - \frac{1}{2}\, t^2 . \] Upon introducing this expression into the next equation, we find \[ \ddot{x}_1 = 2 x_0 = 2t - t^2 , \quad x_1 (0) = \dot{x}_1 (0) = 0 . \] Its solution is \[ x_1 (t) = \frac{1}{3} \, t^3 - \frac{1}{12}\, t^4 . \] Note that only the particular solutions change with every step because the initial conditions for all x_k(t) are homogeneous. Usually, the algebra quickly becomes quite complicated, but today we can make good use of a computer algebra system such as Mathematica for symbolic manipulation. This yields \[ x(t) = t - \frac{t^2}{2} + \varepsilon \left( \frac{t^3}{3} - \frac{1}{12}\,t^4 \right) + \varepsilon^2 \left( - \frac{t^4}{4} + \frac{11}{60}\,t^5 - \frac{11}{360} ,t^6\right) + O\left( \varepsilon^3 \right) . \]

Analytic solution:    Despit that it is not possible to express the full solution of the IVP (1.4), (1.5) via elementary functions, its solution can be expressed integral (quadrature) form.

Upon multiplying equation (1.4) with ẋ, we find \[ \frac{\text d}{{\text d}t} \left( \frac{\dot{x}^2}{2} \right) = \frac{\text d}{{\text d}t} \left( \frac{1}{\varepsilon} \cdot \frac{1}{1 + \varepsilon x} \right) , \] so upon integration we obtain an ordinary differential equation of the first order when independent variable is missing: \[ \frac{\dot{x}^2}{2} - \frac{1}{\varepsilon} \cdot \frac{1}{1 + \varepsilon x} = \mbox{constant}. \]    ■

Example 1A: We consider the same problem discussed in the previous example, but now we solve it using the Adomian decomposition method (ADM for short). \[ \ddot{x} = - \frac{1}{(1 + \varepsilon x)^2} , \quad x(0) = 0, \quad \dot{x}(0) = =1 . \] The ADM is looking for a solution in the form of infinite series \[ x(t) = x_0 (t) + x_1 (t) + x_2 (t) + \cdots , \] where \[ x_0 (t) = t - \frac{1}{2}\, t^2 . \] is the solution of the IVP: \[ \ddot{x}_0 = -1, \qquad x(0) = 0, \quad \dot{x}_0 (0) = 1 . \]    ■

Example 2: In the previous section, we introduced simple pendulum model that consists of weightless rigid rod of fixed length ℓ with a point of mass m at one end; the mass is free to rotate about a vertical axis at the point of suspension. Assuming that pendulum oscillates with a fixed plane, the angular displacement θ^✶ of the bob is a solution of the following initial value problem: \[ \begin{split} \frac{{\text d}^2 \theta^{\ast}}{{\text d} t^{{\ast}2}} + \omega_0^2 \sin\theta^{\ast} = 0 \qquad \left( \omega_0^2 = \frac{g}{\ell} \right) , \\ \theta^{\ast}(0) = a. \qquad \left. \frac{{\text d} \theta^{\ast}}{{\text d}t^{\ast}} \right\vert_{t^{\ast} = 0} = \Omega . \end{split} \] Here θ^✶ = 0 when the pendulum hangs vertically downward. With "usual" assumption of small displacement, one replaces sinθ^✶ by θ^✶ and obtains \[ \theta^{\ast} (t^{\ast}) = a\,\cos\left( \omega_0 t^{\ast} \right) + \frac{\Omega}{\omega_0}\,\sin \left( \omega_0 t^{\ast} \right). \] We would like to improve this approximation in a systematic manner. For simplicity, we restrict our discussion to the case Ω = 0 (pendulum released from the rest).

The first step is to introduce dimensionless variables Θ and t: \[ t = \omega_0 t^{\ast} , \qquad \Theta = \frac{\theta^{\ast}}{a} . \] The problem now takes the form: \[ \frac{{\text d}^2 \Theta}{{\text d}t^2} + \frac{1}{a}\,\sin \left( a\Theta \right) = 0, \qquad \Theta (0) =1, \quad \left. \frac{{\text d}\Theta}{{\text d}t} \right\vert_{t=0} = 0 . \] The new dependent variable Θ measures angular displacement in units of the initial angular displacement 𝑎. The new independent variable t measures time in units of 1/ω₀, where 2π/ω₀ is the period of small oscillations. The introduction of these dimensionless variables has provided us with a problem in which there appears a single parameter 𝑎, and we want to develop an approximate solution under the assumption that 𝑎 is small compared to unity. With these scaling variables, the reference angle 𝑎 is the largest value of θ^✶(t^✶). The reference time is 1/ω₀.

A basic tool in our analysis is the https://mathworld.wolfram.com/MaclaurinSeries.html">Maclaurin expansion \[ F(a) = \sum_{i\ge 0} C_k a^k , \qquad C_k = \frac{1}{k!}\,\frac{{\text d}^k F(0)}{{\text d}a^k} . \] Assuming convergence, the successive partial sums of the series above provide ever more accurate approximations to F for small 𝑎. However, this approximation requires the function to be infinitely differentiable at the origin. On the other hand, if only a few derivatives are known to exist, we can use the Maclauren formula in the form of a finite sum plus the reminder: \[ F(a) = \sum_{i=0}^n C_i a^i + \left[ \frac{{\text d}^{n+1} F}{{\text d} a^{n+1}} \left( \xi \right) \right] \frac{a^{n+1}}{(n+1)!} , \qquad 0 \le \xi \le a . \] No convergence questions arise here and bounds on the (n + 1)st derivative permit estimates of the error.

It is assumed in future exposition that various functions have infinite series of Maclaurin type. If these functions have only a finite number of derivatives, however, there is a tacit understanding that that the series is appropriately terminated. Also, if a function F depends on variables other than 𝑎, partial derivatives must be used in the Maclaurin expansion, and the coefficients are no longer constant. Thus, if F = F(t, 𝑎), then \[ F(t, a) = \sum_{i\ge 0} C_i (t)\,a^i , \qquad C_i (t) = \frac{1}{i!}\,\frac{\partial^i F}{\partial a^i} . \] It is usually only feasible to obtain the first few terms in the various series. To make sure this is done systematically, we must keep track of the terms we have neglected. To do this, we use the big Oh notation, where by O(𝑎^k) we mean a collection of terms containing k-th or higher powers of 𝑎. For example, \[ \sin a = a - \frac{a^3}{3!} + O\left( a^5 \right) , \qquad \cos a = 1 + O\left( a^2 \right) . \]    ■

Example 2A: Assume that the dependent variable can be expanded in a series of powers of the small parameter: \[ \Theta (t, a) = \sum_{i\ge 0} \Theta_i (t)\,a^i . \tag{2A.1} \] Substituting this series into the pendulum equation, we get \[ \sum_{i\ge 0} \ddot{\Theta}_i a^i + a^{-1} \sin \left( a \sum_{i\ge 0} \Theta_i a^i \right) = 0 \qquad \left( \cdot \equiv \frac{\text d}{{\text d}t} \right) , \] where validity of term-by-term differentiation has been taken for granted.

Next step is to use the Faà di Bruno's formula to expand sine function into power series. Faà di Bruno's formula is an identity in mathematics generalizing the chain rule to higher derivatives. It is named after the Italian mathematician Francesco Faà di Bruno (1855, 1857), although he was not the first to state or prove the formula. In 1800, more than 50 years before Faà di Bruno, the French mathematician Louis François Antoine Arbogast had stated the formula in a calculus textbook, which is considered to be the first published reference on the subject. This formula is not for practical usage because the number of terms grows as a snow ball (exponentially). \begin{align*} \frac{1}{a} \sin \left( \sum_{i\ge 0} \Theta_i a^i \right) &= \\ &= \frac{1}{a} \left( a \sum_{i\ge 0} \Theta_i a^i \right) - \frac{a^{-1}}{3!} \left( a \sum_{i\ge 0} \Theta_i a^i \right)^3 + \frac{a^{-1}}{5!} \left( a \sum_{i\ge 0} \Theta_i a^i \right)^5 - \cdots . \end{align*} We aim to retain terms only through O(𝑎²), but during intermediate calculations we regain some higher terms to make sure that nothing crucial has been neglected. So we find \begin{align*} \frac{1}{a} \sin \left( \sum_{i\ge 0} \Theta_i a^i \right) &= \Theta_0 + a \Theta_1 + a^2 \Theta_2 + O\left( a^3 \right) - \frac{a^2}{6} \left[ \Theta_0^3 + O(a) \right] + O\left( a^4 \right) \\ &= \Theta_0 + a \Theta_1 + a^2 \left( \Theta_2 - \frac{\Theta_0^3}{6} \right) + O \left( a^3 \right) . \end{align*}

Our next step is aimed to collect all terms in the equation and equate them to zero the successive coefficients in the series. Using the latter equation, we write \begin{align*} \ddot{\Theta}_0 + \Theta_0 &= 0 , \\ \ddot{\Theta}_1 + \Theta_1 &= 0 , \\ \ddot{\Theta}_2 + \Theta_2 &= \frac{1}{6}\,\Theta_0^3 . \end{align*} Our original nonlinear differential equation (pendulum) has been replaced by a sequence of linear differential equations. This is a characteristic feature of the perturbation method.

since every differential equation in the chain of equations is of order two, we need to impose two initial conditions inherited from the original initial value problem. In present case, Θ(0) = 1 implies that \[ 1 = \sum_{i\ge 0} \Theta (0)\,a^i = 0 \] or \[ 0 = \left\{ \Theta_0 (0) -1 \right\} + a \Theta_1 (0) + a^2 \Theta_2 (0) + O\left( a^3 \right) . \] Hence, we require that \[ \Theta_0 (0) = 1 , \quad \Theta_1 (0) = 0, \quad \Theta_2 (0) = 0 , \quad \mbox{etc.} \] Similarly, to satisfy the homogeneity of velocity, we require that \[ \dot{\Theta}_0 (0) = 0 , \quad \dot{\Theta}_1 (0) = 0 , \quad \dot{\Theta}_2 (0) = 0 , \quad \mbox{ect.} \] Next, we solve successively the sequence of initial value problems. In the present case, we first have to solve \[ \ddot{\Theta}_0 + \Theta_0 = 0, \quad \Theta_0 (0) = 1, \quad \dot{\Theta}_0 (0) = 0 . \] This linear differential equation has the general solution \[ \Theta_0 (t) = A\,\cos t + B\,\sin t , \] with A and B being some constants. The particular solution satisfying the initial conditions becomes \[ \Theta_0 (t) = \cos t , \] which is the usual approximation, in dimensionless form.

The next approximation is identically zero, Θ₁(t) ≡ 0.

Similarly, Θ₂(t) is a solution of the following initial value problem: \[ \ddot{\Theta}_2 + \Theta_2 = \frac{1}{6}\,\Theta_0^3 = \frac{1}{6}\,\cos^3 t , \quad \Theta_2 (0) = 0, \quad \dot{\Theta}_2 (0) = 0 . \] Probably the easiest way to solve this IVP is to use the Laplace transform: \[ \left( \lambda^2 + 1 \right) \Theta^L = \frac{1}{6}\cdot \frac{\lambda \left( \lambda^2 +7 \right)}{\lambda^4 + 10 \lambda^2 + 9} , \]

where \[ \Theta^L = {\cal L} \left[ \Theta_2 (t) \right]_{t\to \lambda}= \int_0^{\infty} e^{-\lambda t} \Theta_2 (t)\,{\text d} t \] is the Laplace transform of function Θ₂(t). Solving algebraic equation and applying the inverse Laplace transform, we obtain \[ \Theta_2 (t) = {\cal L}^{-1} \left[ \Theta^L \right]_{\lambda]to t} = \frac{\sin t}{96} \left[ 6t + \sin 2t \right] H(t) , \] where \[ H(t) = \begin{cases} 1 , & \quad \mbox{if } t > 0 , \\ 1/2, & \quad \mbox{when } t=0 , \\ 0, & \quad \mbox{if } t < 0 , \end{cases} \] is the Heaviside function.

Knowing a correction to the usual approximation cost, we can now determine when this correction is small. We see that that the term proportional to the bounded quantity sint sin(2t) in the correction is small for all time provided only that 𝑎 is small. It is otherwise with the term proportional to t sint; the latter is small compared to unity if and only if 𝑎²t ⪡ 1. It is becoming clear that we have different situations, depending upon whether or not we examine the approximation for a fixed finite time interval [0, T] or for the unbounded interval [0, ^infin;).

Let us consider finite time interval first. For the time interval [0, T], the correction is small and the approximation Θ(t, 𝑎) ≈ cost appears to be consistent provided that 𝑎²T ⪡ 1. We would expect, in fact, that the approximation is asymptotically valid as 𝑎↓0, uniformly for t ∈ [0, T]. By this we mean that we anticipate that the exact solution and the approximation can be made to differ by an arbitrarily small amount for the entire range of the independent variable provided the parameter is taken sufficiently small. Symbolically, in the present case, uniform asymptotic validity of the approximation cost requires that, given any ε > 0, these exists a positive constant 𝑎₀, which may depend on ε and T but must be independent of t, such that \[ \left\vert \Theta (a, t) - \cos t \right\vert < \varepsilon \qquad \mbox{whenever } 0 \le t \le T , \quad 0 \le |a| \le a_0 . \]

Now consider the unbounded interval [0, ∞). For a given initial angular displacement 𝑎, the term 𝑎²t sint oscillates with ever-increasing amplitude. Such behavior can be rejected on physical ground. Hence, our perturbation scheme has apparently not provided an approximation that is asymptotically valid for t ∈ [0, ∞). There are known methods for obtaining improved expansions, but they will be presented later. For now, let us be content with the gact that all appears well in finite interval.

Let us look at our approximations from a slightly different point of view. At the begining of our discussion, we estimated the validity of the usual approximation cost by examining the next term in the perturbation expansion. It is important to realize the advantage of this approach compared to the basic simplification procedure of the opening section. For the pendulum problem, that procedure would entail making the approximation 𝑎⁻¹sin(𝑎Θ) ≈ Θ, finding the approximate solution &Theta₀ = cost, and checking for consistency by determining whether 𝑎⁻¹sin(𝑎Θ₀) ≈ Θ₀. Now \[ \lim_{\varepsilon \to 0} \ \frac{a^{-1} \sin a\Theta_0}{\Theta_0} = 1 \qquad \mbox{uniformly in $t$.} \] since |Θ₀| ≤ 1 for all t, so that the approximation is apparently consistent when 𝑎 is small. However, the true solution also satisfies |Θ(t)| ≤ 1. Hence, the limit above remains true when Θ is substituted for Θ₀ and the approximation is genuinely consistent. When time elapses, the period of oscillation evaluated from the first term, cost, is incorrect brining about grave errors in the prediction of the true angular displacement Θ. The deficiency in the approximation was not signaled by the basic simplification procedure, but it was revealed by the appearance of the next term in the perturbation expansion. This fact serves as an antidote to the examples of wretched consistent approximations.

It is hard to think of an instance in which the failure of an approximation is not indicated by the appearance of the next term. Be that as it may, the antidote is not universal because in complicated problems a forbidding amount of additional work may be required to compute one more term. The simple substitution required by the consistency check in the basic simplification procedure may be all that it is feasible to do.

Higher terms can be obtained by the perturbation method. However, the amount of work in calculating coefficients is hard to achive without the aid of a computer algebra system.

It is helpful to take advantage of the fact that the identical vanishing of Θ₁(t) is not accidental. Indeed, Θ_i(t) = 0 for all odd i. To see this, observe that the governing differential equation is unchanged if −𝑎 is substituted for 𝑎. Assuming that the IVP for the pendulum equation has a unique solution, this means that Θ(t, 𝑎) = Θ(t, −𝑎). In other words, Θ is an even function of its second argument. Hence, coefficients in series representation of Θ have all odd-ordered derivatives of Θ(t, 𝑎) with respect to 𝑎 vanish when 𝑎 = 0. So these coefficients satisfy \[ \Theta_k (t) = \frac{1}{k!} \left[ \frac{\partial^k \Theta (t, a)}{\partial a^k} \right]_{a=0} \equiv 0 \qquad \mbox{for odd $k$}. \] We conclude that there are no odd powers of 𝑎 in the series for Θ(t, 𝑎). To obtain higher corrections, it is best to take advantage of this and assume from the outset that Θ(t, 𝑎) has an expansion in powers of 𝑎².

By calculating more terms, we expect to improve the validity of our approximation. This is true for the bounded interval [0, T]. In fact, it can be shown that the series for Θ(t, 𝑎) is convergent when |𝑎| ≤ 𝑎₀, for some sufficiently small 𝑎₀. The convergence is uniform for t in [0, T]. That is, given any ε > 0 there exists a nonnegative integer N₀ (which is independent of t) such that \[ \left\vert \Theta (t, a) - \sum_{i=0}^N \Theta_i (t)\, a^i \right\vert < \varepsilon \] whenever N > N₀, |𝑎| ≤ 𝑎₀, t ∈ [0, T].

Implied by this convergence is the fact that the series for Θ(t, 𝑎) is asymptotic as 𝑎 → 0, uniformly in t for t ∈ [0, T[. This means that for any ε > 0 and any fixed N, there exists a positive number 𝑎₁ (which is independent of t) such that \[ a^{-N} \left\vert \Theta (t, a) - \sum_{i=0}^N \Theta_i (t)\, a^i \right\vert < \varepsilon \] whenever |𝑎| ≤ 𝑎₁, t ∈ [0, T]. In terms of the notation of asymptotic analysis, \[ \Theta (t, a) = \sum_{i=0}^N \Theta_i (t)\, a^i + o\left( a^N \right) \qquad \mbox{uniformly for } t \in [0, T] . \] Since obtaining more terms of the series is not easy at all, the asymptotic result may be more useful than the convergence proof. Be advised in distinction of these two approaches. In convergent approximation one considers the parameter as fixed and thinks of improving accuracy by taking more terms. In asymptotic approximation one considers a fixed number of terms and thinks of improving accuracy by letting the parameter approach a "favorable" value. Taking more terms does not improve the situation as far as the unbounded interval is concerned.    ■

Example 3: Consider an undamped nonlinear oscillator described by Duffing’s equation \[ y'' + y + \varepsilon y^3 = 0 , \] where the prime denotes a derivative with respect to time t. We look for solutions y(t, ε) that satisfy the initial conditions \[ y(0, \varepsilon ) = 1 , \qquad y'(0, \varepsilon ) = 0 . \] We look for straightforward expansion of an asymptotic solution as ε → 0, \[ y(t, \varepsilon ) = y_0 (t) + \varepsilon \, y_1 (t) + O\left( \varepsilon^2 \right) . \] The leading-order perturbation equations are \[ \begin{split} y''_0 + y_0 &= 0 , \\ y_0 (0) &= 1, \quad y'_0 (0) = 0 . \end{split} \] Its solution is \[ y_0 (t) = \cos t . \] The next-order perturbation equations are \[ \begin{split} y''_1 + y_1 + y_0^3 &=0 , \\ y_1 (0) = 0 , \quad y'_1 (0) &= 0 , \end{split} \] with the solution \[ y_1 (t) = \frac{1}{32} \left[ \cos 3t - \cos t \right] - \frac{3}{8}\,t\,\sin t . \] This solution contains a secular term that grows linearly in t. As a result, the expansion is not uniformly valid in t, and breaks down when t = O(ε) and εy₁ is no longer a small correction to y₀.

The solution is, in fact, a periodic function of t. The straightforward expansion breaks down because it does not account for the dependence of the period of the solution on ε. The following example illustrates the difficulty.

For instance, we have the following Taylor expansion as ε → 0: \[ \cos \left[ \left( 1 + \varepsilon \right) t \right] = \cos t - \varepsilon t\,\sin t + O\left( \varepsilon^2 \right) . \] This asymptotic expansion is valid only when t ⪡ 1/ε.

To construct a uniformly valid solution, we introduced a stretched time variable \[ \tau = \omega (\varepsilon )\,t , \] and write y = y(τ, ε). We require that y is a 2π-periodic function of τ . The choice of 2π here is for convenience; any other constant period — for example 1 — would lead to the same asymptotic solution. The crucial point is that the period of y in τ is independent of ε (unlike the period of y in t).

Since d/dt = ωd/dτ, the function y(τ, ε) satisfies \begin{align*} \omega^2 y'' + y + \varepsilon y^3 &= 0 , \\ y(0, \varepsilon ) = 1, \quad y' (0, \varepsilon ) &= 0, \\ y(\tau + 2 \pi , \varepsilon ) &= y(\tau , \varepsilon ) , \end{align*} where the prime denotes a derivative with respect to τ. We look for an asymptotic expansion of the form \begin{align*} y(\tau , \varepsilon ) &= y_0 (\tau ) + \varepsilon y_1 (\tau ) + O\left( \varepsilon^2 \right) , \\ \0mega (\varepsilon ) &= \omega_0 + \varepsilon \omega_1 + O\left( \varepsilon^2 \right) . \end{align*} Using this expansion in the equation and equating coefficients of ε⁰, we find that \begin{align*} \omega_0^2 y''_0 + y_0 &= 0 , \\ y_0 (0) &= 1 , \quad y'_0 (0) = 0 , \\ y_0 (\tau + 2\pi ) &= y_0 (\tau ) . \end{align*} The solution is \begin{align*} y_0 (\tau ) &= \cos\tau . \\ \omega_0 &= 1 . \end{align*} After setting ω₀ = 1, we find that the next order perturbation equations are \begin{align*} y''_1 + y_1 + 2\omega_1 y''_0 + y_0^3 &= 0, \\ y_1 (0) = 0 , \quad y'_1 (0) &= 0 , \\ y_1 (\tau + 2\pi ) &= y_1 (\tau ) . \end{align*} Using the solution for y₀ in the ODE for y1, we get \begin{align*} y''_1 + y_1 &= 2 \omega_1 \cos\tau - \cos^3 \tau \\ &= \left( 2\omega_1 - \frac{3}{4} \right) \cos\tau - \frac{1}{4}\, \cos 3\tau . \end{align*} We only have a periodic solution if \[ \omega_1 = \frac{3}{8} , \] and then \[ y_1 (t) = \frac{1}{32} \left[ \cos 3\tau - \cos \tau \right] . \] It follows that \begin{align*} y &= \cos t + \frac{1}{32} \,\varepsilon \left[ \cos 3\omega t - \cos \omega t \right] + O\left( \varepsilon^2 \right) , \\ \omega &= 1 + \frac{3}{8}\,\varepsilon + O\left( \varepsilon^2 \right) . \end{align*} This expansion can be continued to arbitrary orders in ε.

The appearance of secular terms in the expansion is a consequence of the nonsolvability of the perturbation equations for periodic solutions.

In the equation for y1, after replacing τ by t, we had \[ f(t) = 2\omega_1 \cos t - \cos 3t . \] This function is orthogonal to sint, and \[ \langle f, \cos t \rangle = 2\pi \left\{ 2\omega_1 \,\overline{\cos^2 t} - \overline{\cos^4 t} \right\} , \] where the overline denotes the average value, \[ \overline{f} = \frac{1}{2\pi} \int_{-\pi}^{\pi} f(t)\, {\text d}t . \] Since \[ \overline{\cos^2 t} = \frac{1}{2} , \qquad \overline{\cos^4 t} = \frac{3}{8} , \] the solvability condition implies that ω₁ = ⅜.    ■

Example 4: We will compute the amplitude of the limit cycle of the van der Pol equation with small damping \[ y'' + \varepsilon \left( y^2 -1 \right) y' + y = 0 . \] This ODE describes a self-excited oscillator, whose energy increases when |y| < 1 and decreases when |y| > 1. It was proposed by the Dutch physicist Balthasar van der Pol as a simple model of a beating heart. The ODE has a single stable periodic orbit, or limit cycle.

We have to determine both the period T(ε) and the amplitude 𝑎(ε) of the limit cycle. Since the ODE is autonomous, we can make a time-shift so that y′(0) = 0. Thus, we want to solve the ODE subject to the conditions that \begin{align*} y(t + T, \varepsilon ) &= y(t, \varepsilon ) , \\ y(0, \varepsilon ) &= a(\varepsilon ), \quad y' (0) = 0 . \end{align*} Using the Poincaré-Lindstedt method, we introduce a strained variable \[ \tau = \omega t , \] and look for a 2π-periodic solution y(τ, ε), where ω = 2π/T. Since d/dt = ωd/dτ, we have \begin{align*} \omega^2 y'' + \varepsilon \omega \left( y^2 - 1 \right) y' + y &= 0 , \\ y(\tau + 2\pi , \varepsilon ) &= y(\tau , \varepsilon ) , \\ y(0, \varepsilon ) &= a , \quad y' (0, \varepsilon ) = 0 , \end{align*} where the prime denotes a derivative with respect to τ We look for asymptotic expansions, \begin{align*} y (\tau , \varepsilon ) &= y_0 (\tau ) + \varepsilon y_1 (\tau ) + O\left( \varepsilon^2 \right) , \\ \omega (\varepsilon ) &= \omega_0 + \varepsilon \omega_1 + O\left( \varepsilon^2 \right) , \\ a(\varepsilon ) &= a_0 + \varepsilon a_1 + O\left( \varepsilon^2 \right) . \end{align*} Using these expansions in the equation and equating coefficients of ε⁰, we find that \begin{align*} \omega_0^2 y''_0 + y_0 &= 0 , \\ y_0 (\tau + 2\pi ) &= y_0 (\tau ) , \\ y_0 (0) &= a_0 , \quad y'_0 (0) = 0 . \end{align*} The solution is \begin{align*} y_0 (\tau ) = a_0 \cos \tau , \\ \omega_0 &= 1 . \end{align*} The next order perturbation equations are \begin{align*} y''_1 + y_1 + 2 \omega_1 y''_0 + \left( y_0^2 -1 \right) y'_0 &= 0 , \\ y_1 (\tau + 2\pi ) &= y_1 (\tau ) , \\ y_1 (0) &= a_1 , \quad y'_1 (0) = 0 . \end{align*} Using the solution for y₀ in the ODE for y₁, we find that \[ y''_1 + y_1 = 2\omega_1 \cos\tau + a_0 \left( a_0^2 \cos^2 \tau -1 \right) \sin\tau . \] The solvability conditions, that the right and side is orthogonal to sin τ and cos τ imply that \[ \frac{1}{8}\,a_0^3 - \frac{1}{2}\,a_0 = 0, \quad \omega_1 = 0 . \] We take 𝑎₀ = 2, the solution 𝑎₀ = 2 2 corresponds to a phase shift in the limit cycle by π, and 𝑎₀ = 0 corresponds to the unstable steady solution y = 0. Then \[ y_1 (\tau ) = - \frac{1}{4}\,\sin 3\tau + \frac{3}{4}\,\sin\tau + \alpha_1 \cos\tau . \] At the next order, in the equation for y;₂, there are two free parameters, ( 𝑎₁, ω;₂), which can be chosen to satisfy the two solvability conditions. The expansion can be continued in the same way to all orders in ε.    ■

1. Kato, T., Perturbation Theory for Linear Operators, Springer, 1980.
2. Lin, C.C. and Segel, L.A., Mathematics Applied to Deterministic Problems in the Natural Sciences, SIAM, Society for Industrial and Applied Mathematics, 1988.