Asymptotic Methods

Mohammed Akram Fellah

► To cite this version:

Mohammed Akram Fellah. Asymptotic Methods. 2021. hal-03338677

HAL Id: hal-03338677

https://hal.science/hal-03338677

Preprint submitted on 9 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Asymptotic Methods

Mohammed Akram Fellah

Department of Applied Mathematics, University of Waterloo

Waterloo, Ontario, Canada

mafellah@uwaterloo.ca

September 6, 2021

In this lecture notes, we will introduce Asymptotics, then we will give a short glimpse on Perturbation theory (regular versus singular), which plays a crucial role especially in theoretical physics.

Our goal is to find Asymptotic series that approximates the values of integrals depending on some parameter or the solutions of differential equations. These methods that lead to obtain more effective algorithms of numerical evaluation are called: Asymptotic Methods.
# Contents

<table><tr><td>1</td><td>Introduction to Asymptotics</td><td>3</td></tr><tr><td>1.1</td><td>A simple example</td><td>3</td></tr><tr><td>2</td><td>Dimensional Analysis</td><td>5</td></tr><tr><td>3</td><td>Series in General</td><td>7</td></tr><tr><td>3.1</td><td>Tests for Convergence (or divergence)</td><td>8</td></tr><tr><td>3.2</td><td>Radius of convergence</td><td>9</td></tr><tr><td>4</td><td>Asymptotic Series</td><td>9</td></tr><tr><td>4.1</td><td>Order relations</td><td>9</td></tr><tr><td>4.2</td><td>Asymptotic expansions</td><td>10</td></tr><tr><td>5</td><td>A First Glimpse of Perturbation Theory</td><td>12</td></tr><tr><td>5.1</td><td>Regular versus singular perturbations</td><td>12</td></tr><tr><td>5.2</td><td>Scaling</td><td>13</td></tr><tr><td>6</td><td>The Regular Perturbation Approach</td><td>14</td></tr><tr><td>7</td><td>Invariant Manifolds in ODEs</td><td>14</td></tr><tr><td>8</td><td>Centre Manifolds</td><td>17</td></tr></table>
<table>
  <tbody>
    <tr>
      <td>9</td>
      <td>Polynomial Equations and the Newton Polygon</td>
      <td>21</td>
    </tr>
    <tr>
      <td>10</td>
      <td>The Curse of the Secular Terms</td>
      <td>24</td>
    </tr>
    <tr>
      <td>11</td>
      <td>Poincaré-Linstedt Method</td>
      <td>26</td>
    </tr>
    <tr>
      <td>11.1</td>
      <td>Shohat's expansion</td>
      <td>28</td>
    </tr>
    <tr>
      <td>12</td>
      <td>The Method of Multiple Scales</td>
      <td>29</td>
    </tr>
    <tr>
      <td>12.1</td>
      <td>Multiple Scales: a BVP</td>
      <td>32</td>
    </tr>
    <tr>
      <td>13</td>
      <td>Matched Asymptotic Expansions</td>
      <td>35</td>
    </tr>
    <tr>
      <td>13.1</td>
      <td>MAE: Higher Order Matching</td>
      <td>37</td>
    </tr>
    <tr>
      <td>13.2</td>
      <td>Scaling boundary layers</td>
      <td>38</td>
    </tr>
    <tr>
      <td>13.3</td>
      <td>Locating boundary layers</td>
      <td>39</td>
    </tr>
    <tr>
      <td>13.4</td>
      <td>A problem with an interior boundary layer</td>
      <td>40</td>
    </tr>
    <tr>
      <td>14</td>
      <td>The WKB Method</td>
      <td>42</td>
    </tr>
    <tr>
      <td>14.1</td>
      <td>The main idea</td>
      <td>42</td>
    </tr>
    <tr>
      <td>14.2</td>
      <td>Eigenvalue problems</td>
      <td>46</td>
    </tr>
    <tr>
      <td>15</td>
      <td>Averaging</td>
      <td>47</td>
    </tr>
    <tr>
      <td>15.1</td>
      <td>Averaging: the 2-dimensional case</td>
      <td>49</td>
    </tr>
    <tr>
      <td>16</td>
      <td>A taste of renormalization</td>
      <td>51</td>
    </tr>
    <tr>
      <td>17</td>
      <td>Optimal Truncations</td>
      <td>53</td>
    </tr>
    <tr>
      <td>18</td>
      <td>Borel summation</td>
      <td>55</td>
    </tr>
    <tr>
      <td>19</td>
      <td>Asymptotics of Integrals</td>
      <td>57</td>
    </tr>
    <tr>
      <td>19.1</td>
      <td>Watson's Lemma</td>
      <td>58</td>
    </tr>
    <tr>
      <td>19.2</td>
      <td>Laplace's method</td>
      <td>61</td>
    </tr>
    <tr>
      <td>19.3</td>
      <td>Method of stationary phase</td>
      <td>65</td>
    </tr>
  </tbody>
</table>
# 1 Introduction to Asymptotics

**Asymptotics** is the art and the science of using the smallness of the small, the largeness of the large, and then trying to interpolate.

Many great names in the history of Mathematics have contributed to this area, including Euler, Laplace, Lagrange, Poincaré, Borel, and recently Ecalle, Ramis, and Berry.

The questions (spanning applied mathematics, mathematical physics, and dynamical systems theory) that Asymptotics deals with can be classified in under four very rough headings:

I. Given the system of ODEs $\dot{x} = f(x)$, $f(0) = 0$, want to understand the phase space picture for $\|x\|$ small.

II. Given the system of ODEs $\dot{x} = f(x, t)$, $x(0) = x_0$ (in initial value problem), want to understand the behaviour of the solution $x(t)$ for $t$ small or for $t$ large ("asymptotically in time").

III. Given $\epsilon \ll 1$, a "small parameter", want to understand the integral $\int_{I(\epsilon)} F(x, \epsilon) dx$.

IV. Given $\epsilon \ll 1$, want to understand $F(x, \epsilon) = 0$, where this can be a scalar algebraic equation, a system thereof, a system of ODEs (both initial and boundary value problems), a PDE, etc.

**Notes:** 1. $F(x, \epsilon) = 0$ can be considered as a "perturbation" of $F(x, 0) = 0$. Perturbation Theory is a major tool in Physics, and so can be considered as a branch of asymptotics.

2. $\epsilon$ small is equivalent to $1/\epsilon$ large.

3. By "understanding" in applied mathematics one classically meant finding an approximate solution (and understanding the nature of the error of approximation). Recently though the interest has shifted to understanding why different approximation methods give the results they do and trying to figure out what this says about the problem.

Asymptotics is the door through which divergent series enter mathematics. Since much of asymptotics is algorithmic, this area lies at the crossroads of applied mathematics, complex analysis, and computer algebra.

## 1.1 A simple example

It is always good to know where the words we use come from. Asymptote comes from Greek a-symptotos, "not-together-thrown", and in elementary mathematics sands for a straight line which is never quite reached by a curve. Let us see a simple asymptotic analysis of asymptotes.

For $x > 1$, let consider $f(x) = x^2/(x - 1)$. Its graph obviously has the vertical asymptote $x = 1$,
and we will show that it also has a slant asymptote.

First, let us see what happens as $x \to 1^+$. Then $x - 1$ is small and positive; let us put $y = x - 1$. In terms of $y$ we have

$$f(y) = \frac{(y+1)^2}{y} = \frac{1}{y} + 2 + y.$$

This is an exact equality; we have written $f$ as a finite sum of terms, each one of which is larger than then next if $y$ is small, with the **leading term** being, not surprisingly, $1/y$.

For large $x$, we find by polynomial division

$$f(x) = \frac{x^2}{x-1} = x + 1 + \frac{1}{x-1}.$$

We could leave it like that, as in the previous case, but let us rewrite this in terms of powers of $x$.

Since for $|z| < 1$, we have the famous result

$$\frac{1}{1-z} = 1 + z + z^2 + \dots, \qquad (1)$$

we write

$$f(x) = x + 1 + \frac{1}{x-1} = x + 1 + \frac{1}{x} \left[ 1 + \frac{1}{x} + \frac{1}{x^2} + \dots \right],$$

and again, the equality here makes sense as the series on the rhs of the above equation is convergent.

**Important!** Please review the notion of a convergent series (in terms of partial sums) and the various (ratio, Leibniz) tests for convergence.

To summarise this example: when $x - 1$ was small and when $x$ was large (that is, $1/x$ was small) we found a representation of the function as a series in powers of the small parameter (in the second case, the leading term was $x = (1/x)^{-1}$), in such a way that each successive term was smaller than the preceding one. Note that in the first case the series was finite and in the second one, infinite but convergent.

This, roughly, is what we do in asymptotics, except that the series we generate, often algorithmically, can fail to converge quite spectacularly. Let us now go back and play with the series (1).

**Exercise:** Prove that

$$\frac{1}{1-z} = 1 + z + z^2 + \dots + z^n + \frac{z^{n+1}}{1-z}.$$

**Exercise:** What the number $n$ of terms that you have to sum if you want to approximate 100 by

$$100 = \frac{1}{1 - 0.99} \approx 1 + 0.99 + \dots + 0.99^n$$

with the error of less than $1/100$?

The first exercise is a useful result we will use later, and the second, (which uses the first) shows that convergent series can converge very slowly and involve a lot of computation.
Now consider (1) beyond its radius of convergence. If $z = 1$ the series does not make much sense as we are dividing by zero; if $z = 2$ we have the very suspect statement that $-1$ "equals" $1 + 2 + \dots$. However, if $z = -1$, the result looks somewhat justifiable:

$$ \frac{1}{2} = 1 - 1 + 1 - 1 + 1 - 1 + \dots $$

Euler believed that the infinite sum in the rhs was $1/2$ by a simple argument: if it is $S$, then it is also $1 - S$, hence $S = 1/2$.

There are theories ("summability theories") which assign values to (some) non-convergent series, such as the one above, while assigning to convergent series their classical sum; Borel and Cesáro are examples of such theories. Much work in asymptotics is done by creating a divergent series representation (of an integral, solution of an ODE, etc.), trying to sum it by some reasonable summability theory, and then thinking what is the meaning of the result.

## 2 Dimensional Analysis

In perturbation theory, we often utilize a small parameter. It is therefore important to understand how to verify that the equations at hand involve one. We shall illustrate the ideas involved using the equation of a driven, damped linear oscillator with rigid amplitude constraints. The equation of motion is

$$ m \frac{d^2 x}{dt^2} + c \frac{dx}{dt} + kx = f \cos \omega t, \quad |x| < x_0. \qquad (2) $$

Here $x$ is the displacement of the center of mass. The equation in this, dimensional, form involves 6 parameters: $m$, the mass, $c$, the damping constant, $k$, the spring constant, $f$, the amplitude of the forcing, $\omega$, the frequency of the forcing, and the maximal displacement, $x_0$. At $|x| = x_0$ the equation is supplemented by an "impact rule" of the form

$$ \frac{dx}{dt}(t^+) = -r \frac{dx}{dt}(t^-), \quad (t^+ - t^-) \to 0, $$

where $0 < r \le 1$ is the coefficient of restitution.

It is not clear which of the terms can be said to be premultiplied by a small parameter. Consider the first term. Whether $m$ is even perceived as "small" or not depends on the units in which it is measured. Changing the units of time or length will also necessarily change the coefficient of $d^2x/dt^2$. In any case, comparing quantities having different dimensions is futile. Thus we need to work with quantities that are independent of the units we use. The process of rewriting the equations in a form involving such, non-dimensional, quantities, is called *non-dimensionalization* and the procedure is as follows

**Step 1.** List all variables and parameters together with their dimensions.

So we have
<table><thead><tr><th>Variable</th><th>Dimension</th></tr></thead><tbody><tr><td>Displacement, <i>x</i></td><td>meters</td></tr><tr><td>Time, <i>t</i></td><td>seconds</td></tr><tr><th>Parameters</th><th>Dimension</th></tr><tr><td>Maximal displacement, <i>x</i><sub>0</sub></td><td>meters</td></tr><tr><td>Mass, <i>m</i></td><td>kg</td></tr><tr><td>Damping coefficient, <i>c</i></td><td>kg/seconds</td></tr><tr><td>Spring constant, <i>k</i></td><td>kg/seconds<sup>2</sup></td></tr><tr><td>Forcing amplitude, <i>f</i></td><td>kg×m/second<sup>2</sup></td></tr><tr><td>Forcing frequency, <i>ω</i></td><td>1/seconds</td></tr></tbody></table>

Note that each term in the equation has the dimensions of force, as it should be, since (2) is just Newton's second law of motion.

**Step 2.** For each variable $v$ in the equation, form a group of parameters $p$ with the same dimensions as $v$ and introduce $v/p$ as a new variable. This step can be done in any number of ways; usually common sense and an ability to find a physical explanation for the resulting *dimensionless* groups of parameters serve as a guide. Clearly, $x$ has the same dimension as $x_0$. Thus we have a new displacement variable,

$$X = \frac{x}{x_0}.$$

Note that this variable always takes values between 0 and 1. As far as time is concerned, we could choose as our group of parameters with dimension of $t$ (seconds), for example, $p = m/c$, but the usual choice is

$$p = \sqrt{m/k}, \qquad (3)$$

so that the new time variable is

$$T = \frac{t}{\sqrt{m/k}}.$$

**Step 3.** Now we express the old variables in terms of the new variables and plug them into the equation, use the chain rule, and obtain

$$kx_0 \frac{d^2 X}{dT^2} + \frac{cx_0 \sqrt{k} dX}{\sqrt{m}dT} + kx_0 X = f \cos(\omega \sqrt{m/k}T), \quad |X| < 1.$$

Dividing all the terms by $kx_0$, we have

$$\frac{d^2 X}{dT^2} + \frac{c}{\sqrt{mk}} \frac{dX}{dT} + X = \frac{f}{kx_0} \cos(\omega \sqrt{m/k} T), \quad |X| < 1.$$

Finally, introducing the new dimensionless parameter groups $2\alpha = c/\sqrt{km}$, $F = f/(kx_0)$, and $\Omega = \omega\sqrt{m/k}$, we obtain the non-dimensionalized form of (2):

$$\frac{d^2 X}{dT^2} + 2\alpha \frac{dX}{dT} + X = F \cos \Omega T, \quad |X| < 1. \qquad (4)$$
**Notes.** 1. The choice of the time-scale (3) is justified by referring to the undamped, undriven harmonic oscillator equation.

2. Notice the three important things about (4):

* It involves only three parameters, as compared with the original six.

* Nothing in it depends on the choice of units of measurement, so that statements such as “F is small” make sense.

* One can talk about “geometrically similar” systems. For example, two systems with different $f$ and $x_0$ will behave in the same way if the ratio $f/x_0$ is the same in both.

Non-dimensionalization provides the justification for the (enlightened) use of scale models.

Non-dimensionalization can be formalised; the result is the branch of applied mathematics known as dimensional analysis. Here is an easy non-trivial application of the kind of results one can get by dimensional analysis. Consider a freely oscillating mass $m$ suspended on a mass-less string of length $l$. Suppose that the arc of the swing does not influence the period of the oscillation. Show by dimension arguments that the period of oscillation is independent of $m$ and derive a formula for it up to a constant.

## 3 Series in General

We start with a sequence of real numbers, {$a_n$}, $n = 0, 1, \dots$. The object of interest is the infinite series,

$$ \sum_{n=0}^{\infty} a_n . $$

Call the result of this summation $S$, the sum of the infinite series. The basic question are, what does $S$ mean, that is, how it is defined, and for what sequences {$a_n$} is $S$ finite. To define $S$ we introduce partial sums $S_N$,

$$ S_N = \sum_{n=0}^{N} a_n . $$

Then we define the infinite series, $S$, to be the limit of the sequence {$S_N$} as $N \to \infty$,

$$ S = \lim_{N \to \infty} S_N. \qquad (5) $$

Straight away we can see examples of sequences, for which finite sums do not converge to a limit. A good one we discussed yesterday is

$$ a_n = (-1)^{n-1} . $$
Then $S_N = 0$ for $N$ even and $S_N = 1$ for $N$ odd. So clearly (57) does not exist.

A much better behaved example is

$$a_n = x^{n-1}, \text{ where } -1 < x < 1.$$

This is the well known geometric series. Then

$$S_N(x) = 1 + x + x^2 + \dots + x^{N-1} = \frac{1-x^N}{1-x},$$

so that $S(x) = 1/(1-x)$. Note that if $x = -1$, we are back to the first, oscillating example. Also note that as $x \to 1^-$, $S \to \infty$. If $|x| > 1$ then $|a_n|$ grow without bound, and one cannot hope to have a finite sum of the series.

## 3.1 Tests for Convergence (or divergence)

So, if $S := \lim_{N \to \infty} S_N$ exists, we say that the series $\sum_{n=0}^{\infty} a_n$ is **convergent**, and that its sum is $S$. If the limit of the sequence of partial sums does not exist, we say that the series is **divergent**.

A simple test for divergence is the following:

**Proposition 1.** If $a_n$ does not tend to zero as $n \to \infty$, the series $\sum_{n=0}^{\infty} a_n$ diverges.

From this result it is easy to check that $\sum_{n=0}^{\infty} x^n$ diverges if $|x| > 1$.

Note that the above proposition does not mean that if $a_n \to 0$ as $n \to \infty$, the series $\sum_{n=0}^{\infty} a_n$ converges. In fact, we will show below that the harmonic series

$$1 + \frac{1}{2} + \frac{1}{3} + \frac{1}{4} + \dots$$

diverges. On the other hand, if the sequence is alternating, that is, if $\operatorname{sign}(a_n \cdot a_{n+1}) = -1$ for all $n$, we have the Leibniz test:

**Proposition 2.** If the sequence $\{a_n\}$ is alternating decreasing sequence such that $a_n \to 0$ as $n \to \infty$, then the series $\sum_{n=0}^{\infty} a_n$ converges.

Thus

$$1 - \frac{1}{2} + \frac{1}{3} - \frac{1}{4} + \dots$$

converges, in fact to $\ln 2$.

Amaze your friends. We know that

$$1 - \frac{1}{2} + \frac{1}{3} - \frac{1}{4} + \dots = \ln 2.$$

Multiply both sides by 2. Thus

$$2 - 1 + \frac{2}{3} - \frac{1}{2} + \frac{2}{5} - \frac{1}{3} + \frac{2}{7} - \frac{1}{4} + \dots = 2 \ln 2.$$
Now rearrange the latter series thusly:

$$ (2-1) - \frac{1}{2} + \left(\frac{2}{3} - \frac{1}{3}\right) - \frac{1}{4} + \left(\frac{2}{5} - \frac{1}{5}\right) - \dots = 2 \ln 2. $$

However, the right-hand side is just the original series, so that $\ln 2 = 2 \ln 2$, or $1 = 2...$. What went wrong here?

The series for which $\sum_{n=0}^{\infty} |a_n|$ converges are called **absolutely convergent**. The series in the example above is not absolutely convergent (thus it is called conditionally convergent).

Next we quote some widely used tests for convergence.

**Proposition 3.** Let $\lim_{n \to \infty} |a_{n+1}/a_n| = L$. If $L < 1$, the series $\sum_{n=0}^{\infty} a_n$ converges.

This is the *ratio test*. Note that for the harmonic series $L = 1$.

Now we state the *integral test*. To that end, suppose there exists a continuous decreasing function $f(x)$ such that $a_n = f(n)$. Then

**Proposition 4.** $\sum_{n=0}^{\infty} a_n$ converges if $\int_0^\infty f(x) dx$ exists and diverges otherwise.

This gives a quick proof of divergence of the harmonic series. Here $f(x) = 1/(x+1)$.

## 3.2 Radius of convergence

Finally, consider the power series $A(x) := \sum_{n=0}^{\infty} a_n x^n$. Let

$$ R = \sup_{|x|} \left\{ \sum_{n=0}^{\infty} a_n x^n \text{ converges} \right\}. $$

We call $R$ the **radius of convergence** of the series $A(x)$. For example, if $a_n x^n = x^n$, the radius of convergence is 1, and if $a_n x^n = x^n/n!$ the radius of convergence is infinite (that can be easily proved by the ratio test).

# 4 Asymptotic Series

## 4.1 Order relations

Suppose we have two functions, $f(x)$ and $g(x)$ defined in a domain $U$ and that a point $x_0 \in \bar{U}$ (the closure of $U$; the definitions below work perfectly well in the complex domain). Then we have

**Definition 5.** $f(x)$ is "large $O$ of $g(x)$," written as $f(x) = O(g(x))$ as $x \to x_0$ in $U$ if there exists a number $A > 0$, such that for some neighbourhood $N$ of $x_0$ we have that $|f(x)| \le A|g(x)|$ for all $x \in N \cap U$.
So, for example, $\sin x = O(x)$ as $x \to 0$, but $\sin x \neq O(x^2)$.

The relation between these two functions is described by

**Definition 6.** $f(x)$ is “little o of $g(x)$,” written as $f(x) = o(g(x))$ as $x \to x_0$ in $U$ if for any number $A > 0$ there exists a neighbourhood $N_A$ of $x_0$ so that we have that $|f(x)| \le A|g(x)|$ for all $x \in N_A \cap U$.

