Perturbation methods

Regular and singular perturbation, boundary layers and WKB.

In this chapter 3

Tick a unit to add it to your PDF basket.

Perturbation theory

Perturbation theory

Method

Regular perturbation expansions

Adapted

Add to PDF

Expands the unknown of a problem carrying a small parameter $\epsilon$ in powers of $\epsilon$ and solves the resulting hierarchy order by order: the unit gives the non-degeneracy condition that makes every coefficient computable, the recipe, and the secular terms that destroy the approximation after a time of order $1/\epsilon$ .

Prerequisites: Taylor's theorem, Implicit and inverse function theorems, Constant-coefficient equations, Variation of parameters, The damped and driven harmonic oscillator, Asymptotic series

Definition(A perturbation family and its regular expansion)AdaptedAdapted from Strogatz (1994), §7.6, Eq. (4), p. 216; Stakgold & Holst (2011), §9.3, Eqs. (9.3.1)-(9.3.2), pp. 584-585

Many problems of physics come as a family $P(\epsilon)$ indexed by a real parameter $\epsilon$ , in which the member $P(0)$ can be solved and the member of interest has $\epsilon$ small but non-zero. Stakgold and Holst write such a family out for a boundary-value problem; Strogatz writes one for an initial-value problem. In either case the unknown is a function $u(\,\cdot\,;\epsilon)$ , and the base problem is $P(0)$ with its known solution $u_0$ .

A regular perturbation expansion of order $k$ is the ansatz

u(\,\cdot\,;\epsilon) \;=\; u_0 + \epsilon\,u_1 + \epsilon^{2}u_2 + \cdots + \epsilon^{k}u_k + o(\epsilon^{k}) ,

in which the coefficients $u_1, u_2, \dots$ do not depend on $\epsilon$ and are to be found from $P(\epsilon)$ itself. Strogatz calls the technique regular perturbation theory. Three things are being assumed, and each has to be checked rather than hoped for.

The limit is of the same type. Setting $\epsilon = 0$ in $P(\epsilon)$ must leave a problem of the same kind — same order of equation, same number of boundary conditions, same number of roots — and not a degenerate one. When it does not, the family is singularly perturbed and this unit's method fails; see Singular perturbations and boundary layers.
The limit solution is non-degenerate. The base solution must be isolated in a sense made precise in the next block, so that $u_0$ is the unique starting point of a branch and not one of several merging ones.
The expansion is uniform on the region of interest. The error terms must be small uniformly in the independent variable over the whole interval or domain where the answer is wanted. For an initial-value problem posed on $0 \le t < \infty$ this is the assumption that fails most often, and the remark below says how.

Nothing in the ansatz requires the series to converge. In many problems the coefficients $u_j$ grow so fast that the series diverges for every $\epsilon \neq 0$ and is useful only as an asymptotic series, truncated at a finite order.

Method(Regular perturbation of a simple root)AdaptedAdapted from Stakgold & Holst (2011), §9.3, pp. 588-589

Let $U \subseteq \mathbb{R}^{2}$ be open, let $m \ge 1$ , and let $F : U \to \mathbb{R}$ be of class $C^{m}$ . Suppose $(x_\ast, 0) \in U$ and

F(x_\ast, 0) = 0 , \qquad \frac{\partial F}{\partial x}(x_\ast,0) \;\neq\; 0 .

Then:

(Existence of the branch.) There are $\delta > 0$ and $\eta > 0$ and exactly one function $x : (-\eta,\eta) \to (x_\ast - \delta, x_\ast + \delta)$ with $x(0) = x_\ast$ and $F(x(\epsilon),\epsilon) = 0$ for all $|\epsilon| < \eta$ ; moreover $x$ is of class $C^{m}$ .
(The expansion exists.) For every $k \le m$ there are constants $c_1, \dots, c_k$ with $x(\epsilon) \;=\; x_\ast + c_1\epsilon + c_2\epsilon^{2} + \cdots + c_k\epsilon^{k} + o(\epsilon^{k}) \qquad (\epsilon \to 0) ,$ namely $c_j = x^{(j)}(0)/j!$ .
(The coefficients are computable in order.) Substitute the truncated expansion into $F$ , expand in powers of $\epsilon$ by Taylor's theorem about $(x_\ast,0)$ , and set the coefficient of $\epsilon^{j}$ to zero. For each $j \le k$ the resulting equation has the form $\frac{\partial F}{\partial x}(x_\ast,0)\,c_j \;+\; R_j(c_1,\dots,c_{j-1}) \;=\; 0 ,$ where $R_j$ depends only on the earlier coefficients and on partial derivatives of $F$ at $(x_\ast,0)$ . Hence $c_1, \dots, c_k$ are determined successively and uniquely.

Stakgold and Holst use exactly this device — a scalar equation $F(\delta,\epsilon) = 0$ with $\partial F/\partial\delta \neq 0$ at the origin, resolved by the implicit function theorem — to show that a simple eigenvalue of a compact self-adjoint operator moves along a well-defined branch when the operator is perturbed. The statement above is the scalar case on its own, stripped of the operator setting.

Idea of the proofOriginal

Order-by-order bookkeeping looks like an act of faith: one writes down a series whose terms are unknown, substitutes it into the problem, and reads off equations. It works for one reason only. The unknown enters $F$ at first order through $\partial F/\partial x$ , so when the powers of $\epsilon$ are collected, the new coefficient $c_j$ appears in the $\epsilon^{j}$ equation linearly, multiplied by the single number $\partial F/\partial x(x_\ast,0)$ , and everything else in that equation is already known. If that number is not zero, each equation is a one-line division. If it is zero, the $\epsilon^{j}$ equation says nothing about $c_j$ , and the whole scheme collapses at the first order where this happens. Everything else — existence of the branch, smoothness, and the fact that the numbers produced really are the Taylor coefficients of a genuine solution — is the implicit function theorem.

ProofOriginal

The branch. Apply the implicit function theorem to $F$ at $(x_\ast,0)$ , solving for $x$ in terms of $\epsilon$ . The hypothesis $\partial F/\partial x(x_\ast,0) \neq 0$ is precisely its non-degeneracy condition, and $F \in C^{m}$ gives a solution function $x(\epsilon)$ of class $C^{m}$ , unique within a neighbourhood. This proves (1).
The expansion. Since $x \in C^{m}$ on $(-\eta,\eta)$ , Taylor's theorem gives, for each $k \le m$ ,
$x(\epsilon) = \sum_{j=0}^{k} \frac{x^{(j)}(0)}{j!}\,\epsilon^{j} + o(\epsilon^{k}) ,$
with $x^{(0)}(0) = x_\ast$ . This is (2), with $c_j = x^{(j)}(0)/j!$ .
The shape of the $\epsilon^{j}$ equation. Write $h = x - x_\ast$ and expand $F$ about $(x_\ast,0)$ by Taylor's theorem in two variables to total order $k$ :
$F(x_\ast + h,\ \epsilon) \;=\; \sum_{\substack{p+q \le k \\ p+q \ge 1}} \frac{1}{p!\,q!}\,\partial_x^{p}\partial_\epsilon^{q}F(x_\ast,0)\; h^{p}\epsilon^{q} \;+\; o\bigl((|h|+|\epsilon|)^{k}\bigr) ,$
the constant term being absent because $F(x_\ast,0) = 0$ . Now put $h = c_1\epsilon + \cdots + c_k\epsilon^{k} + o(\epsilon^{k})$ , which is $O(\epsilon)$ , so that $h^{p}\epsilon^{q} = O(\epsilon^{p+q})$ . Collect the coefficient of $\epsilon^{j}$ for a fixed $j \le k$ . A term with $p + q \le j$ can contribute, and the coefficient $c_j$ can appear in it only through a factor of $h$ carrying the power $\epsilon^{j}$ ; since $h = O(\epsilon)$ , that requires $p = 1$ and $q = 0$ . Every other contribution to the $\epsilon^{j}$ coefficient is built from $c_1, \dots, c_{j-1}$ and from derivatives of $F$ at $(x_\ast,0)$ . Writing those contributions as $R_j(c_1,\dots,c_{j-1})$ , the coefficient of $\epsilon^{j}$ is
$\partial_x F(x_\ast,0)\,c_j + R_j(c_1,\dots,c_{j-1}) .$
Why that coefficient vanishes. By (1), $F(x(\epsilon),\epsilon) = 0$ identically on $(-\eta,\eta)$ . Hence the function $\epsilon \mapsto F(x(\epsilon),\epsilon)$ , which is $C^{m}$ , has all derivatives zero at $\epsilon = 0$ up to order $m$ , so its Taylor polynomial of order $k$ is identically zero and every coefficient of $\epsilon^{j}$ , $j \le k$ , vanishes.
Solving. Combining steps 3 and 4, $\partial_x F(x_\ast,0)\,c_j = -R_j(c_1,\dots,c_{j-1})$ for $j = 1,\dots,k$ . Because $\partial_x F(x_\ast,0) \neq 0$ — the hypothesis used a second time, now not for existence but for solvability — this determines $c_j$ uniquely from its predecessors, by induction on $j$ . This is (3). $\square$

Remark(Where the hypotheses are used, and what fails without them)Original

The condition $\partial_x F(x_\ast,0) \neq 0$ carries the whole statement, and it does so twice: once in the implicit function theorem, to produce a branch at all, and once at every order, to solve for $c_j$ . Both failures are visible in two-line examples.

A double root. $F(x,\epsilon) = x^{2} - \epsilon$ has $F(0,0) = 0$ and $\partial_x F(0,0) = 0$ . The roots are $x = \pm\sqrt{\epsilon}$ : there are two branches for $\epsilon > 0$ and none for $\epsilon < 0$ , and neither branch has an expansion in powers of $\epsilon$ , because $\sqrt{\epsilon}$ is not differentiable at $0$ . The order- $\epsilon$ equation of the method confirms the diagnosis without any knowledge of the answer: it reads $0 \cdot c_1 = 1$ , which has no solution. The cure is to change the gauge, expanding in powers of $\epsilon^{1/2}$ instead of $\epsilon$ .
A root that escapes. $F(x,\epsilon) = \epsilon x^{2} + x - 1$ loses a root at $\epsilon = 0$ : the quadratic has two roots for $\epsilon \neq 0$ , the linear limit only one. The surviving branch is regular, $x = 1 - \epsilon + 2\epsilon^{2} + O(\epsilon^{3})$ . The other root is $x = -1/\epsilon - 1 + O(\epsilon)$ and runs off to infinity as $\epsilon \to 0$ ; no expansion in non-negative powers of $\epsilon$ can reach it. This is the algebraic model of the degeneracy that produces boundary layers in differential equations, where setting $\epsilon = 0$ lowers the order and a solution is lost; see Singular perturbations and boundary layers.
Finite smoothness. $F \in C^{m}$ bounds the number of terms. With $F$ merely $C^{1}$ the branch exists and is differentiable, so $c_1$ is meaningful, but $c_2$ need not exist.

Nothing above claims convergence. Even for $F$ real-analytic, where the branch is analytic and the series converges for small $|\epsilon|$ , the radius of convergence is set by the nearest value of $\epsilon$ at which two branches collide, which is usually not known in advance. For differential equations convergence is the exception; the expansion is normally an asymptotic series and the honest statement is the $o(\epsilon^{k})$ of the theorem, for fixed $k$ as $\epsilon \to 0$ .

Method(Regular perturbation, order by order)AdaptedAdapted from Strogatz (1994), §7.6, pp. 216-217; Stakgold & Holst (2011), §9.3, Eqs. (9.3.1)-(9.3.5), pp. 584-585

