Part C · the analytic toolkit · spine

Perturbation theory: the series that diverges and is still right

Almost nothing has an exact answer. What it has instead is an expansion in a small parameter — often divergent, usually non-uniform, and, handled correctly, accurate to ten digits.

needs 06, 02, 08~3 h slides + 6 h problemsproblems: ./mr problems 21

Unit 02 said it early and plainly: almost no equation of motion has a closed-form solution, and that is not a gap in the literature. Unit 02's answer was the qualitative picture. This unit's answer is the quantitative one. When a problem is close to one you can solve — a spring that is nearly linear, an orbit that is nearly Keplerian, a fluid that is nearly inviscid — you can expand in the smallness and get numbers. Most of the quantitative physics you will ever read is done this way, and most of it rests on three facts that are rarely said aloud:

the series you get usually diverges, and that does not matter, because what you want is not convergence but a controlled error at a fixed number of terms;
the naive expansion is usually non-uniform: correct for a while, or away from some edge, and then badly wrong — and the place where it goes wrong is where the physics is;
the small parameter is not handed to you: you find it by nondimensionalising, which is unit 00's dimensional habit put to work.

How to read this deck

Two calculations run through it, joined on purpose. The first is the Duffing oscillator $\ddot x + x + \varepsilon x^3 = 0$: the naive expansion produces a term $t\sin t$ that wrecks it after a time of order $1/\varepsilon$, and we learn to remove it. The second is the orbit equation $w'' + w = 1 + \varepsilon w^2$ that general relativity hands to unit 37: the same calculation produces the same kind of term — and there it must not be thrown away, because it is the advance of Mercury's perihelion. Sections 1 and 2 are about what an expansion means; 3 and 4 about secular terms; 5 about boundary layers; 6 about WKB; 7 closes the orbit.

Every number quoted was computed. Proofs are on the slides where they are short and in these notes under Full proof where they are not; anything imported is labelled cited with where its proof lives. Integral asymptotics — Laplace's method, stationary phase — belong to unit 20 and are cited here, never used as a step.

21 · overview

Two oscillators, one term

the Duffing oscillator · the term is removed $$\ddot x + x + \varepsilon x^3 = 0, \qquad x(0) = A,\ \dot x(0) = 0 .$$

1. Expand $x = x_0 + \varepsilon x_1 + \cdots$: $x_0 = A\cos t$, and $x_1$ is driven by $-x_0^3$, which contains $-\tfrac34A^3\cos t$ — a push at the oscillator's own frequency.

2. Resonance puts $-\tfrac38A^3\,t\sin t$ into $x_1$: unbounded, though the true $|x| \le A$ for ever.

3. That term is a frequency shift, Taylor-expanded: $A\cos\big((1 + \varepsilon\omega_1)t\big) = A\cos t - \varepsilon\omega_1At\sin t + \cdots$

4. Give the frequency its own expansion and demand periodicity: $\omega = 1 + \tfrac38\varepsilon A^2$.

the relativistic orbit · the term is kept $$w'' + w = 1 + \varepsilon w^2, \qquad w = p/r .$$

$r$ is the planet's distance from the Sun at angle $\theta$, $p$ a fixed length of its orbit, $' = \mathrm{d}/\mathrm{d}\theta$.

1. At $\varepsilon = 0$, $w_0 = 1 + e\cos\theta$: unit 08's ellipse, of eccentricity $e$.

2. $w_0^2$ contains $2e\cos\theta$, again a push at the natural frequency: $w_1$ gets $e\,\theta\sin\theta$.

3. Read as in Duffing's step 3, $w \approx 1 + e\cos\big((1 - \varepsilon)\theta\big)$: perihelion comes round $2\pi\varepsilon$ late each orbit.

4. For Mercury $\varepsilon = 8.0\times10^{-8}$: $43''$ a century, observed. Here the term is the answer.

The same mathematics, opposite verdicts. The plan. 1 finding $\varepsilon$ · 2 asymptotic series · 3 secular terms · 4 removing them · 5 boundary layers · 6 WKB · 7 the orbit. A dotted underline is a glossary term (click, or G); a solid one like the Duffing oscillator is the running example (hover); N shows the notes.

Where each step is done. Duffing: step 1 is Regular perturbation, the naive way; step 2 is Worked: the naive first order; step 3 is What the secular term really is, which draws it as a point that leaves its circle along the tangent; step 4 is Worked: Poincaré–Lindstedt for Duffing, using the condition proved on The solvability condition, and Checking the frequency tests the answer against the exact one. The orbit: the equation and its small parameter are The relativistic orbit equation, nondimensionalised, where $p$ (the semi-latus rectum of the Newtonian ellipse, $a(1 - e^2)$ for semi-major axis $a$) and $\varepsilon = 3GM/c^2p$ are defined; steps 1–2 are Worked: the first-order orbit; step 3 is Reading the secular term: the perihelion advances, which draws the rosette; step 4 is Mercury: 43 seconds of arc a century. Remove it, or keep it? then says why one term goes and the other stays — the question the whole unit is built to answer.

Why these two. The Duffing oscillator is the simplest nonlinear oscillator there is, and the textbook place to see a secular term; physically it is a stiffening (or, with $\varepsilon < 0$, softening) spring, and the pendulum to its first nonlinear order is one of them. The orbit equation is the reason this unit is on the course's spine: unit 37 computes the perihelion advance of Mercury from exactly this equation, and a reader who has only seen secular terms "removed" will find that calculation looks like a trick. It is not a trick, and this unit is where that is made honest.

The course's order matters here. We use unit 00 (dimensions, the Buckingham theorem, quadrature for a conservative force), unit 01 (Grönwall's inequality), unit 02 (conserved energy, centres, the van der Pol limit cycle), unit 03 (partial derivatives and the chain rule), unit 06 (linear constant-coefficient equations, resonance, variation of parameters, the complex exponential) and unit 08 (the Binet equation). Nothing else.

section 1

Small parameters, and where they come from

"Small" means nothing until you say compared with what. Nondimensionalising says it.

21 · small parameters

Nondimensionalise first

A mass on a spring with a cubic stiffening: $m\ddot X + kX + \beta X^3 = 0$, released from rest at $X = a$. Six quantities, $X, t, m, k, \beta, a$, in three dimensions M, L, T.

nondimensionalisation

Choose a scale for each variable from the parameters, so that the scaled variables are of size one: here $x = X/a$ and $t' = t\sqrt{k/m}$. Divide by one coefficient. What is left over is a small parameter or it is not.

$$ka\,x'' + ka\,x + \beta a^3x^3 = 0,\ \text{ so}$$ $$x'' + x + \varepsilon x^3 = 0,\qquad \varepsilon = \frac{\beta a^2}{k}.$$

What $\varepsilon$ says. $\beta a^3/ka$: the cubic force over the linear one, at the amplitude. "Small" now has a meaning.

Buckingham, as a check

Six quantities, dimension matrix of rank 3, so three groups (unit 00): $X/a$, $t\sqrt{k/m}$ and $\beta a^2/k$. The scaled equation contains exactly those. Whatever the answer is, it is a function of $\varepsilon$ alone once time is measured in units of $\sqrt{m/k}$.

the pendulum is one of these

$\ddot\theta + \sin\theta = 0$ with $\sin\theta = \theta - \theta^3/6 + O(\theta^5)$. Put $\theta = A x$: $x'' + x - \tfrac{A^2}{6}x^3 = 0$. So Duffing with $\varepsilon = -A^2/6$, amplitude 1.

The scaled equation, line by line. With $X = a x$ and $t = \sqrt{m/k}\,t'$, the chain rule gives $\ddot X = a\,(k/m)\,x''$ where $'$ is $\mathrm{d}/\mathrm{d}t'$. Substituting, $ka\,x'' + ka\,x + \beta a^3x^3 = 0$; divide by $ka$. The initial conditions become $x(0) = 1$, $x'(0) = 0$. In the running example we keep the amplitude $A$ visible (scale by some fixed length instead of by the amplitude) because it makes the answer's amplitude dependence easy to read: then $\varepsilon = \beta L^2/k$ for the fixed length $L$, and the combination that matters is always $\varepsilon A^2$.

Why the choice of scales matters. The scales must be the ones the solution actually has. If you scale $X$ by a length much bigger than the amplitude, every term looks small and the "small parameter" you find is an artefact. The rule: each scaled variable, and each of its derivatives that appears, should be of size one in the region you care about. That rule is also how boundary layers are found (section 5): a region where the derivatives are not of size one in the original scaling needs its own.

Dimensions. $[\beta] = \mathsf{MT^{-2}L^{-2}}$ (force over length cubed), $[k] = \mathsf{MT^{-2}}$, so $[\beta a^2/k] = 1$.

The pendulum. Unit 00 expanded the exact period to $T \approx 2\pi\sqrt{l/g}\,(1 + \theta_0^2/16)$. Section 4 recovers that factor from the Duffing frequency with $\varepsilon = -A^2/6$ — the first check that the perturbation theory is right.

21 · small parameters

Regular or singular: the test

two kinds of problem

Set $\varepsilon = 0$. If the reduced problem is of the same kind — same degree, same order, all conditions still imposable — and its solution is the limit of the true one, the problem is a regular perturbation: expand in powers of $\varepsilon$ and go. If not, it is a singular perturbation.

the tell

$\varepsilon$ multiplying the highest power or derivative. Setting it to zero lowers the degree or the order: roots or boundary conditions are lost, and something must happen on a scale where the "small" term is not small.

$\varepsilon x^2 + x - 1$ for $\varepsilon = 0.2$ and $0.1$, and the line it becomes at $\varepsilon = 0$. One root stays near $1$; the other sits near $-1/\varepsilon$ and runs off to $-\infty$ as $\varepsilon \to 0$ — the root the reduced problem cannot see.

What to look for in the figure. The parabola's vertex is at $x = -1/2\varepsilon$, depth $-1 - 1/4\varepsilon$; halving $\varepsilon$ doubles the distance to the far root and the parabola flattens towards the straight line $x - 1$. Near $x = 1$ the curves are almost indistinguishable: there $\varepsilon x^2$ really is small, and the root moves smoothly with $\varepsilon$ — regular. Near the far root $\varepsilon x^2$ is as big as $x$, so dropping it is never justified there — singular. The figure is computed from the exact roots, $(-1 \pm \sqrt{1 + 4\varepsilon})/2\varepsilon$.

Why "same kind" is the right test. A quadratic has two roots, a linear equation one; a second-order differential equation takes two boundary conditions, a first-order one only one. If setting $\varepsilon = 0$ changes that count, the $\varepsilon = 0$ problem cannot be the limit of the full one everywhere, and a plain power series in $\varepsilon$ must miss something.

21 · small parameters · the smallest example

Where the lost root went: dominant balance

the smallest example: $\varepsilon x^2 + x - 1 = 0$

At $\varepsilon = 0$: one root, $x = 1$. But a quadratic has two. Put $x = 1 + \varepsilon x_1 + \varepsilon^2x_2$ and collect: $x_1 = -1$, $x_2 = 2$. The other root is lost — it is at $x \approx -1/\varepsilon$, where "$\varepsilon x^2$ is small" is false. Rescale $x = X/\varepsilon$: $X^2 + X - \varepsilon = 0$, regular, $X = -1 - \varepsilon + \dots$

At $\varepsilon = 0.1$: exact $0.916080$ and $-10.916080$; the series through $\varepsilon^2$, $1 - \varepsilon + 2\varepsilon^2$ and $-1/\varepsilon - 1 + \varepsilon - 2\varepsilon^2$ (notes), give $0.920$ and $-10.920$.

Finding the scale of the lost root is dominant balance: try $x = \varepsilon^{-\alpha}X$, and choose $\alpha$ so that two terms are of the same size and the rest are smaller. Here $\alpha = 1$ balances $\varepsilon x^2$ against $x$.

Duffing passes the test — $\varepsilon$ multiplies $x^3$, not $\ddot x$ — and is regular on any fixed time interval. Section 3 shows it is not regular on $[0, \infty)$: a second kind of singularity, in the domain rather than the equation.

The regular root, in full. Substitute $x = 1 + \varepsilon x_1 + \varepsilon^2 x_2 + O(\varepsilon^3)$: $\varepsilon(1 + 2\varepsilon x_1) + 1 + \varepsilon x_1 + \varepsilon^2x_2 - 1 + O(\varepsilon^3) = 0$. Order $\varepsilon$: $1 + x_1 = 0$. Order $\varepsilon^2$: $2x_1 + x_2 = 0$, so $x_2 = 2$. Hence $x = 1 - \varepsilon + 2\varepsilon^2 + O(\varepsilon^3)$. Exact: $x = (-1 + \sqrt{1 + 4\varepsilon})/2\varepsilon$, and the binomial series of $\sqrt{1+4\varepsilon}$ gives $1 - \varepsilon + 2\varepsilon^2 - 5\varepsilon^3 + \dots$ (the Catalan numbers, with signs), agreeing.

The singular root. The two roots sum to $-1/\varepsilon$, so the second is $-1/\varepsilon - 1 + \varepsilon - 2\varepsilon^2 + \dots$ — a series that starts at order $\varepsilon^{-1}$, which no power series in $\varepsilon$ starting at $\varepsilon^0$ can represent. That is what "singular" means here: the expansion needs a different gauge.

Dominant balance, spelled out. With $x = \varepsilon^{-\alpha}X$, the three terms are of sizes $\varepsilon^{1-2\alpha}$, $\varepsilon^{-\alpha}$ and $1$. The possible balances: $\varepsilon^{1-2\alpha} = \varepsilon^{-\alpha}$ gives $\alpha = 1$, with the third term $1 \ll \varepsilon^{-1}$ — consistent, the singular root. $\varepsilon^{-\alpha} = 1$ gives $\alpha = 0$, with $\varepsilon x^2 = O(\varepsilon)$ smaller — consistent, the regular root. $\varepsilon^{1-2\alpha} = 1$ gives $\alpha = 1/2$, but then the middle term is $\varepsilon^{-1/2} \gg 1$, larger than the two that balanced — inconsistent, discarded. Each consistent balance is one family of solutions. The same bookkeeping locates boundary layers in section 5.

section 2

Asymptotic series

What "$\approx$" means, precisely — and why a divergent series can deserve it.

21 · asymptotic series