Thus $x^2 = o(\sin x)$ as $x \to 0$. To get a good grasp on these definitions consider that they imply that if $g(x) \neq 0$ in $U \setminus x_0$, then $f(x) = o(g(x))$ means that $f(x)/g(x)$ goes to zero as $x \to x_0$, while $f(x) = O(g(x))$ means that $f(x)/g(x)$ remains bounded $x \to x_0$. Of course, if $f(x) = o(g(x))$, then $f(x) = O(g(x))$ (for all $A > 0$ in the above definition).

Another point worth noticing is that $f(x) = O(g(x))$ does not imply $f'(x) = O(g'(x))$.

**Exercise** Concoct an example to show that.

A very good example is the following: $f(x) = x \sin(1/x)$ and $g(x) = x$. Since $|\sin x| \le 1$ for all $x$, we clearly have that $f(x) = O(g(x))$ as $x \to 0$. On the other hand, as $x \to 0$ $\lim \sin(1/x)$ does not exist, so that $f(x) \neq o(g(x))$. Moreover, since $|1/\sin(1/x)|$ is not bounded as $x \to 0$, $g(x) \neq O(f(x))$.

Finally, we shall need below the following relations:

$$O(o(f)) = o(f) = o(O(f)) \text{ as } x \to x_0. \qquad (6)$$

The first of these equalities means the following: Let $g = o(f)$, $h = O(g)$. Then the statement is that $h = o(g)$.

**Exercise** Prove these statements.

## 4.2 Asymptotic expansions

Armed with these definitions, we can finally define asymptotic expansions. First consider a sequence of continuous functions $(\phi_n(x))_{n=1}^{\infty}$ defined on $U$.

**Definition 7.** $(\phi_n(x))_{n=1}^{\infty}$ is an asymptotic sequence as $x \to x_0$ in $U$ if

$$\phi_{n+1}(x) = o(\phi_n(x))$$

as $x \to x_0$.

Obvious examples are $\phi_n(x) = x^n$, $\phi_n(x) = \sin^n x$ as $x \to 0$.

Now let $f(x)$ be continuous in $U$ and consider a series $\sum_{n=0}^{\infty} a_n \phi_n(x)$, where $(\phi_n(x))_{n=1}^{\infty}$ is an asymptotic sequence as $x \to x_0$. We shall call such a series “formal,” since no statement concerning its convergence is made for any $x \in U$. Then we have the following definition.
**Definition 8.** The formal series $\sum_{n=0}^{\infty} a_n \phi_n(x)$, where $(\phi_n(x))_{n=1}^{\infty}$ is an asymptotic sequence as $x \to x_0$ is called an (infinite) asymptotic series expansion (of Poincaré type) of $f(x)$ as $x \to x_0$ if

$$f(x) - \sum_{n=0}^{m} a_n \phi_n(x) = O(\phi_{m+1}(x))$$

for all non-negative $m$ as $x \to x_0$.

This relation between the function $f(x)$ and the series is henceforth denoted by

$$f(x) \sim \sum_{n=0}^{\infty} a_n \phi_n(x).$$

Note that from the above definition and (6) we can immediately derive the following two relations:

$$\lim_{x \to x_0} \frac{f(x) - \sum_{n=0}^{m} a_n \phi_n(x)}{\phi_m(x)} = 0, \qquad (7)$$

and

$$a_{m+1} = \lim_{x \to x_0} \frac{f(x) - \sum_{n=0}^{m} a_n \phi_n(x)}{\phi_{m+1}(x)}. \qquad (8)$$

For $a_0$ we have

$$a_0 = \lim_{x \to x_0} \frac{f(x)}{\phi_0(x)}.$$

Note that from the second of these results it follows that (given an asymptotic sequence!) the asymptotic expansion is unique (argue as usual: suppose it is not ...).

**Note:** Consider the function $f(x) = e^{-x}$ and take $\phi_n(x) = x^{-n}$. Clearly, $\phi$'s are an asymptotic sequence as $x \to \infty$. Verify that in the asymptotic expansion of $f(x)$ we have $a_n = 0$ for all $n$. This teaches us two things. First, $f(x)$ is an example of a function that is transcendentally small as $x \to \infty$. Second, that any function $g(x)$ and $g(x) + e^{-x}$ will have the same asymptotic expansion with the asymptotic sequence as above as $x \to \infty$.

In perturbation problems, the asymptotic series generated are of the following special form. For definiteness assume that in $F(x, \epsilon) = 0$ the unknown $x$ is a function of the independent variable $t$. Then we seek an asymptotic representation

$$x(t) \sim \sum_{n=0}^{\infty} a_n(t) \phi_n(\epsilon),$$

where $\{\phi_n(\epsilon)\}$ is an asymptotic sequence as $\epsilon \to 0$ for each $t$ in some prescribed set $I$.

For example for an IVP, $I$ might be all of $\mathbb{R}_+$, while for a BVP it might be an interval, say $I = [0, 1]$, on which the problem is defined.

Suppose

$$x(t) \sim 1 + \epsilon t + \dots$$
This is a fine (start of an) asymptotic expansion on $[0, 1]$, as then $|\epsilon t| \le \epsilon$, and the second term is $o(1)$ for all allowable $t$ as $\epsilon \to 0$, but the same series is not an asymptotic series on $\mathbb{R}_+$ as after time of $O(1/\epsilon)$ the second term is not $o(1)$. In this case we say the asymptotic expansion is not **uniformly valid** on $\mathbb{R}_+$.

Similarly,

$$x(t) \sim 1 + \ln(1 - t)\epsilon + \dots$$

is not a uniformly valid asymptotic expansion on $[0, 1]$ as at $t = 1 - \exp(-1/\epsilon)$, $\ln(1 - t)\epsilon$ is not $o(1)$.

Very often the asymptotic sequence $(\phi_n(\epsilon))_{n=1}^{\infty}$ chosen has $\phi_n(\epsilon) = \epsilon^n$, but that is **not necessarily** the case; sometimes if $\phi_m(\epsilon) = \epsilon^m$ for some $m$, the next term might be not $\phi_{m+1}(\epsilon) = \epsilon^{m+1}$, but for example $\phi_{m+1}(\epsilon) = \ln(\epsilon)\epsilon^{m+1}$; remember the remark about the art of Asymptotics!

So the main task of asymptotics is the construction of asymptotic series. These series do not have to be convergent!

**Example: Let**

$$f(x) \sim \sum_{n=1}^{\infty} n!x^n.$$

Then the series in the rhs is an asymptotic series as $x \to 0$ but it has **zero** radius of convergence: if $a_n = x^n n!$, then

$$\left| \frac{a_{n+1}}{a_n} \right| = (n + 1)|x|,$$

which is larger than 1 for each $x$ for $n$ sufficiently large.

# 5 A First Glimpse of Perturbation Theory

## 5.1 Regular versus singular perturbations

Perturbation problems come in two guises: singular and regular. To see the differences between these two types of problem, let us look at two very simple examples. First, consider the quadratic equation

$$x^2 - 3x + 2 + \epsilon = 0. \qquad (9)$$

Its solutions are easily found to be

$$x_1 = \frac{3}{2} + \frac{\sqrt{1 - 4\epsilon}}{2}$$

and

$$x_2 = \frac{3}{2} - \frac{\sqrt{1 - 4\epsilon}}{2}.$$
The crucial point is that both $x_1$ and $x_2$ are continuous functions of $\epsilon$ at $\epsilon = 0$, that is, solutions for $\epsilon = 0$ ($x_1 = 2$ and $x_2 = 1$) "continue" nicely to solutions for $\epsilon \neq 0$ small.

This is the hallmark of regular perturbation problems: addition of a term involving a small parameter changes solutions only slightly.

**Exercise:** Graph $x_1$ and $x_2$ as functions of $\epsilon$. What happens as $\epsilon$ grows?

One should compare the behaviour of solutions of (9) with behaviour of solutions of the following quadratic equation:

$$ \epsilon y^2 - 3y + 2 = 0. \qquad (10) $$

If we set now $\epsilon = 0$ we obtain only one solution, $y = \frac{2}{3}$. Moreover, if we write down explicitly the dependence of solutions of (10) on $\epsilon$, that is,

$$ y_1 = \frac{3 + \sqrt{9 - 8\epsilon}}{2\epsilon} $$

and

$$ y_2 = \frac{3 - \sqrt{9 - 8\epsilon}}{2\epsilon}. $$

It is easy to see that as $\epsilon \to 0_+$, $y_1 \to \infty$, while $y_2 \to \frac{2}{3}$ (use L'Hôpital's rule to convince yourself of that).

Thus we see that

i. the number of solutions for $\epsilon = 0$ and for $\epsilon \neq 0$ small is not the same;

ii. as $\epsilon \to 0$, one of the solutions "blows up."

Both these phenomena are distinguishing signs of singular perturbation: inclusion of a "small" term has a drastic influence on the structure of the set of solutions.

## 5.2 Scaling

An important tool in treating singular perturbation problems is that by *scaling* the variables we can turn a singular perturbation problem into a regular perturbation one. Here we just illustrate the main idea using (10).

Introduce the new variable $z = \epsilon y$. Then $z$ satisfies the equation

$$ z^2 - 3z + 2\epsilon = 0, \qquad (11) $$

which, clearly, is a regular perturbation problem.

Of course the question now is raised, how did I know to make this particular scaling...
# 6 The Regular Perturbation Approach

Recall: we started with a singularly perturbed problem,

$$ \epsilon y^2 - 3y + 2 = 0. $$

Let us see how to find an asymptotic expansion for the root of this equation that goes to infinity as $\epsilon \to 0$. Recall that the scaling $z = \epsilon y$ led us to the regular perturbation problem (11).

The regular perturbation approach involves expanding the solution in an asymptotic series using an asymptotic sequence of functions of $\epsilon$.

Consider $z$ which solves (11). Write it as

$$ z \sim \sum_{n=0}^{\infty} z_n \epsilon^n. $$

Then we have

$$ z_0^2 - 3z_0 + \epsilon(2z_0z_1 - 3z_1 + 2) + O(\epsilon^2) = 0. \quad (12) $$

From this it follows that

$$ z_0^2 - 3z_0 = 0 $$

by taking the limit as $\epsilon \to 0$. This equation gives us two solutions, $z_0 = 0$ and $z_0 = 3$. We take $z_0 = 3$ (as if we take $z_0 = 0$, we will not get a solution $y$ that blows up as $\epsilon \to 0$). $z_1$ satisfies the equation

$$ 2z_0z_1 - 3z_1 + 2 = 0, $$

obtained by collecting all the terms multiplying $\epsilon^1$ and putting them to zero. [More precisely, we divide (12) by $\epsilon$ and again take the limit as $\epsilon \to 0$.]

Thus we have $z_1 = -2/3$. In the original variable $y$, this corresponds to the solution

$$ y \sim \frac{3}{\epsilon} - \frac{2}{3} + \dots, $$

which blows up as $\epsilon \to 0_+$; this representation also shows the rate of the blow up.

**Remarks:** Note that after the $O(1)$ terms, we are solving linear equations at each step. Note also that the algorithm can easily be implemented in something like Sage.

# 7 Invariant Manifolds in ODEs

Centre manifold theory is a relatively new branch of dynamical systems theory, which has found many uses in addition to the determination of stability in nonhyperbolic rest points of ODEs. My main reason in introducing it is that the computation of approximations is algorithmic and that it immediately confronts us with divergent series. Unfortunately, a bit of preparation is needed.
We start by a short discussion of stable, unstable, and centre subspaces in linear systems. Consider a linear system

$$x' = Ax, \quad x(0) = x^0,$$

where $A$ is a $2 \times 2$ (or in general, $n \times n$) matrix. Its solution can be written in the form

$$x(x^0, t) = e^{At}x^0.$$

Suppose that $v_j$ is an eigenvector of $A$ corresponding to (a real) eigenvalue $\lambda_j$. Then it is not hard to prove that $v_j$ is an eigenvector of $e^{At}$ corresponding to the eigenvalue $e^{\lambda_j t}$ for all $t$.

Now observe

$$x(Cv_j, t) = e^{tA}Cv_j = Ce^{t\lambda_j}v_j,$$

which means that if the initial condition points in the direction of $v_j$, so will the solution $x(Cv_j, t)$ for all time. Thus the set $\text{Sp}(v_j)$ is **invariant** under the flow generated by $x' = Ax$

If $v_i$ is another eigenvector, then clearly $\text{Sp}(v_j, v_i)$ is invariant. In general, for $n$-dimensional linear system with eigenvectors of $A$ being $v_1, \dots, v_n$, the linear subspace $\text{Sp}(v_1, \dots v_k)$, $k \le n$ is invariant. Note that if $\lambda_j$ is a complex eigenvalue, $\lambda_j^*$, its complex conjugate, is also an eigenvalue and the subspace $\text{Sp}(\text{Re}(v_j), \text{Im}(v_j))$ is also invariant.

So if $A$ is $n \times n$, we can decompose the whole space $\mathbb{R}^n$ into a direct sum of three invariant subspaces:

$$E^s = \text{stable subspace} = \text{Sp}\{v_j | \text{Re}(\lambda_j) < 0\},$$

$$E^u = \text{unstable subspace} = \text{Sp}\{v_j | \text{Re}(\lambda_j) > 0\},$$

$$E^c = \text{centre subspace} = \text{Sp}\{v_j | \text{Re}(\lambda_j) = 0\}.$$

$$\text{Let dim } E^s = n^s, \text{ dim } E^u = n^u, \text{ dim } E^c = n^c. \text{ Clearly } n^s + n^u + n^c = n.$$

How do the solutions of the linear system $x' = Ax$ behave on the various subspaces? On $E^s$ (look at its definition!) solutions decay exponentially. On $E^u$ they grow exponentially, while on $E^c$ neither of those two alternatives obtain. On $E^c$ many things can happen: solutions can be constant, they can oscillate, or have, say, polynomial growth in time. Let us do an example. Consider $x' = Ax$ with

$$A = \begin{bmatrix} 0 & 0 \\ 1 & 0 \end{bmatrix}.$$

Then $A$ has two zero eigenvalues and $E^c = \mathbb{R}^2$. What happens to the solutions? It is not hard to solve the system of equations to get

$$x(t) = \begin{bmatrix} 1 & 0 \\ t & 1 \end{bmatrix} x^0.$$

(Note that

$$e^{At} = \begin{bmatrix} 1 & 0 \\ t & 1 \end{bmatrix}.$$
Try to verify that.) Note that this means that solutions through initial conditions of the form $(0, \alpha)$ are constant in time.

With this discussion out of the way, let us consider nonlinear systems

$$x' = f(x) \tag{13}$$

in $\mathbb{R}^n$ (you can always think about $n = 2$). It generates a flow $\phi$, i.e. $\phi(x^0, t)$ maps an initial condition $x^0$ to the solution of (13) $x(x^0, t)$ at time $t \in \mathbb{R}$. First, we need

**Definition 9.** $S \subset \mathbb{R}^n$ is a local invariant manifold for (13) if for all $x^0 \in S$ we have that $\phi(x^0, t) \in S$ for all $|t| < T$ for some number $T > 0$.

Now consider a rest point $\bar{x}$ and a little neighbourhood of it, $U$. Define

$$W_{\text{loc}}^{s}(\bar{x}) = \{x \in U \mid \phi(x, t) \to \bar{x} \text{ as } t \to \infty, \phi(x, t) \in U \ \forall t > 0\}.$$

This is the **local stable manifold** of $\bar{x}$. The local **unstable manifold** is defined similarly:

$$W_{\text{loc}}^{u}(\bar{x}) = \{x \in U \mid \phi(x, t) \to \bar{x} \text{ as } t \to -\infty, \phi(x, t) \in U \ \forall t < 0\}.$$

Convince yourself that indeed these sets are invariant. To show that they are manifolds is not that easy...

The local stable and unstable manifolds of a **hyperbolic** rest point $\bar{x}$ of the nonlinear system (13) are connected with stable and unstable subspaces of the origin of the corresponding linearized system,

$$y' = J(\bar{x})y. \tag{14}$$

by the

**Theorem 10 (Stable manifold theorem).** If the rest point $\bar{x}$ of (13) is hyperbolic, then the local unstable and stable manifolds $W_{\text{loc}}^{u}(\bar{x})$ and $W_{\text{loc}}^{s}(\bar{x})$ have the same dimensions as the stable and unstable manifolds of the linearized system (14), $E^s$ and $E^u$. Furthermore, $W_{\text{loc}}^{s}$ is tangent to $E^s$ and $W_{\text{loc}}^{u}$ is tangent to $E^u$ at $\bar{x}$.

Be sure you know how to make a drawing of this theorem! Finally the rigorous definition of the **global** stable and unstable manifolds is

$$W^s(\bar{x}) = \bigcup_{t \le 0} \phi(W_{\text{loc}}^s(\bar{x}), t),$$

$$W^u(\bar{x}) = \bigcup_{t \ge 0} \phi(W_{\text{loc}}^u(\bar{x}), t),$$

What this means is: to get $W^u$ integrate forward the points on the local unstable manifold $W_{\text{loc}}^u$, etc.
Let us have an example of $E^s, E^u$, and the stable/unstable manifolds. Consider

$$ 
\begin{aligned}
x'_{1} &= x_{1} \\
x'_{2} &= -x_{2} + x_{1}^{2}.
\end{aligned}
 $$

Not hard to see that the only rest point is $\bar{x} = (0, 0)$. Furthermore, the Jacobian matrix evaluated at the rest point is

$$ J(0, 0) = \begin{bmatrix} 1 & 0 \\ 0 & -1 \end{bmatrix} $$

with eigenvalues +1 and -1. $E^s$ of the linearized system is the $x_2$-axis, while $E^u$ is the $x_1$-axis. The nonlinear system can also be solved.

First of all note that if $x^0$ is of the form $(0, a)$ (that is, if the initial condition is on the $x_2$-axis), the solution is just $(0, ae^{-t})$. This means that $W^s(\bar{x})$ is the $x_2$-axis. (In this case it coincides with $E^s$). How about the unstable manifold (which must by the stable manifold theorem be one-dimensional, just like $E^u$)? Dividing the second equation by the first, we have

$$ \frac{dx_2}{dx_1} = -\frac{x_2}{x_1} + x_1, $$

solving which we have

$$ x_2(x_1) = \frac{x_1^2}{3} + \frac{C}{x_1}. $$

These curves are (unions of) orbits in the phase plane. **One of them has to be the unstable manifold.**

It has to pass through the origin and (by the stable manifold theorem again) be tangent to the unstable subspace $E^u$ there. The only candidate is the curve corresponding to $C = 0$. Thus the unstable manifold of the rest point at the origin is the set

$$ \{(x_1, x_2) \in \mathbb{R}^2, | x_2 = \frac{x_1^2}{3}\}. $$

What we have not yet done is to ask the following question: what is tangent to $E^c$ at a nonhyperbolic rest point? But this is precisely what we call a local centre manifold: a **local centre manifold** is a locally invariant manifold tangent to $E^c$ at $\bar{x}$. It is denoted by $W_{\text{loc}}^c(\bar{x})$.

# 8 Centre Manifolds

To consider properties of, and to be able to construct, a centre manifold, we shall first look at systems of equations in a special form. Below I will discuss everything for systems in $\mathbb{R}^n$, but of course you should keep in mind the case of $n = 2$. Suppose that (13) has a rest point at the origin, $\bar{x} = 0$, and that its variables $x$ can be decomposed into two parts, $y \in \mathbb{R}^k$ and $z \in \mathbb{R}^m$, $k + m = n$ (in our
applications we shall always have $k = m = 1$), so that (13) has the form

$$ 
\begin{aligned}
y' &= Ay + g_1(y, z) \\
z' &= Bz + g_2(y, z).
\end{aligned} 
\qquad (15) 
$$

Here $g_1(0, 0) = g_2(0, 0) = 0$; $J_1(0, 0)$ and $J_2(0, 0)$, linearizations of $g_1$ and $g_2$, are zero matrices of suitable order; all eigenvalues of A have zero real parts and all eigenvalues of B have negative real parts. Please note what we have done: we have moved the rest point to the origin, separated out the linear parts and diagonalized.

**Example.** This example will be our guinea pig. Please study it well and compare with the general setup.

$$ 
\begin{aligned}
y' &= yz \\
z' &= -z + \alpha y^2.
\end{aligned} 
\qquad (16) 
$$

Let us work out what is what here. $n = 2$, $k = m = 1$. $(0, 0)$ is a rest point. $A = 0$ (definitely a matrix with a zero eigenvalue); $B = -1$ (definitely a matrix with a negative eigenvalue); $g_1(y, z) = yz$. So $g(0, 0) = 0$, $J_1(y, z) = (z, y)$, a $1 \times 2$ matrix; $J_1(0, 0) = (0, 0)$. Work out $g_2$ and $J_2(y, z)$. This system (16) is thus in a suitable form for us to apply to it the centre manifold theory to understand the stability of the rest point $(0, 0)$, which is, of course, nonhyperbolic.

**Definition 11.** A local centre manifold is a local invariant manifold $z = h(y)$, such that $h(0) = 0$ and the graph of $h$ is tangent to $E^c$ at the origin.

In other words, on it the “stable” variables $z$ are functions of the “centre” variables $y$. The upshot of the theory is that near the rest point at the origin behaviour of the flow (and thus stability) of the rest point is determined by

$$ y' = Ay + g_1(y, h(y)). $$