Non-dimensionalise first, then name $\epsilon$ . Scale the variables so that the small parameter is a dimensionless ratio. Without this step "small" has no meaning and the ordering of terms is arbitrary.
Inspect the limit problem. Set $\epsilon = 0$ and ask whether what remains is a problem of the same type with a solvable, isolated solution. If the order of the differential equation drops, a boundary condition becomes unenforceable, or roots disappear, stop: the problem is singular, and Singular perturbations and boundary layers is the relevant method.
Insert the ansatz. Write $u = u_0 + \epsilon u_1 + \epsilon^{2}u_2 + \cdots$ and substitute. Expand every nonlinear function of $u$ in a Taylor series about $u_0$ , and expand every coefficient that depends on $\epsilon$ as well.
Collect powers of $\epsilon$ and equate each coefficient to zero. This is legitimate because the equation is to hold for a range of $\epsilon$ , not one value: a function of $\epsilon$ that vanishes identically has all Taylor coefficients zero.
Distribute the data. Initial or boundary data that do not depend on $\epsilon$ are carried entirely by $u_0$ ; every higher coefficient gets the corresponding homogeneous data. If the data do depend on $\epsilon$ , expand them too and give the $\epsilon^{j}$ part to $u_j$ .
Solve the hierarchy downwards. The $O(1)$ problem is the base problem. Every later problem has the same linear operator — the linearisation of the base problem about $u_0$ — with a right-hand side assembled from the coefficients already found. One linear solve per order, in order.
Watch for a kernel. If that common operator has a non-trivial kernel, the order- $j$ problem is solvable only if its right-hand side satisfies a solvability condition (orthogonality to the kernel of the adjoint; see The Fredholm alternative and modified Green's functions). That condition is not an obstacle but the point: it determines a parameter left free at lower order, such as an eigenvalue or a frequency. In initial-value problems on an unbounded interval the same kernel shows up instead as resonant forcing and a secular term; see the next remark.
State the range of validity. Compare the first term dropped with the last term kept. The truncation at order $k$ is informative only where $|\epsilon^{k+1}u_{k+1}| \ll |\epsilon^{k}u_k|$ , and the coefficients $u_j$ often grow with the independent variable, which turns this into a restriction on the variable as well as on $\epsilon$ .

Example(Kepler's equation: the eccentric anomaly to third order)Original

A body on an elliptical orbit of eccentricity $e$ has eccentric anomaly $E$ related to the mean anomaly $M$ (which grows linearly in time) by Kepler's equation

E - e\sin E \;=\; M .

This is transcendental in $E$ and has no closed-form inverse, but $e$ is small for most planets, so expand in it. Take $M$ fixed and set $F(E,e) = E - e\sin E - M$ . At $e = 0$ the root is $E_\ast = M$ , and

\frac{\partial F}{\partial E}(M,0) \;=\; 1 - e\cos E\Big|_{e=0} \;=\; 1 \;\neq\; 0 ,

so the statement above applies: there is a unique smooth branch $E(e)$ with $E(0) = M$ . (In fact $\partial_E F = 1 - e\cos E \ge 1 - e > 0$ for all $E$ whenever $0 \le e < 1$ , so the root is unique and the branch is global in $E$ ; only the expansion in powers of $e$ is a local statement.)

Order by order. Write $E = M + eE_1 + e^{2}E_2 + e^{3}E_3 + O(e^{4})$ and expand $\sin E$ about $M$ :

\sin E = \sin M + (E - M)\cos M - \tfrac12 (E-M)^{2}\sin M + O\bigl((E-M)^{3}\bigr) .

Substituting and collecting powers of $e$ in $E - e\sin E - M = 0$ :

$O(e)$ : $\;E_1 - \sin M = 0$ , so $E_1 = \sin M$ .
$O(e^{2})$ : $\;E_2 - E_1\cos M = 0$ , so $E_2 = \sin M\cos M = \tfrac12\sin 2M$ .
$O(e^{3})$ : $\;E_3 - \bigl(E_2\cos M - \tfrac12 E_1^{2}\sin M\bigr) = 0$ , so $E_3 = \sin M\cos^{2}M - \tfrac12\sin^{3}M = \sin M - \tfrac32\sin^{3}M$ .

Note the shape promised by part 3 of the statement: at every order the new coefficient appears alone, multiplied by $\partial_E F(M,0) = 1$ , and everything else is already known. Hence

E \;=\; M + e\sin M + \tfrac12 e^{2}\sin 2M + e^{3}\bigl(\sin M - \tfrac32\sin^{3}M\bigr) + O(e^{4}) .

A number to check. Take $M = 1$ (radian) and $e = 0.2$ , close to Mercury's eccentricity. Solving $E - 0.2\sin E = 1$ numerically gives $E = 1.185324$ to six decimals. The truncations are

order kept	value	error
$0$	$1.000000$	$-1.85 \times 10^{-1}$
$1$	$1.168294$	$-1.70 \times 10^{-2}$
$2$	$1.186480$	$+1.16 \times 10^{-3}$
$3$	$1.186062$	$+7.4 \times 10^{-4}$

The first two corrections each cut the error by about an order of magnitude, as a series in $e = 0.2$ should. The third does less, because the $O(e^{4})$ coefficient $\tfrac13\sin 4M - \tfrac16\sin 2M$ happens to be large at $M = 1$ ; adding it gives $1.185416$ , error $+9.2\times10^{-5}$ . This is normal behaviour and not a sign of failure: the error after keeping $k$ terms is $O(e^{k+1})$ as $e \to 0$ for fixed $k$ , which says nothing about the ordering of individual terms at a fixed $e$ . (Kepler's equation is the rare case where the convergence question is settled: the series in $e$ converges for $e$ below the Laplace limit, about $0.6627$ .)

Remark(Secular terms: why the expansion expires at

t \sim 1/\epsilon

)AdaptedAdapted from Strogatz (1994), §7.6, Eqs. (15)-(17) and the two objections, pp. 217-218

On a bounded interval of the independent variable, a regular expansion with smooth coefficients is as good as its last term. On an unbounded interval it usually is not, and the reason is visible in the hierarchy of step 6 of the method.

Consider an oscillatory base problem, so that the common linear operator of the hierarchy is $d^{2}/dt^{2} + \omega_0^{2}$ and the leading solution $u_0$ oscillates at frequency $\omega_0$ . The right-hand side of the order- $\epsilon$ problem is built out of $u_0$ , so it generally contains a component at that same frequency $\omega_0$ . That is resonant forcing of an undamped oscillator, and its response grows linearly: $u_1$ contains a term proportional to $t\sin\omega_0 t$ or $t\cos\omega_0 t$ . Strogatz calls such a term, one that grows without bound in $t$ , a secular term.

A secular term does not mean the calculation is wrong. The bound in the theorem is an estimate at fixed $t$ as $\epsilon \to 0$ , and it holds. What a secular term destroys is uniformity: the two-term approximation $u_0 + \epsilon u_1$ is a genuine improvement on $u_0$ only while $\epsilon|u_1| \ll |u_0|$ , that is while

\epsilon t \;\ll\; 1 , \qquad\text{equivalently}\qquad t \;\ll\; 1/\epsilon .

Beyond that the "correction" is larger than the thing it corrects, and the formula is worse than useless — it can predict amplitudes the exact solution provably cannot reach, as the next example shows.

Strogatz makes two further observations about what has gone wrong, both worth keeping. First, a secular term is often the first Taylor coefficient of something perfectly bounded: $e^{-\epsilon t} = 1 - \epsilon t + O(\epsilon^{2}t^{2})$ , so a slow exponential decay of amplitude shows up in a truncated $\epsilon$ -expansion as linear growth, with the wrong sign and the wrong qualitative behaviour. Recovering the decay from the series would need infinitely many terms, which defeats the purpose. Second, a small shift in frequency, $\omega = \omega_0 + O(\epsilon)$ , also appears as a secular term, because $\cos(\omega_0 + \epsilon\omega_1)t = \cos\omega_0 t - \epsilon\omega_1 t\sin\omega_0 t + \cdots$ ; the phase error accumulates and becomes an $O(1)$ error after $t \sim 1/\epsilon$ . Both diagnoses say the same thing: the solution has two time scales, a fast one of order $1$ and a slow one of order $1/\epsilon$ , and a series in $\epsilon$ with coefficients depending on $t$ alone cannot represent the slow one. The standard repairs — rescaling the frequency before expanding, or treating $t$ and $\epsilon t$ as independent variables, which Strogatz calls two-timing (§7.6, p. 218) — build the slow scale into the ansatz from the start. They are not developed here.

Example(A weakly nonlinear oscillator, and the failure of the two-term formula)Original

Take the Duffing equation (Strogatz, §7.6, Eq. (3), p. 215) with a soft start from rest:

\ddot x + x + \epsilon x^{3} = 0 , \qquad x(0) = a , \quad \dot x(0) = 0 , \qquad 0 < \epsilon \ll 1 .

The hierarchy. Put $x = x_0 + \epsilon x_1 + O(\epsilon^{2})$ and expand $x^{3} = x_0^{3} + O(\epsilon)$ . The initial data do not depend on $\epsilon$ , so by step 5 they go entirely to $x_0$ :

O(1): \quad \ddot x_0 + x_0 = 0 , \quad x_0(0) = a, \ \dot x_0(0) = 0 ; \qquad O(\epsilon): \quad \ddot x_1 + x_1 = -x_0^{3} , \quad x_1(0) = 0, \ \dot x_1(0) = 0 .

The base solution is $x_0 = a\cos t$ . For the forcing, use $\cos^{3}t = \tfrac34\cos t + \tfrac14\cos 3t$ :

\ddot x_1 + x_1 \;=\; -\tfrac34 a^{3}\cos t \;-\; \tfrac14 a^{3}\cos 3t .

The first term on the right is at the natural frequency of the operator on the left — exactly the resonance described in the remark. Solving by variation of parameters (or by the standard resonant trial $Ct\sin t$ ) and imposing $x_1(0) = \dot x_1(0) = 0$ ,

x_1(t) \;=\; -\frac{3a^{3}}{8}\,t\sin t \;+\; \frac{a^{3}}{32}\bigl(\cos 3t - \cos t\bigr) ,

so the two-term approximation is

x(t) \;\approx\; a\cos t \;+\; \epsilon\left[-\frac{3a^{3}}{8}t\sin t + \frac{a^{3}}{32}(\cos 3t - \cos t)\right] .

The first bracketed term is secular; the second is a bounded third-harmonic correction, the signature of a nonlinearity.

Numbers. Take $a = 1$ and $\epsilon = 0.1$ . Integrating the equation numerically to twelve-digit tolerance and comparing:

$t$	numerical $x(t)$	two-term formula	error
$1$	$0.504662$	$0.503965$	$-7.0\times10^{-4}$
$5$	$0.449763$	$0.460200$	$+1.0\times10^{-2}$
$20$	$-0.305609$	$-0.280878$	$+2.5\times10^{-2}$
$40$	$-0.805957$	$-1.779979$	$-9.7\times10^{-1}$

At $t = 1$ the formula is good to four decimals. At $t = 40$ it is not merely inaccurate but impossible: energy conservation gives

\tfrac12\dot x^{2} + \tfrac12 x^{2} + \tfrac{\epsilon}{4}x^{4} \;=\; \tfrac12 a^{2} + \tfrac{\epsilon}{4}a^{4} ,

and since $\tfrac12\dot x^{2} \ge 0$ while $s \mapsto \tfrac12 s + \tfrac{\epsilon}{4}s^{2}$ is strictly increasing for $s \ge 0$ , this forces $x(t)^{2} \le a^{2}$ , that is $|x(t)| \le a = 1$ , for all $t$ . The two-term formula reports $-1.78$ . The crossover is where the remark predicts it: $\epsilon t = 1$ at $t = 10$ , and the error has grown to a few per cent by $t = 20$ and to order one by $t = 40$ , where $\epsilon t = 4$ .

What the secular term was trying to say. Guess instead that the nonlinearity only shifts the frequency, $x \approx a\cos\omega t$ with $\omega = 1 + \tfrac38\epsilon a^{2}$ . Expanding in $\epsilon$ at fixed $t$ reproduces the secular term exactly:

a\cos\bigl[(1 + \tfrac38\epsilon a^{2})t\bigr] \;=\; a\cos t - \frac{3a^{3}}{8}\epsilon\,t\sin t + O(\epsilon^{2}t^{2}) .

So the growth was the expansion of a phase, not of an amplitude — the second of the two diagnoses in the remark. The frequency-shifted formula, with $\omega = 1.0375$ , gives $-0.790424$ at $t = 40$ against the numerical $-0.805957$ , an error of $1.6\times10^{-2}$ rather than $9.7\times10^{-1}$ , and it stays bounded for all time. Finding the shift $\tfrac38\epsilon a^{2}$ without guessing is what the methods named at the end of the remark do.

Application(First-order shifts: eigenvalues, energy levels, orbits)AdaptedAdapted from Stakgold & Holst (2011), §9.3, Eqs. (9.3.23) and (9.3.26)-(9.3.30), pp. 589-590; Riley, Hobson & Bence (2006), §19.2, p. 671

The single most used consequence of this method in physics is the first-order shift of an eigenvalue, and it is a case where step 7 of the method — the solvability condition — does the work.

Stakgold and Holst take the Sturm-Liouville family

-u'' + (1 + \epsilon x^{2})u = \mu u , \qquad 0 < x < 1 , \qquad u(0) = u(1) = 0 .

At $\epsilon = 0$ the eigenvalues are simple, $\mu_n = n^{2}\pi^{2} + 1$ , with normalised eigenfunctions $e_n(x) = \sqrt{2}\sin n\pi x$ . Differentiating the problem with respect to $\epsilon$ and setting $\epsilon = 0$ gives an inhomogeneous problem for $\partial_\epsilon u$ whose operator is $-d^{2}/dx^{2} + 1 - \mu_n$ , and that operator annihilates $e_n$ . So the problem is solvable only if its right-hand side is orthogonal to $e_n$ — and that condition, which at first looks like an obstruction, is what determines the unknown $\partial_\epsilon\mu(0)$ :

\mu_n(\epsilon) \;\approx\; \mu_n + \epsilon\int_0^{1} x^{2}e_n^{2}(x)\,dx .

The first-order shift is computable without knowing how the eigenfunction moves. For this example $\int_0^1 x^{2}\cdot 2\sin^{2}(n\pi x)\,dx = \tfrac13 - 1/(2n^{2}\pi^{2})$ , so

\mu_n(\epsilon) \;\approx\; n^{2}\pi^{2} + 1 + \epsilon\left(\frac13 - \frac{1}{2n^{2}\pi^{2}}\right) .

For $n = 1$ and $\epsilon = 0.1$ this predicts $\mu_1 \approx 10.897872$ ; a finite-difference computation of the eigenvalue gives $10.897861$ , an agreement to five decimals, and for $\epsilon = 1$ the prediction $11.152277$ against the computed $11.151164$ is still good to three decimals.

In quantum mechanics the same formula is the first thing one reaches for. Riley, Hobson and Bence quote it in the form a physicist uses: if a Hamiltonian acquires an extra small term $H'$ , then to first order the ground state's energy moves by the expectation value $\langle 0|H'|0\rangle$ , evaluated in the unperturbed ground state. That is the display above with $\int_0^1 x^{2}e_n^{2}$ read as $\langle n | x^{2} | n\rangle$ and $\epsilon x^{2}$ as $H'$ : Stark and Zeeman shifts, fine structure and the anharmonic corrections to a vibrational spectrum are all computed this way, and the unperturbed eigenfunctions suffice.

Two further standing uses are worth naming. Celestial mechanics is where the method began: the expansion of Kepler's equation in the eccentricity, as in the example above, together with expansions of planetary orbits in the mass ratios, is how ephemerides were computed for two centuries. And in fluid and plasma problems the same scheme appears as a stability calculation — expand about a known equilibrium in the amplitude of a disturbance and read the growth rate off the order- $\epsilon$ problem — in which the solvability condition of step 7 becomes the dispersion relation.

Method · advanced

Singular perturbations and boundary layers

Adapted

Add to PDF

Treats problems in which setting the small parameter to zero lowers the order of the equation, so that a boundary condition can no longer be imposed: the unit proves that the solution of a linear model problem splits into a smooth outer part plus a correction confined to a layer of width $\epsilon$ , gives the recipe of inner and outer expansions and matching, and shows how dominant balance fixes the width of the layer.

Prerequisites: Regular perturbation expansions, Asymptotic series, Watson's lemma, Constant-coefficient equations, Solving first-order equations

Definition(Singular perturbation, outer region, boundary layer)AdaptedAdapted from Evans (2010), §4.5.1, p. 200; Strogatz (1994), §3.5, A Singular Limit, pp. 68-69

Let $P(\epsilon)$ be a family of differential equations depending on a parameter $\epsilon > 0$ . The feature that earns the name singular perturbation is a failure of continuity in $\epsilon$ : however small the parameter is taken, the solution of $P(\epsilon)$ need not resemble the solution of $P(0)$ . Evans accounts for this by the weight carried by the derivatives of highest order, which are informally the terms that settle how a solution behaves; the perturbations he calls singular are exactly those that attach the factor $\epsilon$ to a term of order higher than any already present, so that the term governing the perturbed equation is the one that vanishes in the limit.

For a two-point problem on an interval the consequence is concrete. If the equation is of second order for every $\epsilon > 0$ but of first order at $\epsilon = 0$ , then the limit problem accepts one boundary condition and not two. Strogatz calls such a limit singular, precisely because the data can no longer be imposed, and names the subject that deals with it singular perturbation theory. The method of Regular perturbation expansions does not apply: its first requirement is that $P(0)$ be a problem of the same type as $P(\epsilon)$ .

What happens instead is that the solution ceases to have a single scale. Fix a point $x_0$ of the closed interval and a positive function $\delta(\epsilon)$ with $\delta(\epsilon) \to 0$ .

The outer region is the part of the interval on which $x - x_0$ is of order $1$ . There the solution is smooth and has a regular expansion in powers of $\epsilon$ , called the outer expansion, whose leading term solves the limit problem $P(0)$ .
The inner region, or boundary layer, is the set where $x - x_0$ is of order $\delta(\epsilon)$ . There the right variable is the inner variable $X = (x - x_0)/\delta(\epsilon)$ , and the expansion of the solution as a function of $X$ is the inner expansion. The number $\delta(\epsilon)$ is the width of the layer.
Matching is the requirement that the two expansions describe the same function in the intermediate range $\delta(\epsilon) \ll |x - x_0| \ll 1$ , where both are supposed to hold.

Strogatz's term boundary layer comes from fluid mechanics, where such a region is a thin sheet of fluid next to a solid surface; he also applies it to a thin interval of time at the start of a motion, and we shall do the same. Nothing in the definition requires $x_0$ to be an endpoint: layers can sit at an interior point, or, as in Evans's dye problem, along a curve inside a two-dimensional region.

Method(A linear model problem)Original

Let $a > 0$ and $b, \alpha, \beta$ be real constants, and for $\epsilon > 0$ consider

\epsilon y'' + a y' + b y \;=\; 0 \quad \text{on } [0,1] , \qquad y(0) = \alpha , \quad y(1) = \beta .

Define the outer solution $y_{\mathrm{out}}$ as the solution of the limit problem $a y' + by = 0$ , $y(1) = \beta$ , namely

y_{\mathrm{out}}(x) \;=\; \beta\, e^{(b/a)(1-x)} ,

and the composite approximation

y_{\mathrm{comp}}(x) \;=\; y_{\mathrm{out}}(x) \;+\; \bigl(\alpha - y_{\mathrm{out}}(0)\bigr)\, e^{-ax/\epsilon} .

Then there are numbers $\epsilon_0 > 0$ and $C$ , both depending only on $a$ , $b$ , $\alpha$ , $\beta$ and not on $\epsilon$ , such that for every $\epsilon \in (0,\epsilon_0)$ the boundary-value problem has exactly one solution $y_\epsilon$ and

\sup_{0 \le x \le 1} \bigl| y_\epsilon(x) - y_{\mathrm{comp}}(x) \bigr| \;\le\; C\,\epsilon .

The approximation is therefore uniform on the whole interval, including inside the layer. Note what each piece does: $y_{\mathrm{out}}$ satisfies the condition at $x = 1$ exactly and the one at $x = 0$ not at all; the added exponential is of size $\alpha - y_{\mathrm{out}}(0)$ at $x = 0$ , repairs the condition there, and is smaller than any power of $\epsilon$ once $x$ is of order $1$ .

Idea of the proofOriginal

A constant-coefficient equation can be solved exactly, so the whole phenomenon is visible in the two roots of $\epsilon m^{2} + am + b = 0$ . Their product is $b/\epsilon$ and their sum is $-a/\epsilon$ , so one of them must be enormous: as $\epsilon \to 0$ one root stays near $-b/a$ , which is the root of the limit equation $am + b = 0$ , while the other runs off to $-a/\epsilon$ . The first gives a smooth function of $x$ ; the second gives $e^{-ax/\epsilon}$ , which equals $1$ at $x = 0$ and is beyond-all-orders small as soon as $x$ exceeds a few multiples of $\epsilon$ . So the exact solution is a smooth outer part plus a spike at the left end whose only duty is to fix the value at $x = 0$ , and the layer width $\epsilon$ is not assumed but read off from the lost root.

The proof has only to make this quantitative. Two estimates do it: each root differs from its stated limit by $O(\epsilon)$ , and the coefficients in front of the two exponentials differ from the ones in $y_{\mathrm{comp}}$ by $O(\epsilon)$ . One subtlety is that the fast exponential is $e^{m_-x}$ rather than $e^{-ax/\epsilon}$ , and these differ by a factor $e^{cx}$ with $c$ of order $1$ — an $O(x)$ relative error, not an $O(\epsilon)$ one. It is harmless because it is multiplied by $e^{-ax/\epsilon}$ , and $x e^{-ax/\epsilon}$ never exceeds $\epsilon/(ea)$ .

ProofOriginal

Throughout, $0 < \epsilon < \epsilon_1 := a^{2}/\bigl(8(1+|b|)\bigr)$ , and $C_1, C_2, \dots$ denote constants built from $a,b,\alpha,\beta$ alone. A second explicit restriction on $\epsilon$ is imposed in step 6; the $\epsilon_0$ of the statement is the smaller of the two bounds.

The two roots. The characteristic equation $\epsilon m^{2} + am + b = 0$ has discriminant $a^{2} - 4\epsilon b$ . Since $4\epsilon|b| \le 4\epsilon_1(1+|b|) = a^{2}/2$ , we have
$\tfrac12 a^{2} \;\le\; a^{2} - 4\epsilon b \;\le\; \tfrac32 a^{2} ,$
so $s := \sqrt{a^{2}-4\epsilon b}$ is real and $a/2 \le s \le 2a$ . The roots $m_\pm = (-a \pm s)/(2\epsilon)$ are real and distinct.
The slow root. Rationalising, $m_+ = \dfrac{s^{2}-a^{2}}{2\epsilon(s+a)} = \dfrac{-4\epsilon b}{2\epsilon(s+a)} = \dfrac{-2b}{s+a}$ . Hence, using $a - s = (a^{2}-s^{2})/(a+s) = 4\epsilon b/(a+s)$ and $s + a \ge a$ ,
$\Bigl| m_+ + \frac{b}{a} \Bigr| \;=\; |b|\,\frac{|a - s|}{a(s+a)} \;=\; \frac{4\epsilon b^{2}}{a(s+a)^{2}} \;\le\; \frac{4b^{2}}{a^{3}}\,\epsilon \;=:\; K_1\epsilon .$
In particular $|m_+| \le |b|/a + K_1\epsilon_1 =: K_2$ , a bound independent of $\epsilon$ .
The fast root. From $m_- = -(a+s)/(2\epsilon)$ and $s \ge a/2$ we get $m_- \le -3a/(4\epsilon) < 0$ . Also
$m_- + \frac{a}{\epsilon} \;=\; \frac{a - s}{2\epsilon} \;=\; \frac{2b}{a+s} \;=:\; c , \qquad |c| \;\le\; \frac{2|b|}{a} \;=:\; c_0 ,$
so that $e^{m_-x} = e^{-ax/\epsilon}e^{cx}$ exactly. This is where $a > 0$ is used: it makes $m_-$ negative, so the fast mode decays as $x$ increases and the layer sits at the left endpoint.
Existence and uniqueness. Because the roots are distinct, every solution is $y = Ae^{m_+x} + Be^{m_-x}$ . The boundary conditions give the linear system $A + B = \alpha$ , $Ae^{m_+} + Be^{m_-} = \beta$ , whose determinant is $e^{m_-} - e^{m_+} \ne 0$ since $m_- \ne m_+$ and the exponential is injective on $\mathbb{R}$ . So there is exactly one solution $y_\epsilon$ .
A name for the exponentially small quantity. Put $E := e^{m_-} \le e^{-3a/(4\epsilon)}$ . For every $k$ there is a constant with $E \le C_k\epsilon^{k}$ , so any term of size $O(E)$ may be absorbed into the $O(\epsilon)$ of the conclusion.
The coefficient of the slow mode. Solving the system, $A = (\beta - \alpha E)/(e^{m_+} - E)$ . Let $A_0 := y_{\mathrm{out}}(0) = \beta e^{b/a}$ , so that $y_{\mathrm{out}}(x) = A_0e^{-(b/a)x}$ . First compare $A$ with $\beta e^{-m_+}$ :
$A - \beta e^{-m_+} \;=\; \frac{\beta - \alpha E}{e^{m_+}-E} - \frac{\beta}{e^{m_+}} \;=\; \frac{E\bigl(\beta - \alpha e^{m_+}\bigr)}{e^{m_+}\bigl(e^{m_+}-E\bigr)} .$
By step 2, $e^{-K_2} \le e^{m_+} \le e^{K_2}$ . By step 5, $E \le e^{-3a/(4\epsilon)}$ , and this is at most $\tfrac12 e^{-K_2}$ as soon as $\epsilon < \epsilon_2 := 3a/\bigl(4(K_2+\ln 2)\bigr)$ . Set $\epsilon_0 := \min(\epsilon_1,\epsilon_2)$ and assume $\epsilon < \epsilon_0$ from here on. Then $|e^{m_+} - E| \ge \tfrac12 e^{-K_2}$ , the numerator is at most $E(|\beta| + |\alpha|e^{K_2})$ , and $|A - \beta e^{-m_+}| \le C_3E$ . Next, by step 2 and $|e^{u}-e^{v}| \le e^{\max(u,v)}|u-v|$ ,
$\bigl| \beta e^{-m_+} - A_0 \bigr| \;=\; |\beta|\,\bigl| e^{-m_+} - e^{b/a} \bigr| \;\le\; |\beta|e^{K_2}K_1\epsilon .$
Adding, $|A - A_0| \le C_4\epsilon$ . Since $A + B = \alpha$ and $B_0 := \alpha - A_0$ , also $|B - B_0| \le C_4\epsilon$ .
The slow part of the error. For $0 \le x \le 1$ ,
$A e^{m_+x} - y_{\mathrm{out}}(x) \;=\; (A - A_0)e^{m_+x} \;+\; A_0\bigl( e^{m_+x} - e^{-(b/a)x} \bigr) .$
The first term is at most $C_4\epsilon\,e^{K_2}$ in absolute value. For the second, $|m_+x + (b/a)x| \le K_1\epsilon$ by step 2 and $x \le 1$ , so again by $|e^{u}-e^{v}| \le e^{\max(u,v)}|u-v|$ the bracket is at most $e^{K_2}K_1\epsilon$ . Hence this part of the error is $\le C_5\epsilon$ , uniformly in $x$ .
The fast part of the error. Using step 3,
$B e^{m_-x} - B_0e^{-ax/\epsilon} \;=\; e^{-ax/\epsilon}\Bigl[ (B - B_0)e^{cx} + B_0\bigl(e^{cx}-1\bigr) \Bigr] .$
Since $e^{-ax/\epsilon} \le 1$ and $|e^{cx}| \le e^{c_0}$ , the first inner term contributes at most $C_4e^{c_0}\epsilon$ . For the second, $|e^{cx}-1| \le |c|x\,e^{c_0} \le c_0e^{c_0}x$ , so its contribution is at most
$|B_0|c_0e^{c_0}\,\sup_{x \ge 0} x\,e^{-ax/\epsilon} \;=\; |B_0|c_0e^{c_0}\,\frac{\epsilon}{ea} ,$
the supremum being attained at $x = \epsilon/a$ . So this part of the error is also $\le C_6\epsilon$ .
Conclusion. Since $y_\epsilon = Ae^{m_+x}+Be^{m_-x}$ and $y_{\mathrm{comp}} = y_{\mathrm{out}} + B_0e^{-ax/\epsilon}$ , steps 7 and 8 give $|y_\epsilon - y_{\mathrm{comp}}| \le (C_5+C_6)\epsilon$ for every $x \in [0,1]$ . Take $C = C_5 + C_6$ . $\square$

Remark(Where the hypotheses are used, and what fails without them)Original

Two boundary conditions against a first-order limit. This is the singularity itself, and it is why Regular perturbation expansions cannot be used. A regular expansion $y = Y_0 + \epsilon Y_1 + \cdots$ with coefficients independent of $\epsilon$ puts $aY_0' + bY_0 = 0$ at leading order, a first-order equation; it can meet $y(0)=\alpha$ or $y(1)=\beta$ , not both, and no later term can help, because every $Y_j$ obeys the same first-order operator with homogeneous data at whichever end was used. The expansion does not fail by being inaccurate; it fails by not existing.
$a > 0$ decides which condition the outer solution keeps. It enters at step 3, where it makes the fast root negative. If $a < 0$ the fast mode is $e^{|a|x/\epsilon}$ , which is negligible except near $x = 1$ : the layer moves to the right endpoint and the outer solution must carry $y(0) = \alpha$ instead. The stated formula is then not merely inaccurate but divergent. Take $a = -1$ , $b = 0$ , $\alpha = 1$ , $\beta = 0$ . The exact solution of $\epsilon y'' - y' = 0$ with these data is $y_\epsilon(x) = \bigl(e^{1/\epsilon} - e^{x/\epsilon}\bigr)/\bigl(e^{1/\epsilon}-1\bigr)$ , which is close to $1$ for $x$ away from $1$ and drops to $0$ in a layer of width $\epsilon$ at $x = 1$ ; the formula of the statement would return $y_{\mathrm{out}} \equiv 0$ plus $e^{x/\epsilon}$ , which at $x=1$ is $e^{1/\epsilon}$ .
$a$ must not vanish, and the width depends on which terms balance. Every bound in steps 2 and 3 is expressed through $1/a$ and becomes vacuous as $a$ shrinks. If $a = 0$ the equation is $\epsilon y'' + by = 0$ and the analysis above says nothing, correctly so: for $b < 0$ the solution has layers of width $\epsilon^{1/2}$ at both ends, and for $b > 0$ it oscillates with wavelength of order $\epsilon^{1/2}$ over the whole interval and the problem is not even uniquely solvable at the values $\epsilon = b/(n^{2}\pi^{2})$ , $n \in \mathbb{N}$ . The lesson is that $\epsilon$ is the width only when the first-derivative term is present; step 4 of the method below computes the width in general.
A coefficient that changes sign. If $a$ is replaced by a function $a(x)$ the statement goes through, with $e^{-ax/\epsilon}$ replaced by $\exp\Bigl(-\frac{1}{\epsilon}\int_0^{x} a(t)\,dt\Bigr)$ , as long as $a(x) \ge a_0 > 0$ throughout; the worked example below is of this kind. Where $a$ changes sign the structure is different again: the layer moves into the interior, and whether it exists at all depends on the direction of the sign change.
The bound is asymptotic, not a guarantee at a given $\epsilon$ . The conclusion is $C\epsilon$ with $C$ unknown, and the example below settles at $C \approx 1.88$ . At $\epsilon = 0.1$ that is an absolute error of $0.13$ at a point where the solution itself is about $1.5$ . Nothing in the theorem forbids that, and nothing in it says when $\epsilon$ is small enough.
$\epsilon > 0$ . For $\epsilon < 0$ the fast root is large and positive and the two-point problem is badly behaved: the sign of the small parameter is part of the data, not a detail.

Method(Matched asymptotic expansions)AdaptedAdapted from Evans (2010), §4.5.1, pp. 200, 203-204; Strogatz (1994), §3.5, pp. 64-68

Steps 1 and 4 follow the two cited works: the nondimensionalisation and two-scale diagnosis are Strogatz's, and the choice of rescaled variables that makes every quantity of order one is the step Evans carries out for the dye plume. Steps 2, 3 and 5 to 8 are the standard practice of the subject and are not taken from either work; they are our responsibility.

Nondimensionalise, then identify $\epsilon$ . Scale the variables so that the small parameter is a dimensionless group. Strogatz's warning is worth repeating. Sending the mass to zero is not a usable limit, because the mass multiplies more than one term and several of them vanish together; only after dividing through by a reference force, and choosing the time scale so that the derivatives are of order one, does a single dimensionless group appear whose smallness means something.
Check that the limit is singular. Set $\epsilon = 0$ . If the order of the equation drops, or a boundary or initial condition becomes unenforceable, continue here; otherwise use Regular perturbation expansions.
Outer expansion. In the region where the independent variable is of order $1$ , write $y \sim Y_0 + \epsilon Y_1 + \cdots$ and solve the reduced problem, imposing the condition at the end away from the layer. Which end that is, is fixed by the sign of the coefficient of the highest surviving derivative: the layer sits where the discarded fast mode decays into the interior. If you guess wrong, step 6 will have no solution — that is the diagnostic.
Fix the position and width by dominant balance. Put $x = x_0 + \delta(\epsilon)X$ and rewrite the equation in $X$ . Choose $\delta$ as the largest gauge for which the highest derivative returns at the same order as some other term; at that $\delta$ the lost term is back in play. For $\epsilon y'' + ay'$ this means $\epsilon/\delta^{2} \sim 1/\delta$ , so $\delta = \epsilon$ . For $\epsilon y'' + by$ , with no first derivative, it means $\epsilon/\delta^{2} \sim 1$ , so $\delta = \epsilon^{1/2}$ . Evans does exactly this for the dye plume, rescaling the transverse coordinate by $\epsilon^{-1/2}$ and the amplitude by $\epsilon^{1/2}$ so that both become quantities of order one; rescaling the unknown as well as the variable is often necessary, and the right powers are those that leave no large and no vanishing coefficient in the leading inner equation.
Inner expansion. Expand $y \sim \mathcal{Y}_0(X) + \delta\,\mathcal{Y}_1(X) + \cdots$ and impose the boundary condition at $x_0$ . The leading inner equation has constant coefficients in the typical case, because the coefficient functions are frozen at their values at $x_0$ , so it can be solved in closed form. It will leave one constant free: the inner problem has lost the boundary condition at the far end.
Match. Require that the two descriptions agree in the overlap. At leading order this is Prandtl's condition $\lim_{X \to \infty} \mathcal{Y}_0(X) \;=\; \lim_{x \to x_0} Y_0(x) ,$ and it is exactly what determines the constant left free in step 5. (If the limit on the left does not exist, the width chosen in step 4 is wrong.)
Composite approximation. Add the two and subtract what they share, so that the shared part is not counted twice: $y_{\mathrm{comp}}(x) \;=\; Y_0(x) \;+\; \mathcal{Y}_0\!\left(\frac{x-x_0}{\delta}\right) \;-\; \lim_{X\to\infty}\mathcal{Y}_0(X) .$ In the outer region the last two terms cancel to leading order and $Y_0$ survives; in the layer the first and last cancel and $\mathcal{Y}_0$ survives.
Check afterwards. Verify that $y_{\mathrm{comp}}$ meets both boundary conditions to the claimed order; that the terms dropped in each region really are smaller there than the ones kept; and that the inner solution settles to its limit over $X$ of order $1$ and not over some larger range, which would mean the width is wrong.

Example(A variable-coefficient problem solved exactly and asymptotically)Original

Take

\epsilon y'' + (1+x)\,y' + y \;=\; 0 \quad \text{on } [0,1] , \qquad y(0) = 0 , \quad y(1) = 1 ,

with $0 < \epsilon \ll 1$ . Here $a(x) = 1 + x \ge 1 > 0$ , so by the remark above the layer is at $x = 0$ .

The method. Outer. The reduced equation is $(1+x)Y_0' + Y_0 = 0$ , that is $\bigl[(1+x)Y_0\bigr]' = 0$ , so $Y_0 = c/(1+x)$ ; the condition at the far end, $Y_0(1) = 1$ , gives $c = 2$ and

Y_0(x) \;=\; \frac{2}{1+x} .

It fails the condition at $x = 0$ , where it takes the value $2$ instead of $0$ .

Width. With $x = \delta X$ the equation becomes $(\epsilon/\delta^{2})\mathcal{Y}'' + (1/\delta)(1+\delta X)\mathcal{Y}' + \mathcal{Y} = 0$ ; the first two coefficients match when $\delta = \epsilon$ , which is step 4.

Inner. With $X = x/\epsilon$ and $y(x) = \mathcal{Y}(X)$ , multiplying through by $\epsilon$ gives exactly

\mathcal{Y}'' + \mathcal{Y}' + \epsilon\bigl(X\mathcal{Y}' + \mathcal{Y}\bigr) \;=\; 0 ,

so at leading order $\mathcal{Y}_0'' + \mathcal{Y}_0' = 0$ , whence $\mathcal{Y}_0 = A + Be^{-X}$ . The boundary condition $\mathcal{Y}_0(0) = 0$ gives $B = -A$ .

Matching and composite. Step 6 requires $\lim_{X\to\infty}\mathcal{Y}_0 = A$ to equal $\lim_{x\to0}Y_0 = 2$ , so $A = 2$ and $\mathcal{Y}_0(X) = 2\bigl(1-e^{-X}\bigr)$ . The shared limit is $2$ , so

y_{\mathrm{comp}}(x) \;=\; \frac{2}{1+x} \;+\; 2\bigl(1-e^{-x/\epsilon}\bigr) \;-\; 2 \;=\; \frac{2}{1+x} \;-\; 2e^{-x/\epsilon} .

It satisfies $y(0) = 0$ exactly and $y(1) = 1 - 2e^{-1/\epsilon}$ , in error by an amount smaller than every power of $\epsilon$ .

An exact solution to test it against. This problem can be integrated in closed form, which is why it was chosen. Because $(1+x)y' + y = \bigl[(1+x)y\bigr]'$ , the equation is $\bigl[\epsilon y' + (1+x)y\bigr]' = 0$ , so $\epsilon y' + (1+x)y$ is a constant $K$ . That first-order linear equation, with $y(0)=0$ and the integrating factor $e^{\varphi(x)/\epsilon}$ where $\varphi(x) = x + \tfrac12x^{2}$ , gives

y(x) \;=\; \frac{K}{\epsilon}\,I(x) , \qquad I(x) \;:=\; \int_0^{x} e^{-(\varphi(x)-\varphi(t))/\epsilon}\,dt ,

and $y(1) = 1$ fixes $K/\epsilon = 1/I(1)$ , so that $y_\epsilon = I(x)/I(1)$ with no unknown constants. (Evaluating $I$ in terms of the imaginary error function is what we used for the numbers below.)

The same formula reproduces both expansions, which is a check on the method rather than on the arithmetic. For $x$ of order $1$ the integrand is concentrated within $O(\epsilon)$ of the upper limit, where $\varphi(x)-\varphi(t) \approx (1+x)(x-t)$ ; Watson's lemma gives $I(x) = \epsilon/(1+x) + O(\epsilon^{2})$ and $I(1) = \epsilon/2 + O(\epsilon^{2})$ , so $y_\epsilon(x) \to 2/(1+x)$ , the outer solution. For $x = \epsilon X$ the exponent is $\varphi(x)-\varphi(t) = (x-t) + O(\epsilon^{2})$ , so $I(\epsilon X) = \epsilon(1-e^{-X}) + O(\epsilon^{2})$ and $y_\epsilon \to 2(1-e^{-X})$ , the inner solution.

Numbers at $\epsilon = 0.01$ . The layer occupies roughly $0 \le x \lesssim 0.05$ , where the solution climbs from $0$ to about $1.9$ ; the outer part then decays slowly to $1$ .

$x$	exact $y_\epsilon(x)$	$y_{\mathrm{comp}}(x)$	$y_{\mathrm{comp}} - y_\epsilon$
$0.005$	$0.784349$	$0.776988$	$-7.36\times10^{-3}$
$0.02$	$1.707852$	$1.690114$	$-1.77\times10^{-2}$
$0.05$	$1.905716$	$1.891286$	$-1.44\times10^{-2}$
$0.10$	$1.828934$	$1.818091$	$-1.08\times10^{-2}$
$0.25$	$1.606397$	$1.600000$	$-6.40\times10^{-3}$
$0.50$	$1.335975$	$1.333333$	$-2.64\times10^{-3}$
$0.75$	$1.143745$	$1.142857$	$-8.88\times10^{-4}$

The rate. Scanning the whole interval:

$\epsilon$	$\max_{[0,1]}\lvert y_{\mathrm{comp}}-y_\epsilon\rvert$	located at $x =$	error $/\epsilon$
$0.1$	$1.328\times10^{-1}$	$0.170$	$1.33$
$0.01$	$1.801\times10^{-2}$	$0.0239$	$1.80$
$0.001$	$1.881\times10^{-3}$	$0.00255$	$1.88$

The error falls by a factor of ten each time $\epsilon$ does, and the ratio settles near $1.88$ : the uniform $O(\epsilon)$ of the model theorem, with an observed constant slightly below $2$ . The worst point moves towards the origin with $\epsilon$ , as it must, since the error lives in the layer.

One further number shows that the layer is genuinely of width $\epsilon$ and not an artefact. Differentiating the exact formula at the origin gives $y_\epsilon'(0) = 1/I(1)$ ; at $\epsilon = 0.01$ this is $199.4975$ , against $y_{\mathrm{comp}}'(0) = 2/\epsilon - 2 = 198$ , a relative discrepancy of $0.75\%$ . A slope of order $1/\epsilon$ across a rise of order $1$ is a layer of width $\epsilon$ .

Application(Strong damping, fast flow, and why centred differences oscillate)AdaptedAdapted from Strogatz (1994), §3.5, pp. 64-69; Evans (2010), §4.5.1, pp. 200-204; Stakgold & Holst (2011), §6.1 Example 5, p. 353; Stakgold & Holst (2011), §9.4, p. 612 and Remark 1, p. 616

An initial layer: the overdamped limit. Strogatz's bead of mass $m$ sliding on a hoop of radius $r$ spinning at rate $\omega$ , with damping constant $b$ , obeys $mr\ddot\phi + b\dot\phi = -mg\sin\phi + mr\omega^{2}\sin\phi\cos\phi$ . Measuring time in units of $T = b/mg$ and dividing by $mg$ turns this into

\epsilon\,\frac{d^{2}\phi}{d\tau^{2}} + \frac{d\phi}{d\tau} \;=\; f(\phi) , \qquad f(\phi) = \sin\phi\,(\gamma\cos\phi - 1) , \qquad \epsilon = \frac{m^{2}gr}{b^{2}} , \quad \gamma = \frac{r\omega^{2}}{g} ,

so "overdamped" has the precise meaning $b^{2} \gg m^{2}gr$ . This is step 1 of the method, and it is Strogatz's point that it cannot be skipped: the naive limit $m \to 0$ removes the gravitational and centrifugal terms along with the inertia.

The limit $\epsilon = 0$ leaves $d\phi/d\tau = f(\phi)$ , of first order, which cannot accept both the initial position and the initial velocity. The resolution is a layer at $\tau = 0$ of width $\epsilon$ . In the inner variable $\Theta = \tau/\epsilon$ , writing $\Omega = d\phi/d\tau$ , the system becomes $d\phi/d\Theta = \epsilon\Omega$ and $d\Omega/d\Theta = f(\phi) - \Omega$ , so at leading order the position is frozen at $\phi_0$ while the velocity relaxes,

\Omega(\Theta) \;=\; f(\phi_0) + \bigl(\Omega_0 - f(\phi_0)\bigr)e^{-\Theta} .

The velocity therefore forgets its initial value over a time $\tau \sim \epsilon$ and arrives at the value the reduced equation demands; after that the slow drift $d\phi/d\tau = f(\phi)$ is correct. This is precisely Strogatz's phase-plane picture: a trajectory starting off the curve $\Omega = f(\phi)$ feels an acceleration of order $1/\epsilon$ and is thrown onto that curve almost at once, then creeps along it to a fixed point. He calls the transient a boundary layer in time, a thin interval near $\tau = 0$ . The same calculation is the justification of what is called adiabatic elimination of a fast variable in chemical kinetics and in laser physics.

Fluid mechanics, where the name comes from. Strogatz names the high-Reynolds-number limit of fluid flow as the singular limit that produces the thin layers over an aircraft wing. The structure is the one above: dropping the viscous term lowers the order of the Navier–Stokes equations, so the no-slip condition at the surface can no longer be imposed, and it is restored inside a layer. Step 4 gives the thickness. With $U$ a free-stream speed, $L$ a chord length and $\nu$ the kinematic viscosity, the balance inside the layer is viscous diffusion across it against advection along it, $\nu U/\delta^{2} \sim U^{2}/L$ , so $\delta \sim \sqrt{\nu L/U} = L\,\mathrm{Re}^{-1/2}$ . The small parameter is $\epsilon = 1/\mathrm{Re}$ and the width is $\epsilon^{1/2}$ , not $\epsilon$ , because the advection that opposes diffusion runs along the wall rather than across it.

Transport with weak diffusion. Stakgold and Holst's diffusion-absorption problem $-u'' + q(x)u = f(x)$ on $0 < x < 1$ , with prescribed concentrations at two walls, is their standard illustration: small $q$ is a regular perturbation, whereas very large $q$ — the same thing as a small diffusion constant — requires singular perturbation. Take $q$ constant for the purpose of the scaling and divide by it: the equation becomes $\epsilon u'' - u = -f/q$ with $\epsilon = 1/q$ , which has no first-derivative term, so step 4 predicts layers of width $q^{-1/2}$ at both walls, with the interior settled at $u \approx f/q$ . Their Remark 1 on the nonlinear absorption problem $-\Delta U = -\lambda U^{p}$ on a bounded convex domain gives the corresponding statement independently, for absorption exponents $p<1$ only: in that range large $\lambda$ produces a dead core, a set on which $U$ vanishes identically, and the core grows until it covers all of the domain apart from a thin collar next to the boundary, whose thickness they give as $\lambda^{-1/2}$ . That the same power comes out of a one-line dominant balance is a useful check that step 4 is doing real work. The restriction on $p$ is not a technicality, and Stakgold and Holst make the point with the case $p = 1$ on $|x| < a$ , whose solution $U = \cosh(\sqrt{\lambda}x)/\cosh(\sqrt{\lambda}a)$ is strictly positive for every $\lambda$ : no dead core, and therefore no collar, although $U$ is still exponentially small away from the walls in a layer of the same width. What the two cases differ in is the Lipschitz behaviour of the absorption term at $U = 0$ , which they identify as the source of the phenomenon: the derivative of $U^{p}$ is infinite at $U = 0+$ only when $p < 1$ , and it is that infinite slope which lets a solution reach zero at a finite distance instead of merely approaching it.

Evans's example is of the same family but with the layer along a curve rather than at a wall. Dye injected at the origin into a steady plane flow $\mathbf{b}$ satisfies $\operatorname{div}(u\mathbf{b}) = \delta_0$ , which concentrates it on the streamline through the origin; adding a small diffusion, $-\epsilon\Delta u^{\epsilon} + \operatorname{div}(u^{\epsilon}\mathbf{b}) = \delta_0$ , spreads it into a plume. Evans estimates the plume's width at time $t$ along the streamline as $O\bigl((\epsilon t)^{1/2}\bigr)$ , and hence, the total mass being fixed, its height as $O(\epsilon^{-1/2})$ ; rescaling the normal coordinate by $\epsilon^{-1/2}$ and the density by $\epsilon^{1/2}$ so that both are of order one turns the elliptic equation into a parabolic one in the inner variables. The width is again $\epsilon^{1/2}$ , for the reason given above: across the plume the flow has no component to oppose the diffusion.

A numerical consequence. Discretising $\epsilon y'' + ay' = 0$ with centred differences on a uniform mesh of spacing $h$ gives $(1+P)y_{j+1} - 2y_j + (1-P)y_{j-1} = 0$ with $P = ah/(2\epsilon)$ , half the cell Péclet number $\mathrm{Pe} = ah/\epsilon$ , whose decaying mode is $r = (1-P)/(1+P)$ in place of the true $e^{-ah/\epsilon}$ . For $P > 1$ , that is for $\mathrm{Pe} > 2$ — any mesh that does not resolve the layer, $h > 2\epsilon/a$ — this is negative, and the computed solution oscillates from node to node instead of decaying. The remedies are the standard ones: upwind the first-derivative term, or refine the mesh inside the layer where the asymptotics say the layer is. The analysis of this unit is what tells one where to refine, and the composite approximation is a good initial guess for the iteration. See Finite-difference approximations, Stability and stiff equations and Advection, diffusion and reaction; the companion technique for problems whose small parameter multiplies the highest derivative but whose solutions oscillate rather than decay is The WKB approximation.

Method

The WKB approximation

Adapted

Add to PDF

Builds two approximate solutions of the equation $y'' + q(x)y = 0$ when $q$ varies slowly, by letting the phase accumulate as $\int\sqrt{q}$ and the amplitude follow as $q^{-1/4}$ ; gives the exact residual that measures the error, the precise validity criterion, the connection formulas that carry a solution through a turning point, and the quantization condition they produce.

Prerequisites: Asymptotic series, Structure of linear equations and the Wronskian, Airy functions

Definition(Slowly varying coefficients, the WKB ansatz, turning points)AdaptedAdapted from Riley, Hobson & Bence (2006), §25.7, Eq. (25.42), p. 895; Hassani (2013), §14.8, Eq. (14.43), p. 450

The subject is the single equation

\frac{d^{2}y}{dx^{2}} + q(x)\,y \;=\; 0 ,

which Hassani writes in exactly this form and Riley, Hobson and Bence write as $y'' + k_0^{2}f(x)y = 0$ , separating out a large constant wavenumber $k_0$ and setting $n(x) = \sqrt{f(x)}$ by analogy with a refractive index. The two notations are related by $q = k_0^{2}f$ , so that $\sqrt{q} = k_0 n$ ; both are used below, and translations are given where it matters. The function $q$ may be complex, which is what allows the method to describe an absorbing medium.

When $q$ is a positive constant the solutions are $e^{\pm i\sqrt{q}\,x}$ : a wave of fixed wavenumber $\sqrt{q}$ and fixed amplitude. The local wavelength of such a wave is $2\pi/\sqrt{q}$ , and $1/\sqrt{q}$ is a wavelength divided by $2\pi$ . Riley, Hobson and Bence call $f$ slowly varying when it changes little over a range of $x$ of order $k_0^{-1}$ , that is over about one wavelength; Hassani puts the same demand on $q$ and reads it as the change in $q$ across one wavelength being small compared with $|q|$ itself. The precise inequality is stated below; it is a condition on $q'$ relative to $q^{3/2}$ , not a condition on $q'$ alone, so no function is "slowly varying" in the abstract — it is slowly varying where $|q|$ is large enough.

The WKB ansatz is to look for a solution of the form

y(x) \;=\; A\,e^{i\varphi(x)}

with $A$ constant and $\varphi$ a possibly complex phase, and to determine $\varphi$ approximately on the assumption that $\varphi''$ is small. The letters are those of Wentzel, Kramers and Brillouin, who used the method on wave-mechanical problems in 1926; Riley, Hobson and Bence note that Jeffreys had studied it earlier and Green had used it in the first half of the nineteenth century.

A point $x_0$ with

q(x_0) \;=\; 0

is a turning point. The name comes from the quantum-mechanical case $q = 2m(E - V)/\hbar^{2}$ , where $q$ vanishes exactly where a classical particle of energy $E$ would reverse direction. Turning points are where the construction below breaks down, and separate statements — the connection formulas — are needed to carry a solution past one.

Method(The WKB solutions and the exact size of their error)AdaptedAdapted from Riley, Hobson & Bence (2006), §25.7.2, Eq. (25.49), p. 898; Hassani (2013), §14.8, Eqs. (14.45)-(14.46), pp. 450-451

Let $I$ be an open interval on which $q$ is twice continuously differentiable and nowhere zero, fix $x_0 \in I$ , and fix a branch of $\sqrt{q}$ on $I$ (possible since $q$ is continuous and non-vanishing). Define the two WKB functions

y_\pm(x) \;=\; \frac{1}{q(x)^{1/4}}\,\exp\left(\pm\,i\int_{x_0}^{x}\sqrt{q(u)}\,du\right) ,

which in the notation of Riley, Hobson and Bence read $y_\pm = n^{-1/2}\exp(\pm ik_0\int n)$ , with the amplitude $f^{-1/4}$ that they emphasise is the only change from the pure phase-memory guess. Then:

They are exact solutions of a nearby equation. For every $x \in I$ , $y_\pm'' + \bigl(q - \delta\bigr)\,y_\pm \;=\; 0 , \qquad\text{where}\qquad \delta(x) \;=\; \frac{5}{16}\frac{q'(x)^{2}}{q(x)^{2}} - \frac{1}{4}\frac{q''(x)}{q(x)} .$
They are independent. The Wronskian of $y_+$ and $y_-$ is the non-zero constant $-2i$ .
Validity criterion. $y_\pm$ approximate solutions of $y'' + qy = 0$ exactly to the extent that $|\delta| \;\ll\; |q| .$ Where $q''$ is negligible beside $q'^{2}/q$ this reduces to $|q'| \ll \sqrt{16/5}\,|q|^{3/2}$ , which up to the numerical constant is Hassani's condition $|q'/(2\sqrt{q})| \ll |q|$ and, after substituting $q = k_0^{2}n^{2}$ , is exactly the condition $|k_0^{-1}n'| \ll |n|^{2}$ of Riley, Hobson and Bence.

Part 1 is an identity, not an approximation, and parts 2 and 3 follow from it; the derivation of the form of $y_\pm$ is the iteration of the two source texts, reproduced in the proof block.

Idea of the proofOriginal

Two ideas, taken in order.

The first is phase memory, which is Riley, Hobson and Bence's name for it. If $q$ were constant the phase would advance at the fixed rate $\sqrt{q}$ per unit length. If $q$ varies slowly, let the phase advance at the local rate instead and accumulate it: $\varphi = \int^{x}\sqrt{q}$ . This already captures how a wave stretches and compresses in a varying medium, and it is exact in the limit of no variation.

The second idea is that the leftover error from the first is itself computable, and absorbing it changes the amplitude rather than the phase. Substituting the pure phase-memory guess leaves a residue proportional to $q'$ ; feeding that residue back into the equation for $\varphi$ adds an imaginary piece $\tfrac{i}{4}\log q$ to the phase, and $\exp(i \cdot \tfrac{i}{4}\log q) = q^{-1/4}$ . So the correction is a real amplitude factor $q^{-1/4}$ , which is exactly what conservation of flux demands: the wave slows where $q$ is small, so it piles up there. One further round of the same iteration would change the phase again, and so on; the point of stopping after two is that the identity in part 1 of the statement then measures precisely what has been left out.

Proof(Construction of the WKB solutions, and the error identity)AdaptedAdapted from Riley, Hobson & Bence (2006), §25.7.1-25.7.2, Eqs. (25.45)-(25.49), pp. 896-898; Hassani (2013), §14.8, Eqs. (14.44)-(14.46), pp. 450-451

Step 1: the equation for the phase. Put $y = Ae^{i\varphi}$ . Then $y' = i\varphi' y$ and $y'' = (i\varphi'' - \varphi'^{2})y$ , so $y'' + qy = 0$ becomes, after dividing by $y \neq 0$ ,

\varphi'^{2} - i\varphi'' - q \;=\; 0 .

This is exact and nonlinear, and has no solution in closed form for general $q$ . The sign convention matters here and is discussed in the remark below.

Step 2: first iterate. Drop $\varphi''$ on the assumption that the phase curves slowly. Then $\varphi' = \pm\sqrt{q}$ and

\varphi(x) \;=\; \pm\int_{x_0}^{x}\sqrt{q(u)}\,du ,

which is the phase-memory solution. It is not exact: differentiating gives $\varphi'' = \pm q'/(2\sqrt{q})$ , which is small only when $|q'| \ll |q|^{3/2}$ , and that is where Hassani's validity condition comes from.

Step 3: second iterate. Put that estimate of $\varphi''$ back into the exact equation of step 1, so that $\varphi'^{2} = q + i\varphi'' \approx q \pm iq'/(2\sqrt{q})$ . Take the square root and expand the bracket, which is legitimate because the second term inside it is small by step 2:

\varphi' \;\approx\; \pm\sqrt{q}\left(1 \pm \frac{i\,q'}{2\,q^{3/2}}\right)^{1/2} \;\approx\; \pm\sqrt{q}\left(1 \pm \frac{i\,q'}{4\,q^{3/2}}\right) \;=\; \pm\sqrt{q} + \frac{i}{4}\frac{q'}{q} .

The two signs multiply to give $+$ in the last term, so the amplitude correction is the same for both solutions — a point worth noticing, since it is why both WKB functions carry the same $q^{-1/4}$ .

Step 4: integrate and exponentiate. Integrating,

\varphi(x) \;=\; \pm\int_{x_0}^{x}\sqrt{q(u)}\,du + \frac{i}{4}\log q(x) + \text{const} ,

and since $\exp\bigl(i\cdot\tfrac{i}{4}\log q\bigr) = \exp\bigl(-\tfrac14\log q\bigr) = q^{-1/4}$ , the ansatz $y = Ae^{i\varphi}$ gives the two functions $y_\pm$ of the statement. This is the route of both sources: Riley, Hobson and Bence carry a factor $k_0$ throughout and arrive at $n^{-1/2}\exp(\pm ik_0\int n)$ , Hassani works with $q$ and arrives at $q^{-1/4}\exp(\pm i\int\sqrt{q})$ , and the two agree because $n^{-1/2} = (q/k_0^{2})^{-1/4}$ differs from $q^{-1/4}$ only by the constant $k_0^{1/2}$ , which is absorbed in $A$ .

Step 5: the error identity (part 1 of the statement). Nothing above says how good $y_\pm$ is, so compute the residual exactly. Write $S(x) = \int_{x_0}^{x}\sqrt{q}$ , so $S' = \sqrt{q}$ and $S'' = q'/(2\sqrt{q})$ , and $y_\pm = q^{-1/4}e^{\pm iS}$ . Differentiating twice by the product rule,

y_\pm'' \;=\; \Bigl[\bigl(q^{-1/4}\bigr)'' \;\pm\; 2i\,\bigl(q^{-1/4}\bigr)'S' \;\pm\; i\,q^{-1/4}S'' \;-\; q^{-1/4}(S')^{2}\Bigr]e^{\pm iS} .

The last term is $-q\,y_\pm$ , which cancels $+qy_\pm$ . Of the three that remain, the two carrying $\pm i$ cancel each other: $\bigl(q^{-1/4}\bigr)' = -\tfrac14 q^{-5/4}q'$ , so $2(q^{-1/4})'S' = -\tfrac12 q^{-5/4}q'\cdot q^{1/2} = -\tfrac12 q^{-3/4}q'$ , while $q^{-1/4}S'' = q^{-1/4}q'/(2\sqrt{q}) = \tfrac12 q^{-3/4}q'$ . (This cancellation is the content of step 3: the amplitude $q^{-1/4}$ was chosen to kill the first-order residue, for either sign.) What is left is

y_\pm'' + q\,y_\pm \;=\; \bigl(q^{-1/4}\bigr)''\,q^{1/4}\,y_\pm .

Differentiating $q^{-1/4}$ twice gives $\bigl(q^{-1/4}\bigr)'' = \tfrac{5}{16}q^{-9/4}q'^{2} - \tfrac14 q^{-5/4}q''$ , so multiplying by $q^{1/4}$ ,

y_\pm'' + q\,y_\pm \;=\; \left(\frac{5}{16}\frac{q'^{2}}{q^{2}} - \frac{1}{4}\frac{q''}{q}\right)y_\pm \;=\; \delta\,y_\pm ,

which rearranges to $y_\pm'' + (q - \delta)y_\pm = 0$ . Note that $\delta$ is the same for both signs, because it came only from the amplitude.

Step 6: independence. $y_+y_-' - y_-y_+'$ is computed directly: $y_+y_- = q^{-1/2}$ and $y_\pm' = \bigl(-\tfrac14 q^{-5/4}q' \pm i q^{-1/4}\sqrt{q}\bigr)e^{\pm iS}$ , so the terms with $q'$ cancel in the difference and

W(y_+,y_-) \;=\; q^{-1/2}\bigl(-i\sqrt{q} - i\sqrt{q}\bigr) \;=\; -2i ,

a non-zero constant, so the two functions are linearly independent. (Constancy is also forced by Structure of linear equations and the Wronskian, since both solve the same second-order equation $y'' + (q-\delta)y = 0$ , which has no first-derivative term.)

Step 7: the criterion. By step 5, $y_\pm$ solves the intended equation with $q$ replaced by $q - \delta$ . Hence the fractional error made in the coefficient of the equation is exactly $|\delta/q|$ , and $y_\pm$ is a good approximation precisely when that is small. Expanding $\delta$ , the first term dominates unless $q''$ is anomalously large, giving $\tfrac{5}{16}q'^{2}/q^{2} \ll |q|$ , that is $|q'| \ll \sqrt{16/5}\,|q|^{3/2}$ . Substituting $q = k_0^{2}n^{2}$ , so $q' = 2k_0^{2}nn'$ and $|q|^{3/2} = k_0^{3}|n|^{3}$ , this becomes $2k_0^{2}|nn'| \ll \sqrt{16/5}\,k_0^{3}|n|^{3}$ , that is $|n'| \ll \sqrt{4/5}\,k_0|n|^{2}$ , which is the criterion $|k_0^{-1}n'| \ll |n|^{2}$ quoted by Riley, Hobson and Bence up to the constant. $\square$

Remark(Where the hypotheses are used, and two differences between the sources)Original

$q \ne 0$ is used twice, and it is the real restriction. The formula contains $q^{-1/4}$ , which is undefined at a zero of $q$ , and it contains $\sqrt{q}$ , which is not differentiable there. Worse, the criterion itself fails: as $q \to 0$ with $q'$ bounded away from zero, $\delta \sim \tfrac{5}{16}q'^{2}/q^{2}$ blows up while $q$ shrinks, so $|\delta/q| \to \infty$ . The WKB functions are therefore not merely inaccurate near a turning point, they are unbounded where the true solution is perfectly finite. Hassani names both ways the construction can break down — too fast a variation of $q$ , or a zero of $q$ inside the region — and says that the zero is the harder of the two to live with, precisely because the usual task is to carry one solution from a region where $q$ is positive into a region where it is negative. That is what the connection formulas below are for.

Twice-differentiability is used only in the error identity, where $q''$ appears; the functions $y_\pm$ themselves need only $q \in C^{1}$ . If $q$ has an isolated corner, $q''$ contains a point mass and the identity fails at that one point, but the phase integral $\int\sqrt{q}$ is unaffected, since its integrand is still continuous. The worked example below is of exactly this kind.

Difference 1: the sign convention in the phase equation. Riley, Hobson and Bence's Eq. (25.47) reads, with their $k_0^{2}n^{2}$ for $q$ ,

-\varphi'^{2} + i\varphi'' + q \;=\; 0 ,

which is the identity derived in step 1 above. Hassani's Eq. (14.44) is printed as $\varphi'^{2} + i\varphi'' - q = 0$ , which differs in the sign of the $i\varphi''$ term. This unit follows Riley, Hobson and Bence. The difference has no effect on the answer: Hassani's own next step inserts $\varphi'' = +q'/(2\sqrt{q})$ with the sign that the Riley-Hobson-Bence convention produces, and his Eq. (14.46) is the same two-term solution as theirs. A reader deriving the iteration from his Eq. (14.44) as printed would get the amplitude $q^{+1/4}$ instead of $q^{-1/4}$ .

Difference 2: how the quantization condition is argued. Both sources reach $\int\sqrt{q}\,dx = (n+\tfrac12)\pi$ across a well, but by different arguments, and the unit follows the Riley-Hobson-Bence one, given below, because it produces every $n$ at once. Hassani's Example 14.8.1 equates the two cosine forms continued inward from the two turning points and concludes that their amplitudes are equal. Read literally that is impossible — one cosine's argument increases with $x$ while the other's decreases, so the two functions cannot agree pointwise — and the correct statement is that they agree up to sign, since an eigenfunction is determined only up to a constant multiple. Taking the two signs in turn gives $\int\sqrt{q} = (2k+\tfrac12)\pi$ and $\int\sqrt{q} = (2k+\tfrac32)\pi$ , whose union is his stated $(n+\tfrac12)\pi$ for $n = 0,1,2,\dots$

The expansion does not converge. Iterating steps 2-4 further produces more terms, and Arfken, Weber and Harris record the standing fact about the result: the WKB expansion, one of the earliest approximations of quantum mechanics still in use, is an asymptotic series. More terms help for a while and then stop helping; Riley, Hobson and Bence's §25.7.3 makes the same point, adding the practical observation that the validity criterion can always be met in principle by making $k_0$ or $|z|$ large, but that physical parameters cannot be varied at will, so this is often of no use.

Method(Connection formulas at a simple turning point)AdaptedAdapted from Hassani (2013), §14.8, Eqs. (14.47)-(14.48), p. 451; Riley, Hobson & Bence (2006), §25.7, Eqs. (25.52)-(25.53) and pp. 899-900

Let $x_0$ be a turning point with $q(x_0) = 0$ and $q > 0$ to the right of $x_0$ , $q < 0$ to the left, and suppose the validity criterion holds in regions on both sides, away from $x_0$ . Suppose further that the solution of interest decays exponentially as one moves left from $x_0$ . Then the solution is, to the left of $x_0$ ,

y(x) \;\approx\; \frac{C}{\bigl(-q(x)\bigr)^{1/4}}\,\exp\left(-\int_{x}^{x_0}\sqrt{-q(u)}\,du\right) ,

and the same solution is, to the right of $x_0$ ,

y(x) \;\approx\; \frac{2C}{q(x)^{1/4}}\,\cos\left(\int_{x_0}^{x}\sqrt{q(u)}\,du \;-\; \frac{\pi}{4}\right) ,

with one and the same constant $C$ . These are the formulas Hassani quotes, declining to derive them; reflecting the $x$ -axis turns these into the mirror-image pair, which applies when the signs of $q$ on the two sides of the turning point are interchanged.

Two features are the whole content. The factor $2$ in the oscillatory amplitude and the phase $-\pi/4$ measured from the turning point are not conventions: they are forced, and they are what makes the quantization condition come out with a half-integer.

Riley, Hobson and Bence supply the reason, which is worth recording even though the derivation is not reproduced here. Near a simple turning point, write $q(x) \approx q'(x_0)(x - x_0)$ ; the equation becomes Stokes' equation $d^{2}y/dz^{2} = zy$ after the stretch $x - x_0 = \mu z$ with $|\mu| = |q'(x_0)|^{-1/3}$ , the sign of $\mu$ chosen so that $z > 0$ lies in the region where $q < 0$ (in the quantum case $|\mu| = [\hbar^{2}/2m|V'(x_0)|]^{1/3}$ ). Its decaying solution is the Airy function $\mathrm{Ai}(z)$ , which is finite at $z = 0$ where the WKB functions are not. The two WKB functions $z^{-1/4}\exp(\mp\tfrac23 z^{3/2})$ are each asymptotic to $\mathrm{Ai}$ in one sector of the complex $z$ -plane only, and the combination of them that equals $\mathrm{Ai}$ changes across radial lines — the Stokes lines, at $\arg z = 0, 2\pi/3, 4\pi/3$ for Stokes' equation. The factor $2$ and the phase $-\pi/4$ are the recorded values of that change: they are read off the two known asymptotic forms of $\mathrm{Ai}(z)$ , decaying like $\tfrac12\pi^{-1/2}z^{-1/4}e^{-2z^{3/2}/3}$ for $z$ real and positive and oscillating like $\pi^{-1/2}(-z)^{-1/4}\sin\bigl(\tfrac23(-z)^{3/2} + \tfrac{\pi}{4}\bigr)$ for $z$ real and negative. The general phenomenon — that the constants needed to build a solution change as the argument crosses certain lines in the complex plane — is the Stokes phenomenon.

Method(The quantization condition for a well)AdaptedAdapted from Riley, Hobson & Bence (2006), §25.7, pp. 900-901; Hassani (2013), §14.8, Example 14.8.1, pp. 451-452

Let $q$ be positive on an interval $(a,b)$ and negative outside it, with simple turning points at $a$ and $b$ , and let the validity criterion hold inside the well away from the two ends. A solution that decays exponentially in both outer regions exists only if

\int_{a}^{b}\sqrt{q(x)}\,dx \;=\; \left(n + \tfrac12\right)\pi , \qquad n = 0, 1, 2, \dots

For the one-dimensional Schrödinger equation $-\tfrac{\hbar^{2}}{2m}\psi'' + V\psi = E\psi$ , where $q = 2m(E-V)/\hbar^{2}$ and $a, b$ are the classical turning points $V(a) = V(b) = E$ , this reads

\int_{a}^{b}\sqrt{2m\bigl(E - V(x)\bigr)}\,dx \;=\; \left(n + \tfrac12\right)\pi\hbar ,

the Bohr-Sommerfeld quantization condition of pre-1925 quantum mechanics, in Hassani's words; it is the condition that the classical action across the well be a half-integer number of units of $\pi\hbar$ .

Why. Apply the connection formulas at each end. At $b$ the decaying outer solution continues inward as $\tfrac{2C}{q^{1/4}}\cos\bigl(\int_{x}^{b}\sqrt{q} - \tfrac\pi4\bigr)$ , and at $a$ as $\tfrac{2D}{q^{1/4}}\cos\bigl(\int_{a}^{x}\sqrt{q} - \tfrac\pi4\bigr)$ . These must be the same solution of a second-order equation, so the two expressions and their derivatives must agree at any interior point. Riley, Hobson and Bence write the two matching equations and observe that requiring them to give a consistent value of the amplitude ratio forces

\sin\left[\left(\int_{a}^{x}\sqrt{q} - \frac\pi4\right) + \left(\int_{x}^{b}\sqrt{q} - \frac\pi4\right)\right] \;=\; 0 ,

in which the interior point $x$ has dropped out because the two integrals add to $\int_a^b\sqrt{q}$ . Hence $\int_a^b\sqrt{q} - \tfrac\pi2 = m\pi$ for an integer $m$ , and since $\sqrt q > 0$ the left side is positive, giving $\int_a^b\sqrt q = (n+\tfrac12)\pi$ with $n = m \ge 0$ . Their calculation also notes that the terms involving $q'$ in the derivative matching cancel on their own, so this particular conclusion does not even need the validity criterion a second time.

Method(Using WKB on a second-order equation)AdaptedAdapted from Riley, Hobson & Bence (2006), §25.7, pp. 895-901; Hassani (2013), §14.8, pp. 450-452

Put the equation in the form $y'' + q(x)y = 0$ . Divide out the coefficient of $y''$ . If the result $y'' + p(x)y' + r(x)y = 0$ still has a first-derivative term, remove it by $y = w\exp\bigl(-\tfrac12\int p\bigr)$ ; then $w'' + qw = 0$ with $q = r - \tfrac12 p' - \tfrac14 p^{2}$ .
Locate the turning points, the zeros of $q$ , and separate the axis into regions where $q > 0$ (oscillatory) and $q < 0$ (exponential).
Test the criterion region by region. Compute $\delta = \tfrac{5}{16}q'^{2}/q^{2} - \tfrac14 q''/q$ and check $|\delta| \ll |q|$ . Record where it fails; it always fails in a neighbourhood of every turning point, and the width of that neighbourhood is the honest measure of how much of the problem WKB cannot see.
Write the solutions. In a region with $q > 0$ use $\frac{C_\pm}{q^{1/4}}\exp\left(\pm i\int\sqrt{q}\right) \quad\text{or}\quad \frac{1}{q^{1/4}}\left\{\cos,\ \sin\right\}\left(\int\sqrt{q}\right) ;$ in a region with $q < 0$ use $(-q)^{-1/4}\exp\bigl(\pm\int\sqrt{-q}\bigr)$ . Choose the decaying combination wherever a boundary condition at infinity requires it.
Join regions with the connection formulas, not by matching the WKB forms directly at the turning point, where both are singular. Keep the factor $2$ and the phase $-\pi/4$ .
Read off the eigenvalue condition. For a well bounded by two turning points this is the quantization condition above; the turning points themselves depend on the eigenvalue, so it is an implicit equation for $E$ .
Judge the answer by step 3, not by the formula's appearance. The approximation improves as the phase integral across the well grows — in a bound-state problem that means high quantum numbers — and is at its worst for the ground state.

Example(A triangular well: WKB levels against the exact Airy spectrum)Original

Take a particle of mass $m$ in the one-dimensional potential $V(x) = F|x|$ with $F > 0$ : a linear restoring force of constant magnitude, pointing towards the origin. This is the potential of a ball bouncing under gravity, reflected to make it symmetric, and it is a good test because its exact spectrum is known in closed form.

WKB. With $q(x) = 2m\bigl(E - F|x|\bigr)/\hbar^{2}$ the turning points are at $\pm a$ with $a = E/F$ , and

\int_{-a}^{a}\sqrt{q}\,dx \;=\; \frac{2\sqrt{2m}}{\hbar}\int_{0}^{E/F}\bigl(E - Fx\bigr)^{1/2}dx \;=\; \frac{2\sqrt{2m}}{\hbar}\cdot\frac{2E^{3/2}}{3F} \;=\; \frac{4\sqrt{2m}\,E^{3/2}}{3\hbar F} .

Setting this equal to $(n+\tfrac12)\pi$ and solving for $E$ ,

E_n \;=\; \left[\frac{3\pi\hbar F\left(n+\tfrac12\right)}{4\sqrt{2m}}\right]^{2/3} .

In units with $\hbar = m = F = 1$ this is $E_n = \bigl[3\pi(n+\tfrac12)/(4\sqrt{2})\bigr]^{2/3}$ .

The exact answer. In the same units the equation is $\psi'' + 2(E - |x|)\psi = 0$ . For $x > 0$ the substitution $u = 2^{1/3}(x - E)$ turns it into $d^{2}\psi/du^{2} = u\psi$ , so the solution decaying as $x \to +\infty$ is $\psi = \mathrm{Ai}\bigl(2^{1/3}(x-E)\bigr)$ ; see Airy functions. Because $V$ is even, the eigenfunctions are even or odd, and the matching condition at the origin is $\psi'(0) = 0$ or $\psi(0) = 0$ respectively. Hence

\text{even } n: \ \mathrm{Ai}'\bigl(-2^{1/3}E\bigr) = 0 , \qquad\qquad \text{odd } n: \ \mathrm{Ai}\bigl(-2^{1/3}E\bigr) = 0 ,

so $E = 2^{-1/3}|a_k|$ with $a_k$ the zeros of $\mathrm{Ai}'$ and of $\mathrm{Ai}$ in turn.

Comparison (units $\hbar = m = F = 1$ ):

$n$	WKB $E_n$	exact $E_n$	relative error, per cent
$0$	$0.885341$	$0.808617$	$+9.5$
$1$	$1.841584$	$1.855757$	$-0.76$
$2$	$2.588754$	$2.578096$	$+0.41$
$3$	$3.239735$	$3.244608$	$-0.15$
$5$	$4.378976$	$4.381671$	$-0.06$

The criterion explains the table. Here $q = 2(E - |x|)$ , so $q' = \mp 2$ and $q'' = 0$ away from the origin, giving $\delta = \tfrac{5}{16}\cdot\tfrac{4}{q^{2}} = \tfrac{5}{4q^{2}}$ , and the criterion $|\delta| \ll q$ becomes

q^{3} \;\gg\; \tfrac54 .

The best place in the well is the bottom, $x = 0$ , where $q = 2E$ . For $n = 0$ , $q(0) = 1.771$ and $|\delta/q| = 0.225$ — a 22 per cent error in the coefficient of the equation at the most favourable point, and the criterion fails outright ( $|\delta| > q$ ) for $|x| > 0.347$ , which is the outer 61 per cent of the well. A 9.5 per cent error in the level is a mild punishment for that. For $n = 5$ , $q(0) = 8.758$ and $|\delta/q| = 0.0019$ , with the criterion failing only in the outer 12 per cent of the well, and the level is right to six parts in ten thousand. This is the behaviour promised in step 7 of the method: WKB is an approximation for the highly excited states that becomes, almost by accident, tolerable for the ground state.

One honest caveat. $V = F|x|$ has a corner at the origin, where $q'$ jumps and $q''$ has a point mass, so the error identity does not hold at that single point. It does not matter: $\sqrt{q}$ is continuous there, so the phase integral and hence the quantization condition are unaffected, and $\delta$ is finite everywhere else in the well. A smooth potential with the same shape away from the origin gives the same table to the accuracy shown.

Application(Waves in slowly varying media, and barrier penetration)AdaptedAdapted from Riley, Hobson & Bence (2006), §25.7.1, pp. 895-897; Hassani (2013), §14.8, Example 14.8.1, pp. 451-452

Optics and radio propagation. The reason Riley, Hobson and Bence set the problem up with a refractive index $n(x)$ is that this is the method's oldest use. A wave in a medium whose properties change over distances long compared with a wavelength keeps a phase that is the accumulated local phase, $k_0\int n$ , and an amplitude that scales as $n^{-1/2}$ . The condition for this to be valid — that $n$ change little over a few wavelengths — is comfortably satisfied by visible light in glass and by radio waves in the ionosphere, which is why ray tracing works at all; the places where it fails are precisely the interesting ones, such as the height at which a radio wave is reflected, where the effective $q$ passes through zero and the turning-point analysis above is needed. Allowing $n$ to be complex, which costs nothing in the derivation, describes a wave in an absorbing medium: the imaginary part of $n$ turns the oscillation into a decaying one, and $\mathrm{Im}\int n$ is the accumulated attenuation.

Quantum spectra. The quantization condition is how energy levels were estimated before the Schrödinger equation could be solved, and it is still the fastest way to get the level density of a complicated potential. Hassani's reading of it is the physical one: the integral is the classical action accumulated across the classically allowed region, and the half-integer on the right is the accumulated $\pi/4$ of phase from each of the two turning points. Two checks are standard. For a harmonic oscillator the condition is exact, returning $E_n = (n+\tfrac12)\hbar\omega$ including the zero-point energy; and for a well $V = V_0|x|^{2s}$ Riley, Hobson and Bence obtain $E_n \propto (n+\tfrac12)^{2s/(s+1)}$ , which reproduces the equal spacing of the oscillator at $s = 1$ and the $n^{2}$ spacing of a square well in the limit of large $s$ .

Barrier penetration. The exponential branch of the construction is the other standard use. Where $q < 0$ the WKB solutions are $(-q)^{-1/4}\exp\bigl(\pm\int\sqrt{-q}\bigr)$ , so a wave crossing a classically forbidden region between two turning points is attenuated by a factor whose logarithm is $\int\sqrt{-q}\,dx$ across the barrier. Squaring to get a probability gives a transmission that is exponentially small in $2\int\sqrt{-q}$ . This is the estimate behind the Gamow factor for alpha decay, behind field emission of electrons from a metal surface, and behind the exponential dependence of a scanning tunnelling microscope's current on tip height; in each case the enormous dynamic range of the measured quantity comes from a quantity sitting in an exponent, and WKB is accurate enough for the exponent even where it would be poor for a prefactor.