Order notation

order notation

Let $f, g$ be defined near $x_0$ (possibly $x_0 = \infty$), $g \neq 0$ near $x_0$. As $x \to x_0$:

$f = O(g)$ if $|f/g|$ is bounded on some punctured neighbourhood of $x_0$;
$f = o(g)$ if $f/g \to 0$;
$f \sim g$, asymptotically equivalent, if $f/g \to 1$, that is, $f - g = o(g)$.

as $x \to 0$

$\sin x = x + O(x^3)$; $\ \sin x \sim x$; $\ 1 - \cos x = o(x)$; $\ x^{10} = o(x^9)$.

$e^{-1/x} = o(x^n)$ as $x \to 0^+$, for every $n$: smaller than every power.

as $x \to \infty$

$\ln x = o(x^\delta)$ for every $\delta \gt 0$; $\ x^2 + x \sim x^2$; $\ x^2 + x \sim x^2 + 1$, although their difference is unbounded.

Every statement is about a limit, never about a particular $x$. That is both the power and the trap: "$O(\varepsilon^2)$" says nothing about the constant, and nothing about $\varepsilon = 0.3$.

Proof that $e^{-1/x} = o(x^n)$ as $x \to 0^+$. Put $s = 1/x \to \infty$. Then $e^{-1/x}/x^n = s^n e^{-s}$, and $e^s \ge s^{n+1}/(n+1)!$ (one term of its series, all terms positive), so $s^ne^{-s} \le (n+1)!/s \to 0$. $\square$

Proof that $\ln x = o(x^\delta)$. Put $x = e^u$: $\ln x/x^\delta = u e^{-\delta u} \le 2/(\delta^2 u) \to 0$ by the same one-term bound $e^{\delta u} \ge (\delta u)^2/2$. $\square$

The last example. $(x^2 + x)/(x^2 + 1) \to 1$, so $x^2 + x \sim x^2 + 1$, although $(x^2 + x) - (x^2 + 1) = x - 1$ is unbounded. $\sim$ controls the relative error only. When an absolute error matters — a phase, which enters as $\cos(\text{phase})$ — you need more, and that distinction is the whole of section 3.

Notation warning. $f = O(g)$ is a one-way statement dressed as an equation: $x^3 = O(x^2)$ as $x \to 0$ does not make $O(x^2) = x^3$. Read "$=$" as "is in the class".

21 · asymptotic series

Poincaré's definition

asymptotic sequence and asymptotic expansion

$(\phi_n)_{n \ge 0}$ is an asymptotic sequence as $x \to x_0$ if $\phi_{n+1} = o(\phi_n)$ for every $n$. Then $f \sim \sum_{n\ge0} a_n\phi_n$ means

$$f - \sum_{n=0}^{N} a_n\phi_n = o(\phi_N)\ \text{ for every } N .$$

Poincaré, 1886. The series need not converge: the definition is about a fixed $N$ as $x \to x_0$, not about $N \to \infty$ at fixed $x$.

uniqueness of the coefficients

If $f$ has an expansion in $(\phi_n)$, its coefficients are $a_N = \lim_{x\to x_0}\big(f - \sum_{n\lt N} a_n\phi_n\big)/\phi_N$.

Why: divide the definition at $N$ by $\phi_N$. $\square$

but not of the function

$e^{-1/x} \sim 0 + 0\cdot x + 0\cdot x^2 + \cdots$ as $x \to 0^+$. So $f$ and $f + e^{-1/x}$ have the same expansion in powers. An expansion determines the coefficients, never the function.

Duffing's frequency, as an expansion

$\lambda = \varepsilon A^2$	$\omega$ exact	$\omega - (1 + \tfrac38\lambda)$	$\div\lambda^2$
0.2	1.072000	$-3.0\times10^{-3}$	$-0.0750$
0.1	1.036717	$-7.8\times10^{-4}$	$-0.0783$
0.01	1.003742	$-8.2\times10^{-6}$	$-0.0816$

The error is $O(\lambda^2)$, so $1 + \tfrac38\lambda$ is the two-term expansion — and the last column is converging to the next coefficient, $-21/256 = -0.0820$, which the problem sheet derives by carrying section 4's method one order further.

Where the "exact" column comes from. Energy is conserved (unit 00): $\tfrac12\dot x^2 + \tfrac12x^2 + \tfrac14\varepsilon x^4 = \tfrac12A^2 + \tfrac14\varepsilon A^4$. So the period is a quadrature, $T = 4\int_0^A \mathrm{d}x/\sqrt{(A^2 - x^2) + \tfrac12\varepsilon(A^4 - x^4)}$. Substituting $x = A\sin\varphi$ and using $A^4 - x^4 = A^4\cos^2\varphi\,(1 + \sin^2\varphi)$,

$$T = 4\int_0^{\pi/2}\frac{\mathrm{d}\varphi}{\sqrt{1 + \tfrac12\lambda(1 + \sin^2\varphi)}},\qquad \lambda = \varepsilon A^2 ,$$

a proper integral of a smooth function, evaluated numerically to twelve digits; $\omega = 2\pi/T$. It depends on $\varepsilon$ and $A$ only through $\lambda$, as the scaling slide predicted.

The second coefficient $-21/256$ is derived on the problem sheet ("Duffing to second order") by carrying Poincaré–Lindstedt one order further. The table is a check on that derivation, not a substitute for it: $-0.0816$ at $\lambda = 0.01$ differs from $-0.0820$ by an amount of order $\lambda$, as the next term predicts.

Different sequences, different coefficients. Uniqueness is relative to the sequence. As $x \to 0$, $1/(1-x) = \sum x^n$, and also $= \sum (1+x)x^{2n}$, an expansion in the asymptotic sequence $\phi_n = (1+x)x^{2n}$. The problem sheet proves both uniqueness statements.

History. Poincaré gave the definition in 1886 (Acta Mathematica 8), motivated by exactly this unit's subject: the series of celestial mechanics, which astronomers used with great success and which he could prove diverge. Stieltjes, the same year, studied the series of the next slide.

21 · asymptotic series

Convergent is not the same as useful

convergent, and useless

$e^{-20} = \sum_n (-20)^n/n!$ converges. Its terms climb to $4.31\times10^{7}$ (at $n = 19, 20$) before they fall, and they must cancel to $2.06\times10^{-9}$.

You need terms up to $n = 74$ for $0.1\%$ accuracy. Summed in double precision, 121 terms give anything from $0.7\times10^{-9}$ to $6\times10^{-9}$, depending on how the terms are rounded — never the right answer, however many terms.

divergent, and excellent

$\displaystyle\sum_n (-1)^n n!\,x^n$ diverges for every $x \ne 0$: the ratio of terms is $(n+1)x \to \infty$.

Yet at $x = 0.04$ its first twenty-five terms give the function it expands to an error of $8.6\times10^{-11}$ — ten correct digits.

The next slide puts the two side by side, term by term.

The $e^{-20}$ numbers. The term ratio is $|t_{n+1}/t_n| = 20/(n+1)$, which exceeds 1 until $n = 19$, so the largest terms are $t_{19}$ and $t_{20}$, of size $20^{20}/20! = 4.310\times10^{7}$. In exact arithmetic the partial sums first get within $0.1\%$ of $e^{-20} = 2.0612\times10^{-9}$ at $n = 74$. In IEEE double precision each term near the peak carries a rounding error of order $10^{7}\times 10^{-16} = 10^{-9}$ — the size of the answer — so the computed sum is decided by the rounding. Three natural ways, all with $n = 0, \dots, 120$: each term $(-20)^n/n!$ correctly rounded and then summed exactly gives $7.2\times10^{-10}$; the same terms added left to right give $1.48\times10^{-9}$; terms built by the recurrence $t_n = t_{n-1}\cdot(-20/n)$ and added left to right give $6.1\times10^{-9}$. The right way is $1/e^{20}$, a series of positive terms: the problem is cancellation, not convergence.

The divergence. For $x \ne 0$ the terms $n!\,|x|^n$ do not tend to zero (their ratio tends to infinity), so the series diverges by the term test.

The general lesson is the one Poincaré drew: convergent series of celestial mechanics could need hundreds of terms, while divergent ones gave the answer in three. Physicists had known this empirically for a century; the definition on the previous slide made it mathematics.

21 · asymptotic series

The two series, term by term

$|$term $n|$ on a log scale. The convergent series must cancel terms $10^{16}$ times larger than its answer (dashed). The divergent one reaches $10^{-10}$ before it turns: stop there.

Convergence is a statement about $N \to \infty$ at fixed $x$. Usefulness is a statement about small $N$ at the $x$ you have. Neither implies the other.

An asymptotic series is one whose first few terms are good when $x$ is near $x_0$, whatever the tail does. In the picture: what matters is how low the curve gets before it turns, not where it ends.

How the figure was computed. $\log_{10}$ of $20^n/n!$ and of $n!\,(0.04)^n$, from $\ln n!$ by the log-gamma function, for $n = 0, \dots, 80$; every second point is marked. The first peaks at $n = 19, 20$ at $4.31\times10^7$; the second bottoms out at $n = 25$ at $1.75\times10^{-10}$, and keeps rising after $n = 25$ because each new term is $(n+1)\times0.04$ times the last.

Reading it. For the convergent series, the answer $2.06\times10^{-9}$ (dashed) lies sixteen decades below the largest term, so it emerges only after nearly all of the sum has cancelled — which is why finite precision ruins it. For the divergent series, the error after stopping at the smallest term is about half that term (Optimal truncation makes this precise), so the lowest point of the curve is the accuracy you can have.

21 · asymptotic series

The Stieltjes series, and its remainder

the Stieltjes function $$S(x) = \int_0^\infty \frac{e^{-t}}{1 + xt}\,\mathrm{d}t, \qquad x \ge 0 .$$

For $y \ne -1$, the finite geometric sum is exact: $\dfrac1{1+y} = \displaystyle\sum_{n=0}^{N-1}(-y)^n + \frac{(-y)^N}{1+y}$. Put $y = xt$, multiply by $e^{-t}$, integrate, and use $\int_0^\infty t^ne^{-t}\,\mathrm{d}t = n!$:

$$S(x) = \sum_{n=0}^{N-1}(-1)^n n!\,x^n + R_N(x),$$ $$R_N(x) = (-x)^N\!\int_0^\infty\!\frac{t^Ne^{-t}}{1+xt}\,\mathrm{d}t .$$

the error is smaller than the first omitted term

For $x \gt 0$: $R_N(x)$ has the sign of $(-1)^N$ and $0 \lt |R_N(x)| \lt N!\,x^N$.

Why: $0 \lt 1/(1+xt) \lt 1$ for $t \gt 0$, so $0 \lt \int_0^\infty t^Ne^{-t}/(1+xt) \lt N!$. $\square$

So $|R_N| \lt N!\,x^N = o(x^{N-1})$: the series is the asymptotic expansion of $S$ as $x \to 0^+$, with an explicit error bound — a rare luxury.

$S$ is not exotic: $xS(x)$ solves Euler's equation $x^2y' + y = x$, whose only formal power-series solution is $x$ times this divergent series. The problem sheet shows it.

Full proof of the remainder formula. For every real $y \ne -1$ and $N \ge 1$, $(1 + y)\sum_{n=0}^{N-1}(-y)^n = 1 - (-y)^N$ (telescoping), which rearranges to the identity on the slide. With $y = xt \ge 0$, multiply by $e^{-t}$ and integrate over $[0, \infty)$. Each term is a convergent improper integral (unit 00): $\int_0^\infty t^ne^{-t}\,\mathrm{d}t = n!$ by induction, integrating by parts once, $\int t^ne^{-t} = n\int t^{n-1}e^{-t}$, with the boundary terms vanishing because $t^ne^{-t} \to 0$. The remainder integral converges because its integrand is at most $t^Ne^{-t}$. Nothing was interchanged: the sum is finite. That is the whole difference between this derivation and the illegitimate one — expand $1/(1+xt)$ in its full geometric series and integrate term by term — which fails because the geometric series diverges for $t \gt 1/x$.

Dimensions and a limit. Everything is dimensionless. At $x = 0$, $S(0) = \int_0^\infty e^{-t} = 1$, matching the first term.

Why the sign matters. The error has the sign of the first omitted term and is smaller than it. So consecutive partial sums bracket $S$ — you get a rigorous interval, not just an estimate. This is special to series of Stieltjes type; in general an asymptotic series gives no such guarantee, and the only honest error estimate is "about the size of the first omitted term", which is a heuristic.

21 · asymptotic series · example

Optimal truncation

worked example 1 — $S(0.2)$, term by term

$n$	term $(-1)^nn!\,0.2^n$	sum to $n$	minus $S$
0	1	1	$+0.148$
1	$-0.2$	0.8	$-0.052$
2	0.08	0.88	$+0.028$
3	$-0.048$	0.832	$-0.020$
4	0.0384	0.8704	$\mathbf{+0.018}$
5	$-0.0384$	0.832	$-0.020$
6	0.04608	0.87808	$+0.026$
7	$-0.064512$	0.813568	$-0.039$
8	0.10322	0.91679	$+0.065$

$S(0.2) = 0.852111$. Terms shrink while $(n+1)x \lt 1$, then grow. Stop at the smallest term: through $n = 4$ the sum is off by $0.018$, less than the first omitted term $0.0384$, and the correction still needed, $-0.018$, has that term's sign.

optimal truncation, and its error

Truncate at $N = \lfloor 1/x\rfloor$, where the terms are smallest. Then

$$|R_N(x)| \le \frac{e^2}{x}\,e^{-1/x}.$$

The error is exponentially small — smaller than every power of $x$ — though no finite truncation has error zero.

Reading the table. The partial sum through $n$ is $S_{n+1}$ in the remainder notation (it has $n+1$ terms), so "minus $S$" is $-R_{n+1}$. The theorem says $R_{n+1}$ has the sign of the first omitted term $t_{n+1}$ and is smaller than it. Row $n = 4$: $S - (\text{sum}) = -0.0183$, and the first omitted term is $t_5 = -0.0384$: same sign, smaller, exactly as proved. Every row obeys the same rule; check two more before moving on.