This sometimes entails an enormous reduction of the dimension of the problem. In our case, say, once we know the centre manifold $z = h(y)$, all we have to do is just to investigate

$$ y' = yh(y), $$

a one-dimensional problem. Centre manifolds can be computed!

There are three theorems in the centre manifold theory.

**Theorem 12 (Existence).** If $f(\cdot)$ in (13) is $C^r$ ($r$ times differentiable in all its variables) then there exists a local centre manifold $z = h(y)$, $\|y\| \le \delta$; $h(\cdot)$ is $C^{r-1}$.

On the centre manifold, the flow is given by

$$ y' = Ay + g_1(y, h(y)). \qquad (17) $$
**Theorem 13 (Equivalence and convergence).** (a) Suppose that the zero solution of (17) is stable/asymp-
totically stable/unstable. Then the zero solution of (13) is stable/asymptotically stable/unstable.

(b) Suppose that the zero solution of (17) is locally asymptotically stable. Let $x(t) = (y(t), z(t))$ be a solution of (13) with $y(0), z(0)$ small enough. Then there exists a solution of (17), say $u(t)$, such that as $t \to \infty \lim \|y(t) - u(t)\| = 0$ and $\lim \|z(t) - h(u(t))\| = 0$.

(b) shows that if the origin is stable, the centre manifold attracts nearby solutions. (a) is the reason
we are doing the centre manifold reduction in the first place. And here is how we do it:

First, we need some notation. If we have a function $z = h(y)$, then its linearization at $y$ will be denoted $Dh(y)$. Now, consider (15). If $z = h(y)$, then $z' = Dh(y)y'$ (just the chain rule) and so we have

$$
Dh(y)y' = Bh(y) + g_2(y, h(y)).
$$

We also have

$$
y' = Ay + g_{1}(y, z).
$$

Thus

$$
Dh(y)(Ay + g_1(y, h(y))) = Bh(y) + g_2(y, h(y)). \quad (18)
$$

This equation, (18) and the additional conditions $h(0) = 0$ (the centre manifold must pass through the origin) and the fact it must be tangent to a subspace at the origin, define the centre manifold. That condition is $Dh(0) = 0_{m \times k}$, by which I mean a zero matrix. Remember that if $y \in \mathbb{R}^k, z \in \mathbb{R}^m$, $Dh(y)$ is a $m \times k$ matrix; if that bothers you, think about $k = m = 1$, as in our guinea pig example; then $Dh(y)$ is nothing else but $h'(y)$. Of course in general (18) cannot be solved. However, there is a fix. Take any function $u : \mathbb{R}^k \to \mathbb{R}^m$ defined in a neighbourhood of the origin and set

$$
(Mu)(y) := Du(y)(Ay + g_1(y, u(y))) - Bu(y) - g_2(y, u(y)). \quad (19)
$$

Note that $(M(h))(y) = 0$.

**Theorem 14 (Approximation).** Suppose $u$ is a $C^1$ function from $\mathbb{R}^k$ to $\mathbb{R}^m$, such that $u(0) = 0$, $Du(0) = 0_{m \times k}$. Suppose that as $\|y\| \to 0$, $(Mu)(y) = O(\|y\|^q)$, $q > 1$, then as $\|y\| \to 0$, $\|u(y) - h(y)\| = O(\|y\|^q)$.

This allows us to approximate $h(y)$ to any degree of accuracy. What this theorem says is this: if we manage to find a function $u(y)$, such that $(Mu)(y) = O(\|y\|^q)$, then this function approximates $h(y)$ to $O(\|y\|^q)$. For if $\|h(y) - u(y)\| = O(\|y\|^q)$, then we have that $h(y) = u(y) + O(\|y\|^q)$ (check the definition of $O(\cdot)!$). In a moment we shall see precisely how this works; you should realize that it is an algorithm that can be efficiently implemented in a computer algebra system.

Let us go back to our example,

$$
\begin{align}
y' &= yz \\
z' &= -z + \alpha y^2. \tag{20}
\end{align}
$$
Here the defining equation (18) is nothing but

$$h'(y)[yh(y)] + h(y) - \alpha y^2 = 0. \qquad (21)$$

Now we shall see how the approximation algorithm works. Assume that our approximation to the centre manifold is $u(y) = ay^2 + by^3 + cy^4$ and put this ansatz into (21). Then $u'(y) = 2ay + 3by^2 + 4cy^3$. Then we have from (19)

$$(Mu)(y) = ay^2 - \alpha y^2 + O(y^3).$$

Therefore, appealing to the approximation theorem, to approximate the centre manifold to $O(y^3)$ we take $a = \alpha$. That is all there is to it. Though this is sufficient for stability, let us go further on. If $a = \alpha$, we have from (19)

$$(Mu)(y) = by^3 + O(y^4).$$

Thus to approximate the centre manifold to $O(y^4)$ we take $b = 0$. Going on, we have

$$(Mu)(y) = 2\alpha^2 y^4 + cy^4 + O(y^5).$$

Thus we choose $c = -2\alpha^2$ and the approximation to the centre manifold obtained is

$$h(y) = \alpha y^2 - 2\alpha^2 y^4 + O(y^5).$$

The stability of (0, 0) of (20) is the same as the stability of $y = 0$ of

$$y' = yh(y) = \alpha y^3 + O(y^5).$$

Clearly, if $\alpha > 0$ the origin is unstable, while if $\alpha < 0$, it is locally asymptotically stable.

**Question:** What happen if $\alpha = 0$?

To summarize: we use the expression (19) and look for an approximate centre manifold of the form $u(y) = ay^2 +$ higher order terms, since we must have that $h(0) = 0$ and $h'(0) = 0$. Then we successively equate orders of $y$ to get better and better approximations to the centre manifold.

Let us do another example.

**Example.** Consider

$$y' = yz + \alpha y^3 + \beta z^2 y$$

$$z' = -z + \gamma y^2 + \delta y^2 z.$$

Let us investigate stability of the origin. There is a one-dimensional centre manifold $z = h(y)$ which satisfies

$$h'(y)[yh(y) + \alpha y^3 + \beta h^2(y)y] + h(y) - \gamma y^2 - \delta y^2 h(y) = 0. \qquad (22)$$

$h(0) = h'(0) = 0$. Assume the form $u(y) = ay^2 + by^3 + O(y^4)$ and plug into (19). Start equating orders of $y$:

$$(Mu)(y) = (a - \gamma)y^2 + O(y^3).$$
Thus if $u(y) = \gamma y^2$, then $(Mu)(y) = O(y^3)$ and thus $h(y) = \gamma y^2 + O(y^3)$. Again, this might be enough to establish stability: the flow on the centre manifold is

$$y' = (\alpha + \gamma)y^3 + O(y^4).$$

So if $\alpha + \gamma < 0$, the origin is stable; if $\alpha + \gamma > 0$, it is unstable, while if $\alpha + \gamma = 0$, we must go to higher orders of approximation (left as exercise).

**Exercise:** Consider

$$y' = y^2$$

$$z' = -z.$$

By direct integration, show that the centre manifold is not unique.

# 9 Polynomial Equations and the Newton Polygon

In this section we return to consider problems of the form

$$f(x, \epsilon) = 0,$$

where $f$ is a polynomial in $x$ and $\epsilon$. In short, the technique we use to solve such problems involves scaling to “regularize” the problem, and then expanding the solution in “correct” powers of $\epsilon$. To scale correctly, we use the device known as the *Newton polygon*. The scaled equation shows us exactly in what (fractional) powers of $\epsilon$ to expand. For more information on the Newton polygon, additional examples and a proof of the theorem below, see the book of Chow and Hale.

Thus, suppose we have a polynomial

$$f(x, \epsilon) = \sum_{k \in K} \alpha_k x^{i_k} \epsilon^{j_k}, \qquad (23)$$

where $K$ is a finite set of indices, $i_k, j_k$ are nonnegative integers and $\alpha_k$ are non-zero complex numbers for all $k \in K$.

The Newton polygon is constructed as follows. In the positive quadrant of the $(i, j)$ plane mark the points $P_k = (i_k, j_k)$ for $k \in K$. Then the Newton polygon of $f$ is the convex polygonal (piecewise linear) curve, such that all the points $P_k$ lie either on, or above it. Draw the Newton polygon of of the polynomial

$$f(x, \epsilon) = x^3 - \epsilon x + 2\epsilon^2 x^2 + \epsilon^4. \qquad (24)$$

Clearly, some of the points $P_k$ necessarily lie on the Newton polygon. These points are the more “important ones” in the following sense:

Consider a piece of the polygonal boundary defined by the line connecting any two points, $P_1 = (i_1, j_1)$ and $P_2 = (i_2, j_2)$, $i_2 > i_1$. Suppose that the slope of this line is $m$, that $i_2 - i_1 = n$ and that $P_1$ and $P_2$ are the only points on that segment of the boundary. Then we have
**Theorem 15.** 1. If $m = 0$ there are $n$ solutions of (23) of the form

$$x(\epsilon) = x_0 + y(\epsilon),$$

where $y(\epsilon)$ is analytic in $\epsilon$ in a neighbourhood of zero, $y(\epsilon) = O(1)$ as $\epsilon \to 0$, and $x_0$ is some solution of $f(x, 0) = 0$.

2. If $m \neq 0$ there are $n$ solutions of (23) of the form

$$x(\epsilon) = \epsilon^{-m}y(\epsilon),$$

where $y(\epsilon)$ is analytic in some fractional power of $\epsilon$ (to be determined from the scaled equation) in a neighbourhood of zero and $y(\epsilon) = O(1)$ as $\epsilon \to 0$.

The proof is based on the implicit function theorem. The assumption that $P_1$ and $P_2$ are the only points on the segment is sufficient (but not necessary) to ensure that the roots are simple and that the implicit function theorem can be used. It is not hard to see why this theorem should be true: if $m = 0$, we have a regular perturbation problem (so that the result is just an application of the implicit function theorem), otherwise such scaling of $x$ results in at least two terms (corresponding to $P_1$ and $P_2$) in the equation for $y$ having the same minimal power of $\epsilon$. Upon division by that power of $\epsilon$ we obtain a regular perturbation problem.

The case of multiple roots requires special treatment.

Let us consider (24). There are two segments on the Newton polygon: one connecting (0, 4) to (1, 1) and the other connecting (1, 1) to (3, 0). For the first segment we have $n = 1$ and $m = -3$.

Thus the correct scaling to obtain a solution is $x = \epsilon^3 y$, which gives us the equation

$$y^3 \epsilon^9 + 2 \epsilon^8 y^2 - \epsilon^4 y + \epsilon^4 = 0.$$

As expected, two last terms match. We divide by $\epsilon^4$ to get the final scaled equation

$$y^3 \epsilon^5 + 2 \epsilon^4 y^2 - y + 1 = 0.$$

This we can now expand in a regular perturbation expansion in powers of $\epsilon$:

$$y = 1 + \epsilon y_1 + \epsilon^2 y_2 + \dots$$

It is easy to see that $y_1 = y_2 = y_3 = 0$, so that we have obtained a representation for one of the roots (out of three) of (24):

$$x_1 = \epsilon^3 + O(\epsilon^7).$$

Now let us deal with the other segment. Here $n = 2$ and $m = -1/2$, So the correct scaling is $x = \epsilon^{1/2} y$, which gives

$$y^3 \epsilon^{3/2} + 2 \epsilon^3 y^2 - \epsilon^{3/2} y + \epsilon^4 = 0.$$
As expected, we can match the first and the third term. Dividing by $\epsilon^{3/2}$, we have

$$y^3 + 2\epsilon^{3/2}y^2 - y + \epsilon^{5/2} = 0.$$

which provides us with two more roots. Call them $v$ and $w$. They can be now determined from the expansions

$$v = 1 + \epsilon^{1/2}v_1 + \epsilon v_2 + \dots \text{ and}$$

$$w = -1 + \epsilon^{1/2}w_1 + \epsilon w_2 + \dots$$

Clearly $v_1 = v_2 = 0$ and same is true for $w_1, w_2$, so that we have the representation for the two other roots of (24)

$$x_2 = \epsilon^{1/2} + O(\epsilon^2) \text{ and}$$

$$x_3 = -\epsilon^{1/2} + O(\epsilon^2).$$

Note that if the coefficient of $\epsilon x$ in (24) were 1 instead of -1, $x_1$ would have been the only real root.

Now it is a good time to go back to the example of Lecture 3:

$$\epsilon y^2 - 3y + 2 = 0.$$

Here the Newton polygon has two segments, one connecting (0, 0) to (1, 0) and the other connecting (1, 0) to (2, 1). The first one has $m = 0$ and corresponds to the solution $y = 2/3 + O(\epsilon)$ and the other one has $m = 1$ and hence tells us to scale $y = z/\epsilon$, which is precisely what I did in Lecture 3.

Let us consider now the following quadratic polynomial:

$$(1 - \epsilon)x^2 - 2x + 1 = 0. \qquad (25)$$

Hinch devotes to it quite a lot of space. The problem here is that if $\epsilon = 0$, the equation has a double root $x = 1$.

Let us draw its Newton polygon. It is on the $x$ axis and there are three points on a straight line segment. The treatment in the case of a double root (which in general is uncovered after scaling) is very simple: translate to it and restart.

Thus in our case we introduce a new variable, $w = x - 1$. In the new variable the equation (25) is

$$w^2 - \epsilon w^2 - 2\epsilon w - \epsilon = 0,$$

whose Newton polygon has only one segment, $n = 2$, $m = -1/2$, so that we scale by $w = \epsilon^{1/2}y$ and obtain, on division by $\epsilon$ the scaled equation

$$y^2 - \epsilon y^2 - 2\epsilon^{1/2}y - 1 = 0,$$

the solutions of which can be expanded in powers of $\epsilon^{1/2}$. Thus we have that solutions of (25) have the representation

$$x_{1,2} = 1 \pm \epsilon^{1/2} + O(\epsilon).$$
# 10 The Curse of the Secular Terms

So far the only perturbation type problem we have looked at have been of a polynomial type and the Newton polygon gave us a way of de-singularising it and then using regular perturbation expansion. Let us start our treatment of perturbation techniques in **differential equations** and see spectacular ways in which regular perturbation expansion fails.

Before we do that, let us consider the difference between

$$y'' + y = \sin(2t) \quad \text{and} \quad y'' + y = \sin t.$$

The general solution of the first is

$$y(t) = A \sin(t) + B \cos(t) - 1/3 \sin(2t),$$

while that of the second is

$$y(t) = A \sin(t) + B \cos(t) - 1/2 t \cos(t).$$

Please note that the solution of the first one is $2\pi$-periodic, while that of the second one is not.

Let us consider the van der Pol equation,

$$y'' + \epsilon(y^2 - 1)y' + y = 0. \quad (26)$$

It is possible to show that for any positive $\epsilon$ this equation possesses a periodic solution, which is unique (up to shifts of phase) and globally attracting. This result must seem very puzzling: for $\epsilon = 0$, (26) reduces to the linear oscillator equation, which possess an infinite number of periodic solutions, none of which is asymptotically stable. The moment the non-linearity is “switched on” only one of them “survives” and acquires very strong stability properties.

We are interested in approximating this limit cycle for small $\epsilon$. Thus, let us choose our initial point somewhere on the limit cycle, taking, for example, $y'(0) = 0$ (note that we are not prescribing $y(0)$ at the same time), and crank the regular perturbation engine, that is, substitute the representation

$$y(\epsilon, t) \sim y_0(t) + \epsilon y_1(t) + \epsilon^2 y_2(t) + \dots \quad (27)$$

into (26). This results in the following sequence of equations (up to $O(\epsilon^3)$):

$$y_0'' + y_0 = 0,$$

$$y_1'' + y_1 = y_0'(1 - y_0^2),$$

$$y_2'' + y_2 = y_1'(1 - y_0^2) - 2y_0 y_0' y_1.$$

Solving the $O(1)$ equation for $y_0$, we obtain

$$y_0(t) = a_0 \cos t,$$
for some so far unknown $a_0$. Substituting this expression into the $O(\epsilon)$ equation, we see that $y_1$ satisfies

$$y_1'' + y_1 = \left(\frac{a_0^3}{4} - a_0\right) \sin t + \frac{a_0^3}{4} \sin 3t.$$

Here we have used the relation

$$\sin^3 t = \frac{3}{4}\sin t - \frac{1}{4}\sin 3t.$$

Consider this last equation. The problem with it is that it has a $\sin t$ term. Such a term, a multiple of one of the two linearly independent solutions of the homogeneous equation, will give rise to a term in the solution, which involves a power of $t$.

Indeed, its solution is

$$y_1(t) = -\left(\frac{a_0^3}{4} - a_0\right) \frac{t}{2} \cos t - \frac{a_0^3}{32} \sin 3t + a_1 \cos t + b_1 \sin t.$$

The secular term $t \cos t$ becomes unbounded as $t \to \infty$. The approximation

$$y(\epsilon, t) \approx y_0(t) + \epsilon y_1(t)$$

is definitely not a good one for large times if the coefficient of the $\epsilon t \cos t$ term is not zero. This term will become $O(1)$ when $t$ is $O(\frac{1}{\epsilon})$, so that the expansion (27) is no longer an asymptotic expansion at all. As we are looking for a **periodic** solution, we must require that all secular terms be removed. This leads directly to

$$a_0 = 2,$$

Thus our lowest order approximation to the limit cycle is

$$y_0(t) = 2 \cos t,$$

which is actually not bad, as we shall see later on.

The problem is that it just cannot be made better, and that in fact including more terms will only make matters worse.

Continuing, we use the initial condition $y'(0) = 0$, which means that $y_i'(0) = 0$ for all $i$, to obtain that $b_1 = \frac{3}{4}$, so that

$$y_1(t) = -\frac{1}{4}\sin 3t + \frac{3}{4}\sin t + a_1 \cos t.$$

Thus at this stage $a_1$ is not yet determined. Our hope would be that the solution of the equation for $y_2(t)$ will involve one secular term, which a judicious choice of $a_1$ will remove. Unfortunately, it is not so.

After some trigonometry, the $O(\epsilon^2)$ equation reads

$$y_2'' + y_2 = \frac{1}{4} \cos t + 2a_1 \sin t - \frac{3}{2} \cos 3t + 3a_1 \sin 3t + \frac{5}{4} \cos 5t,$$

from which we see that no choice of $a_1$ will get rid of both (resonant) $\cos t$ and $\sin t$ terms and hence secular terms will appear in $y_2(t)$.
As a different example of the failure of the regular perturbation technique, let us consider the simple
**linear** initial value problem

$$
y'' + 2\epsilon y' + (1 + \epsilon^2)y = 0, \quad y(0) = 0, \quad y'(0) = 1. \tag{28}
$$

This problem can be solved explicitly, and its solution is

$$
y(\epsilon, t) = e^{-\epsilon t} \sin t.
$$

Suppose we apply the regular perturbation procedure to this problem, Solving the resulting equa-
tions, we obtain

$$
y(\epsilon, t) \sim \sin t - \epsilon t \sin t + \frac{1}{2} \epsilon^2 t^2 \sin t + \dots \quad (29)
$$

Note again the appearance of secular terms. Do you notice something strange here? We will return
to this example later.

11 Poincaré-Linstedt Method

Above we saw that the regular perturbation expansion failed miserably for the van der Pol equation

$$
y'' + \epsilon(y^2 - 1)y' + y = 0.
$$

A way to overcome the difficulties encountered in our treatment of the van der Pol equation (26) is
the method due to Poincaré and further developed by Linstedt. It is also known as the *method of
strained coordinates*. To see the reasoning behind this method, we must understand why the regular
perturbation expansion technique failed for (26). Formally, that is, from a purely algebraic point of
view, we failed because there just were not enough arbitrary constants available to kill off all the
secular terms.

The deeper, physical, reason is uncovered as follows. Suppose we succeeded in applying the regular
perturbation technique. Then the resulting asymptotic solution **would have the same frequency**
ω = 1 for all values of ε small. This is the root of the problem: we should expect the non-linearity
to change the frequency of the oscillation. The basis of the Poincaré-Linstedt method is to allow
from the beginning for this variation in frequency.

Thus, we take the frequency ω also to have an asymptotic expansion in powers of ε,

$$
\omega \sim 1 + \epsilon k_1 + \epsilon^2 k_2 + \dots \tag{30}
$$

It is these undetermined coefficients $k_i$ that will allow us to get rid of the secular terms. We expand
the solution of (26) as

$$
y(\epsilon, t) \sim y_0(\omega t) + \epsilon y_1(\omega t) + \epsilon^2 y_2(\omega t) + \dots \quad (31)
$$
If we set $\tau = \omega t$ (this is the "strained coordinate"), (26) becomes

$$ \omega^2 \frac{d^2 y}{d\tau^2} + \epsilon\omega(y^2 - 1)\frac{dy}{d\tau} + y = 0. \quad (32) $$

Below we use now primes to denote differentiation with respect to $\tau$. If we substitute (30) and (31) into (32) and collect terms, we obtain (up to $O(\epsilon^3)$):

$$ y_0'' + y_0 = 0, \quad (33) $$

$$ y_1'' + y_1 = y_0'(1 - y_0^2) - 2k_1 y_0'' \quad (34) $$

$$ y_2'' + y_2 = y_1'(1 - y_0^2) - 2y_0 y_1' - 2k_1 y_1'' - (2k_2 + k_1^2)y_0'' + k_1(1 - y_0^2)y_0'. \quad (35) $$

Again, we must discuss initial conditions before we start solving the above equations. We take $y'(\epsilon, 0) = 0$ from which it follows that $y_i'(0) = 0, i = 0, 1, \dots$

Solving the $O(1)$ equation (33) for $y_0$ with $y_0'(0) = 0$, we obtain as before

$$ y_0(\tau) = a_0 \cos \tau, $$

for some so far unknown $a_0$.

Now, with this information the $O(\epsilon)$ equation (34) becomes

$$ y_1'' + y_1 = \left(-a_0 + \frac{1}{4}a_0^3\right) \sin \tau + 2k_1 a_0 \cos \tau + \frac{1}{4}a_0^3 \sin(3\tau). $$

For $y_1(\tau)$ to have no secular terms, we must, as before, put $a_0 = 2$ (to get rid of the $\sin \tau$ term), and also have in addition $k_1 = 0$ (to get rid of the $\cos \tau$ term). As in the regular perturbation procedure, this leads to

$$ y_1(\tau) = -\frac{1}{4}\sin 3\tau + \frac{3}{4}\sin \tau + a_1 \cos \tau. $$

A computation shows that the $O(\epsilon^2)$ equation (35) now reads

$$ y_2'' + y_2 = \left(4k_2 + \frac{1}{4}\right) \cos \tau + 2a_1 \sin \tau - \frac{3}{2} \cos 3\tau + 3a_1 \sin 3\tau + \frac{5}{4} \cos 5\tau, $$

from which it is clear that the secular terms can be removed if we take $a_1 = 0$ and $k_2 = -\frac{1}{16}$.

Thus, the dependence of the frequency on $\epsilon$ has the leading term behaviour of

$$ \omega \sim 1 - \frac{1}{16}\epsilon^2 + \dots $$

Thus, at this stage we have obtained the following (asymptotic) approximation to our solution:

$$ y(\epsilon, t) = 2 \cos \tau + \epsilon \left( -\frac{1}{4} \sin 3\tau + \frac{3}{4} \sin \tau \right) + O(\epsilon^2), $$

where $\tau = t\omega$ and

$$ \omega = 1 - \frac{1}{16}\epsilon^2 + O(\epsilon^3). $$
One can continue like this for as long as one likes. Of course, all this can be automated in Sage.

Let us think for a moment about the Poincaré-Linstedt method. It is a constructive method, the result of which is an asymptotic expansion, that we hope approximates a periodic solution. By itself it does not (a) prove that there is a periodic solution to approximate in the first place and (b) tell us anything about its stability. These are two topics we shall deal with briefly when we discuss averaging.

## 11.1 Shohat's expansion

Considering the definition of an asymptotic expansion, it is clear that the fixed number of terms approximations become progressively worse (though sometimes surprisingly slowly so) as $\epsilon$ grows. This leads to the question how the above technique can be augmented and adapted to deal with larger values of the parameter. In this section we present an idea in this direction.

The idea of Shohat (see the book by Bellman) is very simple: it is to introduce a different expansion parameter which is $O(\epsilon)$ for small $\epsilon$ and $O(1)$ as $\epsilon \to \infty$ and then to use the Poincaré-Linstedt method. We shall demonstrate the procedure using the van der Pol equation (26).

First we introduce our new expansion variable,

$$r = \frac{\epsilon}{1 + \epsilon}.$$

Now we assume the dependence of the new time variable $\tau$ on $t$ to be of the form

$$\tau = \omega t,$$

with

$$\omega \sim 1 + a_1 \epsilon + a_2 \epsilon^2 + \dots$$

As before, (26) takes the form

$$\omega^2 \frac{d^2 y}{d\tau^2} + \omega \epsilon (y^2 - 1) \frac{dy}{d\tau} + y = 0.$$

Now we crank again the perturbation engine, this time by setting

$$y(t, \epsilon) \sim y_0(\tau) + r y_1(\tau) + r^2 y_2(\tau) + \dots,$$

$$\epsilon = \frac{r}{1-r} \sim r + r^2 + r^3 + \dots,$$

so that

$$\epsilon\omega \sim r + c_2 r^2 + c_3 r^3 + \dots,$$

where $c_i$ are defined in terms of the coefficients $a_j$ of the asymptotic expansion of $\omega$ in $\epsilon$. Continuing as in the Poincaré-Linstedt technique, we obtain the following values for the first three $c_i$'s:

$$c_2 = 1, \quad c_3 = \frac{15}{16}, \quad c_4 = \frac{13}{16}.$$

The very good fit of the 4-term asymptotic solution with numerical simulations for various values of $\epsilon$, prompted Bellman to say that “it is tempting to conjecture that the Shohat series converges for all $\epsilon \ge 0$”.
# 12 The Method of Multiple Scales

Let us try to understand what went wrong in the example (28),

$$y'' + 2\epsilon y' + (1 + \epsilon^2)y = 0, \quad y(0) = 0, \quad y'(0) = 1,$$

whose exact solution was

$$y(\epsilon, t) = e^{-\epsilon t} \sin t.$$

The regular perturbation approach gave

$$y(\epsilon, t) \sim \sin t - \epsilon t \sin t + \frac{1}{2} \epsilon^2 t^2 \sin t + \dots \quad (36)$$

There was a question in class, whether we could know beforehand that the regular perturbation expansion will fail. A partial answer is like this: since the damping term (the coefficient of $y'$) is $O(\epsilon)$, the leading term of the regular perturbation expansion will be a $\sin t$ function, which completely misses the fact that there is a decay towards 0, and since all the consequent terms can only be either sines and cosines, or sines and cosines multiplied by powers of $t$, this decay will not be expressed in any finite series.

Note that (36) is nothing else but $\sin t$ multiplied by the Maclaurin series for $e^{-\epsilon t}$.

In (28) we have fast oscillation coexisting with the slow decay, so two different time scales, $t$ and $\epsilon t$. For equations with such structure, Cole introduced the method of multiple scales (also known as two-timing, two variable expansion method, though of course the number of time scales can be larger than two). To illustrate this idea, let us reconsider (28). We assume that our solution $y(t, \epsilon)$ is of the form

$$y(t, \epsilon) = Y(t, \eta),$$

where $\eta = \epsilon t$ is the “slow” time.

The chain rule gives us that

$$\frac{dy}{dt} = \frac{\partial Y}{\partial t} + \frac{\partial Y}{\partial \eta} \frac{d\eta}{dt} = \frac{\partial Y}{\partial t} + \epsilon \frac{\partial Y}{\partial \eta}$$

and similarly

$$\frac{d^2 y}{dt^2} = \frac{\partial^2 Y}{\partial t^2} + 2\epsilon \frac{\partial^2 Y}{\partial t \partial \eta} + \epsilon^2 \frac{\partial^2 Y}{\partial \eta^2}.$$

Note that implementing the multiple scales method involves dealing with *partial differential equations*, albeit simple ones.

With this information, our equation becomes

$$\frac{\partial^2 Y}{\partial t^2} + 2\epsilon \frac{\partial^2 Y}{\partial t \partial \eta} + \epsilon^2 \frac{\partial^2 Y}{\partial \eta^2} + 2\epsilon \left( \frac{\partial Y}{\partial t} + \epsilon \frac{\partial Y}{\partial \eta} \right) + (1 + \epsilon^2)Y = 0.$$
This can be rewritten as

$$ \frac{\partial^2 Y}{\partial t^2} + Y + 2\epsilon \left( \frac{\partial^2 Y}{\partial t \partial \eta} + \frac{\partial Y}{\partial t} \right) + \epsilon^2 \left( \frac{\partial^2 Y}{\partial \eta^2} + 2 \frac{\partial Y}{\partial \eta} + Y \right) = 0. $$

We also assume that our solution $Y(t, \eta)$ admits the asymptotic representation

$$ Y(t, \eta) \sim Y_0(t, \eta) + \epsilon Y_1(t, \eta) + \dots $$

In this, and in other examples below we shall determine only the leading term behaviour, as obtaining more terms means either introducing a strained coordinate or more time scales, both of which are cumbersome to implement. Of course, we shall see that in (28) the leading term of the asymptotic expansion is the exact solution!

We plug the asymptotic representation into the equation (28) and have that the $O(1)$ term is

$$ \frac{\partial^2 Y_0}{\partial t^2} + Y_0 = 0, \qquad (37) $$

which has as its general solution

$$ Y_0(t, \eta) = A(\eta) \sin t + B(\eta) \cos t. $$

Let us consider the initial conditions $y(0) = 0$, $y'(0) = 1$ of (28) in this setup. The first of these means that $Y_i(0, 0) = 0$ for all $i$, while the second means that

$$ \frac{\partial Y_0}{\partial t}(0, 0) = 1 $$

and

$$ \frac{\partial Y_i}{\partial \eta}(0, 0) + \frac{\partial Y_{i+1}}{\partial t}(0, 0) = 0 $$

for all $i \ge 0$.

Applying these initial conditions, we get that $B(0) = 0$, and that $A(0) = 1$. The $O(\epsilon)$ equation reads now

$$ \begin{aligned} \frac{\partial^2 Y_1}{\partial t^2} + Y_1 &= -2 \left( \frac{\partial^2 Y_0}{\partial t \partial \eta} + \frac{\partial Y_0}{\partial t} \right) \\ &= -2 \left( \frac{dA(\eta)}{d\eta} + A(\eta) \right) \cos t + 2 \left( \frac{dB(\eta)}{d\eta} + B(\eta) \right) \sin t. \end{aligned} \qquad (38) $$

To have solutions of (38) bounded for all time $t$, we must require the coefficients of $\sin t$ and $\cos t$ to be zero. This gives us the slow flow equations

$$ \frac{dA(\eta)}{d\eta} + A(\eta) = 0, \quad \frac{dB(\eta)}{d\eta} + B(\eta) = 0, $$

the solutions of which, subject to $A(0) = 1$ and $B(0) = 0$, are $A(\eta) = e^{-\eta}$ and $B(\eta) = 0$, so that $Y_0(t, \eta) = e^{-\eta} \sin t = e^{-\epsilon t} \sin t$, which is the exact solution. Note that as before we used solvability conditions in the $O(\epsilon)$ equation to determine parts of the solution of the $O(1)$ equation.
Let us apply this procedure (again to $O(\epsilon)$) to the van der Pol equation,

$$y'' + \epsilon(y^2 - 1)y' + y = 0.$$

Using the same notation as above, we get the same $O(1)$ equation,

$$\frac{\partial^2 Y_0}{\partial t^2} + Y_0 = 0,$$

the solution of which is

$$Y_0(t, \eta) = A(\eta) \cos t + B(\eta) \sin t.$$

The $O(\epsilon)$ equation is

$$\frac{\partial^2 Y_1}{\partial t^2} + Y_1 = -2 \frac{\partial^2 Y_0}{\partial t \partial \eta} + (1 - Y_0^2) \frac{\partial Y_0}{\partial t},$$

or, using the above form of $Y_0(t, \eta)$, we have

$$\frac{\partial^2 Y_1}{\partial t^2} + Y_1 = c_1 \sin t + c_2 \cos t + c_3 \sin 3t + c_4 \cos 3t,$$

where

$$c_1 = 2 \frac{dA}{d\eta} - A + \frac{A}{4}(A^2 + B^2),$$

$$c_2 = -2 \frac{dB}{d\eta} + B - \frac{B}{4}(A^2 + B^2),$$

$$c_3 = \frac{A}{4}(A^2 - 3B^2),$$

$$c_4 = \frac{B}{4}(B^2 - 3A^2).$$

Clearly, to get rid of the secular terms, we must have $c_1 = 0$, $c_2 = 0$, which is a system of two coupled ordinary differential equations, the slow flow equations.

In general, they are not easy to solve, but in this particular case passing to polar coordinates $A(\eta) = R(\eta) \cos \theta(\eta)$, $B(\eta) = R(\eta) \sin \theta(\eta)$, we obtain

$$\frac{dR}{d\eta} = \frac{R}{2} \left(1 - \frac{R^2}{4}\right),$$

$$\frac{d\theta}{d\eta} = 0,$$

so that the equations decouple. The second of these has the solution $\theta(\eta) = \theta_0$, and represents an arbitrary shift in phase, while the first one is readily shown to be

$$R(\eta) = \frac{2}{\sqrt{1 + \left(\frac{4}{R(0)^2} - 1\right) e^{-\eta}}}. \quad (39)$$

So

$$\begin{aligned} y(t, \epsilon) \sim & Y_0(t, \eta) + O(\epsilon) = R(\epsilon t) \cos(\theta_0) \cos t \\ & + R(\epsilon t) \sin(\theta_0) \sin t + O(\epsilon) = R(\epsilon t) \cos(t - \theta_0) + O(\epsilon), \end{aligned} \quad (40)$$
where $R(\cdot)$ is given by (39).

From (39) we can derive a lot of information. First, $\lim_{\eta \to \infty} R(\eta) = 2$ for any $R(0)$, which gives us a strong indication that the limit cycle of the van der Pol equation is stable. Also we have the rate of approach to the limit cycle.

(40) indicates that the period of the limit cycle is $2\pi$, so while this form of the method of multiple scales gives us the dynamics of the transient approach to the limit cycle, it is clear that we still have to stretch the fast time $t$ to obtain the correction to the frequency as we did in the Poincaré-Linstedt method.

## 12.1 Multiple Scales: a BVP

Now we will see how to use the idea of multiple scales in a boundary value problem.

Recall: we used the following approach to IVPs: we somehow deduce that our solution $y(t)$ can be written in the form

$$y(t, \epsilon) = Y(t, \eta),$$

where last time $\eta = \epsilon t$ was a “slow” time scale.

Then using the chain rule, we wrote

$$\frac{dy}{dt} = \frac{\partial Y}{\partial t} + \frac{\partial Y}{\partial \eta} \frac{d\eta}{dt} = \frac{\partial Y}{\partial t} + \epsilon \frac{\partial Y}{\partial \eta}$$

and

$$\frac{d^2 y}{dt^2} = \frac{\partial^2 Y}{\partial t^2} + 2\epsilon \frac{\partial^2 Y}{\partial t \partial \eta} + \epsilon^2 \frac{\partial^2 Y}{\partial \eta^2}.$$

We also assumed that our solution $Y(t, \eta)$ admits the asymptotic representation

$$Y(t, \eta) \sim Y_0(t, \eta) + \epsilon Y_1(t, \eta) + \dots$$

and we looked for the leading term of the expansion, $Y_0$.

Next we will deal with a BVP,

$$\epsilon y'' + y' + xy^2 = 0, \quad (41)$$

with the boundary conditions $y(0) = 0$ and $y(1) = 1$.

Here I am using $x$ as an independent variable, and primes denote differentiation with respect to $x$.

Note that this is a singular perturbation problem: putting $\epsilon = 0$ reduces the order of the equation and will make it impossible in general to meet both boundary conditions.

Before we can apply the method, we need some preparatory thinking. Let us solve the equation when $\epsilon = 0$. Then we have

$$y' + xy^2 = 0,$$
which is separable and gives as a general solution

$$y_o(x) = \frac{2}{x^2 + C}, \qquad (42)$$

a monotone decreasing function. If is clear that it impossible to choose C so that the condition $y(0) = 0$ is satisfied. On the other hand, the condition $y(1) = 1$ is satisfied if we choose $C = 1$.

The conclusion of this computation is that something drastic is happening close to $x = 0$; far from $x = 0$ $y_o(x)$ should give a good approximation for small $\epsilon > 0$.

So it looks like the additional (spatial) scale is a **fast scale**. That means that we want $\eta = x/\epsilon^\alpha$, with $\alpha > 0$. The easiest but not the quickest way to work out $\alpha$ is by trial and error.

We will discuss the question of scaling more when we talk about matched asymptotic expansions.

Some more notation: if there is a drastic change on a fast time scale we say that there is a **boundary layer**. If the layer is close to the boundary, as in our case, we simply say there is a **boundary layer**; if it is in the interior of the domain, it is an **interior boundary layer**.

Let us choose $\eta = \frac{x}{\epsilon}$. I will explain later why this is the right choice.

Thus we have

$$y(x, \epsilon) = Y(x, \eta), \quad \eta = x/\epsilon.$$

Using the chain rule, we have

$$\frac{dy}{dx} = \frac{\partial Y}{\partial x} + \frac{\partial Y}{\partial \eta} \frac{d\eta}{dx} = \frac{\partial Y}{\partial x} + \frac{1}{\epsilon} \frac{\partial Y}{\partial \eta}$$

and

$$\frac{d^2 y}{dx^2} = \frac{\partial^2 Y}{\partial x^2} + \frac{2}{\epsilon} \frac{\partial^2 Y}{\partial x \partial \eta} + \frac{1}{\epsilon^2} \frac{\partial^2 Y}{\partial \eta^2}.$$

Finally, we assume the expansion

$$Y(x, \eta) \sim Y_0(x, \eta) + \epsilon Y_1(x, \eta) + \dots, \qquad (43)$$

and again we will only find $Y_0(x, \eta)$.

In the new variables, the equation itself becomes

$$\frac{1}{\epsilon} \left( \frac{\partial^2 Y}{\partial \eta^2} + \frac{\partial Y}{\partial \eta} \right) + \left( 2 \frac{\partial^2 Y}{\partial x \partial \eta} + \frac{\partial Y}{\partial x} + x Y^2 \right) + \epsilon \frac{\partial^2 Y}{\partial x^2} = 0.$$

Using expansion (43), we have

$$\frac{1}{\epsilon} \left( \frac{\partial^2 Y_0}{\partial \eta^2} + \frac{\partial Y_0}{\partial \eta} \right) + \left( \frac{\partial^2 Y_1}{\partial \eta^2} + \frac{\partial Y_1}{\partial \eta} + 2 \frac{\partial^2 Y_0}{\partial x \partial \eta} + \frac{\partial Y_0}{\partial x} + x Y_0^2 \right) = O(\epsilon).$$

So now we consecutively put the coefficients of different orders of $\epsilon$ to zero, starting with the $O(1/\epsilon)$ term.
At $O(1/\epsilon)$, we have

$$ \frac{\partial^2 Y_0}{\partial \eta^2} + \frac{\partial Y_0}{\partial \eta} = 0, $$

from which we conclude that

$$ Y_0(x, \eta) = A(x) + B(x)e^{-\eta}. $$

Note that even though we are solving a partial differential equation, it involves only one variable, so all “constants” will be really functions of the slow variable $x$.

Moving on to the $O(1)$ terms, we need to compute $\frac{\partial Y_0}{\partial x}$ and $\frac{\partial^2 Y_0}{\partial x \partial \eta}$. These are

$$ \frac{\partial Y_0}{\partial x} = A'(x) + B'(x)e^{-\eta} \quad \text{and} \quad \frac{\partial^2 Y_0}{\partial x \partial \eta} = -B'(x)e^{-\eta}. $$

So the $Y_1$ equation becomes

$$ \frac{\partial^2 Y_1}{\partial \eta^2} + \frac{\partial Y_1}{\partial \eta} = -xB^2e^{-2\eta} - (2xAB - B')e^{-\eta} - (A' + xA^2). $$

The resonant terms (in $\eta$) are the “constants” $A' + xA^2$ and $(2xAB - B')e^{-\eta}$. If the first is not zero, it will give rise in $Y_1$ to a term proportional to $\eta$, i.e. of $O(1/\epsilon)$; the second term will give rise to a term proportional to $\eta e^{-\eta}$. This does not look too bad, but if it so happens that $A(x) \equiv 0$, then we could have $Y_1 = C(x)\eta e^{-\eta}$ and hence in

$$ Y(x, \eta) = B(x)e^{-\eta} + \epsilon C(x)\eta e^{-\eta} + O(\epsilon^2), $$

when $x = O(1)$, the second term is of the same order as the first, which destroys asymptoticity.

The upshot of this not completely satisfactory discussion is that we are advised to put the resonant terms to zero. Thus we have the two differential equations,

$$ \frac{dA}{dx} = -xA^2 \quad \text{and} \quad \frac{dB}{dx} = 2xAB. $$

From this we get that $A(x) = \frac{2}{a+x^2}$, arguing as we did for $y_0$, and since $dA/dB = -A/2B$, $B = b/A^2$ for some constants $a, b$. Hence

$$ B(x) = c(a + x^2)^2, \quad c = b/4. $$

So

$$ y(x, \epsilon) \sim \frac{2}{a+x^2} + c(a+x^2)^2 e^{-\eta} + O(\epsilon). $$

Now we need to work out the constants $a$ and $c$ using the boundary conditions.

First of all, $y(0, \epsilon) = 0$ leads to $c = -2a^{-3}$. Secondly, at $x = 1$ $e^{-\eta}$ is vanishingly small, so we put that term to zero. Hence at $x = 1$, $y(1, \epsilon) \approx \frac{2}{a+1} = 1$, so that $a = 1$, and our final result is that

$$ y(x, \epsilon) \sim \frac{2}{1+x^2} - 2(1+x^2)e^{-x/\epsilon} + O(\epsilon). $$
It is worthwhile to see how well we have satisfied the boundary conditions for $x = 1$. If $\epsilon = 0.1$, $\frac{2}{1+x^2} - 2(1+x^2)e^{-x/\epsilon}$ at $x = 1$ is 0.9996368006, and if $\epsilon = 0.01$, it is 1 to 8 digit precision.