Full proof of the exponentially small bound. For $N \ge 1$, $\ln k \le \int_k^{k+1}\ln s\,\mathrm{d}s$ for $k = 1, \dots, N-1$, so $\ln N! \le \int_1^N \ln s\,\mathrm{d}s + \ln N = N\ln N - N + 1 + \ln N$, that is, $N! \le e\,N^{N+1}e^{-N}$. With $N = \lfloor 1/x\rfloor$ we have $Nx \le 1$ and $N \gt 1/x - 1$, so

$$|R_N| \lt N!\,x^N \le e\,N\,(Nx)^N e^{-N} \le e\cdot\tfrac1x\cdot e^{-(1/x - 1)} = \tfrac{e^2}{x}e^{-1/x}.\qquad\square$$

At $x = 0.1$ this bound is $3.4\times10^{-3}$ against a true error of $1.8\times10^{-4}$: crude in the constant, right in the exponential, which is what matters.

21 · asymptotic series

Optimal truncation, plotted

$|S_N - S(x)|$ after $N$ terms of the Stieltjes series. Each curve is a V with its point at $N = 1/x$ (ringed): the smaller $x$, the later the turn and the deeper the V.

The best errors. $x = 0.2$: $1.8\times10^{-2}$ at $N = 5$. $\ x = 0.1$: $1.8\times10^{-4}$ at $N = 10$. $\ x = 0.04$: $8.6\times10^{-11}$ at $N = 25$.

Each is about half the smallest term, which Stirling's formula (cited, unit 20) puts at $\sqrt{2\pi/x}\,e^{-1/x}$. The depth of the V goes like $e^{-1/x}$: exponentially small, as the theorem said.

Right of the point, every added term makes the answer worse.

How the figure was computed. $S(x)$ by numerical quadrature of its integral to fifteen digits, the partial sums exactly, $N = 1, \dots, 40$. The error alternates in sign (the remainder theorem), so the log of its absolute value is plotted. At $x = 0.2$ the minimum is at $N = 5$ ($0.0183$), with $N = 6$ close behind ($0.0201$), since the terms $n = 4$ and $n = 5$ are equal in size ($0.0384$); the table on the previous slide is the left end of the $x = 0.2$ curve.

The sharp constant, cited. Stirling's formula $N! \sim \sqrt{2\pi N}\,N^Ne^{-N}$ is proved in unit 20 by Laplace's method; it gives the smallest term as $\sqrt{2\pi/x}\,e^{-1/x}$. The "about half" is an observation from these computed numbers — at $x = 0.04$ the best error is $8.65\times10^{-11}$ against $\sqrt{\pi/2x}\,e^{-1/x} = 8.70\times10^{-11}$ — not a theorem proved here.

The lesson for the whole unit. For a given $x$ there is a best accuracy the series can deliver, and adding terms beyond the smallest makes things worse. In physics the small parameter is fixed by nature — $\varepsilon = 7.99\times10^{-8}$ for Mercury — so that best accuracy is what you actually get.

21 · asymptotic series

Borel summation, named

Borel summation

Given $\sum a_nx^n$, form $B(t) = \sum a_n t^n/n!$ — often convergent where the original is not. If $B$ extends to $[0, \infty)$ and the integral exists, define

$$\text{sum} := \int_0^\infty e^{-t}\,B(xt)\,\mathrm{d}t .$$

The motive: $n! = \int_0^\infty t^ne^{-t}\,\mathrm{d}t$, so formally $\sum a_n x^n = \int_0^\infty e^{-t}\sum a_n(xt)^n/n!\,\mathrm{d}t$ — the swap being exactly the step that is not allowed.

on the Stieltjes series

$a_n = (-1)^nn!$ gives $B(t) = \sum(-t)^n = 1/(1+t)$ for $|t| \lt 1$, and the formula $1/(1+t)$ is defined on all of $[0, \infty)$. Then $\int_0^\infty e^{-t}/(1 + xt)\,\mathrm{d}t = S(x)$ — the function, recovered exactly.

named, not developed

When Borel's sum exists and when it equals "the" function is a theorem with hypotheses (Watson's, cited: Bender and Orszag, ch. 8; Hardy, Divergent Series, for the theory). This course uses optimal truncation instead, which needs nothing.

Why name it at all. Because the question it answers — can a divergent series determine its function? — has the answer "not by itself" (the $e^{-1/x}$ slide) and "yes, with extra analytic information" (Borel). Quantum field theory's perturbation series are believed to be divergent and of this factorial type, and Borel summation is the language in which their meaning is discussed. You will see the name; now you know what it is.

Why not develop it. Making "$B$ extends to $[0,\infty)$" precise wants analytic continuation, which is unit 20's complex analysis. For the Stieltjes series we side-stepped that: the formula $1/(1+t)$ is simply defined on $[0, \infty)$, and the integral is the function we started from.

A check on the definition. For a convergent geometric series $\sum x^n$ ($|x| \lt 1$), $B(t) = e^t$ and $\int_0^\infty e^{-t}e^{xt}\,\mathrm{d}t = 1/(1 - x)$, the ordinary sum. Borel summation agrees with ordinary summation where both apply — a property called regularity, which holds in general (cited: Hardy, Divergent Series).

section 3

Secular terms

The naive expansion of the Duffing oscillator, and the term that kills it after a time $1/\varepsilon$.

21 · secular terms

Regular perturbation, the naive way

Treat Duffing as a regular problem: substitute

$$x(t) = x_0(t) + \varepsilon x_1(t) + \varepsilon^2x_2(t) + \cdots$$

into $\ddot x + x + \varepsilon x^3 = 0$ and into the initial conditions, and set each power of $\varepsilon$ to zero separately. This is the perturbation hierarchy:

$$\begin{aligned} \varepsilon^0:&\quad \ddot x_0 + x_0 = 0,\\ \varepsilon^1:&\quad \ddot x_1 + x_1 = -x_0^3,\\ \varepsilon^2:&\quad \ddot x_2 + x_2 = -3x_0^2x_1, \end{aligned}$$

with $x_0(0) = A$, $\dot x_0(0) = 0$, and zero data for every $x_n$, $n \ge 1$.

the structure to notice

Every order is the same linear operator $L = \mathrm{d}^2/\mathrm{d}t^2 + 1$, driven by a known right-hand side built from lower orders. One nonlinear problem has become a sequence of forced linear oscillators: unit 06, over and over.

The one identity we need, from $\cos^3 t = \big(\tfrac{e^{it} + e^{-it}}{2}\big)^3$:

$$\cos^3 t = \tfrac34\cos t + \tfrac14\cos 3t .$$

So the $\varepsilon^1$ equation is an undamped oscillator driven at its own frequency, $-\tfrac34A^3\cos t$, plus a harmless $-\tfrac14A^3\cos 3t$.

Why the hierarchy is legitimate. The solution $x(t;\varepsilon)$ depends smoothly on $\varepsilon$ for each fixed $t$, because the right-hand side is a polynomial in $(x, \varepsilon)$. Differentiating the equation in $\varepsilon$ at $\varepsilon = 0$ gives precisely the equations above, with $x_n = \frac{1}{n!}\partial_\varepsilon^n x|_{\varepsilon=0}$. So the $x_n$ are the Taylor coefficients of $x$ in $\varepsilon$ at fixed $t$. Everything about to go wrong goes wrong in the dependence on $t$, not in the existence of the coefficients. (Smooth dependence on a parameter is the one piece here not proved in unit 01, which proved continuous dependence; it is cited: Hartman, Ordinary Differential Equations, ch. V, on dependence on initial conditions and parameters. What we need is the weaker statement proved directly, by Grönwall, on "How long the naive answer is good".)

The identity. Expanding the cube, $(e^{it} + e^{-it})^3 = e^{3it} + 3e^{it} + 3e^{-it} + e^{-3it} = 2\cos 3t + 6\cos t$; divide by 8. Unit 06's complex exponential, used as algebra.

Dimensions. In scaled variables everything is dimensionless; in physical ones each order carries the same units as $X$, and $\varepsilon$ carries none — which is why the hierarchy can add them.

21 · secular terms · example

Worked: the naive first order

worked example 2 — solve the $\varepsilon^1$ equation

Step 1 — the order-zero solution. $x_0 = A\cos t$.

Step 2 — the forcing. $-x_0^3 = -\tfrac34A^3\cos t - \tfrac14A^3\cos 3t$.

Step 3 — the off-resonant part. Try $c\cos 3t$: $(-9 + 1)c = -\tfrac14A^3$, so $c = \tfrac1{32}A^3$.

Step 4 — the resonant part. $c\cos t$ fails ($L\cos t = 0$): unit 06's resonance. Try $b\,t\sin t$: $L(t\sin t) = 2\cos t$, so $2b = -\tfrac34A^3$, $b = -\tfrac38A^3$.

Step 5 — the initial conditions. Add the homogeneous $d\cos t$ (the $\sin t$ part is zero since $\dot x_1(0) = 0$ already holds). $x_1(0) = \tfrac1{32}A^3 + d = 0$, so $d = -\tfrac1{32}A^3$.

$$x_1 = \underbrace{-\tfrac38A^3\,t\sin t}_{\text{secular}} + \tfrac1{32}A^3(\cos3t - \cos t).$$

Step 6 — check. $\dot x_1(0) = 0 + \tfrac1{32}A^3(0 - 0) = 0$ ✓, and $L x_1 = -\tfrac34A^3\cos t + \tfrac1{32}A^3(-8\cos 3t) = -x_0^3$ ✓.

The term $t\sin t$ is a secular term: it grows without bound, and it comes from forcing a linear oscillator at its own frequency.

Step 4 in full. $\frac{\mathrm{d}}{\mathrm{d}t}(t\sin t) = \sin t + t\cos t$, and $\frac{\mathrm{d}^2}{\mathrm{d}t^2}(t\sin t) = 2\cos t - t\sin t$. Adding $t\sin t$ gives $L(t\sin t) = 2\cos t$. This is exactly unit 06's "no damping, driven at $\omega_0$: the ansatz fails" slide; there the amplitude grew linearly because energy was pumped in at resonance. Here nothing pumps energy in — the system is conservative — so the growth cannot be physical. What the secular term really is, two slides on, shows what it is instead.

Step 5 in full. $\dot x_1 = -\tfrac38A^3(\sin t + t\cos t) + \tfrac1{32}A^3(-3\sin 3t + \sin t) + (\text{homogeneous})'$; every term vanishes at $t = 0$, so the homogeneous part needs no $\sin t$.

Why "secular". Latin saeculum, an age or a century. Eighteenth-century astronomers — Lagrange and Laplace above all — found terms growing like $t$ in their expansions of planetary orbits and called them secular inequalities, because they showed up over centuries rather than within a single orbit. The word has stuck to every term in a perturbation series that grows without bound in the independent variable. Section 7 returns the word to the astronomers.

The ✓ in step 6. Always substitute back. It costs a line, and the resonant coefficient is where factors of 2 go missing.

21 · secular terms

Why the naive answer dies at t ~ 1/ε

$\varepsilon = 0.1$, $A = 1$. The exact solution (numerical) never leaves the shaded band $|x| \le 1$; the naive $x_0 + \varepsilon x_1$ drifts out of phase and then out of the band, starting near $t = 1/\varepsilon$.

the true motion is bounded

For $\varepsilon \ge 0$, $|x(t)| \le A$ for all $t$.

Why: energy $\tfrac12\dot x^2 + \tfrac12x^2 + \tfrac14\varepsilon x^4$ is conserved, and $\tfrac12x^2 + \tfrac14\varepsilon x^4$ is increasing in $|x|$. $\square$

$t$	exact	$x_0 + \varepsilon x_1$
1	0.5047	0.5040
10	$-0.5834$	$-0.6320$
50	0.0004	1.4561
100	$-1.0000$	2.7584

The numbers. Exact values from a numerical integration with relative tolerance $10^{-12}$; energy conservation was checked to hold to that accuracy, and $\max|x| = 1.0000$ over $[0, 100]$, confirming the bound. At $t = 1$ the naive answer is good to $7\times10^{-4}$, of order $\varepsilon^2 = 10^{-2}$ times a modest constant. At $t = 10 = 1/\varepsilon$ it is off by $0.05$. By $t = 100$ it is $2.76$, outside the band $|x| \le 1$ the exact motion can never leave.