Finally, it is a good idea to see a graph of this approximate solution.

<figure><img src="image_2.png" /></figure>

It shows that the solution is a splicing together of two different bits, an “inner solution” that satisfies the boundary condition at $x = 0$, changes $O(1)$ on the scale of $\epsilon$ (using the “inner variable” $x/\epsilon$) and an “outer solution” that is a good approximation outside of the boundary layer. These ideas are the building blocks of matched asymptotic expansions that we are going to explore next.

## 13 Matched Asymptotic Expansions

Let us start with an example. It still should be a mystery how we locate the boundary layer or how we know how to choose “inner variables”, but it is best first to see the method at work. Again, we will only be computing the first term of the asymptotic expansion, though unlike the case of the method of multiple scales, this method can be iterated without much difficulty.

We will do an example due to Nayfeh; it is very popular!

Consider

$$ \epsilon y'' + y' = 2x + 1, \quad x \in (0, 1), \quad y(0) = 1, \quad y(1) = 4. \tag{44} $$

First of all let us work out its exact solution. It is

$$ y(x) = x^2 + x - 2x\epsilon + \frac{(1 + 2\epsilon) \exp(-x/\epsilon)}{\exp(-1/\epsilon) - 1} - \frac{(2 + 2\epsilon - \exp(-x/\epsilon))}{\exp(-1/\epsilon) - 1}. $$

which for small $\epsilon$ (so that $\exp(-1/\epsilon)$ is negligible), simplifies to

$$ y(x) = x^2 + x + 2 - \exp(-x/\epsilon) + \epsilon(2 - 2x - 2\exp(-x/\epsilon)) + \text{EST}, $$

where EST stands for exponentially small terms. Note that boundary conditions are satisfied up to EST.

Now let us do a regular perturbation expansion. Let us “assume” that there is a boundary layer at $x = 0$. So we will take our regular expansion solution to satisfy the condition $y(1) = 4$. Since we are expanding $y(x) \sim y_0(x) + \epsilon y_1(x) + \dots$, this means that $y_0(1) = 4$ and $y_j(1) = 0$ for all $j \ge 1$.
So we have $y_0' = 2x + 1$, i.e. $y_0(x) = x^2 + x + a$, and the boundary condition gives us that $a = 2$. Continuing in the same way, we get that $y_1' = -y_0'' = -2$, $y_1(1) = 0$, so that $y_1(x) = 2(1-x)$. Hence

$$y_0(x) \sim x^2 + x + 2 + 2\epsilon(1-x) + O(\epsilon^2),$$

where I called this expansion the “outer expansion”, hence the subscript o. It is clear that I can continue like that as long as I wish (though in this particular problem there is no need).

Now let us introduce the “inner” variable, $X = x/\epsilon$. I will explain later how I know how to scale correctly. Let us set $y(x) = Y(X)$. By the chain rule, (44) becomes

$$\frac{1}{\epsilon} \left( \frac{d^2 Y}{dX^2} + \frac{dY}{dX} \right) = 2\epsilon X + 1. \qquad (45)$$

Now let us expand $Y(X)$ in powers of $\epsilon$, $Y(X) \sim Y_0(X) + \epsilon Y_1(X) + \dots$. The boundary conditions are $Y_0(0) = 1$ and $Y_j(0) = 0$ for $j \ge 1$. So for $Y_0$ we have

$$\frac{d^2 Y_0}{dX^2} + \frac{dY_0}{dX} = 0, \quad Y_0(0) = 1.$$

which is a second order equation, but we only have one boundary condition, which turns out to be magic. We obtain

$$Y_0(X) = A + (1 - A)e^{-X}.$$

Let us compute the next term. $Y_1$ satisfies the equation

$$\frac{d^2 Y_1}{dX^2} + \frac{dY_1}{dX} = 1, \quad Y_1(0) = 0,$$

so that $Y_1(X) = B(1 - e^{-X}) + X$.

Hence the “inner expansion”, $y_i(x) = Y(X)$ is given by

$$y_i(x) \sim A + (1 - A)e^{-X} + \epsilon(B(1 - e^{-X}) + X) + O(\epsilon^2).$$

[It is best to write $y_i(x)$ at this stage in terms of $X$ as we will see next.]

The outer expansion is valid for $0 < x < 1$ and the inner one is valid for $X \ge 0$. Now we have a bit of witchcraft called matching conditions. They were originally developed by Prandtl and then much discussed by Lagerstrom. So let us do the **Prandtl matching condition** which is only good for leading expansion terms.

So we have: $y_i(x) \sim A + (1 - A)e^{-X} + O(\epsilon) = Y_0(X) + O(\epsilon)$ and $y_0(x) \sim x^2 + x + 2 + O(\epsilon) = y_0(x) + O(\epsilon)$. The Prandtl matching condition says that

$$\lim_{X \to \infty} Y_0(X) = \lim_{x \to 0} y_0(x). \qquad (46)$$

So from matching we have $A = 2$.
We call this limit “the common part”. And then the recipe for a combined leading term which is valid for all $x \in [0, 1]$ is

$$y_c(x) \sim y_0(x) + Y_0(X) - \text{common part} + O(\epsilon) = x^2 + x + 2 - e^{-x/\epsilon} + O(\epsilon).$$

This trick will not work if we want higher order matching. We will see next time how to matching in a more reasoned way. Meanwhile let us check: for $\epsilon = 0.05$, $y_c(0) = 1$, $y_c(1) = 3.999999998$ and let us plot the various different leading terms:

<figure><img src="image_3.png" /></figure>

## 13.1 MAE: Higher Order Matching

To summarise: we considered the singularly perturbed boundary value problem,

$$\epsilon y'' + y' = 2x + 1, x \in (0, 1), y(0) = 1, y(1) = 4.$$

We found the outer expansion,

$$y_o(x) \sim x^2 + x + 2 + 2\epsilon(1-x) + O(\epsilon^2),$$

and the inner expansion

$$y_i(x) \sim A + (1-A)e^{-X} + \epsilon(B(1-e^{-X}) + X) + O(\epsilon^2).$$

Finally, we used Prandtl's matching principle to match the leading terms of the two expansions to produce the leading term of the combined (matched) expansion,

$$y_c(x) \sim x^2 + x + 2 - e^{-x/\epsilon} + O(\epsilon).$$

The (Lagerstrom) idea of higher order matching is more sophisticated. Let us see how it works. Let us rewrite the outer expansion:

$$y_o(x) \sim x^2 + x + 2 + 2\epsilon(1-x) + O(\epsilon^2),$$

and the inner expansion

$$y_i(x) \sim A + (1-A)e^{-X} + \epsilon(B(1-e^{-X}) + X) + O(\epsilon^2).$$
Now set $x = t\epsilon^\alpha$, $\alpha \in (0,1)$. Let us rewrite the outer expansion in the new variable $t$:

$$y_o \sim \epsilon^{2\alpha}t^2 + \epsilon^\alpha t + 2 + 2\epsilon - 2t\epsilon^{1+\alpha} + O(\epsilon^2).$$

If $x = t\epsilon^\alpha$, then $X = x/\epsilon = t/\epsilon^{1-\alpha}$. Let us rewrite the inner expansion $y_i$ in terms of $t$:

$$y_i \sim A + (1-A)e^{-t/\epsilon^{1-\alpha}} + \epsilon B - \epsilon B e^{-t/\epsilon^{1-\alpha}} + t\epsilon^\alpha + O(\epsilon^2).$$

Let us now choose $t = O(1)$. Then the exponential terms are vanishingly small, and what remains is

$$y_i \sim A + \epsilon B + t\epsilon^\alpha + O(\epsilon^2).$$

Let us choose $\alpha > 1/2$, so that $2\alpha > 1$. We want the inner and the outer expansions to match at this intermediate scale. Now we are ready to match.

$$y_o \sim 2 + \epsilon^\alpha t + 2\epsilon + o(\epsilon)$$

should match

$$y_i \sim A + \epsilon^\alpha t + B\epsilon + o(\epsilon).$$

So we recover $A = 2$ and $B = 2$. The common part now is $2 + x + 2\epsilon$. As before, we form the combined expansion (now to $O(\epsilon^2)$):

$$\begin{aligned} y_c(x) & \sim y_o(x) + y_i(x) - \text{common part} + O(\epsilon^2) \\ &= x^2 + x + 2 + 2\epsilon(1-x) + 2 - e^{-x/\epsilon} + \epsilon(2 - 2e^{-x/\epsilon} \\ & \quad + x/\epsilon) - (2+x+2\epsilon) + O(\epsilon^2) \\ &= x^2 + x + 2 - e^{-x/\epsilon} + \epsilon(2 - 2x - 2e^{-x/\epsilon}) + O(\epsilon^2), \end{aligned}$$

which matches the exact solution.

## 13.2 Scaling boundary layers

Let us go back to our original problem,

$$\epsilon y'' + y' = 2x + 1, \quad x \in (0,1), \quad y(0) = 1, \quad y(1) = 4.$$

Let us assume that the inner variable for the boundary layer close to $x = 0$ is $X = x/\epsilon^\nu$. Then in the new variables, putting $y(x) = Y(X)$, we have

$$\epsilon^{1-2\nu} \frac{d^2 Y}{dX^2} + \epsilon^{-\nu} \frac{dY}{dX} = 2\epsilon^{\nu} X + 1.$$

We want to choose $\nu$ is such a way, that after algebraic simplification, the term involving the highest derivative term can be matched with another term, such that all the rest of the terms are of higher
order in $\epsilon$. Choosing $\nu = 1/2$, which would match the second derivative with 1, does not work, because then the leading term will be the $\epsilon^{-\nu} \frac{dY}{dX}$ term.

Choosing $\nu = 1/3$, matching the second derivative with the $2\epsilon^\nu X$ term does not work either, as the first derivative term is still the dominant one. On the other hand, making the choice $\nu = 1$ and multiplying the equation by $\epsilon$ gives

$$ \frac{d^2 Y}{dX^2} + \frac{dY}{dX} = \epsilon(2\epsilon X + 1), $$

to which the regular perturbation procedure can be applied.

Let us consider a very nice problem due to N. Trefethen.

$$ \epsilon y'' + y(2 - x^2) = 1, \quad x \in (-1, 1), \quad y(-1) = y(1) = 0. $$

Then you should convince yourself that there are boundary layers at $-1$ and $1$; in fact you can draw the solution for small $\epsilon$ without solving anything!

Let us work close to $x = -1$. We assume that the inner variable for the boundary layer close to $x = -1$ is $X = (1+x)/\epsilon^\nu$. So $x = -1 + X\epsilon^\nu$. Then we have

$$ \epsilon^{1-2\nu} \frac{d^2 Y}{dX^2} + Y - 1 + 2YX\epsilon^\nu - X^2 Y\epsilon^{2\nu} = 0. $$

Clearly, here the correct scaling is $\nu = 1/2$, as then the first 3 terms are matched and the rest are of higher order. The same scaling is also the correct one for the boundary layer at $x = 1$, though of course the correct definition of the inner variable now is $X = (1-x)/\epsilon^{1/2}$.

## 13.3 Locating boundary layers

Next we need to deal with locating boundary layers. In general it is not an easy problem, but for problems of the form

$$ \epsilon y'' + f(x)y' + g(x)y = 0, $$

with $x \in I = (a, b)$, $x(a) = \alpha$, $x(b) = \beta$, and $f(x)$ and $g(x)$ continuous, there are well-established rules. I will state them and then will discuss in detail one case, to explain the reasoning behind them.

The rules all depend on the behaviour of the function $f(x)$, the reason for which will become clear.

### Boundary Layer Location Rules

1. If $f(x) > 0$ on $I$, the boundary layer is at $x = a$.

2. If $f(x) < 0$ on $I$, then boundary layer is at $x = b$.

3. If there is a (unique) point $x_0 \in I$ such that $f(x_0) = 0$, then there are 2 cases: (a) if $f(a) < 0$ and $f(b) > 0$, then there is only an interior boundary layer at $x = x_0$, and (b) if $f(a) > 0$ and $f(b) < 0$, then there are boundary layers at $x = a$ and $x = b$ and a possible internal boundary layer at $x = x_0$.
Assuming for a second that rules 1 and 2 hold, let me explain why 3 makes sense. I will just give the reasoning for 3(a). If $f(x)$ only has a zero at $x_0$ and $f(a) < 0$, it means that in the interval $(a, x_0)$, $f(x) < 0$. Let us consider our problem on that interval. Rule 2 says that there will be no boundary layer at $a$ and a boundary layer at $x_0$. Applying the same logic to the interval $(x_0, b)$ and using Rule 1, we see that there is no boundary layer at $b$, but a boundary layer at $x_0$, i.e. an internal boundary layer.

Let us see an example showing why Rule 2 is correct. Consider on $I = (0, 1)$

$$ \epsilon y'' + (x-2)y' + y = 0, \quad y(0) = 3, \quad y(1) = 2. $$

First of all there must be a boundary layer somewhere as putting $\epsilon$ to zero reduces the order of the equation and in general (and in this case) would not allow us to satisfy both boundary conditions. By Rule 2, there will be a boundary layer at $x = 1$, since $x - 2 < 0$ on $I$. Let us show that there cannot be boundary layer at $x = 0$.

Suppose there is, so that the outer solution will satisfy the boundary condition at $x = 1$. Hence its leading term is the solution of

$$ (x-2)y'_{0} + y_{0} = 0, \quad y_{0}(1) = 2, $$

so that $y_0(x) = \frac{2}{2-x}$. Therefore $y_0 = \frac{2}{2-x} + O(\epsilon)$.

As you can check, the correct boundary layer scaling gives the inner variable $X = x/\epsilon$. In the inner variables, setting $y(x) = Y(X)$, we get

$$ \frac{d^2 Y}{dX^2} + (\epsilon X - 2) \frac{dY}{dX} + \epsilon Y = 0, $$

and we have to satisfy $Y(0) = 3$. The leading term of the inner solution thus is $Y_0(X) = 3 - A + Ae^{2X}$, so that the inner solution is $Y(X) = Y_0(X) + O(\epsilon)$. Now apply the Prandtl matching principle: to match we must have

$$ \lim_{x \to 0} \frac{2}{2-x} = \lim_{X \to \infty} (3 - A + Ae^{2X}). $$

But the right-hand limit only exists if $A = 0$, in which case we have $1 = 3$, which shows that matching is impossible. It is the negative sign of $f(x) = (x-2)$ that led to the exponential ($e^{2X}$ here) term that made matching impossible.

## 13.4 A problem with an interior boundary layer

We will finish this section by doing an example with an internal boundary layer to see what is involved.

Consider the boundary value problem

$$ \epsilon y'' + xy' - xy = 0, \quad x \in (-1, 1), \quad y(-1) = 1, \quad y(1) = 1. $$
By Rule 3(a) there are no boundary layers but there will be an internal boundary layer at $x = 0$. There will be two outer expansions, one for $x < 0$ and one for $x > 0$. For $x > 0$ the leading term of the outer expansion satisfies

$$y' - y = 0, y(1) = 1,$$

which gives $y_0^+(x) = e^{x-1}$. Similarly, the leading term of the outer expansion satisfies the same equation with the condition $y(-1) = 1$, which leads to $y_0^-(x) = e^{x+1}$.

Now it is time to define the inner expansion. If we set as usual $X = x/\epsilon^\nu$ and let $y(x) = Y(X)$, the equation becomes

$$\epsilon^{1-2\nu} \frac{d^2 Y}{dX^2} + X \frac{dY}{dX} - \epsilon^\nu XY = 0.$$

Hence the right scaling will be $\nu = 1/2$, i.e. the inner variable is $X = x/\sqrt{\epsilon}$, and the leading term of the inner expansion satisfies

$$\frac{d^2 Y_0}{dX^2} + X \frac{dY_0}{dX} = 0.$$

This a separable first order equation for $dY_0/dX$, which gives

$$\frac{dY_0}{dX} = a_0 e^{-X^2/2},$$

so that

$$Y_0(X) = b_0 + a_0 \int_0^X e^{-u^2/2} du.$$

Note that there are no boundary conditions on $Y_0(X)$, and that we have generated two constants of integration which will be used to match with the two outer expansions using the Prandtl matching principle.

We have

$$\lim_{x \to 0^+} y_0^+(x) = e^{-1} = b_0 + a_0 \int_0^\infty e^{-u^2/2} du = \lim_{X \to \infty} Y_0(X).$$

and similarly

$$\lim_{x \to 0^-} y_0^-(x) = e^1 = b_0 + a_0 \int_0^{-\infty} e^{-u^2/2} du = \lim_{X \to -\infty} Y_0(X).$$

Note that

$$\int_0^\infty e^{-u^2/2} du = \sqrt{\frac{\pi}{2}} = -\int_0^{-\infty} e^{-u^2/2} du,$$

This means that we have a system of two linear equations for $a_0$ and $b_0$:

$$b_0 + a_0 \sqrt{\frac{\pi}{2}} = e^{-1}$$

$$b_0 - a_0 \sqrt{\frac{\pi}{2}} = e^1.$$

Solving them, we have that

$$Y_0(X) = \frac{1}{2}(e^{-1} + e^1) + \frac{1}{\sqrt{2\pi}}(e^{-1} - e^1) \int_0^X e^{-u^2/2} du.$$
Finally, we combine the outer and the inner expansions to get the leading term of the expansion,

$$y_c(x) = \begin{cases} y_0^-(x) + Y_0(X) - e^1 & \text{if } x < 0 \\ y_0^+(x) + Y_0(X) - e^{-1} & \text{if } x > 0, \end{cases}$$

where I have subtracted the respective common parts from the expansions.

Check that this function is continuous; you might want to see whether it is $C^1$. The positive part of it looks like this:

<figure><img src="image_4.png" /></figure>

# 14 The WKB Method

## 14.1 The main idea

Let us consider another very spectacular failure of the regular perturbation procedure. Consider the very simple boundary value problem,

$$\epsilon y'' + y = 0, \quad y(0) = 0, y(1) = 1.$$

If you try the RPE, you will get identically zero, which is not very informative. On the other hand, the exact solution is

$$y_{ex}(x) = \frac{\sin(x/\sqrt{\epsilon})}{\sin(1/\sqrt{\epsilon})}.$$

So not only is there no solution when $\epsilon = 1/(n^2\pi^2)$ (try to work out why), but for small $\epsilon$ it oscillates very fast. No amount of scaling will help here, and a new idea is needed.

A very powerful approach was developed by Wentzel, Kramers and Brillouin in 1920s, hence WKB, but with contributions from Jeffreys, hence sometimes WKBJ.

It is applicable to boundary value problems of the form

$$y'' + \lambda^2 f(x)y = 0, \quad x \in I = (a, b), \quad y(a) = \alpha, y(b) = \beta, \quad (47)$$

where $\lambda$ is very large. So $\lambda^2 = 1/\epsilon$. We will construct asymptotic expansions in $1/\lambda$, really only the $O(\lambda)$ and the $O(1)$ terms.

It might seem that the form (47) is very restrictive. First of all, the method can handle equations of the form

$$y'' + f(x, \lambda)y = 0, \quad (48)$$
and equations in the more general form

$$y'' + g(x, \lambda)y' + f(x, \lambda)y = 0$$

can be reduced to the above form as follows:

Set $y = uz$. Then

$$y' = u'z + z'u; \quad y'' = u''z + 2u'z' + uz''$$

Plug this into the equation and divide by $u$. We get

$$z'' + \left(\frac{2u'}{u} + g\right) z' + \left(\frac{u''}{u} + g\frac{u'}{u} + f\right) z = 0.$$

Now choose $u$ so that $2u'/u + g = 0$, i.e. $\log u = -1/2 \int g dx$, and we have the equation for $z$,

$$z'' + \left(\frac{u''}{u} - \frac{g^2}{2} + f\right) z = 0,$$

from which $y$ can be recovered.

The main idea of the method is very simple. The way we solve second order linear boundary value problems is first to find two linearly independent solutions and then to match the boundary conditions. The way the WKB method goes about finding these linearly independent solutions is best motivated by considering the linear problem,

$$y'' + K\lambda^2 y = 0,$$

where $K > 0$ is a constant. Then in complex notation the general solution is

$$y = A e^{i\lambda\sqrt{K}x} + B e^{-i\lambda\sqrt{K}x}.$$

For negative $K$ we have

$$y = A e^{\lambda\sqrt{-K}x} + B e^{-\lambda\sqrt{-K}x}.$$

In both cases (the first one oscillatory), we have that the two linearly independent solutions are of the form $e^{\lambda 0}$, which is what led W, K, and B to suggest that we look for solutions of (47) or (48) in the form

$$y = e^{\lambda Y(x, \lambda)}.$$

So here we reduce the problem to something we know how to handle by a nonlinear change of variable. Let us rewrite (47) using $Y$.

$$y' = \lambda Y' e^{\lambda Y(x, \lambda)}, \quad y'' = \lambda^2 (Y')^2 e^{\lambda Y(x, \lambda)} + \lambda Y'' e^{\lambda Y(x, \lambda)}.$$

Substitute this into (47), divide by $\lambda^2$, cancel out the exponents, and we have

$$\frac{1}{\lambda} Y'' + (Y')^2 + f(x) = 0. \qquad (49)$$
We are interested in large $\lambda$, so in small $1/\lambda$. Now strangely enough (due to the square!) we are able to recover the two required linearly independent solutions simply by applying the regular perturbation procedure to (49).

We write

$$Y = Y_0 + \frac{1}{\lambda}Y_1 + O\left(\frac{1}{\lambda^2}\right)$$

Then we have that $Y_0$ satisfies the relation

$$(Y_0')^2 = -f(x).$$

From now on we stipulate that $f(x)$ is of a definite sign on $I$, either positive or negative. The problems in which “turning points”, points where $f$ changes sign, occur are hard and need matched asymptotics; they are the next stage if you are interested in this topic. So if $f(x)$ is positive, we have two possibilities:

$$Y_0' = i(f(x))^{1/2} \text{ and } Y_0' = -i(f(x))^{1/2};$$

while for $f(x) < 0$ we have

$$Y_0' = (-f(x))^{1/2} \text{ and } Y_0' = -(-f(x))^{1/2};$$

So the two possibilities for $Y_0$ are:

$$Y_0 = \begin{cases} \pm i \int (f(x))^{1/2} dx & \text{if } f(x) > 0 \\ \pm \int (-f(x))^{1/2} dx & \text{if } f(x) < 0, \end{cases}$$

and in the first case we have as the leading term approximation

$$y(x) = A e^{i\lambda \int (f(x))^{1/2} dx} + B e^{-i\lambda \int (f(x))^{1/2} dx},$$

while in the second one we have

$$y(x) = A e^{\lambda \int (-f(x))^{1/2} dx} + B e^{-\lambda \int (-f(x))^{1/2} dx},$$

Note that all integration constants that one generates will be absorbed into $A$ and $B$, and we will not mention them again.

It turns out that $Y_1$ is quite easy to compute too. It satisfies the equation

$$Y_0'' + 2Y_0'Y_1' = 0, \text{ so that } Y_1' = -\frac{Y_0''}{2Y_0'}.$$

Let us consider the case of positive and negative $f$ separately (I will no longer indicate its dependence on $x$). So if $f > 0$, $Y_0' = \pm if^{1/2}$. Hence

$$Y_0'' = \pm \frac{i}{2} f^{-1/2} f',$$

so that

$$Y_1' = -\frac{Y_0''}{2Y_0'} = -\frac{1}{4} \frac{f'}{f}.$$
Hence for both possibilities of $Y_0$, we have that

$$Y_1 = \log f^{-1/4}.$$

The same analysis for $f < 0$ gives for both plus and minus versions of $Y_0$ that

$$Y_1 = \log(-f)^{-1/4}.$$

Hence since $y = e^{\lambda Y_0 + Y_1 + O(1/\lambda)}$, where $e^{O(1/\lambda)} = 1 + O(1/\lambda)$, we have that

1. For $f > 0$,

$$y = \left[ \frac{A}{f^{1/4}(x)} \exp \left( i\lambda \int f^{1/2}(x) dx \right) + \frac{B}{f^{1/4}(x)} \exp \left( -i\lambda \int f^{1/2}(x) dx \right) \right] (1 + O(1/\lambda)).$$

2. For $f < 0$,

$$y = \left[ \frac{A}{(-f)^{1/4}(x)} \exp \left( \lambda \int (-f)^{1/2}(x) dx \right) + \frac{B}{(-f)^{1/4}(x)} \exp \left( -\lambda \int (-f)^{1/2}(x) dx \right) \right] (1 + O(1/\lambda)).$$

These are called the **WKB approximations** and they can sometimes be used “off the shelf”. It is clear that they are only valid if $f$ is non-zero. Also note that for $f > 0$ we can write

$$y = \frac{A}{f^{1/4}(x)} \cos \left( \lambda \int f^{1/2}(x) dx \right) + \frac{B}{f^{1/4}(x)} \sin \left( \lambda \int f^{1/2}(x) dx \right) + O(1/\lambda),$$

since the functions involved are $O(1)$. An example of an off-the shelf use is this: approximate for large $\lambda$ the solution of the boundary value problem

$$y'' + \lambda^2 x^4 y = 0, \quad y(1) = 1, \quad y(2) = 2.$$

Here $f(x) = x^4 > 0$,

$$\int f(x)^{1/2} dx = \int x^2 dx = \frac{x^3}{3},$$

so that

$$y(x) = \frac{A}{x} \cos \left( \frac{\lambda x^3}{3} \right) + \frac{B}{x} \sin \left( \frac{\lambda x^3}{3} \right) + O(1/\lambda),$$

and $A$ and $B$ are found from the boundary conditions which give the linear system

$$1 = A \cos \left( \frac{\lambda}{3} \right) + B \sin \left( \frac{\lambda}{3} \right),$$

$$2 = A \cos \left( \frac{8\lambda}{3} \right) + B \sin \left( \frac{8\lambda}{3} \right).$$
So the off-the shelf formula are good for $y'' + \lambda^2 f(x)y = 0$. Now let us do an example when these formulae do not work, but the approach can still be used. Consider

$$y'' - (\lambda^2 x^4 + \lambda x^3)y = 0, \quad x > 0$$

and let us find the general form of the solution.

Using the WKB ansatz $y = e^{\lambda y}$, writing $Y = Y_0 + 1/\lambda Y_1 + O(1/\lambda^2)$, and dividing by $\lambda^2$, we have

$$\frac{1}{\lambda} \left( Y_0'' + \frac{1}{\lambda} Y_1'' \right) + \left( (Y_0')^2 + 2 \frac{1}{\lambda} Y_0' Y_1' \right) - x^4 - \frac{1}{\lambda} x^3 = O(1/\lambda^2).$$

So $Y_0' = \pm x^2$ and hence $Y_0 = \pm x^3/3$.

For $Y_1$ we have

$$\pm 2x \pm 2x^2 Y_1' - x^3 = 0,$$

from which we conclude that

$$Y_1' = \pm \frac{x}{2} - \frac{1}{x},$$

so that $Y_1 = -\log x \pm x^2/4$, and the final result is

$$y(x) = \left( \frac{A}{x} e^{\lambda x^3/3 + x^2/4} + \frac{B}{x} e^{-\lambda x^3/3 - x^2/4} \right) (1 + O(1/\lambda)).$$

## 14.2 Eigenvalue problems

One of the main applications of WKB method is to produce asymptotics of eigenvalues for elliptic eigenvalue problems. Recall that for a problem of the form

$$y'' + \lambda^2 f(x)y = 0, \quad y(0) = y(1) = 0,$$

$y \equiv 0$ is always a solution, but there is an unbounded set of numbers $\lambda_n \in \mathbb{R}$, for which there are non-zero solutions $y_n$ to the problem (eigenfunctions). Also note that if $y_n$ is a solution to the problem, so is $\alpha y_n$ for any $\alpha$. So an “eigenvalue problem” is the problem of finding the pairs $(\lambda_n, y_n)$ of eigenvalues and eigenfunctions. Here is how we use the off-the shelf formulae to find the asymptotics of eigenvalues in a typical problem (of course the method is not restricted only to such eigenvalue problems):

Consider the eigenvalue problem

$$y'' + \lambda^2 (1+x)^2 y = 0, \quad y(0) = y(1) = 0.$$

Here $f(x) = (1+x)^2 > 0$, $(f(x))^{1/2} = 1+x$, $\int f^{1/2}(x) dx = x + x^2/2$, $(f(x))^{1/4} = \sqrt{1+x}$. So we immediately have that the WKB approximation is

$$y(x) = \frac{A}{\sqrt{1+x}} \cos(\lambda(x + x^2/2)) \\ + \frac{B}{\sqrt{1+x}} \sin(\lambda(x + x^2/2)) + O(1/\lambda).$$
The boundary condition at $x = 0$ means that $A = 0$. Therefore the boundary condition at $x = 1$ reads

$$0 = \frac{B}{\sqrt{2}} \sin(3\lambda/2).$$

This will be possible for non-zero $B$ if the sine term is 0, but that means that $\lambda$ has to take one of the values

$$\lambda_n = \frac{2}{3}n\pi.$$

These are our (approximate) eigenvalues, and the approximation gets better and better as $n$ becomes larger. I cannot say more than that. As for the eigenfunctions, it is easier to give an error estimate:

$$y_n(x) = \frac{B}{\sqrt{1+x}} \sin\left(\frac{2}{3}n\pi\left(x + \frac{x^2}{2}\right)\right) + O(1/n),$$

because $\lambda_n = O(n)$.

# 15 Averaging

The averaging method was invented by Lagrange in the XVIIIth century and developed by Mitropol'ski, Krylov and Bogolyubov in the thirties. We deal with $n$-dimensional systems of the form

$$x' = \epsilon f(x, t, \epsilon), \qquad (50)$$

where $f$ is periodic in $t$ with period $T$ and $\epsilon$ is small. The method can also be applied to autonomous systems, as we shall show. Note that below, by a change of variable, we obtain an asymptotic expansion of the equation, not of the solution, as in preceding chapters.

To the system (50) we associate the *averaged* (autonomous) system

$$y' = \epsilon \frac{1}{T} \int_{0}^{T} f(y, t, 0) dt \equiv \epsilon \bar{f}(y). \qquad (51)$$

Then we have

**Theorem 16.** Consider equation (50), where $f$ is $C^2$ and of period $T$ in $t$. Then

1. There is a change of variable

$$x(t) = y(t) + \epsilon w(y, t, \epsilon), \qquad (52)$$

such that $y(t)$ satisfies

$$y' = \epsilon \bar{f}(y) + O(\epsilon^2).$$

2. If $y(t)$ solves (51), $x(0) = x_0$, $y(0) = y_0$ and $\|x_0 - y_0\| = O(\epsilon)$, then $\|x(t) - y(t)\| = O(\epsilon)$ for $t = O(1/\epsilon)$.
3. If $p$ is a hyperbolic rest point of the (51), then there exists $\epsilon_0 > 0$, such that for all $\epsilon \in (0, \epsilon_0)$ (50) has a hyperbolic $T$-periodic solution $\gamma(t)$, such that $\|\gamma(t) - p\| = O(\epsilon)$ for all time, having the same stability as $p$.

Let us understand what this theorem does for us. 2. tells us that the averaged equation gives for bounded times a good approximation to the original equation, while the most important part, 3., tells us that finding periodic solutions and investigating their stability reduces to the question of finding stationary points of the autonomous averaged equation and investigating their stability.

To make the procedure of deriving the $y$ equation less mysterious and also to show what kind of results can be obtained, let us consider the simple example

$$x' = \epsilon x \sin^2 t. \tag{53}$$

Here

$$\bar{f}(x) = \frac{1}{\pi} \int_{0}^{\pi} x \sin^2 t \, dt = \frac{x}{2}.$$

To define $w$, put

$$\tilde{f}(x, t, \epsilon) := x \sin^2 t - \frac{x}{2} = -\frac{x}{2} \cos 2t,$$

and now set

$$\frac{\partial w}{\partial t} = -\frac{y}{2} \cos 2t,$$

so that

$$w = -\frac{y}{4} \sin 2t.$$

If we now introduce (52) into (53), we obtain

$$y' + \epsilon \left( \frac{\partial w}{\partial y} y' + \frac{\partial w}{\partial t} \right) = \epsilon \frac{y + \epsilon w}{2} - \epsilon \frac{y + \epsilon w}{2} \cos 2t,$$

or, by remembering what $w$ is,

$$y'(1 - \frac{\epsilon}{4} \sin 2t) = \epsilon \frac{y}{2} - \frac{\epsilon^2}{8} y \sin 2t (1 - \cos 2t). \tag{54}$$

So far no approximation has been made. We must, however, replace the reciprocal of the term in brackets in the left hand-side of the above equation by a polynomial in $\epsilon$. This is done by using the Maclaurin series in $\epsilon$, that is,

$$\frac{1}{1 - \frac{\epsilon}{4} \sin 2t} = 1 + \frac{\epsilon}{4} \sin 2t + \frac{\epsilon^2}{16} \sin^2 2t + \dots,$$

so that (54) becomes

$$y' = \frac{\epsilon}{2} y + \frac{\epsilon^2}{16} y \sin 4t + O(\epsilon^3).$$

So the averaged autonomous system is

$$y' = \frac{\epsilon}{2} y,$$
the solution of which is $y_0 e^{\epsilon t/2}$ and you are advised to compare it to the exact solution of the original equation (53), which is

$$x(t) = x_0 e^{\epsilon((t/2)-\sin(2t))/4},$$

for $x_0$ and $y_0$ close.

As an example of the use of the averaging theorem, consider

$$x' = \epsilon(-x + \cos^2 t).$$

It is in the right form for averaging. Here $f(x, t, \epsilon) = (-x + \cos^2 t)$. Thus

$$\bar{f}(y) = \frac{1}{\pi} \int_{0}^{2\pi} (-y + \cos^2 t) dt = -y + 1/2.$$

Therefore the averaged equation is

$$y' = \epsilon(-y + 1/2).$$

It clearly has an asymptotically stable rest point $y = 1/2$, so the original equation has an asymptotically stable $\pi$-periodic solution close to $x = 1/2$.

## 15.1 Averaging: the 2-dimensional case

What Krylov and Bogolyubov realised was that the same idea could be applied to autonomous equations of the form $x'' + \epsilon h(x, x') + x = 0$. This is not obvious, so let us see their argument. In the usual way, write $x' = y$, so that we have a system,

$$\begin{cases} x' = y, \\ y' = -\epsilon h(x, y) - x. \end{cases}$$

Now put $x(t) = a(t) \cos \theta(t)$, $y(t) = a(t) \sin \theta(t)$, and our system of equations becomes

$$\begin{cases} a' = -\epsilon h(a \cos \theta, a \sin \theta) \sin \theta, \\ \theta' = -1 - \epsilon a^{-1} h(a \cos \theta, a \sin \theta) \cos \theta. \end{cases} \quad (55)$$

Now put $\phi = \theta + t$ and the equations become

$$\begin{cases} a' = -\epsilon h(a \cos(-t + \phi), a \sin(-t + \phi)) \sin(-t + \phi), \\ \phi' = -\epsilon a^{-1} h(a \cos(-t + \phi), a \sin(-t + \phi)) \cos(-t + \phi). \end{cases}$$

So now we can apply the averaging theorem. Note that

$$\begin{align*} & \int_{0}^{2\pi} h(a \cos(-t + \phi), a \sin(-t + \phi)) \sin(-t + \phi) dt \\ &= - \int_{\phi}^{-2\pi+\phi} h(a \cos z, a \sin z) \sin z dz \\ &= \int_{0}^{2\pi} h(a \cos z, a \sin z) \sin z dz. \end{align*}$$
It is customary to use the same letters for the averaged variables and I hope this causes no confusion.

Hence the averaged system is

$$ egin{cases} a' = -\epsilon p_0(a) \\ \theta' = -1 - \epsilon a^{-1} r_0(a), \end{cases} $$

where

$$ p_0(a) = \frac{1}{2\pi} \int_{0}^{2\pi} h(a \cos z, a \sin z) \sin z \, dz $$

and

$$ r_0(a) = \frac{1}{2\pi} \int_{0}^{2\pi} h(a \cos z, a \sin z) \cos z \, dz $$

Note the amazing simplification we have achieved! The first equation is an autonomous 1st order equation. So the averaging theorem tells us that the limit cycle of the van der Pol equation (for small $\epsilon$) is unique and asymptotically stable: In this case

$$ p_0(a) = \frac{1}{2\pi} \int_{0}^{2\pi} a \sin z (a^2 \cos^2 z - 1) \sin z \, dz = \frac{a(a^2 - 4)}{8}. $$

So just from one integral we know that not only is the limit cycle unique and asymptotically stable, but also that its amplitude is $2 + O(\epsilon)$. Not bad. But we can do more. We can estimate its period as well to $O(\epsilon^2)$.

From (55) we have that

$$ \frac{da}{d\theta} = \frac{\epsilon h(a \cos \theta, a \sin \theta) \sin \theta}{1 + \epsilon h(a \cos \theta, a \sin \theta) \cos \theta}, \quad (56) $$

This equation means that in one period of $a$, $\theta$ changes by $2\pi$.

Suppose that $p_0(a)$ has a simple zero at $a = a_0$. So from the averaging theorem, if $a(t)$ is a periodic solution of (56), $a(t) = a_0 + O(\epsilon)$. Hence from

$$ \theta' = -1 - \epsilon a^{-1} h(a \cos \theta, a \sin \theta) \cos \theta, $$

we have that

$$ \begin{align*} T &= \int_{0}^{T} dt = \int_{2\pi}^{0} \frac{d\theta}{-1 - \epsilon a^{-1} h(a \cos \theta, a \sin \theta)) \cos \theta} \\ &= \int_{0}^{2\pi} \frac{d\theta}{1 + \epsilon a^{-1} h(a \cos \theta, a \sin \theta) \cos \theta} \\ &= \int_{0}^{2\pi} d\theta (1 - \epsilon a^{-1} h(a \cos \theta, a \sin \theta) \cos \theta) + O(\epsilon^2) \\ &= 2\pi - \epsilon 2\pi a_0^{-1} r_0(a_0) + O(\epsilon^2). \end{align*} $$

In the case of the van der Pol equation we have

$$ r_0(a) = \frac{1}{2\pi} \int_{0}^{2\pi} a \sin z (a^2 \cos^2 z - 1) \cos z \, dz = 0, $$

so that $T = 2\pi + O(\epsilon^2)$.
16 A taste of renormalization

Renormalization group theory has been used to “tame” divergent integrals in QFT for a long time and found many applications in statistical mechanics, various natural sciences, such as geophysics, and even the theory of dynamical systems and PDEs. In late 1980s Goldenfeld, Oono and others noticed that RG ideas can be applied to perturbation problems and since that time, these ideas have been refined by many people in different ways. A very transparent presentation of what this approach is about is in the SIAM Review paper by E. Kirkinis (54, 374–388 (2012)) and I will discuss the first (linear) example from there, adding steps he misses.

As in averaging where the approach miraculously uncoupled the amplitude equation, there is some-
thing mysterious happening here. In most of our investigations, we started with a failure and then
adapted the method to take care of whatever was it that caused the failure. In RG methods, the
philosophy is different: failure is confronted head on, from conviction that as the problem has a
solution, failure simply has to be handled correctly to lead us to it.

Let us for the last time in this course consider the linear damped oscillator,

$$
y'' + \epsilon y' + y = 0, \quad y(0) = 3, \quad y'(0) = 1.
$$

We know that RPE is going to lead to secular terms and that is why we used the method of multiple
scales. The RG approach is to use the structure of the secular terms cleverly. We also need a different
representation of solutions.

So we start with the RPE, $y \sim y_0 + \epsilon y_1 + \epsilon^2 y_2 + O(\epsilon^3)$ and solve equations one by one. However we will write our $y_0$ differently and we will defer using the initial conditions.

So $y_0'' + y_0 = 0$, and since our solution is real, we will write it in the general form

$$
y_0(t) = A e^{it} + A^* e^{-it} := A e^{it} + cc,
$$

where *cc* stand for complex conjugate.

Let us compute the next term in this representation.

$$
y_1'' + y_1 = -y_0' = -iAe^{it} + cc,
$$

from which we get the particular integral in the form $(\alpha t)e^{it}$ as $y_1 = -1/2 Ate^{it} + cc$. Note how neat this notation is. Since $y_n'' + y_n = -y_{n-1}'$, it is not hard to get as many terms as you wish (always putting the general solution of the homogeneous equation to 0).

Just for fun, let me quote two more:

$$
y_2 = \frac{1}{8} A(t^2 - it)e^{it} + cc, \quad y_3 = -\frac{1}{16} A(t^3/3 - it^2)e^{it} + cc,
$$

and so on. The point to notice is that all these terms starting with $y_1$ are secular. Also worth noting is that they all depend linearly on $A$, which would not be the case in a nonlinear problem. Thus the
solution can be written in the form

$$y^"=A(1+\epsilon y_{1A}+\epsilon^2 y_{2A}+\dots)e^{it}+cc,$$

where I put

$$y_{1A}=-\frac{1}{2}t, \quad y_{2A}=\frac{1}{8}A(t^2-it),$$

and so on. Of course this is a formal equality, not even asymptotic.

Let us call the “term” in brackets, Y. So

$$y^"=AYe^{it}+cc.$$

So far, no new ideas. The new idea is now. We say, OK, let us call $AY S(t)$ and see what equation it satisfies. Perhaps we can solve it even though it makes no rigorous sense. So $A = Y^{-1}S$, $A$ is a constant, and hence by implicit differentiation with respect to time,

$$\frac{S'}{Y}=\frac{SY'}{Y^2}.$$

But this means that

$$S'=S\frac{Y'}{Y}.$$

So far so good, but it is getting better:

$$Y'=\epsilon y'_{1A}+\epsilon^2 y'_{2A}+\epsilon^3 y'_{3A}+\dots,$$

where the equality is absolutely formal. Also,

$$\frac{1}{Y}=1-\epsilon y_{1A}+\epsilon^2(y_{1A}^2-y_{2A})+\dots$$

(Just think what we have done here!) Hence

$$\begin{aligned} \frac{Y'}{Y} = & \epsilon y'_{1A} + \epsilon^2 (y'_{2A} - y'_{1A} y_{1A}) \\ & + \epsilon^3 (y'_{3A} - y_{1A} y'_{2A} + y'_{1A} y_{1A}^2 - y'_{1A} y_{2A}) + \dots \end{aligned}$$

I hope you are ready. Let us compute, as Leibniz said.

$$y'_{1A}=-1/2, \quad y'_{2A}-y'_{1A}y_{1A}=-\frac{i}{8}, \quad y'_{3A}-y_{1A}y'_{2A}+y'_{1A}y_{1A}^2-y'_{1A}y_{2A}=0,$$

i.e. all the terms are independent of $t$! I suggest you do one more step. And so the equation for $S$ is

$$S'=-1/2\epsilon S+i\left(-\frac{\epsilon^2}{8}-\frac{\epsilon^4}{128}+O(\epsilon^6)\right).$$

Now the equality is at least in the sense of asymptotic series, but the series looks even convergent for small enough $\epsilon$.

Clearly, you will generate a constant of integration when finding the bounded and decaying function $S$. Finally write $y = Se^{it} + cc$ and plug in the initial condition. Whew.
# 17 Optimal Truncations

We have learnt a number of ways of generating series expansions in polynomial and differential equations. Some of these expansions were convergent, and some were divergent.

As we saw in Question 1(b) in Assignment 1, the series can have zero radius of convergence. There we constructed an approximation for the centre manifold of

$$x' = -x^2,$$
$$y' = -y + x^2,$$

and the expansion was $h(x) = \sum_{n=0}^{\infty} (n-1)!x^n$. But note that the centre manifold $h(x)$ satisfies a differential equation, $-h'(x)x^2 = -h(x) + x^2$. So we could try to solve it, and see how the solution is related to the infinite series we generate. Unfortunately the manipulations though feasible would obstruct the ideas required. I will leave this example as an exercise, noting that solution of the ODE is only defined for $x \le 0$. We will do a simpler example, following a 2004 paper of C. Rousseau.

So let us instead consider Euler's equation $y'x^2 = -y + x$ for $x > 0$ small. If we expand $y$ as a power series in $x$, we obtain

$$y \sim \sum_{n=0}^{\infty} (-1)^n n! x^{n+1}, \qquad (57)$$

which clearly has zero radius of convergence. On the other hand, using the method of integrating factors, we have

$$(ye^{-1/x})' = \frac{e^{-1/x}}{x}.$$

Hence

$$y = Ce^{1/x} + \int_{0}^{x} \frac{e^{1/x-1/s}}{s} ds.$$

If we want a solution that remains bounded as $x \to 0_+$, we must put $C = 0$. The integral is better understood if we make the change of variable $z/x = 1/s - 1/x$, so that in terms of $z$, our solution becomes

$$y = I(x) := \int_{0}^{\infty} \frac{e^{-z/x}}{1+z} dz, \qquad (58)$$

from which it is clear that it exists for all $x \ge 0$.

So the questions we have to address are: (i) How well does (57) approximate (58) and (ii) What is the exact relation between the two?

Before we discuss the first question consider the question of approximating $0 = \sin \pi$ by its convergent Maclaurin series. Here the value $\pi$ is fixed, we just have to take more and more terms to get closer and closer to 0. For example taking 40 terms gives an error of approx. $10^{-9}$.

For asymptotic series, taking a fixed $x$ and increasing the number of terms generally does not result in better approximation: For each $x$ there usually is an optimal number of terms we should take. Taking that number of terms is what is called **optimal truncation** of the asymptotic series and the
approximation achieved at optimal truncation is called **superasymptotic**. Any approximation that improves on the superasymptotic approximation has been called (by Berry) **hyperasymptotic**. Let us consider these concepts using our example. Let us put

$$S_k(x) = \sum_{n=0}^{k-1} (-1)^n n! x^{n+1},$$

Now write

$$\frac{1}{1+z} = \sum_{n=0}^{k-1} (-1)^n z^n + (-1)^k \frac{z^k}{1+z}.$$

Putting this relation into the integral and integrating term by term, we have

$$I(x) = \sum_{n=0}^{k-1} (-1)^n \int_{0}^{\infty} e^{-z/x} z^n dz + (-1)^k \int_{0}^{\infty} \frac{e^{-z/x} z^k}{1+z} dz.$$

Now remember that $\int_{0}^{\infty} z^n e^{-z} dz = n!$. Hence

$$\int_{0}^{\infty} e^{-z/x} z^n dz = n! x^{n+1}.$$

Therefore $|I(x) - S_k(x)| = |R_k(x)|$, where the remainder $R_k(x)$ is given by

$$R_k = (-1)^k \int_0^\infty \frac{e^{-z/x} z^k}{1+z} dz.$$

But then

$$|R_k(x)| \le \int_0^\infty z^k e^{-z/x} dz = k! x^{k+1}.$$

So that is a bound on the error.

Now the question arises, for a **given** (small, positive) $x$ how many terms should we take to minimise the error. There are two strategies.

One of them does not use the estimate on the remainder at all. The logic behind it is: in an asymptotic expansion $\sum_{n=0}^{\infty} a_n \phi_n(x)$, take the largest number $N$ of terms for which

$$\frac{|a_N \phi_N(x)|}{|a_{N+1} \phi_{N+1}(x)|} \ge 1.$$

Let us use this in our example. We must have

$$\frac{N! x^{N+1}}{(N+1)! x^{N+2}} \ge 1,$$

so for small $x$, $N \sim 1/x$. Thus the optimal truncation by this recipe is at $[1/x]$ terms.
As another example, $I_0(x)$, the Bessel function of order 0 has as $x \to \infty$ the asymptotic expansion $\sum_{n=0}^{\infty} a_n \phi_n(x)$, where $a_0 = 1$, $a_n$s are computed recursively from

$$a_{n+1} = \frac{(2n+1)^2}{8(n+1)} a_n,$$

and $\phi_n(x) = e^x x^{-n-1/2}$. Then we must have

$$\frac{a_N x^{N+1}}{a_{N+1} x^N} \ge 1,$$

from which it follows that

$$\frac{8x(N+1)}{(2N+1)^2} \ge 1,$$

from which $N$ can be found. So for example, for $x = 10$, $N = 19$.

Though the two are largely equivalent (when?), a more principled approach uses the estimate of the remainder: truncate when the remainder is at a minimum. Let us go back to our example. We had $|R_k(x)| \le k! x^{k+1}$.

For $x$ small, we expect the optimal $k = N$ to be large. But then Stirling's formula gives a good approximation to the factorial,

$$k! \sim \sqrt{2\pi} k^{k+1/2} e^{-k},$$

which already for $k = 5$ is 118! So put $f(k) = \sqrt{2\pi} k^{k+1/2} e^{-k} x^{k+1}$ and let us find the minimum of $f(k)$. Taking logs and differentiating, we have that

$$\frac{f'(k)}{f(k)} = \ln k + \frac{k+1/2}{k} - 1 + \ln x,$$

which again gives $N = [1/x]$.

However, having the remainder is good for estimating how good our approximation is. If we have $k = N \sim 1/x$, the remainder satisfies

$$|R_N(x)| \le \sqrt{2\pi} (1/x)^{(1/x+1/2)} e^{-1/x} x^{1/x+1} = \sqrt{2\pi x} e^{-1/x}.$$

So for $x = 0.2$, truncating after 5 terms, gives a relative error of 4%.

For various hyperasymptotic strategies (he names 8), please see the paper by Boyd, *Acta Applicandae Mathematica* **56**, 1–98 (1999).

We are now ready to confront the second question, which can be formulated as follows: is there a transformation that takes the divergent series (57) into the integral (58).

# 18 Borel summation

A good theory for summation of series is one that assigns to convergent series their classical sum, and makes sense of arithmetic operations and term by term differentiation. In other words, if $S$ is
the sum it assigns to a series (we say that the series is summable in the theory), the theory has to
have the following properties:

* If $\sum_{n=0}^{\infty} a_n$ converges, $S$ is its usual sum;

* If $\sum_{n=0}^{\infty} a_n = S_1$ in the theory, $\sum_{n=0}^{\infty} b_n = S_2$, then $\sum_{n=0}^{\infty} (a_n + b_n) = S_1 + S_2$;

* $\sum_{n=0}^{\infty} a_n$, $\sum_{n=0}^{\infty} b_n$ are absolutely summable, then $\sum_{n=0}^{\infty} a_n b_n = S_1 S_2$;

* If $\sum_{n=0}^{\infty} a_n = S$, $\sum_{n=1}^{\infty} a_n = S - a_0$;

* If $\sum_{n=0}^{\infty} a_n x^n = f(x)$, then $\sum_{n=0}^{\infty} na_n x^{n-1} = f'(x)$.

**Definition 17.** A series $\sum_{n=0}^{\infty} a_n$ is **Borel-summable** if the series $\sum_{n=0}^{\infty} a_n z^n/n!$ has a non-zero radius of convergence and if the integral

$$
\int_0^\infty \left( \sum_{n=0}^\infty \frac{a_n}{n!} z^n \right) e^{-z} dz
$$

converges. Then this integral is the **Borel sum** of the series $\sum_{n=0}^{\infty} a_n$.

Consider $1 - 1 + 1 - 1 + \dots$. Then clearly

$$
\sum_{n=0}^{\infty} \frac{(-1)^n}{n!} z^n = e^{-z},
$$

the series having an infinite radius of convergence. So the alternating series of ones is Borel-
summable and its sum is

$$
\int_0^\infty e^{-2z} dz = \frac{1}{2}.
$$

Let us just check that the Borel sum of a convergent series is the classical sum. So suppose $\sum_{n=0}^{\infty} a_n$ exists in the classical sense and is equal to $S$. Then clearly $\sum_{n=0}^{\infty} a_n$ is Borel-summable. We have

$$
\begin{align*}
S &= \sum_{n=0}^{\infty} a_n = \sum_{n=0}^{\infty} a_n n! / n! = \sum_{n=0}^{\infty} \frac{a_n}{n!} \int_0^{\infty} z^n e^{-z} dz \\
  &= \int_0^{\infty} \left( \sum_{n=0}^{\infty} a_n z^n / n! \right) e^{-z} dz.
\end{align*}
$$

So the usual sum of a convergent series is also its Borel sum.

Here is our pièce de résistance: What is the Borel sum of our series (57)? The series is $\sum_{n=0}^{\infty} (-1)^n n! x^{n+1}$.
Hence

$$
\sum_{n=0}^{\infty} a_n z^n / n! = x \sum_{n=0}^{\infty} (-1)^n (xz)^n = \frac{x}{1+xz}
$$

for all $z$ such that $|z| < 1/x$. So the series is Borel-summable and its Borel sum is

$$
\int_0^\infty \frac{x}{1+xz} e^{-z} dz = \int_0^\infty \frac{e^{-z/x}}{1+z} dz.
$$

Thus Borel summation of the divergent series recovers the solution of the ODE.
19 Asymptotics of Integrals

Now we start on a quick tour of methods for finding asymptotic expansions of integrals. These methods are very extensive and sophisticated and there are whole books written about them, for example by Wong. We will only cover Watson's lemma, Laplace's method, and the method of stationary phase. My main reference is the book by Murray [16]. But let us start with a simple example.

Consider the error function,

$$
\mathrm{Erf}(x) = \frac{2}{\sqrt{\pi}} \int_{0}^{x} e^{-t^2} dt. \qquad (59)
$$

Recall that $\lim_{x \to \infty} \mathrm{Erf}(x) = 1$. We are interested in approximating the error function for large values of $x$. There are simple possible approaches. One is to note that the integrand is a real analytic function for all $t$, expand it in a Taylor series around 0 (which thus converges for all $t$), and interchange summation and integration.

Then we get

$$
\mathrm{Erf}(x) = \frac{2}{\sqrt{\pi}} \sum_{n=0}^{\infty} \frac{(-1)^n x^{2n+1}}{(2n+1)n!}, \quad (60)
$$

so that

$$
\mathrm{Erf}(x) = \frac{2}{\sqrt{\pi}} \left( x - \frac{1}{3}x^3 + \frac{1}{10}x^5 - \dots \right).
$$

This (convergent) expansion involves large numbers, and the final result, which has to be between
zero and 1, is obtained by cancellation between large numbers of opposite signs. Since computers
work to finite precision, it is practically impossible to compute Erf(x) with a given accuracy for a
large enough value of x.

A different approach is to write

$$
\mathrm{Erf}(x) = 1 - \frac{2}{\sqrt{\pi}} \int_{x}^{\infty} e^{-t^2} dt
$$

integrate by parts as many times as desired and estimate the remainder.

Thus

$$
\int_x^\infty e^{-t^2} dt = - \int_x^\infty \frac{1}{2t} \frac{d}{dt}(e^{-t^2}) = \frac{e^{-x^2}}{2x} - \int_x^\infty \frac{e^{-t^2}}{2t^2} dt.
$$

Continuing in this vein, we obtain the following representation

$$
\mathrm{Erf}(x) = 1 - \frac{2}{\sqrt{\pi}} \left( \sum_{n=0}^{n=N} a_n \phi_n(x) + R_{N+1}(x) \right).
$$

Below we use the notation $n? = 1 \cdot 3 \cdots (2n-1)$, $0? = 1$, $1? = 1$. In the expression above

$$
a_n = (-1)^n n?, \quad \phi_n(x) = \frac{e^{-x^2}}{2x(2x^2)^n}
$$
and

$$R_N(x) = (-1)^N \int_x^\infty \frac{N?e^{-t^2}}{(2t^2)^N} dt.$$

A simple computation shows that $\phi_{n+1}(x) = o(\phi_n(x))$ as $x \to \infty$ and that $R_{N+1}(x) = o(\phi_N(x))$ as $x \to \infty$. To see this, note that

$$|R_{N+1}(x)| = \left| \int_x^\infty \frac{(N+1)?}{2t(2t^2)^{N+1}} \frac{d}{dt}(e^{-t^2}) dt \right| \le |a_{N+1}\phi_{N+1}(x)|.$$

Therefore we can write

$$\text{Erf}(x) \sim 1 - \frac{2}{\sqrt{\pi}} \sum_{n=0}^{\infty} a_n \phi_n(x). \qquad (61)$$

and (61) is an asymptotic expansion as $x \to \infty$. Also note that (61) again diverges for all $x$ by the Leibniz test.

## 19.1 Watson's Lemma

In the asymptotic expansion of integrals, methods have to be fitted to the structure of the integrand. Watson's lemma gives off-the shelf formula for asymptotic expansion of integrals of a particular form,

$$I(x) = \int_{0}^{T} e^{-xt} \phi(t) dt,$$

as $x \to \infty$. We will discuss restrictions on $\phi$ below; certainly the integral must make sense! If $T = \infty$ and $x$ is allowed to be complex, Watson's lemma is useful in understanding asymptotics of Laplace transforms, and that is already useful enough! Many integrals not in that form can be brought to that form, as we will see. We shall take

$$\phi(t) = t^{\lambda} g(t),$$

where $\lambda > -1$ to make the function integrable near 0, and we assume that $g(0) \neq 0$, and that $g(t)$ can be expanded in a Taylor expansion around 0. Why around 0?

... Simply because if $x$ is large $e^{-xt}$ is very peaked around zero and all the contributions come from around $t = 0$. The lemma involves $\Gamma(x)$, so something needs to be said about it. $\Gamma(n) = (n-1)!$ for all positive integers. Furthermore, for any positive $\alpha \in \mathbb{R}$,

$$\Gamma(\alpha) = \int_{0}^{\infty} z^{\alpha-1} e^{-z} dz.$$

We will also need the result that for all $a, u > 0$, $(1 + u)^a \le e^{au}$, which is true since raising to a positive degree is a monotone function, and since

$$e^u = 1 + u + \frac{u^2}{2} \cdots \ge 1 + u.$$

Now we are ready to state the Lemma.
**Lemma 18 (Watson's Lemma).** *Under the above conditions*

$$ \int_0^T e^{-xt} t^\lambda g(t) dt \sim \sum_{n=0}^\infty \frac{g^{(n)}(0) \Gamma(\lambda + n + 1)}{n! x^{n+\lambda+1}} \quad \text{as } x \to \infty. $$

Note that $T$ does not appear anywhere! Let us do an “off-the-shelf” example. Let

$$ I(x) = \int_0^{\pi/2} e^{-xt} (1 + \cos t)^{1/2} dt, \quad x \to \infty. $$

Then by the Lemma with $\lambda = 0$, we have

$$ I(x) \sim \frac{2^{1/2}}{x} + O\left(\frac{1}{x^3}\right), $$

as $g'(0) = 0$. Just in case you are wondering about the $\pi/2$. MAPLE gives for that integral with $x = 100 \ 0.1414178208$, while our approximation to $O(1/x^2)$ gives 0.1414213562. Not bad. Let us do a slightly more challenging example. Consider

$$ I(x) = \int_{0}^{\infty} e^{-xt} \log(1 + t^2) dt \quad \text{as } x \to \infty. $$

The problem here is that $\phi(0) = \log(1) = 0$. So we need to expand $\log(1+t^2)$ in a Maclaurin series,

$$ \log(1 + t^2) = t^2 - \frac{1}{2}t^4 + \frac{1}{3}t^6 + O(t^8) = t^2\left(1 - \frac{1}{2}t^2 + \frac{1}{3}t^4\right) + O(t^8). $$

Hence $\lambda = 2$ and the Maclaurin expansion of $g(t)$ is

$$ g(t) = 1 - \frac{1}{2}t^2 + \frac{1}{3}t^4 + O(t^6). $$

So

$$ I(x) \sim \frac{1 \cdot \Gamma(3)}{0! x^3} + 0 - \frac{1 \cdot \Gamma(5)}{2! x^5} + O(1/x^6) = \frac{2}{x^3} - \frac{12}{x^5} + O(1/x^6). $$

**Proof of Watson's Lemma:** We will prove it under the assumption that $T < R$, where $R$ is the radius of convergence the Maclaurin series of $g(t)$. Let us write

$$ g(t) = \sum_{n=0}^{\infty} \frac{g^{(n)}(0)}{n!} t^n := \sum_{n=0}^{\infty} a_n t^n. $$

for $|t| < R$. So

$$ g(t) = \sum_{n=0}^{N} a_n t^n + R_N(t) $$

for any $N$, where for all $t \in [0, T]$, since we assumed $T < R$, we have that

$$ |R_N(t)| \leq Lt^{N+1} $$

for some $L > 0$. (what is it really?)
Hence we can write

$$I(x) = \int_{0}^{T} e^{-xt} t^{\lambda} \sum_{n=0}^{N} a_{n} t^{n} dt + \int_{0}^{T} e^{-xt} t^{\lambda} R_{N}(t) dt := I_{1} + I_{2}.$$

Let us start with $I_2$. We have

$$|I_2| \le L \int_0^T e^{-xt} t^{\lambda+N+1} dt = L x^{-(\lambda+N+2)} \int_0^{xT} e^{-\tau} \tau^{\lambda+N+1} d\tau := LI_3.$$

where I have put $\tau = xt$. Now,

$$\begin{aligned} I_3 &= x^{-(\lambda+N+2)} \left( \int_0^\infty - \int_{xT}^\infty \right) e^{-\tau} \tau^{\lambda+N+1} d\tau \\ &= x^{-(\lambda+N+2)} \left( \Gamma(\lambda + N + 2) - \int_{xT}^\infty e^{-\tau} \tau^{\lambda+N+1} d\tau \right) \\ &= O(x^{-(\lambda+N+2)}) - x^{-(\lambda+N+2)} \int_{xT}^\infty e^{-\tau} \tau^{\lambda+N+1} d\tau. \end{aligned}$$

We need to estimate $I_4 := x^{-(\lambda+N+2)} \int_{xT}^\infty e^{-\tau} \tau^{\lambda+N+1} d\tau$. Let us use our inequality $(1+u)^a \le e^{au}$ to get some traction.

Put $\tau = xT(1+u)$ and we have

$$\begin{aligned} I_4 &= T^{\lambda+N+2} e^{-xT} \int_0^\infty e^{-xTu} (1+u)^{\lambda+N+1} du \\ &\le T^{\lambda+N+2} e^{-xT} \int_0^\infty e^{(-xT+\lambda+N+1)u} du \\ &= \frac{T^{\lambda+N+2}}{x} e^{-xT} \frac{1}{xT - (\lambda + N + 1)} \\ &= T^{\lambda+N+1} e^{-xT} \left( 1 + \frac{\lambda + N + 1}{xT} + O(1/x^2) \right) \\ &= O(e^{-xT}) = o(x^{-(\lambda+N+2)}). \end{aligned}$$

So $I_4$ is $o(x^{-(\lambda+N+2)})$. Hence $I_3$, and therefore $I_2$ are $O(x^{-(\lambda+N+2)})$. It remains to massage $I_1$ into the right shape.

$$\begin{aligned} I_1 &= \int_0^T e^{-xt} t^{\lambda} \sum_{n=0}^{N} a_n t^n dt = \sum_{n=0}^{N} a_n \left( \int_0^{\infty} - \int_T^{\infty} \right) e^{-xt} t^{\lambda+n} dt \\ &= \sum_{n=0}^{N} a_n \Gamma(\lambda + n + 1) x^{-(\lambda+n+1)} + O(e^{-xT}), \end{aligned}$$

using the same argument as before. Putting it all together, we have that

$$I(x) \sim \sum_{n=0}^{N} a_n \Gamma(\lambda + n + 1) x^{-(n+\lambda+1)} + O(x^{-(N+\lambda+2)}).$$
Since this is true for every $N$, the result follows.

To finish off this section, let us do another example. Let us find the asymptotic expansion (so asks Murray) of

$$I(x) = \int_{0}^{1} e^{-xt^{3}} dt \text{ as } x \to \infty.$$

This is not in the right form to use Watson's lemma, but that is easy to fix. Set $z = t^3$, so that we have

$$I(x) = \frac{1}{3} \int_{0}^{1} e^{-xz} z^{-2/3} dz.$$

Now we can use Watson's lemma with $\lambda = -2/3$ and $g(z) \equiv 1$. But then

$$I(x) \sim \frac{\Gamma(1/3)}{3x^{1/3}} + EST.$$

There are no other terms. If we take again $x = 100$, the approximate answer is 0.19238666038, while evaluating the integral numerically, MAPLE gives 0.19238666037. The difference is in the 10th digit, but if you want the error to be $O(10^{-11})$, you are out of luck, there is no optimal truncation to help you.

## 19.2 Laplace's method

Laplace's method relies for its justification on Watson's Lemma, so we will just give the ideas and the construction which leads again to an off-the-shelf result which one needs surprisingly often. Before that, using the same idea as in the last example, we get by Watson's Lemma

$$\int_{0}^{T} e^{-xt^2} dt \sim \frac{1}{2} \left(\frac{\pi}{x}\right)^{1/2} + EST \text{ as } x \to \infty \quad (62)$$

as $\Gamma(1/2) = \sqrt{\pi}$. This will be useful below.

You might want to check also that

$$\int_{0}^{T} te^{-xt^{2}} dt = O(1/x)$$

by explicit integration (or Watson's Lemma if you have nothing better to do).

In Laplace's method we deal with integrands of the form $g(t)e^{xh(t)}$. As $x$ the main contribution to the integral will come from the points $t$ where $h(t)$ is at a maximum. If over the domain of integration there are many maxima, since as we show the contribution of each point of maximum $t_m$ is proportional to $e^{xh(t_m)}$, we need only to worry about the dominant maximum. To develop Laplace's formula, we might for now assume that the maximum of $h(t)$ is at the left end of the interval of integration.

So let us suppose we are integrating over $[0, T]$ and $h(0) > h(t)$ for all $t \in [0, T]$. Then there are two cases. Either (i) $h'(0) = 0$ and $h''(0) < 0$ or (ii) $h'(0) < 0$. The first case is more interesting, and I will develop it in detail, later coming back to give the formulae for case (ii).
In case (i)

$$h(t) = h(0) + 1/2 t^2 h''(0) + 1/6 t^3 h'''(0) + \dots \quad (63)$$

We start by defining a new variable instead of $t$, $s$; in the $s$ variable we will be able to use Watson's Lemma.

We set $h(t) - h(0) = -s^2$. To find $t$ in terms of $s$ we use (63). So we put $t = \alpha_1 s + \alpha_2 s^2 + \dots$ and plug into (63) to compute as many terms of the relation between $t$ and $s$ as we want terms in the asymptotics. So Laplace's method can be taken to any number of terms, but in practice people only use the leading term.

Let me concentrate on finding the leading term. From (63), matching the $s^2$ coefficient gives

$$t = \left( \frac{-2}{h''(0)} \right)^{1/2} s + O(s^2).$$

Next we need to deal with $g(t)$. Since all the contributions are from a neighbourhood of 0, we will write

$$g(t) = g(0) + g'(0)t + 1/2 g''(0)t^2 = g(0) + g'(0) \left( \frac{-2}{h''(0)} \right)^{1/2} s + O(s^2).$$

Now we ready:

Set $A = \sqrt{h(0) - h(T)} > 0$. Then

$$
\begin{align*}
I(x) &= \int_{0}^{T} g(t)e^{xh(t)} dt \\
&= \int_{0}^{A} e^{x(h(0)-s^2)} \left[ g(0) + \left( \frac{-2}{h''(0)} \right)^{1/2} g'(0)s + O(s^2) \right] ds \\
&\sim g(0) \left( \frac{-2}{h''(0)} \right)^{1/2} e^{xh(0)} \int_{0}^{A} e^{-xs^2} ds \\
&\quad + e^{xh(0)} O \left( \int_{0}^{A} s e^{-xs^2} ds \right) \\
&\sim g(0) \left( \frac{-2}{h''(0)} \right)^{1/2} e^{xh(0)} \frac{1}{2} \left( \frac{\pi}{x} \right)^{1/2} + e^{xh(0)} O(1/x).
\end{align*}
$$

Let us rewrite this nicer:

$$I(x) \sim g(0) \left( \frac{-\pi}{2xh''(0)} \right)^{1/2} e^{xh(0)} + e^{xh(0)} O \left( \frac{1}{x} \right). \quad (64)$$

This was in the case when the maximum was 0 and at the left boundary of the interval. If it is at 0, and the domain of integration is $[-\alpha, \beta]$ with $\alpha, \beta$ positive, we need to add the contributions of the
maximum to the integrals over [-$\alpha$, 0] and [0, $\beta$]. Hence if $h(t)$ has a maximum at zero,

$$ \int_{-\alpha}^{\beta} g(t)e^{xh(t)} dt \sim g(0) \left(\frac{-2\pi}{xh''(0)}\right)^{1/2} e^{xh(0)} + e^{xh(0)}O\left(\frac{1}{x^{3/2}}\right). \quad (65) $$

**Exercise:** Somewhat mysteriously $O(1/x)$ has become $O(1/x^{3/2})$. Please show that this indeed is the case!

Finally, there is nothing special about zero. If $h(t)$ has a maximum at a point $t = a$ inside the domain of integration $I$, then we have

$$ \int_I g(t)e^{xh(t)} dt \sim g(a) \left( \frac{-2\pi}{xh''(a)} \right)^{1/2} e^{xh(a)} + e^{xh(a)}O\left(\frac{1}{x^{3/2}}\right). \quad (66) $$

For case (ii), the equivalent of (64) of the relevant formula for integrals over $[0, T]$ is

$$ I(x) \sim \left( \frac{-g(0)}{2xh'(0)} \right) e^{xh(0)} + e^{xh(0)}O\left(\frac{1}{x^2}\right). \quad (67) $$

And now it is time for applications. For future use, let us remember that

$$ \int_0^T e^{-xt^2} dt \sim \frac{1}{2} \left(\frac{\pi}{x}\right)^{1/2} \text{ as } x \to \infty. $$

It is also true by symmetry that $\int_{-a}^a te^{-xt^2} dt = 0$,

Finally, note that

$$ \begin{aligned} \int_0^T t^2 e^{-xt^2} dt &= -\frac{d}{dx} \int_0^T e^{-xt^2} dt \sim -\frac{d}{dx} \left[ \frac{1}{2} \left(\frac{\pi}{x}\right)^{1/2} \right] \\ &\sim \frac{1}{4} \frac{\sqrt{\pi}}{x^{3/2}} + \text{EST as } x \to \infty. \end{aligned} $$

Let us derive Stirling's formula using the Laplace method "off-the-shelf" expression (66), which is valid for an integrand having a maximum inside the interval $I$.

We start with

$$ \begin{aligned} \Gamma(x+1) &= \int_0^\infty t^x e^{-t} dt = \int_0^\infty e^{-t+x \log t} dt \\ &= x \int_0^\infty e^{-x\tau+x \log(x\tau)} d\tau = x^{x+1} \int_0^\infty e^{x(\log(\tau)-\tau)} d\tau. \end{aligned} $$

Here the function $h(\tau)$ is $h(\tau) = \log(\tau) - \tau$, and it is not hard to see that it has a maximum at $\tau = 1$, $h(1) = -1$, $h''(1) = 1$. In this case $g(\tau) \equiv 1$ and so Laplace's formula tells us that

$$ \Gamma(x+1) \sim x^{x+1} e^{-x} \left(\frac{2\pi}{x}\right)^{1/2} + \dots = \sqrt{2\pi x^{x+1/2}} e^{-x} + \dots $$
which since $\Gamma(n + 1) = n!$ for positive integers $n$, gives Stirling's formula

$$n! \approx \sqrt{2\pi n^{n+1/2}}e^{-n}.$$

We saw that the approximation is pretty good even for $n = 5$: Stirling's formula gives 118.0192 as an approximation to 120.

However, it is tempting to see whether we have the know-how to find the next term in the approximation...

So instead of using the off-the-shelf product, let us go through Laplace's way of computing things. First of all let us move the maximum point to the origin by setting $\tau = 1+u$. So we want to compute

$$x^{x+1}e^{-x} \int_{-1}^{\infty} e^{x(\log(u+1)-u)} du := J_1.$$

We will change variables by setting $\log(1+u) - u = -s^2$, and will have to expand $u$ in terms of $s$ up to $O(s^3)$, i.e.

$$u = \alpha_1 s + \alpha_2 s^2 + \alpha_3 s^3 + \dots,$$

so that the integral $J_1$ above is transformed as follows:

$$J_1 = \alpha_1 x^{x+1} e^{-x} \int_{-\infty}^{\infty} e^{-xs^2} ds + 2\alpha_2 x^{x+1} e^{-x} \int_{-\infty}^{\infty} s e^{-xs^2} ds \\ + 3\alpha_3 x^{x+1} e^{-x} \int_{-\infty}^{\infty} s^2 e^{-xs^2} ds + \dots$$

The first term is just Stirling's formula; the second one is zero by symmetry; the third one is the one we want. So we need to compute

$$x^{x+1}e^{-x} \int_{-\infty}^{\infty} s^2 e^{-xs^2} ds = \frac{3}{2} \alpha_3 x^{x+1} e^{-x} \frac{\sqrt{\pi}}{x^{3/2}}.$$

So the whole problem has been reduced to finding $\alpha_3$!

$$\log(u+1) - u = -\frac{1}{2}u^2 + \frac{1}{3}u^3 - \frac{1}{4}u^4 + \dots = -s^2.$$

We plug in $u = \alpha_1 s + \alpha_2 s^2 + \alpha_3 s^3 + \dots$ into the above expression and collect terms in $s^2$. We obtain $\alpha_1 = \sqrt{2}$, $\alpha_2 = 2/3$ and $\alpha_3 = 1/(9\sqrt{2})$. Hence the first correction to Stirling's formula is

$$\frac{\sqrt{2\pi}}{12} e^{-x} x^{x-1/2}.$$

To summarise,

$$n! \approx \sqrt{2\pi n^{n+1/2}} e^{-n} + \frac{\sqrt{2\pi}}{12} e^{-n} n^{n-1/2} + O(e^{-n} n^{n-3/2}).$$

Let us see how good this is for $n = 5$.

$$\frac{\sqrt{2\pi}}{12} e^{-5} 5^{5-1/2} = 1.9667,$$
so we have as a two-term approximation to 120, 118.0192 + 1.9667 = 119.9769, which is quite impressive.

Let us do another example. Let us consider

$$I(x) = \int_{0}^{\pi/2} e^{x \cos t} \log(\lambda + \sin t) dt.$$

here $h(t) = \cos t$ and it has a maximum at the left end-point of the interval. If $\lambda \neq 1$, we can use Laplace's formula in the form suitable for a real maximum being at an endpoint,

$$I(x) \sim g(0) \left( \frac{-\pi}{2xh''(0)} \right)^{1/2} e^{xh(0)} + e^{xh(0)} O(1/x),$$

to give

$$I(x) \sim \log \lambda \left( \frac{-\pi}{2x} \right)^{1/2} e^x + e^x O(1/x).$$

However, what if $\lambda = 1$? In that case we need the expansion of

$$g(t) = \log(1 + \sin t).$$

We have

$$g(t) = 0 + \frac{\cos 0}{1 + \sin 0} t + O(t^2).$$

Since $h(t) - h(0) = -s^2$ leads to $-1/2t^2 + O(t^4) = -s^2$, we have that

$$t = \sqrt{2}s + O(s^2),$$

and hence $g(t) = \sqrt{2}s + O(s^2)$. Therefore we have

$$I(x) \sim \int_{0}^{1} \sqrt{2s} e^{x(1-s^2)} \sqrt{2} ds \sim \frac{e^x}{x} + o(1/x).$$

As an exercise, please check that the expansion is in fact continuous in $\lambda$ as it should be.

## 19.3 Method of stationary phase

In very many problems of Fourier analysis and wave propagation one needs to estimate integrals of the form

$$J(x) = \int_{I} g(t) e^{ixh(t)} dt \text{ as } x \to \infty,$$

where $I$ is an interval.

Note that $\text{Re } J(x) = \int_I g(t) \cos(xh(t)) dt$. As $(\cos(xh(t)))' = -xh'(t) \sin(xh(t))$, it is clear the points where $h'(t) = 0$ (turning points) are important in understanding $J(x)$. In fact, away from turning points, the fast oscillations cancel positive and negative contributions to the integral and so, whether corresponding to minima or maxima of $h(t)$ it is only the neighbourhoods of the turning
points that would make a contribution. Boundary points can also make a contribution, but in the
presence of turning points, their contribution is negligible and we will not consider it.

For orientation, consider the graph of the function $h(t) = \cos(100t(2-t))$ on the interval $[0, 2]$:

<figure><img src="image_5.png" /></figure>

You can see the “stationary phase” around $t = 1$.

A method to deal with such problems had been devised by Stokes and Lord Kelvin in the 1880s. It
is not dissimilar in spirit to Laplace’s method. Before I explain it, let us do a complex integral by
contour integration.

We want to find $\int_{0}^{\infty} e^{-ir^2} dr$. Let us integrate $e^{-z^2}$ along the contour $C$ below, where the angle of the opening of the sector is $\pi/4$.

<figure><img src="image_6.png" /></figure>

As $e^{-z^2}$ is holomorphic, the integral along the contour is zero.

But it can be written as a sum of three integrals:

$$0 = \int_C e^{-z^2} dz = \int_0^R e^{-r^2} dr + iR \int_0^{\pi/4} e^{-R^2 e^{2i\theta}} e^{i\theta} d\theta - e^{i\pi/4} \int_0^R e^{-ir^2} dr.$$

Now take the limit as $R \to \infty$. The first integral goes to $\sqrt{\pi}/2$, the second goes to zero as we show below, and the third one goes to $-e^{i\pi/4} \int_0^\infty e^{-ir^2} dr$. Hence

$$\int_0^\infty e^{-ir^2} dr = \frac{\sqrt{\pi}}{2} e^{-i\pi/4}.$$
The same result holds for the integral from $-\infty$ to $0$, and by taking complex conjugates, we obtain

$$ \int_{-\infty}^{\infty} e^{\pm ir^2} dr = \sqrt{\pi} e^{\pm i\pi/4}. $$

This is all the complex analysis we will need.

Before we continue, let us see indeed that the second integral above goes to zero as $R \to \infty$. Call that integral $I_2(R)$. So we want to show that $|I_2(R)|$ goes to zero, where $|z|$ is the modulus of a complex number $z$. Then

$$ |I_2(R)| = \left| iR \int_0^{\pi/4} e^{-R^2 e^{2i\theta}} e^{i\theta} d\theta \right| \le R \int_0^{\pi/4} e^{-R^2 \cos(2\theta)} d\theta $$

as $|e^{ia}| = 1$ for all $a \in \mathbb{R}$. So estimating the modulus of $I_2(R)$ reduces to Laplace's method!

In fact we are in case (ii): putting $2\theta = t$ we have

$$ |I_2(R)| \le \frac{R}{2} \int_0^{\pi/2} e^{R^2(-\cos t)} dt, $$

so here $x = R^2$, $h(t) = -\cos t$. $h(t)$ reaches a maximum at $t = \pi/2$, $h(\pi/2) = 0$ and $h'(\pi/2) = 1$, so we are in case (ii). Using formula (67) with suitable modifications (what are they?), we have

$$ |I_2(R)| \le \frac{R}{2} \left( \frac{1}{2R^2} + O(1/R^4) \right) $$

so $|I_2(R)| = O(1/R)$ and does go to zero as $R \to \infty$ as promised.

Now we proceed as in Laplace's method: we change variables in such a way that the resulting integral can be computed. Let us assume that $h(t)$ has a unique turning point at $t_0$. In the simplest case, $h''(t_0)$ is non-zero. As in Laplace's method, we introduce a new variable $s$ by

$$ h(t) - h(t_0) = \pm s^2, $$

where we take the plus sign if $h''(t_0) > 0$ (minimum) and minus otherwise. Hence expanding $h(t)$ in Taylor series around $t_0$, we have

$$ t - t_0 = \left( \frac{2}{|h''(t_0)|} \right)^{1/2} s + O(s^2). $$

We also expand $g$:

$$ g(t) = g(t_0) + (t - t_0)g'(t_0) + O((t - t_0)^2) = g(t_0) + O(s). $$

Our interval $I$ only has contributions from a neighbourhood $I_1$ of the turning point, so we might as well write

$$ J(x) = \left( \frac{2}{|h''(t_0)|} \right)^{1/2} g(t_0) e^{ixh(t_0)} \int_{I_1} e^{\pm ixs^2} (1 + O(s)) ds. $$
Now set $r = s/\sqrt{x}$ and we have

$$J(x) \sim \left( \frac{2}{x|h''(t_0)|} \right)^{1/2} g(t_0) e^{ixh(t_0)} \int_{\mathbb{R}} e^{\pm ir^2} \left(1 + \frac{1}{\sqrt{x}} O(r)\right) dr.$$

But we already know how to compute these integrals.

So we have the method of stationary phase formula,

$$J(x) \sim \left( \frac{2\pi}{x|h''(t_0)|} \right)^{1/2} g(t_0) e^{i(xh(t_0) \pm \pi/4)} + O(1/x). \quad (68)$$

Let us see a quick application. Let us compute $J_1(x) = \int_0^2 \cos(xt(2-t)) dt$. This is the real part of $J(x)$ with $h(t) = t(2-t)$, which has a maximum at $t=1$ and $g(t) = 1$. The maximum means we need to take the minus sign, and our estimate is the real part of the above formula, i.e.

$$J_1(x) \sim \sqrt{\frac{\pi}{x}} \cos(x - \pi/4) + O(1/x).$$

Please have a good look at the Method of Steepest Descent too.

## 20 Bibliography

Here is the bibliography for the Asymptotics Course. I am mainly using [7] (for centre manifold theory), and [10], [13] and [17] (for asymptotic methods proper), but you can find all the other books and articles useful. Please contact me if you need guidance. Barenblatt [2] gives the beautiful proof of Pythagoras' theorem by dimensional analysis; the article of Berry [3] is where the concepts of super- and hyper-asymptotics are introduced. Hinch [12] is a very good book on asymptotic methods by a master, with many tricky examples. The book [18] by Rand and Ambruster is based on MACSYMA, which is the same as Maxima, which is included in Sage and available via sage:!maxima.

## References

[1] C. M. Anderson and J. F. Geer, "Power series expansions for the frequency of the limit cycle of the van der Pol equation," SIAM J. Appl. Math., **43**, 678-693 (1982).

[2] G. I. Barenblatt, *Scaling, Self-similarity, and Intermediate Asymptotics: Dimensional Analysis and Intermediate Asymptotics*, CUP, Cambridge 1996.

[3] M. V. Berry, "Asymptotics, superasymptotics, hyperasymptotics," in: *Proc. Workshop on Asymptotics Beyond all Orders*, ed. S. Tanveer, Plenum, New York (1991).
[4] R. Bellman, *Perturbation techniques in Mathematics, Physics and Engineering*, Holt, Rinehart and Winston, New York 1964.

[5] N. Bleistein and R. A. Handelsman, *Asymptotic Expansions of Integrals*, Dover Publications, New York 1986.

[6] T. J. I'a. Bromwich, *An Introduction to the Theory of Infinite Series*, MacMillan and Co., London 1908.

[7] J. Carr, *Applications of Centre Manifold Theory*, Springer-Verlag, New York 1981.

[8] S. N. Chow and J. K. Hale, *Methods of Bifurcation Theory*, Springer-Verlag, Berlin 1982.

[9] J. C. Gibbings, *Dimensional Analysis*, Springer-Verlag, New York 2011.

[10] J. Guckenheimer and P. Holmes, *Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields*, 2nd edition, Springer-Verlag, New York 1986.

[11] G. H. Hardy, *Divergent Series*, Clarendon Press, Oxford 1949.

[12] E. J. Hinch, *Perturbation Methods*, CUP, Cambridge 1991.

[13] D. W. Jordan and P. Smith, *Nonlinear Differential Equations*, OUP, Oxford 2007.

[14] J. P. Keener, *Principles of Applied Mathematics*, Addison-Wesley, Redwood City, CA. 1988.

[15] C. C. Lin and L. A. Segel, *Mathematics Applied to Deterministic Problems in the Natural Sciences*, McMillan Publishing Co., New York 1974.

[16] J. D. Murray, *Asymptotic Analysis*, Springer-Verlag, New York 1984.

[17] A. H. Nayfeh, *Problems in Perturbation*, John Wiley & Sons, New York 1985.

[18] R. H. Rand and D. Ambruster, *Perturbation Methods, Bifurcation Theory, and Computer Algebra*, Applied Mathematical Sciences 65, Springer-Verlag, New York 1987.

[19] J. G. Simmonds and J. E. Mann, Jr. *A First look at Perturbation Theory*, R. E. Krieger Publishing Company, Malabar, Florida 1986.

[20] F. Verhulst, *Nonlinear Differential Equations and Dynamical Systems*, Springer-Verlag, Berlin 1990.