Full proof of the bound. Multiply $\ddot x + x + \varepsilon x^3 = 0$ by $\dot x$: $\frac{\mathrm{d}}{\mathrm{d}t}E = 0$ with $E = \tfrac12\dot x^2 + V(x)$, $V(x) = \tfrac12x^2 + \tfrac14\varepsilon x^4$ (unit 00's first integral). So $V(x(t)) \le E = V(A)$. For $\varepsilon \ge 0$, $V$ is even and strictly increasing on $[0, \infty)$, hence $|x(t)| \le A$. In phase space (unit 02) the orbit is a closed level curve of $E$ around a centre: the motion is periodic. For $\varepsilon \lt 0$ the same holds provided $A^2 \lt 1/|\varepsilon|$, so that the orbit stays inside the well.

Why the table's exact value at $t = 100$ is $-1.0000$. Coincidence, not a bug: the exact frequency is $1.036717$, and $1.036717\times100 = 103.67 \approx 33\pi = 103.67$, so $x \approx A\cos 33\pi = -1$.

21 · secular terms

What the secular term really is

$x = A\cos(\text{phase})$ is the shadow of a point going round a circle. A frequency error $\delta$ puts the true point an arc $A\delta t$ further on; Taylor-expanding moves it the same length along the tangent instead, off the circle — and its shadow leaves $|x| \le A$. Schematic, $\delta t = 1.3$.

a phase error, Taylor-expanded

Expand a slightly wrong frequency: $A\cos\big((1 + \varepsilon\omega_1)t\big) = A\cos t - \varepsilon\omega_1A\,t\sin t + O(\varepsilon^2t^2)$. With $\omega_1 = \tfrac38A^2$ that is exactly $x_0 + $ the secular part of $\varepsilon x_1$. The series is Taylor-expanding a phase error, which is fine while $\varepsilon t \ll 1$ and nonsense after.

uniformly valid

$x_0 + \varepsilon x_1$ has error $O(\varepsilon^2)$ at each fixed $t$, but not uniformly in $t \in [0, \infty)$: no single constant works for all $t$.

The figure, as algebra. Write the motion as the point $(A\cos\phi, A\sin\phi)$ with phase $\phi = (1 + \delta)t$, $\delta = \varepsilon\omega_1$; its first coordinate is $x$. (It is the phase-space point $(x, -\dot x)$ of unit 02, at leading order.) Expanding in $\delta$ at fixed $t$, $A\cos(t + \delta t) = A\cos t - \delta t\,A\sin t + O(\delta^2t^2)$, which is the first coordinate of $P + A\delta t\,(-\sin t, \cos t)$: the point $P$ at angle $t$, moved a distance $A\delta t$ along the unit tangent $(-\sin t, \cos t)$. The arc and the tangent agree to first order in $\delta t$ and part company after. In the drawing, $t = 90^\circ$ and $\delta t = 1.3$: the exact shadow is $-0.96A$, the naive one $-1.30A$, outside the band. The distance from the centre along the tangent is $A\sqrt{1 + (\delta t)^2}$, which is why the naive solution's envelope grows like $\sqrt{1 + (\tfrac38\varepsilon A^2t)^2}$.

The general point. An asymptotic expansion in $\varepsilon$ is a statement about $\varepsilon \to 0$ at fixed everything else. If "everything else" includes a variable ranging over an unbounded set, the expansion may be valid pointwise and fail uniformly. The failure occurs where $\varepsilon t$ is no longer small: the problem has a second time scale, $1/\varepsilon$, which the naive expansion cannot see.

21 · secular terms

How long the naive answer is good: a Grönwall bound

theorem

Let $\varepsilon \ge 0$ and $r = x - x_0 - \varepsilon x_1$. Then for $0 \le t \le T$,

$$|r(t)| \le \tfrac{3}{16}\,\varepsilon^2A^5\,\big(3T^2 + T\big)\,e^{3\varepsilon A^2T}.$$

proof, in outline

$r$ solves $\ddot r + r = -\varepsilon(x^3 - x_0^3)$ with zero data. Since $|x|, |x_0| \le A$, $|x^3 - x_0^3| \le 3A^2|x - x_0| \le 3A^2(\varepsilon|x_1| + |r|)$. Variation of parameters (unit 06) writes $r$ as an integral against $\sin(t - s)$, whose size is at most 1; then unit 01's Grönwall inequality closes the estimate.

full proof in the notes (N)

What it says. On any fixed interval the error is $O(\varepsilon^2)$: the naive expansion is asymptotic there. The bound is useful while $\varepsilon T \lesssim 1$ and useless for $T \gg 1/\varepsilon$, exactly where the table failed.

Check. $\varepsilon = 0.1$, $A = 1$, $T = 1$: bound $0.0101$; the true maximum error on $[0, 1]$ is $7.0\times10^{-4}$. Pessimistic, as Grönwall bounds are, but of the right order in $\varepsilon$.

Full proof. Write $u(t) = \sqrt{r^2 + \dot r^2}$. From $\ddot x + x = -\varepsilon x^3$, $\ddot x_0 + x_0 = 0$ and $\ddot x_1 + x_1 = -x_0^3$ we get $\ddot r + r = -\varepsilon x^3 + \varepsilon x_0^3 =: g(t)$, with $r(0) = \dot r(0) = 0$. Variation of parameters (unit 06) gives the unique solution $r(t) = \int_0^t \sin(t - s)\,g(s)\,\mathrm{d}s$ and $\dot r(t) = \int_0^t\cos(t - s)\,g(s)\,\mathrm{d}s$. So $(r, \dot r)$ is the integral of the vector $(\sin(t - s), \cos(t - s))\,g(s)$, whose length is $|g(s)|$, and the triangle inequality for integrals of vectors gives $u(t) \le \int_0^t |g(s)|\,\mathrm{d}s$.

Now bound $g$. $x^3 - x_0^3 = (x - x_0)(x^2 + xx_0 + x_0^2)$, and $|x| \le A$ (the previous slide) and $|x_0| \le A$, so $|g| \le 3\varepsilon A^2|x - x_0| \le 3\varepsilon A^2(\varepsilon|x_1| + u)$. From the worked example, $|x_1(s)| \le \tfrac38A^3s + \tfrac1{16}A^3$. So for $t \le T$,

$$u(t) \le \underbrace{3\varepsilon^2A^2\!\int_0^T\!\big(\tfrac38A^3s + \tfrac1{16}A^3\big)\mathrm{d}s}_{=\ \frac3{16}\varepsilon^2A^5(3T^2 + T)} + 3\varepsilon A^2\!\int_0^t u(s)\,\mathrm{d}s .$$

This is unit 01's Grönwall hypothesis on $[0, T]$ with constant $A_{\rm G} = \tfrac3{16}\varepsilon^2A^5(3T^2 + T)$ and $L = 3\varepsilon A^2$ (the first integral was extended from $[0,t]$ to $[0,T]$ to make it a constant). Hence $u(t) \le A_{\rm G}e^{3\varepsilon A^2 t}$ on $[0,T]$; put $t = T$, and $|r| \le u$. $\square$

Dimensions. In scaled variables $A$ is dimensionless; restoring units, $\varepsilon A^2$ is the dimensionless group $\beta a^2/k$ and the bound carries one power of length, as $r$ does.

The check's number. $\tfrac3{16}(0.01)(3 + 1)e^{0.3} = 0.0075\times 1.3499 = 0.0101$.

section 4

Removing secular terms

If the secular term is a frequency error, put the frequency in the ansatz.

21 · removing secular terms

The solvability condition

the solvability condition

Let $f$ be continuous and $2\pi$-periodic. The equation $u'' + u = f$ has a $2\pi$-periodic solution if and only if

$$\int_0^{2\pi} f(\tau)\cos\tau\,\mathrm{d}\tau = 0 \quad\text{and}\quad \int_0^{2\pi} f(\tau)\sin\tau\,\mathrm{d}\tau = 0 .$$

only if, in one line

Multiply by $\cos\tau$ and integrate over a period; integrating by parts twice moves both derivatives onto $\cos\tau$ with no boundary terms (periodicity), and $\cos'' + \cos = 0$. Same for $\sin$. $\square$

If: the particular solution $u = \int_0^\tau \sin(\tau - s)f(s)\,\mathrm{d}s$ from variation of parameters satisfies $u(\tau + 2\pi) - u(\tau) = \sin\tau\int_0^{2\pi}\!f\cos - \cos\tau\int_0^{2\pi}\!f\sin = 0$. Details in the notes.

In words: the forcing must have no component along the kernel of $\mathrm{d}^2/\mathrm{d}\tau^2 + 1$. A component along $\cos\tau$ or $\sin\tau$ is resonant, and produces $\tau\sin\tau$ or $\tau\cos\tau$.

Full proof. Only if. Let $u$ be a $2\pi$-periodic $C^2$ solution. Then $\int_0^{2\pi}u''\cos = [u'\cos]_0^{2\pi} + \int_0^{2\pi}u'\sin = 0 + [u\sin]_0^{2\pi} - \int_0^{2\pi}u\cos$, the bracketed terms vanishing by periodicity. So $\int_0^{2\pi}(u'' + u)\cos = 0$, which is $\int f\cos = 0$. The same computation with $\sin$.

If. Variation of parameters (unit 06) gives the solution with zero data at $0$: $u(\tau) = \int_0^\tau\sin(\tau - s)f(s)\,\mathrm{d}s$. Then $u(\tau + 2\pi) - u(\tau) = \int_\tau^{\tau+2\pi}\sin(\tau - s)f(s)\,\mathrm{d}s$ (the integrand is $2\pi$-periodic in $s$ once $\tau$ is fixed, because $f$ and $\sin$ are), and expanding $\sin(\tau - s) = \sin\tau\cos s - \cos\tau\sin s$ turns this into $\sin\tau\int f\cos - \cos\tau\int f\sin$ over one period, which is $0$. The same argument for $u'(\tau) = \int_0^\tau\cos(\tau - s)f(s)\,\mathrm{d}s$ shows $u'$ is periodic too. A function with $u(\tau + 2\pi) = u(\tau)$ for all $\tau$ is periodic. $\square$

The linear-algebra reading. $\mathcal L = \mathrm{d}^2/\mathrm{d}\tau^2 + 1$ on $2\pi$-periodic functions is symmetric for the inner product $\int_0^{2\pi}fg$ (the integration by parts above), and its kernel is spanned by $\cos, \sin$. The condition says $f \perp \ker\mathcal L$ — the same statement as "$Ax = b$ is solvable iff $b \perp \ker A^{\top}$" for a matrix. It is called the Fredholm alternative in that generality; here it was proved by hand.

21 · removing secular terms · example

Worked: Poincaré–Lindstedt for Duffing

worked example 3 — the Poincaré–Lindstedt method on the running example

Step 1 — strain the time. The motion is periodic with an unknown frequency $\omega(\varepsilon)$. Put $\tau = \omega t$, so $x$ is $2\pi$-periodic in $\tau$, and $\omega^2x'' + x + \varepsilon x^3 = 0$.

Step 2 — expand both. $x = x_0(\tau) + \varepsilon x_1(\tau) + \cdots$, $\ \omega = 1 + \varepsilon\omega_1 + \cdots$, so $\omega^2 = 1 + 2\varepsilon\omega_1 + O(\varepsilon^2)$.

Step 3 — order $\varepsilon^0$. $x_0'' + x_0 = 0$, $x_0 = A\cos\tau$.

Step 4 — order $\varepsilon^1$.

$$\begin{aligned}x_1'' + x_1 &= -2\omega_1x_0'' - x_0^3\\ &= \big(2\omega_1A - \tfrac34A^3\big)\cos\tau - \tfrac14A^3\cos3\tau .\end{aligned}$$

Step 5 — demand periodicity. The solvability condition: the $\cos\tau$ coefficient must vanish (the $\sin\tau$ one already does).

$$\boxed{\ \omega_1 = \tfrac38A^2,\qquad \omega = 1 + \tfrac38\varepsilon A^2 + O(\varepsilon^2)\ }$$

Step 6 — the correction. $x_1 = \tfrac1{32}A^3(\cos3\tau - \cos\tau)$, fixing $x_1(0) = 0$. No secular term.

Step 7 — undo the straining. $x \approx A\cos\omega t + \tfrac{\varepsilon A^3}{32}(\cos3\omega t - \cos\omega t)$ — bounded, periodic, and uniformly valid on times of order $1/\varepsilon$.

Step 4, the right-hand side. Substituting into $\omega^2x'' + x + \varepsilon x^3 = 0$ and collecting order $\varepsilon$: $x_1'' + 2\omega_1x_0'' + x_1 + x_0^3 = 0$. With $x_0'' = -A\cos\tau$ and the cube identity, $-2\omega_1x_0'' - x_0^3 = 2\omega_1A\cos\tau - \tfrac34A^3\cos\tau - \tfrac14A^3\cos 3\tau$.

Step 6. Particular solution $\tfrac1{32}A^3\cos3\tau$ (as in worked example 2), plus $d\cos\tau + e\sin\tau$. $x_1'(0) = 0$ gives $e = 0$, $x_1(0) = 0$ gives $d = -\tfrac1{32}A^3$.

Why this works and the naive method did not. The naive method insisted the frequency was 1 at every order and then had to express a frequency shift as a polynomial in $t$ times trigonometric functions. Poincaré–Lindstedt gives the frequency its own expansion, and uses the one freedom this introduces, $\omega_1$, to satisfy the one condition — solvability — that periodicity imposes. At every order there is exactly one new unknown $\omega_n$ and one condition, which is why the method can be continued indefinitely; the problem sheet does order two and gets $\omega_2 = -\tfrac{21}{256}A^4$.

When it applies. Only when the true solution is periodic — here guaranteed by the closed energy level curves of the previous section. For a damped or driven oscillator whose amplitude drifts there is no single frequency to strain, and the method has nothing to hold on to. That is what multiple scales is for.

History. Lindstedt (1882) introduced the device for the three-body problem; Poincaré (Les méthodes nouvelles de la mécanique céleste, 1892–99) proved the resulting series are asymptotic — and that for the three-body problem they generally diverge.

21 · removing secular terms

Checking the frequency

dimensions, restored

$\omega = \sqrt{k/m}\,\big(1 + \tfrac{3\beta a^2}{8k}\big)$. $[\beta a^2/k] = 1$ ✓, and the prefactor is a frequency ✓.

limit: $\varepsilon \to 0$

$\omega \to 1$, the linear oscillator ✓. $\varepsilon \gt 0$ (a stiffening spring) raises the frequency, $\varepsilon \lt 0$ lowers it — as it must, since a stiffer spring is faster ✓.

special case: the pendulum

$\varepsilon = -A^2/6$ at unit amplitude: $\omega = 1 - \tfrac1{16}A^2$, so $T = 2\pi(1 + \tfrac1{16}A^2 + \cdots)$ — unit 00's $\theta_0^2/16$, obtained there from an elliptic integral ✓.

against the quadrature, on Duffing

$\varepsilon A^2$	exact $\omega$	$1 + \tfrac38\varepsilon A^2$	error
0.01	1.003742	1.003750	$8.2\times10^{-6}$
0.1	1.036717	1.037500	$7.8\times10^{-4}$
0.5	1.170781	1.187500	$1.7\times10^{-2}$
1	1.317776	1.375000	$5.7\times10^{-2}$

The error falls a hundredfold when $\varepsilon A^2$ falls tenfold: $O(\varepsilon^2)$, as promised. At $\varepsilon A^2 = 1$ it is still only $4\%$ wrong.

And the time series: at $\varepsilon = 0.1$, $t = 100$, the Poincaré–Lindstedt answer is $-0.9969$ against the exact $-1.0000$, where the naive one said $2.758$.

The pendulum, in numbers. For $\theta_0 = 23^\circ = 0.4014$ rad, the exact pendulum frequency is $\pi/2K(\sin(\theta_0/2))= 0.989937$ (unit 00's elliptic integral, $K$ with modulus $\sin(\theta_0/2)$), and $1 - \theta_0^2/16 = 0.989929$. Unit 00's "one per cent wrong at about $23^\circ$" is the statement $1/0.98994 = 1.0102$. The Duffing truncation of the pendulum (dropping $\theta^5$) is a different equation from the pendulum, and agrees with it only to this order; at order $\varepsilon^2$ the $\theta^5$ term enters too.

The remaining error at $t = 100$. Poincaré–Lindstedt to first order has a frequency error of $7.8\times10^{-4}$, so by $t = 100$ the phase is off by $0.078$ radians, and the error in $x$ is about $0.078$ times $|\sin(\text{phase})|$. At $t = 100$ that sine happens to be near zero, which is why the error there is only $0.003$; at $t = 50$, where the phase error is $0.039$ and the sine is near one, the error is $0.039$. Poincaré–Lindstedt to order $n$ is uniformly valid on times of order $\varepsilon^{-n}$, not for ever. Nothing is uniformly valid for ever unless it is exact.

How the "exact" column was obtained is in the notes to "Poincaré's definition": a quadrature from energy conservation, evaluated numerically.

21 · removing secular terms

Multiple scales

the method of multiple scales

Pretend the solution depends on two times, a fast $t$ and a slow $T = \varepsilon t$, independently: $x(t) = X(t, \varepsilon t)$. By the chain rule (unit 03)

$$\frac{\mathrm{d}}{\mathrm{d}t} = \partial_t + \varepsilon\,\partial_T,$$ $$\frac{\mathrm{d}^2}{\mathrm{d}t^2} = \partial_t^2 + 2\varepsilon\,\partial_t\partial_T + \varepsilon^2\partial_T^2 .$$

Expand $X = X_0 + \varepsilon X_1 + \cdots$ and use the new freedom — how $X_0$ depends on $T$ — to kill secular terms.

On Duffing. Order 1: $X_0 = a(T)\cos\big(t + \phi(T)\big)$ — amplitude and phase are free to drift. Order $\varepsilon$:

$$\begin{aligned}\partial_t^2X_1 + X_1 = {}&2a'\sin\psi + \big(2a\phi' - \tfrac34a^3\big)\cos\psi\\ &- \tfrac14a^3\cos3\psi,\end{aligned}$$

with $\psi = t + \phi$. No secular term iff the amplitude equations hold:

$$a' = 0,\qquad \phi' = \tfrac38a^2 .$$

So $x \approx A\cos\big((1 + \tfrac38\varepsilon A^2)t\big)$ — the same answer, with no assumption of periodicity.

Order $\varepsilon$, in full. Substituting, the $\varepsilon^1$ terms are $\partial_t^2X_1 + X_1 + 2\partial_t\partial_TX_0 + X_0^3 = 0$. With $X_0 = a\cos\psi$, $\partial_tX_0 = -a\sin\psi$ and $\partial_T\partial_tX_0 = -a'\sin\psi - a\phi'\cos\psi$. So $-2\partial_t\partial_TX_0 = 2a'\sin\psi + 2a\phi'\cos\psi$, and $-X_0^3 = -\tfrac34a^3\cos\psi - \tfrac14a^3\cos3\psi$. For $X_1$ to stay bounded in $t$ at fixed $T$, the solvability condition kills the $\sin\psi$ and $\cos\psi$ coefficients: $a' = 0$ and $\phi' = \tfrac38a^2$. With $a(0) = A$, $\phi(0) = 0$: $\phi = \tfrac38A^2T = \tfrac38\varepsilon A^2t$.

Is "pretend there are two independent times" legitimate? As a derivation, no: $x$ is a function of one variable. As a way of generating an ansatz, yes — and the ansatz is then justified the way the naive one was, by a Grönwall estimate on the remainder, which now holds on intervals of length $O(1/\varepsilon)$ (cited: Kevorkian and Cole, Multiple Scale and Singular Perturbation Methods, the chapter on multiple-scale expansions for ordinary differential equations; or Verhulst, Nonlinear Differential Equations and Dynamical Systems, the chapter on averaging, where the equivalent averaging method is proved in full).

21 · removing secular terms

Two clocks, plotted

$\ddot x + 2\varepsilon\mu\dot x + x + \varepsilon x^3 = 0$, $\varepsilon = 0.1$, $\mu = 0.5$, $A = 1$: numerical (solid), multiple scales $a\cos(t + \phi)$ (dashed), and its envelope $\pm a(T)$ (dotted). The fast clock sets the wiggles, the slow one the envelope.

Add weak damping to Duffing. Nothing is periodic any more, so Poincaré–Lindstedt cannot start. Multiple scales gains one term, and the amplitude equations become

$$a' = -\mu a,\qquad \phi' = \tfrac38a^2,$$

so $a = Ae^{-\mu T}$: the amplitude decays on the slow clock $T = \varepsilon t$, and the frequency $1 + \tfrac38\varepsilon a^2$ drifts back down towards 1 as it does.

Largest error over $0 \le t \le 60$: $0.044$, of order $\varepsilon$ — what a first-order expansion promises.

What it buys over Poincaré–Lindstedt. Add weak damping, $\ddot x + 2\varepsilon\mu\dot x + x + \varepsilon x^3 = 0$. There is no periodic solution, so Poincaré–Lindstedt cannot start. Multiple scales gets one extra term, $-2\mu\partial_tX_0 = 2\mu a\sin\psi$, and the amplitude equations become $a' = -\mu a$, $\phi' = \tfrac38a^2$: the amplitude decays as $Ae^{-\mu\varepsilon t}$ and the frequency drifts down towards 1 as it does, $\omega(t) = 1 + \tfrac38\varepsilon A^2e^{-2\mu\varepsilon t}$. The problem sheet does the same for the van der Pol oscillator of unit 02 and finds its limit cycle's amplitude, which Poincaré–Bendixson could only prove exists.

How the figure was computed. The damped equation integrated numerically (relative tolerance $10^{-11}$); the approximation is $x \approx a(T)\cos(t + \phi(T))$ with $a = Ae^{-\mu T}$ and $\phi = \tfrac38A^2\,(1 - e^{-2\mu T})/2\mu$, the solution of the amplitude equations with $a(0) = A$, $\phi(0) = 0$. Over $t \le 60$ the two differ by at most $0.044$.

21 · removing secular terms

Secular terms are a boundary layer in time

Look at what the three methods did to Duffing:

naive: one time variable $t$; valid for $t = O(1)$;
Poincaré–Lindstedt: one strained variable $\omega t$; valid for $t = O(1/\varepsilon)$ provided the motion is periodic;
multiple scales: two variables $t$ and $\varepsilon t$; valid for $t = O(1/\varepsilon)$ in general.

The naive expansion is good on the "inner" time scale and fails on the "outer" one, $T = \varepsilon t = O(1)$. Multiple scales writes down an expansion that is right on both at once.

the unifying idea

Every non-uniformity has two scales, and a correct approximation respects both. Secular growth is two time scales. The next section's boundary layer is two length scales. WKB is a fast phase and a slow amplitude.

choosing a method

Periodic solution, want its frequency: Poincaré–Lindstedt, fastest. Drifting amplitude, damping, transients, slow variation of a parameter: multiple scales. Thin region of rapid change in space: matched expansions.

The analogy made exact. In section 5 we shall have an "inner" variable $X = x/\varepsilon$ near a boundary and an "outer" variable $x$. Here, the roles are: the fast time $t$ plays the inner variable, the slow time $T = \varepsilon t$ plays the outer, and the "boundary" is $t = 0$. The naive expansion is an inner expansion. The amplitude equations for $a(T)$, $\phi(T)$ are the outer equations, and multiple scales solves them simultaneously rather than matching them afterwards. Some texts do treat secular problems by matched expansions in time; the answers agree.

An honest limit. Two scales are enough for first order. At second order new secular terms can appear on the time scale $1/\varepsilon^2$, and a third time $\varepsilon^2t$ is needed. The pattern continues, and for Hamiltonian systems with more than one degree of freedom it can break down entirely — resonances between frequencies produce small divisors, which is the door into unit 17's chaos, not opened here.

section 5

Boundary layers

A small parameter on the highest derivative, and the thin region where it is not small.

21 · boundary layers

A small parameter on the highest derivative

the model problem $$\varepsilon y'' + (1 + x)\,y' + y = 0$$

on $[0, 1]$, with $y(0) = 0$ and $y(1) = 1$.

Set $\varepsilon = 0$: $(1 + x)y' + y = 0$, i.e. $\big((1 + x)y\big)' = 0$, so $y = C/(1 + x)$. A first-order equation takes one boundary condition. $y(1) = 1$ gives $C = 2$; then $y(0) = 2 \ne 0$.

The reduced problem cannot meet both conditions. Near one end, $y$ must change by $O(1)$ over a short distance, where $\varepsilon y''$ is not small: a boundary layer.

this one is exactly solvable

$\varepsilon y'' + (1 + x)y' + y = \big(\varepsilon y' + (1 + x)y\big)'$, so $\varepsilon y' + (1+x)y = K$, constant. A first-order linear equation:

$$y(x) = \frac K\varepsilon\int_0^x e^{[(1+s)^2 - (1+x)^2]/2\varepsilon}\,\mathrm{d}s ,$$

with $K$ fixed by $y(1) = 1$. We shall not need it to find the answer — only to check it.

Physics: Prandtl, 1904. Viscosity multiplies the highest derivative in the equation of fluid motion, and the flow past a wing is inviscid except in a thin layer at its surface — where the drag is decided.

Checking the exact solution. Integrating factor for $\varepsilon y' + (1+x)y = K$: multiply by $\varepsilon^{-1}e^{(1+x)^2/2\varepsilon}$, getting $\big(e^{(1+x)^2/2\varepsilon}y\big)' = (K/\varepsilon)e^{(1+x)^2/2\varepsilon}$; integrate from $0$ using $y(0) = 0$. That is the formula. At $\varepsilon = 0.05$, $K = 1.97433$; at $\varepsilon = 0.01$, $K = 1.99497$. Note that $K = (1+x)y + \varepsilon y'$ is exactly the combination the outer solution makes constant, and the outer value is $K = 2$: the exact constant tends to the outer one as $\varepsilon \to 0$. That is already a first matching statement.

Why the layer is not at $x = 1$. Deferred to "Where the layer goes"; the short version is that a layer at $x = 1$ would have to grow exponentially into the interior.

Prandtl. The Navier–Stokes equations of viscous flow have viscosity multiplying the second derivatives of the velocity, the highest derivatives in them. Setting viscosity to zero gives Euler's equations, which cannot satisfy no-slip at a wall and predict zero drag on any body (d'Alembert's paradox, 1752). Prandtl's 1904 paper, eight pages long, resolved the paradox with exactly this section's idea. Van Dyke later made the method systematic, and the matching rule two slides on is his.

21 · boundary layers

Inner and outer solutions

outer solution

The expansion valid away from the layer, obeying the boundary condition the reduced equation can keep: $y_{\rm out} = 2/(1 + x)$.

inner solution

Stretch the variable near $x = 0$: $x = \varepsilon^\alpha X$. The terms $\varepsilon y''$, $(1+x)y'$, $y$ become of sizes $\varepsilon^{1 - 2\alpha}$, $\varepsilon^{-\alpha}$, $1$. The balance that keeps $y''$ is $\alpha = 1$ — the distinguished limit. With $Y(X) = y(\varepsilon X)$:

$$Y'' + (1 + \varepsilon X)Y' + \varepsilon Y = 0,$$

so at leading order $Y_0'' + Y_0' = 0$.

With $Y_0(0) = 0$: $Y_0 = C\,(1 - e^{-X})$. One constant is still unknown.

$\varepsilon = 0.05$: the exact solution (solid) and the composite (short dashes) climb steeply through the shaded layer, then join the outer solution $2/(1+x)$ (long dashes) — right away from $x = 0$, wrong in the layer.

The inner equation, in full. With $x = \varepsilon X$ and $Y(X) = y(\varepsilon X)$, the chain rule gives $y' = \varepsilon^{-1}Y'$ and $y'' = \varepsilon^{-2}Y''$. Substituting: $\varepsilon\cdot\varepsilon^{-2}Y'' + (1 + \varepsilon X)\varepsilon^{-1}Y' + Y = 0$; multiply by $\varepsilon$. At leading order the $\varepsilon XY'$ and $\varepsilon Y$ terms drop, leaving $Y_0'' + Y_0' = 0$ with general solution $C_1 + C_2e^{-X}$, and $Y_0(0) = 0$ forces $C_1 = -C_2$.

Dominant balance, once more. The other balances: $\alpha = 0$ ($y'$ against $y$) is the outer problem; $\varepsilon^{1-2\alpha} = 1$ gives $\alpha = \tfrac12$, but then $y'$, of size $\varepsilon^{-1/2}$, is larger than both — inconsistent. So there are exactly two regions, and the layer thickness is $\varepsilon$, not $\sqrt\varepsilon$. When the coefficient of $y'$ vanishes at the boundary the balance changes and so does the thickness; the problem sheet's problems all have a non-vanishing coefficient.

The figure is computed: exact from the quadrature on the previous slide; composite from the next slides. The shaded strip is $0 \le x \le \varepsilon$. The exact solution never reaches $2$: it peaks at $1.71$ near $x = 0.16$, because by the time the layer has finished climbing, the outer solution it is joining has already fallen from $2$.

21 · boundary layers

Matching: the intermediate variable and Van Dyke's rule

The inner solution must turn into the outer one somewhere in between. Matching makes that precise, in either of two ways.

through an intermediate variable

Put $x = \varepsilon^\beta\xi$ with $0 \lt \beta \lt 1$ and $\xi$ fixed, so $\varepsilon \ll x \ll 1$. Then

$$y_{\rm out} = \frac2{1 + \varepsilon^\beta\xi} \to 2,$$ $$Y_0 = C\big(1 - e^{-\varepsilon^{\beta-1}\xi}\big) \to C .$$

Both are approximations to the same $y$ in the overlap, so $C = 2$.

Van Dyke's matching rule

The $m$-term inner expansion of the $n$-term outer solution equals the $n$-term outer expansion of the $m$-term inner solution.

With $m = n = 1$: the outer $2/(1 + \varepsilon X)$ in inner variables, to one term, is $2$. The inner $C(1 - e^{-x/\varepsilon})$ in outer variables, to one term, is $C$ (the exponential is beyond all orders). So $C = 2$ again.

a rule, not a theorem

Van Dyke's rule is a procedure that works for power-law expansions and can fail when logarithms appear. The intermediate-variable argument is the one to trust when they disagree.

Why the overlap exists. The outer expansion is valid for $x \gg \varepsilon$ (away from the layer) and the inner for $x \ll 1$ (inside the region where $1 + x \approx 1$). An intermediate variable sits in both. The matching principle — that the two expansions agree there to the order both are valid — is the hypothesis of the method. It can be proved for a given problem by comparing with the exact solution or by a maximum-principle estimate, but in general it is an assumption justified afterwards. For the model problem the exact solution is available, and the next slide checks it.

Van Dyke's rule with logarithms. The standard failures involve terms like $\varepsilon\ln\varepsilon$, which do not sort cleanly into "orders"; the fix is to count $\ln\varepsilon$ with the power it multiplies. Cited: Hinch, Perturbation Methods (CUP, 1991), the chapter on matched asymptotic expansions, which gives both the rule and its failures with examples; Van Dyke, Perturbation Methods in Fluid Mechanics (1964/1975), for the original.

Higher orders. With $m = n = 2$ the rule matches the $O(\varepsilon)$ corrections as well and yields a composite with $O(\varepsilon^2)$ error. The procedure is mechanical; the bookkeeping is the whole difficulty. One order is enough to see the method.

21 · boundary layers

The overlap, on a logarithmic axis

$\varepsilon = 0.001$, $x$ on a log scale so that both scales fit. In the middle band, $0.003 \lt x \lt 0.05$, inner and outer are both within $0.1$ of the exact solution, and both are close to $2$.

On a linear axis the layer is a sliver at the wall. On a log axis every decade gets the same room, and three regions appear:

inner, $x = O(\varepsilon)$: $y$ climbs; only $Y_0$ is right;
overlap, $\varepsilon \ll x \ll 1$: both are right, and both are $\approx 2$ — this is what fixes $C = 2$;
outer, $x = O(1)$: $y$ falls to $1$; only $y_{\rm out}$ is right.

The intermediate variable $x = \varepsilon^\beta\xi$ is a point that stays in the middle band as $\varepsilon \to 0$.

How the figure was computed. Exact solution by the quadrature on A small parameter on the highest derivative, at $\varepsilon = 0.001$; $Y_0 = 2(1 - e^{-x/\varepsilon})$ and $y_{\rm out} = 2/(1+x)$ as derived. The bands are computed, not drawn by eye: the middle one is where both approximations are within $0.1$ of the exact $y$, which is $0.0030 \lt x \lt 0.053$. Sample values: at $x = 0.005$ the errors are $0.008$ (inner) and $0.012$ (outer); at $x = 0.1$ the inner is off by $0.18$ and the outer by $0.001$; at $x = 0.002$ the outer is off by $0.27$ and the inner by $0.002$.

Why the overlap widens as $\varepsilon \to 0$. The inner solution's error grows like $x$ (it ignores the slope of $2/(1+x)$), and the outer one's like $e^{-x/\varepsilon}$, so for a fixed tolerance the overlap runs from a few $\varepsilon$ to a fixed $x$: on a log axis it spans about $\log_{10}(1/\varepsilon)$ decades, less a margin. At $\varepsilon = 0.05$, the value of the previous figures, there is barely room for it, which is why this slide uses $\varepsilon = 0.001$. An intermediate variable with $0 \lt \beta \lt 1$ is a point whose logarithm stays a fixed fraction of the way from $\log\varepsilon$ to $0$: it lands in the middle band for every small enough $\varepsilon$.

21 · boundary layers · example

Worked: the composite solution, against the exact one

worked example 4 — the model problem, assembled

Step 1 — outer. $y_{\rm out} = 2/(1 + x)$, from $y(1) = 1$.

Step 2 — inner. $Y_0 = C(1 - e^{-X})$, $X = x/\varepsilon$, from $y(0) = 0$.

Step 3 — match. $C = 2$; the common limit is $y_{\rm match} = 2$.

Step 4 — combine. The composite expansion is outer $+$ inner $-$ common part:

$$\boxed{\ y_c = \frac2{1 + x} - 2e^{-x/\varepsilon}\ }$$

Step 5 — check the ends. $y_c(0) = 0$ exactly; $y_c(1) = 1 - 2e^{-1/\varepsilon}$, off by an exponentially small amount.

	$x = \varepsilon$	$x = 0.1$	$x = 0.5$	max error
$\varepsilon = 0.05$, exact	1.2300	1.6246	1.3477
$\varepsilon = 0.05$, composite	1.1690	1.5475	1.3332	0.077
$\varepsilon = 0.01$, exact	1.2574	1.8289	1.3360
$\varepsilon = 0.01$, composite	1.2444	1.8181	1.3333	0.018

the checks

Order: the maximum error drops by 4.3 when $\varepsilon$ drops by 5: $O(\varepsilon)$, uniformly on $[0,1]$. Limit: as $\varepsilon \to 0$ at fixed $x \gt 0$, $y_c \to 2/(1+x)$, the outer solution ✓. Dimensions: $x$ and $y$ are dimensionless, and $x/\varepsilon$ says the layer's width is $\varepsilon$ in units of the interval ✓.

Why subtract the common part. In the overlap both $y_{\rm out}$ and $Y_0$ are $\approx 2$; adding them counts that value twice. Subtracting the common limit leaves a function that reduces to $y_{\rm out}$ outside the layer (the inner part minus the common part is exponentially small there) and to $Y_0$ inside it (the outer part minus the common part is $2/(1+x) - 2 = O(x) = O(\varepsilon)$ there).

Where the error comes from. The maximum error sits just outside the layer, at $x \approx 0.10$ for $\varepsilon = 0.05$ and $x \approx 0.024$ for $\varepsilon = 0.01$: it is the neglected $O(\varepsilon)$ term of the outer expansion, plus the inner correction. Its size, $0.077 \approx 1.5\varepsilon$ and $0.018 \approx 1.8\varepsilon$, is $O(\varepsilon)$ with a modest constant, which is what a one-term composite promises. The exact values come from the quadrature formula on "A small parameter on the highest derivative".

The problem sheet does the same for $\varepsilon y'' + y' + y = 0$ — whose leading-order composite turns out to be $e^{1-x} - e^{1-x/\varepsilon}$ — and compares at $\varepsilon = 0.01$ against the exact solution, which there is a sum of two exponentials.

21 · boundary layers

Where the layer goes

the sign of the first-derivative term decides

For $\varepsilon y'' + a(x)y' + b(x)y = 0$ with $a$ of one sign on $[0,1]$: at a layer at $x_0$, the inner equation is $Y'' \pm a(x_0)Y' = 0$, with solutions $1$ and $e^{\mp a(x_0)X}$. The layer must decay into the interior.

$a \gt 0$: layer at the left end, $x = 0$.
$a \lt 0$: layer at the right end, $x = 1$.
$a$ changes sign: interior layers or turning points — a subject of its own, named here.

Test it. A layer at $x = 1$ for the model problem: $X = (1 - x)/\varepsilon$ gives $Y'' - 2Y' = 0$, solutions $1$ and $e^{2X}$. The second grows as you leave the wall — it cannot match any outer solution. No layer there.

$\varepsilon = 0.05$, exact solutions (solid) and outer solutions (dashed). The outer solution is carried in from the end whose condition it keeps (arrow); the layer forms where it arrives.

back to Duffing

Duffing's $\varepsilon$ does not touch $\ddot x$, so there is no layer in the equation. Its non-uniformity lives at large $t$: the naive expansion is the "inner" solution, good for $t = O(1)$, and the amplitude equations, on the slow scale $T = \varepsilon t$, are the "outer" ones. Regular in the equation, singular in the domain.

The inner equation at a general endpoint. At $x_0 = 0$ put $X = x/\varepsilon$: leading order $Y'' + a(0)Y' = 0$, solutions $1, e^{-a(0)X}$; the exponential decays into the interior $X \to +\infty$ iff $a(0) \gt 0$. At $x_0 = 1$ put $X = (1 - x)/\varepsilon$, so $\mathrm{d}/\mathrm{d}x = -\varepsilon^{-1}\mathrm{d}/\mathrm{d}X$: leading order $Y'' - a(1)Y' = 0$, solutions $1, e^{a(1)X}$, decaying iff $a(1) \lt 0$. Exactly one end qualifies when $a$ has a fixed sign, and that is where the layer goes.

The figure. Constant coefficients, so both are exact: $y = c_1e^{r_1x} + c_2e^{r_2x}$ with $\varepsilon r^2 \pm r + 1 = 0$. Left, $\varepsilon y'' + y' + y = 0$, $y(0) = 0$, $y(1) = 1$: outer $e^{1-x}$, keeping $y(1) = 1$, and a layer at $x = 0$. Right, $\varepsilon y'' - y' + y = 0$, $y(0) = 1$, $y(1) = 0$: outer $e^{x}$, keeping $y(0) = 1$, and a layer at $x = 1$. The boundary values were chosen so that neither outer solution is identically zero.

Physically. $a\,y'$ is transport — the analogue of flow carrying something downstream — and $\varepsilon y''$ is diffusion. Diffusion is only needed where the flow runs into a wall it cannot pass: the downstream end. For $a \gt 0$ in the equation as written, "downstream" is the left. (Rewrite with $x \to -x$ and it becomes the familiar picture.)

Interior layers and turning points (where $a$ vanishes) are treated in Bender and Orszag, Advanced Mathematical Methods for Scientists and Engineers, ch. 9. They are not needed until far beyond this course.

section 6

WKB

A wave in a medium that changes slowly: fast phase, slow amplitude.

21 · WKB

The WKB approximation

$\varepsilon^2y'' + q(x)\,y = 0$ with $q \gt 0$: an oscillator whose "stiffness" varies slowly compared with the oscillation. $\varepsilon$ is on the highest derivative — singular — but there is no layer: the solution oscillates on the short scale $\varepsilon$ everywhere.

Try $y = e^{S(x)/\varepsilon}$ with $S = S_0 + \varepsilon S_1 + \cdots$ complex. Then $S'^2 + \varepsilon S'' + q = 0$, and

$$\begin{aligned}\varepsilon^0:&\quad S_0'^2 = -q,\\ \varepsilon^1:&\quad 2S_0'S_1' + S_0'' = 0.\end{aligned}$$

the WKB approximation

$S_0 = \pm i\int\sqrt q$, $\ S_1 = -\tfrac14\ln q$, so

$$y \approx q^{-1/4}\big[C_1\cos\Phi(x) + C_2\sin\Phi(x)\big],$$ $$\Phi(x) = \frac1\varepsilon\int^x\!\sqrt{q(s)}\,\mathrm{d}s .$$

Check: $q$ constant gives $\cos(\sqrt q\,x/\varepsilon)$, the exact solution ✓. Dimensions: with $\varepsilon$ a length and $q$ dimensionless, the phase $\int\sqrt q\,\mathrm{d}x/\varepsilon$ is a pure number ✓.

Names: the eikonal equation fixes the phase, the transport equation the amplitude. If you have done unit 15, $S_0'^2 = -q$ is its eikonal again, the Hamilton–Jacobi equation of a particle whose momentum is $\sqrt q$.

The derivation in full. $y = e^{S/\varepsilon}$ gives $y' = (S'/\varepsilon)y$ and $y'' = (S''/\varepsilon + S'^2/\varepsilon^2)y$. So $\varepsilon^2y'' + qy = (\varepsilon S'' + S'^2 + q)y = 0$. With $S = S_0 + \varepsilon S_1 + O(\varepsilon^2)$: $S'^2 = S_0'^2 + 2\varepsilon S_0'S_1' + O(\varepsilon^2)$, and $\varepsilon S'' = \varepsilon S_0'' + O(\varepsilon^2)$. Order 1: $S_0'^2 + q = 0$, so $S_0' = \pm i\sqrt q$. Order $\varepsilon$: $2S_0'S_1' + S_0'' = 0$. Since $S_0' = \pm i\sqrt q$, the ratio $S_0''/S_0' = (\sqrt q)'/\sqrt q$ is real, so $S_1' = -(\sqrt q)'/2\sqrt q$ and $S_1 = -\tfrac12\ln\sqrt q + \text{const} = -\tfrac14\ln q + \text{const}$. Hence $e^{S/\varepsilon} \approx q^{-1/4}\exp(\pm\tfrac{i}{\varepsilon}\int\sqrt q)$, and the two signs combine (unit 06: the complex exponential as bookkeeping) into the real cos/sin form.

Why $q^{-1/4}$, a second way. The equation has no $y'$ term, so by Abel's identity (unit 06) the Wronskian of any two solutions is constant. For the pair $y_\pm = \mathcal A\,e^{\pm i\Phi}$ with slowly varying amplitude $\mathcal A$, the Wronskian is $y_+y_-' - y_-y_+' = -2i\mathcal A^2\Phi' + (\text{terms in }\mathcal A\mathcal A') $, and the $\mathcal A\mathcal A'$ terms cancel exactly, leaving $-2i\mathcal A^2\sqrt q/\varepsilon$. Constancy, to the accuracy with which the WKB pair solves the equation, forces $\mathcal A^2\sqrt q$ constant: amplitude squared times local wavenumber is conserved, which is $\mathcal A \propto q^{-1/4}$ — and, in time, the adiabatic invariant of the next slide.

The names. WKB is Wentzel, Kramers and Brillouin (1926), who used it on the Schrödinger equation; Jeffreys had it in 1924; Liouville and Green had it in 1837, and "Liouville–Green approximation" is the name careful mathematicians use. Unit 15 remark: writing a wave's phase as $S/\varepsilon$ and finding a Hamilton–Jacobi equation for $S$ at leading order is the "mechanics was optics all along" of unit 15, run backwards. Nothing in this unit depends on it.

21 · WKB

When WKB is valid, and an adiabatic invariant

validity condition

The expansion is consistent while the correction is small beside the leading term: $|\varepsilon S_1'| \ll |S_0'|$, that is

$$\frac{\varepsilon\,|q'|}{4\,q^{3/2}} \ll 1 .$$

In words: the wavelength $2\pi\varepsilon/\sqrt q$ changes by a small fraction of itself over one wavelength. It fails where $q \to 0$.

In time: $\ddot x + \omega(\varepsilon t)^2x = 0$ is the same equation ($s = \varepsilon t$, $q = \omega^2$). WKB says $x \approx a\cos\int\omega\,\mathrm{d}t$ with amplitude $a \propto \omega^{-1/2}$, so the energy $E \approx \tfrac12\omega^2a^2$ is proportional to $\omega$:

$$E/\omega \approx \text{const}.$$

Slowly doubling the frequency doubles the energy: the next slide shows it happen.

This adiabatic invariant is the Duffing lesson again: the multiple-scales amplitude equation for a slowly varying $\omega$ is $(\omega a^2)' = 0$ — the same law, from the same solvability condition.

The validity condition, derived. $S_0' = i\sqrt q$ and $S_1' = -q'/4q$. The ratio of the second term to the first in $S'/\varepsilon$ is $\varepsilon|S_1'|/|S_0'| = \varepsilon|q'|/4q^{3/2}$. This is a consistency condition: necessary for the expansion to make sense, not by itself a bound on the error. Rigorous error bounds for the Liouville–Green approximation, in terms of the total variation of a control function, are cited: Olver, Asymptotics and Special Functions, ch. 6.

The multiple-scales route, sketched in full. Let the fast phase $\theta$ satisfy $\dot\theta = \omega(T)$, $T = \varepsilon t$, and $x = X(\theta, T)$. Then $\ddot x = \omega^2X_{\theta\theta} + \varepsilon(2\omega X_{\theta T} + \omega'X_\theta) + O(\varepsilon^2)$, and the equation becomes $\omega^2(X_{\theta\theta} + X) + \varepsilon(2\omega X_{\theta T} + \omega'X_\theta) = 0$. Order 1: $X_0 = a(T)\cos(\theta + \phi(T))$. Order $\varepsilon$: the $\sin$ coefficient of the forcing is $2\omega a' + \omega'a$, and solvability sets it to zero: $(\omega a^2)' = 0$. The energy is $\tfrac12\omega^2a^2$, so $E/\omega = \tfrac12\omega a^2$ is constant. The same as WKB's $a\propto\omega^{-1/2}$.

History. At the 1911 Solvay conference Lorentz asked what happens to a pendulum whose string is slowly shortened; Einstein answered that $E/\nu$ stays constant. The question mattered because Planck's quanta have $E = h\nu$: only an adiabatic invariant could be quantised consistently. For a real pendulum the equation is not quite $\ddot\theta + (g/l)\theta = 0$, because shortening the string also changes the angular momentum; the problem sheet's WKB problem uses a cleaner setting.

21 · WKB

The adiabatic invariant, plotted

$\ddot x + \omega(t)^2x = 0$ with $\omega = 1 + t/100$, integrated numerically. The frequency doubles; the amplitude follows $\omega^{-1/2}$ (dashed) down to $1/\sqrt2$, so $E \approx \tfrac12\omega^2a^2$ doubles while $E/\omega$ stays put.

a numerical check

$\ddot x + (1 + t/100)^2x = 0$, $x(0) = 1$, $\dot x(0) = 0$, run to $t = 100$: the frequency doubles. Integrated numerically:

energy: $\times1.998$ (WKB: $\times2$);
$E/\omega$: changed by $0.12\%$;
amplitude: $0.7067$ (WKB: $1/\sqrt2 = 0.7071$);
$x(100)$: $0.4926$ (WKB: $0.4944$).

What to look for. Twenty-four oscillations, each slightly faster and slightly smaller than the last; the envelope is not fitted, it is $\omega^{-1/2}$ with the starting amplitude $1$. The WKB phase $\int_0^{t}\omega\,\mathrm{d}t$ reaches $150$ radians at $t = 100$ — the "fast" variable — while $\omega$ itself changes by a factor of only two over the whole run — the "slow" one. That separation, a hundred to one here, is what $\varepsilon$ measures.

The numbers. Integrated with relative tolerance $10^{-12}$. Initially $E = \tfrac12$, $\omega = 1$, so $E/\omega = 0.5$; finally $E/\omega = 0.49941$. The small change is dominated by the abrupt start and stop of the ramp at $t = 0$ and $t = 100$, where $\dot\omega$ jumps; the WKB phase at $t = 100$ is $\int_0^{100}(1 + t/100)\,\mathrm{d}t = 150$.

21 · WKB

Turning points and the connection formulae

Where $q$ changes sign, at a turning point $x = a$, the approximation fails: $q^{-1/4}$ blows up and the validity ratio $\varepsilon|q'|/4q^{3/2}$ does too. On one side the solution oscillates; on the other it grows or decays like $|q|^{-1/4}e^{\pm\frac1\varepsilon\int\sqrt{|q|}}$.

Near $a$, $q \approx q'(a)(x - a)$, and the local equation is Airy's, $y'' = zy$ in a stretched variable $z$ — a boundary layer of width $\varepsilon^{2/3}$.

the connection formulae (cited)

Let $q \gt 0$ for $x \lt a$, $q \lt 0$ for $x \gt a$, with $q'(a) \ne 0$. The solution that decays for $x \gt a$,

$$y \approx |q|^{-1/4}\exp\Big(-\frac1\varepsilon\int_a^x\!\sqrt{|q|}\Big),$$

is, for $x \lt a$,

$$y \approx 2\,q^{-1/4}\cos\Big(\frac1\varepsilon\int_x^a\!\sqrt q - \frac\pi4\Big).$$

Proof: match both sides to the Airy function through its asymptotics, which need unit 20's integral methods. Bender and Orszag, ch. 10.

What is stated and what is not. The formula is cited, not proved: the proof needs the large-argument behaviour of the Airy function, which comes from its integral representation by the method of steepest descent (unit 20). What can be checked here: both sides are WKB solutions on their own side (by the previous slides), and the width of the turning-point region follows from dominant balance. With $q \approx -q_1(x - a)$, $q_1 = -q'(a) \gt 0$, put $x - a = \varepsilon^\gamma z$: the terms $\varepsilon^2y''$ and $q y$ have sizes $\varepsilon^{2 - 2\gamma}$ and $\varepsilon^\gamma$, balanced when $\gamma = 2/3$. So the inner region has width $\varepsilon^{2/3}$, and the inner equation is $y_{zz} = q_1 z\,y$, Airy's up to a rescaling of $z$.

The direction warning. The formula is safe in the direction written: from the decaying solution to the oscillatory one. Read backwards it is unstable, because a tiny error in the phase of the cosine introduces a growing exponential on the other side that swamps the decaying one.

A consequence, named. Between two turning points $a \lt b$, a solution decaying on both outer sides exists only if the two cosines agree, which forces $\frac1\varepsilon\int_a^b\sqrt q = (n + \tfrac12)\pi$. In quantum mechanics, with $\varepsilon = \hbar$ and $q = 2m(E - V)$, this is the Bohr–Sommerfeld rule with the $\tfrac12$ that the old quantum theory missed. For $V = \tfrac12x^2$ (units $m = \hbar = 1$) it gives $E = n + \tfrac12$ — exact, a coincidence special to the harmonic oscillator. Named; not used again before quantum mechanics, which this course does not reach.

Reading. Bender and Orszag ch. 10 up to (not including) the Airy-function asymptotics, which use unit 20's methods; then return after unit 20.

21 · WKB

The turning point, pictured

$\varepsilon^2y'' + q\,y = 0$ with $q = a - x$, $\varepsilon = 0.1$: exact (an Airy function) against the two WKB forms of the connection formula, decaying for $x \gt a$ and $2q^{-1/4}\cos(\cdots - \pi/4)$ for $x \lt a$ (same constant).

Away from $a$, each WKB form lies on the exact curve: oscillation where $q \gt 0$, with wavelength $2\pi\varepsilon/\sqrt q$ and amplitude $q^{-1/4}$ both growing towards $a$; decay where $q \lt 0$.

Near $a$, inside the strip of width $\sim\varepsilon^{2/3}$, both blow up like $|x - a|^{-1/4}$ while the exact solution passes smoothly through. That strip is the turning point's boundary layer.

The connection formula is the statement that the two sides belong to one solution — the factor $2$ and the phase $-\pi/4$ are what make them line up.

How the figure was computed. With $q = a - x$, the equation $\varepsilon^2y'' = (x - a)y$ is solved exactly by $y = \mathrm{Ai}\big((x - a)/\varepsilon^{2/3}\big)$, the Airy function that decays to the right, evaluated here by a library routine. The WKB curves are the two sides of the connection formula, multiplied by the one constant $\varepsilon^{1/6}/2\sqrt\pi$ that matches the decaying side to the Airy function's known large-$x$ behaviour (cited, as on the previous slide). With that single constant the oscillatory side lines up too, which is the connection formula's content, seen rather than proved. The shaded strip is $|x - a| \lt \varepsilon^{2/3} = 0.215$.

What to look for. At $x = a - 2$ the WKB and exact curves are indistinguishable at this scale; closer in, the WKB amplitude rises faster than the exact one and goes to infinity at $a$. The decaying side is fitted to the exact curve's constant; the oscillating side is not fitted at all.

section 7

The orbit: a secular term that is the answer

The same calculation, pointed at Mercury.

21 · the orbit

The relativistic orbit equation, nondimensionalised

Unit 08's Binet equation for the inverse-square law, with $u = 1/r$ and $h = r^2\dot\theta = L/\mu$:

$$u'' + u = \frac{GM}{h^2} .$$

General relativity (cited: unit 37 derives it from the Schwarzschild geodesics) adds one term:

$$u'' + u = \frac{GM}{h^2} + \frac{3GM}{c^2}\,u^2 .$$

Nondimensionalise with the Newtonian semi-latus rectum $p = h^2/GM$: put $w = pu$.

$$\boxed{\ w'' + w = 1 + \varepsilon w^2,\qquad \varepsilon = \frac{3GM}{c^2p}\ }$$

the size of the small parameter

$GM/c^2$ is a length: $1.476$ km for the Sun (half its Schwarzschild radius $2.953$ km). For Mercury, $p = a(1 - e^2) = 5.5460\times10^{10}$ m, so

$$\varepsilon = 7.987\times10^{-8}.$$

Classify it. $\varepsilon$ multiplies $w^2$, not $w''$: regular in the equation, like Duffing. And $\theta$ runs over $[0, \infty)$ as the planet goes round: potentially singular in the domain, like Duffing.

At $\varepsilon = 0$: $w_0 = 1 + e\cos\theta$, unit 08's ellipse, perihelion at $\theta = 0$.

The scaling, line by line. Multiply $u'' + u = GM/h^2 + (3GM/c^2)u^2$ by $p = h^2/GM$: $(pu)'' + pu = 1 + (3GM/c^2)\,p\,u^2 = 1 + \frac{3GM}{c^2p}(pu)^2$. So $w = pu$ obeys the boxed equation.

Dimensions. $[GM] = \mathsf{L^3T^{-2}}$ and $[c^2] = \mathsf{L^2T^{-2}}$, so $[GM/c^2] = \mathsf L$ and $\varepsilon$ is a ratio of lengths ✓. Limits. $c \to \infty$ gives $\varepsilon \to 0$ and Newton back ✓; a wider orbit (larger $p$) gives smaller $\varepsilon$ — relativistic effects are strongest close in ✓. Equivalently $\varepsilon = 3(GM/hc)^2$: three times the square of the orbital speed scale $GM/h$ over $c$. For Mercury $GM/hc = 1.632\times10^{-4}$.

The numbers. $GM_\odot = 1.32712440018\times10^{20}\ \mathrm{m^3\,s^{-2}}$, $c = 299\,792\,458$ m/s, Mercury's $a = 5.7909050\times10^{10}$ m and $e = 0.205630$. Then $p = 5.54604\times10^{10}$ m and $\varepsilon = 3\times1476.6/5.54604\times10^{10} = 7.9875\times10^{-8}$.

What is being assumed. The equation itself. This unit's job is not to derive it — that needs the Schwarzschild metric of unit 36 — but to say, correctly, what its solution is. Everything below is unit 21's mathematics applied to a given ODE.

21 · the orbit · example

Worked: the first-order orbit

worked example 5 — expand the relativistic orbit naively

Step 1 — hierarchy. $w = w_0 + \varepsilon w_1 + \cdots$: $\ w_0'' + w_0 = 1$, and $w_1'' + w_1 = w_0^2$.

Step 2 — order zero. $w_0 = 1 + e\cos\theta$.

Step 3 — the forcing. With $\cos^2\theta = \tfrac12(1 + \cos2\theta)$:

$$w_0^2 = \big(1 + \tfrac12e^2\big) + 2e\cos\theta + \tfrac12e^2\cos2\theta .$$

Step 4 — solve term by term. Constant: $1 + \tfrac12e^2$. $\cos2\theta$: $\tfrac{1}{1-4}\cdot\tfrac12e^2 = -\tfrac16e^2$. $\cos\theta$: resonant, and $L(\theta\sin\theta) = 2\cos\theta$ gives $e\,\theta\sin\theta$.

$$w_1 = 1 + \tfrac12e^2 - \tfrac16e^2\cos2\theta + \underbrace{e\,\theta\sin\theta}_{\text{secular}}$$

(plus any homogeneous solution, which only re-labels $e$ and the perihelion direction).

Step 5 — check. $L\big(e\,\theta\sin\theta\big) = 2e\cos\theta$ ✓; $L\big(-\tfrac16e^2\cos2\theta\big) = -\tfrac16e^2(-4 + 1)\cos2\theta = \tfrac12e^2\cos2\theta$ ✓.

Step 6 — the same structure as Duffing. The constant and $\cos2\theta$ pieces change the orbit's shape by a bounded, tiny amount. The secular piece is worked example 2's term again, in $\theta$ instead of $t$.

Step 1. Substituting and collecting: order 1, $w_0'' + w_0 = 1$; order $\varepsilon$, $w_1'' + w_1 = w_0^2$. The general solution of the first is $1 + e\cos(\theta - \theta_0)$; choosing the axis so that perihelion (largest $w$, smallest $r$) is at $\theta = 0$ gives $\theta_0 = 0$.

Step 3. $(1 + e\cos\theta)^2 = 1 + 2e\cos\theta + e^2\cos^2\theta = 1 + 2e\cos\theta + \tfrac12e^2 + \tfrac12e^2\cos2\theta$.

Step 4, the $\cos2\theta$ piece. Try $c\cos2\theta$: $L(c\cos2\theta) = (-4 + 1)c\cos2\theta = -3c\cos2\theta$; set equal to $\tfrac12e^2\cos2\theta$: $c = -\tfrac16e^2$.

Mercury's non-secular corrections. $\varepsilon(1 + \tfrac12e^2) = 8.16\times10^{-8}$: a fractional change in the orbit's size of that order, far below anything measurable against the uncertainty in $p$ itself. The periodic $\cos2\theta$ piece is smaller still. Only the secular piece accumulates — which is why it is the only one anybody measures.

21 · the orbit

Reading the secular term: the perihelion advances

Exactly as for Duffing, $\cos\big((1 - \varepsilon)\theta\big) = \cos\theta + \varepsilon\theta\sin\theta + O(\varepsilon^2\theta^2)$. So, to this order,

$$w \approx 1 + e\cos\big((1 - \varepsilon)\theta\big) + (\text{bounded, }O(\varepsilon)) .$$

Perihelion — largest $w$ — is where $(1 - \varepsilon)\theta = 2\pi k$, at $\theta = 2\pi k/(1 - \varepsilon)$. Each orbit the perihelion's angle $\varpi$ moves on by

$$\boxed{\ \Delta\varpi = 2\pi\varepsilon = \frac{6\pi GM}{c^2\,a(1 - e^2)}\ \text{per orbit}\ }$$

Poincaré–Lindstedt, with $\varphi = \nu\theta$, gives $\nu = 1 - \varepsilon + O(\varepsilon^2)$: the same $\Delta\varpi$, uniformly valid.

Three turns of $w'' + w = 1 + \varepsilon w^2$, integrated numerically, with $e = 0.5$ and $\varepsilon = 0.05$ — exaggerated six hundred thousand times over Mercury's. Successive perihelia (1–4) step round by $\Delta\varpi$ each orbit; the Newtonian ellipse (dashed) closes.

Poincaré–Lindstedt on the orbit, in full. Put $\varphi = \nu\theta$, $\nu = 1 + \varepsilon\nu_1 + \cdots$, so $\nu^2w_{\varphi\varphi} + w = 1 + \varepsilon w^2$ with $w$ $2\pi$-periodic in $\varphi$. Order 1: $w_0 = 1 + e\cos\varphi$. Order $\varepsilon$: $w_1'' + w_1 = w_0^2 - 2\nu_1w_0'' = \dots + (2e + 2\nu_1e)\cos\varphi + \dots$ Solvability: $\nu_1 = -1$. So the radial oscillation has frequency $\nu = 1 - \varepsilon$ in $\theta$; the orbit returns to perihelion after $\Delta\theta = 2\pi/\nu = 2\pi(1 + \varepsilon + O(\varepsilon^2))$.

The figure. At $\varepsilon = 0.05$ the measured step between perihelia is $0.362$ rad $= 20.7^\circ$, against $2\pi\varepsilon = 0.314$ rad: $15\%$ more, because at so large an $\varepsilon$ the next order is visible (the next slide shows that excess growing in proportion to $\varepsilon$). The advance is in the direction of motion — the planet goes round counterclockwise here, and so does the perihelion. Distances are $r/p$ with the Sun at the origin.

21 · the orbit

Mercury: 43 seconds of arc a century

Mercury

$2\pi\varepsilon = 5.019\times10^{-7}$ rad $= 0.1035''$ per orbit. With a period of $87.969$ days, $415.20$ orbits per century:

$$\Delta\varpi = 42.98''\ \text{per century}.$$

The residual Le Verrier (1859, $38''$) and Newcomb (1882, $43''$) could not explain with any Newtonian perturbation.

checks

Dimensionless ✓. $c \to \infty$: no advance ✓. Numerically, integrating $w'' + w = 1 + \varepsilon w^2$ at $\varepsilon = 10^{-3}$, $e = 0.2$: advance per orbit $0.0062990$ against $2\pi\varepsilon = 0.0062832$ — agreement to $0.25\% \approx 2.5\varepsilon$, the size of the next order ✓.

How wrong $2\pi\varepsilon$ is: the integrated advance over $2\pi\varepsilon$, minus 1. The excess grows like $2.5\varepsilon$ — the next order; the ring is the check in the box, $0.25\%$ at $\varepsilon = 10^{-3}$. For Mercury it is $2\times10^{-7}$ of the $43''$.

Mercury's numbers. $2\pi\times7.9875\times10^{-8} = 5.0187\times10^{-7}$ rad $= 0.10352''$; the sidereal period $87.9691$ d gives $36525/87.9691 = 415.203$ orbits per Julian century, and $415.203\times0.10352'' = 42.98''$.

The numerical check. Integrated from $w = 1.2$, $w' = 0$ with relative tolerance $10^{-12}$, locating successive maxima of $w$. At $\varepsilon = 10^{-2}$ the advance is $0.064461$ against $2\pi\varepsilon = 0.062832$, $2.6\%$ high — ten times the discrepancy at $\varepsilon = 10^{-3}$, so the relative error is $O(\varepsilon)$ as it should be. For a near-circular orbit the exact frequency is $\nu = \sqrt{1 - 2\varepsilon w_c}$ with $w_c = (1 - \sqrt{1 - 4\varepsilon})/2\varepsilon$, which at $\varepsilon = 10^{-2}$ gives an advance of $0.06445$, matching.

The figure. Dots: the orbit equation integrated at $e = 0.2$ for nine values of $\varepsilon$, the advance measured between successive maxima of $w$. Solid: the exact advance of a near-circular orbit, $2\pi/\nu - 2\pi$ with $\nu = \sqrt{1 - 2\varepsilon w_c}$, whose expansion is $2\pi\varepsilon(1 + \tfrac52\varepsilon + \cdots)$; dashed, its first term $2.5\varepsilon$. The eccentricity changes the excess only slightly. At Mercury's $\varepsilon = 8\times10^{-8}$ the correction is $2.5\varepsilon = 2\times10^{-7}$ of the advance, about $10^{-5}$ seconds of arc per century.

History. Le Verrier announced the anomaly in 1859 as $38''$ per century; Newcomb's 1882 re-analysis put it at $43''$. Proposed explanations — an unseen planet (Vulcan), a dust ring, a modified inverse-square law — each failed another test. Einstein's November 1915 calculation — a first-order expansion in the same small parameter, though he evaluated the orbit integral rather than this differential equation — gave $43''$ with no adjustable parameter. He wrote that he was beside himself with joy for days.

21 · the orbit

Remove it, or keep it?

	Duffing	the orbit
secular term	$-\tfrac38\varepsilon A^3t\sin t$	$\varepsilon e\,\theta\sin\theta$
what it encodes	frequency $1 + \tfrac38\varepsilon A^2$	frequency $1 - \varepsilon$
its growth	false: $\|x\| \le A$	false: the orbit does not swell
its phase	not observed	accumulates, and is observed

In both, the secular term is a frequency shift written badly. "Removing" it means re-writing it correctly — never discarding it.

why it must be kept

A Newtonian bound orbit closes: $w$ is $2\pi$-periodic in $\theta$. Demand that of $w_1$ and the solvability condition requires the $\cos\theta$ coefficient of $w_0^2$, which is $2e$, to vanish — false for any non-circular orbit. The equation forbids a closed orbit. The only way to satisfy it is to let the radial frequency differ from 1, and that difference is the precession.

and the naive form is enough

The naive series fails when $\varepsilon\theta \sim 1$: for Mercury, after $2.0$ million orbits, $480\,000$ years. Observations span centuries. Over them, $\varepsilon e\,\theta\sin\theta$ is uniformly accurate — and its coefficient is the $43''$.

The distinction, stated once and carefully. The phrase "remove secular terms" is misleading and is the source of the confusion this section exists to clear up. What Poincaré–Lindstedt and multiple scales remove is the non-uniformity: the false algebraic growth that comes from Taylor-expanding $\cos((1 + \delta)t)$ in $\delta$. What they keep, and compute, is the frequency shift $\delta$ that the secular term's coefficient encodes. For Duffing that shift is interesting but nobody measures the accumulated phase; for the orbit the accumulated phase is the longitude of perihelion, and astronomers have measured it for centuries. Discarding the secular term outright — setting its coefficient to zero because "secular terms are unphysical" — would give zero precession, and it is not even consistent: the solvability argument on the slide shows the equation has no $2\pi$-periodic solution at first order unless $e = 0$.

The circular orbit. If $e = 0$ there is no secular term at all: a circular orbit has no perihelion to advance. The equation then just shifts the radius, $w = 1 + \varepsilon + O(\varepsilon^2)$. ✓

How long the naive answer lasts, in numbers. $1/\varepsilon = 1.252\times10^7$ radians, which is $1.99\times10^6$ orbits, which at $87.969$ days each is $4.80\times10^5$ years. Over a century the naive first-order error is of order $(\varepsilon\theta)^2$ with $\varepsilon\theta = 7.99\times10^{-8}\times2609 = 2.1\times10^{-4}$ — invisible. So unit 37 can, honestly, use the naive expansion and read the precession off the secular coefficient, provided it says what this slide says.

The same trick, elsewhere. The binary pulsar PSR B1913+16 has $\varepsilon = 1.04\times10^{-5}$, 130 times Mercury's; the problem sheet computes its periastron advance, about $4.2^\circ$ per year.

21 · summary

What this unit bought

the toolkit

Nondimensionalise to find $\varepsilon$; classify: regular, singular in the equation, or singular in the domain.
Asymptotic, not convergent: error $o(\phi_N)$ at fixed $N$. Truncate at the smallest term; the error is exponentially small.
Secular terms are frequency shifts expanded badly. Poincaré–Lindstedt for periodic motion, multiple scales otherwise, both through the solvability condition.
Boundary layers: stretch by dominant balance, solve inner and outer, match, composite.
WKB: phase from the eikonal equation, amplitude $q^{-1/4}$ from transport, valid while the medium varies slowly; connection formulae at turning points.

Reading

Hinch, Perturbation Methods (CUP, 1991): the chapters on algebraic equations, asymptotic approximations, matched expansions and multiple scales. Skip the one on integrals until after unit 20. Short, and exactly at this level.

Bender and Orszag, Advanced Mathematical Methods for Scientists and Engineers: ch. 7 (perturbation series), 9 (boundary layers), 10 (WKB, stopping before the Airy asymptotics), 11 (multiple scales). Skip ch. 6 until after unit 20: it uses complex analysis.

22 continuum mechanics · 37 the classical tests, where the orbit is derived and closed · 38 gravitational waves, first-order perturbation theory of spacetime.

What was left out, deliberately. Integral asymptotics — Laplace's method, stationary phase, steepest descent — are unit 20's, and appear here only as citations (Stirling's formula, the Airy asymptotics behind the connection formulae). Averaging theory as a theorem (Bogoliubov–Krylov), KAM, and small divisors belong to a course on Hamiltonian dynamics. Renormalisation-group methods for singular perturbations, which unify all of the above, are a pleasant read afterwards: Chen, Goldenfeld and Oono, Physical Review E 54 (1996).

Reading notes. Hinch's chapter on asymptotic approximations is the best short account of asymptotic series anywhere; his chapter on matched expansions has the Van Dyke failures promised in the matching slide. In Bender and Orszag, the problem sets are graded by difficulty and the hardest are research-level; the easiest are exactly the right practice after this unit's sheet.

21 · checkpoint

Done when

checkpoint

Given a perturbed problem you can, on paper and without notes:

nondimensionalise it, name $\varepsilon$, and say whether it is regular or singular, in the equation or in the domain;
remove a secular term by Poincaré–Lindstedt or multiple scales, and get the frequency correction right including its sign;
construct a boundary layer, match it to the outer solution, and write the composite;
for a divergent asymptotic series, find the optimal truncation point and state the accuracy it achieves;
say, for the relativistic orbit equation, why the secular term must be kept.

Failing it is unambiguous

A layer at the wrong end; a stiffening spring that slows down; terms added past the smallest one; or "secular terms are unphysical, so the precession is zero".

The sheet: ./mr problems 21. When the checkpoint is true: ./mr done 21.

Try the checkpoint on something not in the deck before opening the sheet. A good one: the softening oscillator $\ddot x + x - \varepsilon x^3 = 0$ (a pendulum to its first nonlinear order, with $\varepsilon = 1/6$ in units where the amplitude is $A$). Poincaré–Lindstedt gives $\omega = 1 - \tfrac38\varepsilon A^2$, lower than 1 as a softening spring must be; at $\varepsilon = 1/6$ and $A = 0.4$ that is $\omega = 1 - 0.01 = 0.99$, and the exact pendulum frequency at $\theta_0 = 0.4$ rad is $0.990008$ — agreement to $10^{-5}$. Then decide, without computing, where the layer sits in $\varepsilon y'' - y' + y = 0$ on $[0,1]$ (at $x = 1$: the coefficient of $y'$ is negative), and how many terms of $\sum(-1)^nn!\,(0.05)^n$ to keep (twenty, up to $n = 19$, error $1.1\times10^{-8}$).

Every number in that paragraph was computed. If yours disagree, check your arithmetic before the deck.

Perturbation theory: the series that diverges and is still right

Two oscillators, one term

Small parameters, and where they come from

Nondimensionalise first

Regular or singular: the test

Where the lost root went: dominant balance

Asymptotic series

Order notation

Poincaré's definition

Convergent is not the same as useful

The two series, term by term

The Stieltjes series, and its remainder

Optimal truncation

Optimal truncation, plotted

Borel summation, named

Secular terms

Regular perturbation, the naive way

Worked: the naive first order

Why the naive answer dies at t ~ 1/ε

What the secular term really is

How long the naive answer is good: a Grönwall bound

Removing secular terms

The solvability condition

Worked: Poincaré–Lindstedt for Duffing

Checking the frequency

Multiple scales

Two clocks, plotted

Secular terms are a boundary layer in time

Boundary layers

A small parameter on the highest derivative

Inner and outer solutions

Matching: the intermediate variable and Van Dyke's rule

The overlap, on a logarithmic axis

Worked: the composite solution, against the exact one

Where the layer goes

WKB

The WKB approximation

When WKB is valid, and an adiabatic invariant

The adiabatic invariant, plotted

Turning points and the connection formulae

The turning point, pictured

The orbit: a secular term that is the answer

The relativistic orbit equation, nondimensionalised

Worked: the first-order orbit

Reading the secular term: the perihelion advances

Mercury: 43 seconds of arc a century

Remove it, or keep it?

What this unit bought

Reading

Next

Done when

Failing it is unambiguous

Next