Linear equations

Structure of solutions of linear ODEs and the standard solution methods.

In this chapter 9

Tick a unit to add it to your PDF basket.

Linear equations of order n
Series solutions

Linear equations of order n

Theorem

Structure of linear equations and the Wronskian

Adapted

Add to PDF

For a linear equation of order n whose leading coefficient does not vanish, the solutions of the homogeneous equation form a vector space of dimension exactly n, any n solutions are independent precisely when their Wronskian is nonzero at one point, and Abel's formula gives that Wronskian up to a constant without solving the equation.

Prerequisites: Existence and uniqueness of solutions, Determinants

Definition(Linear equations of order

n

)AdaptedAdapted from Coddington (1989), Ch. 3, §1, p. 103; Stakgold & Holst (2011), §3.1, p. 186

A linear differential equation of order $n$ on a real interval $I$ is an equation

a_0(x)\,y^{(n)} + a_1(x)\,y^{(n-1)} + \cdots + a_n(x)\,y \;=\; b(x) ,

with coefficients $a_0, \dots, a_n$ and right-hand side $b$ complex-valued functions on $I$ . A point where $a_0$ vanishes is a singular point of the equation and needs separate treatment (Classifying singular points, including infinity). Throughout this unit $a_0(x) \neq 0$ on $I$ , so we may divide by $a_0$ and work with the normalised equation. Writing

L(y) \;=\; y^{(n)} + a_1(x)\,y^{(n-1)} + \cdots + a_n(x)\,y

(the coefficients renamed after the division), the equation is $L(y) = b(x)$ . It is homogeneous if $b(x) = 0$ for all $x \in I$ and inhomogeneous otherwise. A solution on $I$ is a function $\varphi$ that has $n$ derivatives at every point of $I$ and satisfies $L(\varphi) = b$ there. We assume once and for all that $a_1, \dots, a_n$ and $b$ are continuous on $I$ ; $I$ may be open, closed, half-open, bounded or unbounded.

$L$ is a linear operator: $L(c_1\varphi_1 + c_2\varphi_2) = c_1 L(\varphi_1) + c_2 L(\varphi_2)$ for constants $c_1, c_2$ and any functions with $n$ derivatives. This single fact is the principle of superposition, and everything below is a consequence of it together with unique solvability of the initial-value problem.

Definition(Linear independence and the Wronskian)AdaptedAdapted from Stakgold & Holst (2011), §3.1, pp. 186-187; Coddington (1989), Ch. 3, §3-§4, pp. 106, 111; Arfken, Weber & Harris (2013), §7.6, Eqs. (7.54)-(7.55), pp. 359-360

Functions $f_1, \dots, f_n$ on $I$ are linearly dependent on $I$ when some constants $c_1, \dots, c_n$ , at least one of them nonzero, satisfy

c_1 f_1(x) + c_2 f_2(x) + \cdots + c_n f_n(x) \;=\; 0 \qquad \text{for every } x \in I ,

and linearly independent on $I$ if they are not. The constants are the same for all $x$ : one relation must hold identically on $I$ , and a relation that happens to hold at isolated points or on part of $I$ does not make the functions dependent.

If $f_1, \dots, f_n$ are continuous on $I$ and have continuous derivatives of orders up to $n-1$ there — the hypothesis under which Stakgold and Holst (Stakgold & Holst (2011), p. 187) define it — their Wronskian is the function

W(f_1, \dots, f_n)(x) \;=\; \det \begin{pmatrix} f_1 & f_2 & \cdots & f_n \\ f_1' & f_2' & \cdots & f_n' \\ \vdots & \vdots & & \vdots \\ f_1^{(n-1)} & f_2^{(n-1)} & \cdots & f_n^{(n-1)} \end{pmatrix} (x) .

For $n = 2$ this is $W(f_1, f_2) = f_1 f_2' - f_1' f_2$ . Solutions of the normalised equation meet this requirement automatically: from $y^{(n)} = -a_1y^{(n-1)} - \cdots - a_ny$ and the continuity of the $a_k$ , a solution has $n$ derivatives of which the first $n-1$ are continuous (indeed the $n$ -th is continuous too). A set of $n$ solutions of $L(y) = 0$ that is linearly independent on $I$ is called a fundamental set, or a basis of solutions.

Theorem(Superposition and the dimension of the solution space)AdaptedAdapted from Coddington (1989), Ch. 3, §3, Thms 4 and 5, pp. 106-108, and §4, Thm 7, pp. 112-113

Let $a_1, \dots, a_n$ be continuous on $I$ and let $\mathcal{S}$ be the set of solutions of $L(y) = 0$ on $I$ .

Superposition. $\mathcal{S}$ is a complex vector space: if $\varphi_1, \dots, \varphi_m \in \mathcal{S}$ and $c_1, \dots, c_m$ are constants, then $c_1\varphi_1 + \cdots + c_m\varphi_m \in \mathcal{S}$ .
Initial data. Fix $x_0 \in I$ . For every choice of constants $\alpha_1, \dots, \alpha_n$ there is exactly one $\varphi \in \mathcal{S}$ with
$\varphi(x_0) = \alpha_1, \quad \varphi'(x_0) = \alpha_2, \quad \dots, \quad \varphi^{(n-1)}(x_0) = \alpha_n ,$
and this solution exists on the whole of $I$ .
Dimension $n$ . There exist $n$ linearly independent solutions on $I$ ; and if $\varphi_1, \dots, \varphi_n$ is any set of $n$ linearly independent solutions on $I$ , then every solution $\varphi$ of $L(y) = 0$ on $I$ can be written as
$\varphi \;=\; c_1\varphi_1 + \cdots + c_n\varphi_n$
with uniquely determined constants $c_1, \dots, c_n$ . Thus $\dim \mathcal{S} = n$ , and every fundamental set is a basis of $\mathcal{S}$ .

Theorem(The Wronskian criterion)AdaptedAdapted from Coddington (1989), Ch. 3, §4, Thm 6, pp. 111-112 and the consequence of Thm 8, p. 116; Stakgold & Holst (2011), §3.1, Thm 3.1.2, pp. 187-188

Let $\varphi_1, \dots, \varphi_n$ be $n$ solutions of $L(y) = 0$ on $I$ , with $a_1, \dots, a_n$ continuous on $I$ , and let $x_0$ be any point of $I$ . The following are equivalent.

$\varphi_1, \dots, \varphi_n$ are linearly independent on $I$ .
$W(\varphi_1, \dots, \varphi_n)(x) \neq 0$ for every $x \in I$ .
$W(\varphi_1, \dots, \varphi_n)(x_0) \neq 0$ .

Equivalently: $n$ solutions are linearly dependent on $I$ if and only if their Wronskian vanishes at a single point of $I$ . In particular the Wronskian of $n$ solutions either never vanishes on $I$ or vanishes identically, and it may be tested at whichever point is most convenient.

Theorem(Abel's formula)AdaptedAdapted from Coddington (1989), Ch. 3, §4, Thm 8 and its corollary, pp. 113-116; Stakgold & Holst (2011), §3.1, Eqs. (3.1.7)-(3.1.8), p. 187

Let $\varphi_1, \dots, \varphi_n$ be any $n$ solutions (independent or not) of $L(y) = 0$ on $I$ and let $x_0 \in I$ . Then, with $W = W(\varphi_1, \dots, \varphi_n)$ ,

W(x) \;=\; \exp\left[-\int_{x_0}^{x} a_1(t)\,dt\right] W(x_0) , \qquad x \in I .

Only the coefficient $a_1$ of $y^{(n-1)}$ enters. Three consequences are worth recording.

If the $a_k$ are constants, $W(x) = e^{-a_1(x - x_0)}\,W(x_0)$ ; if $a_1 = 0$ , $W$ is constant on $I$ .
For the unnormalised equation $a_0 y^{(n)} + a_1 y^{(n-1)} + \cdots + a_n y = 0$ with $a_0 \ne 0$ on $I$ , the same formula reads $W = C e^{-m(x)}$ , where $m$ is any particular solution of $m' = a_1/a_0$ and the constant $C$ depends on the solutions used to form $W$ .
For the second-order equation in self-adjoint form, $(p\,y')' + q\,y = 0$ with $p \ne 0$ (so that $a_0 = p$ and $a_1 = p'$ ), this gives $W(x) = C/p(x)$ .

Idea of the proofOriginal

Everything follows from one observation: reading off the initial data at a fixed point $x_0$ is a linear map

E_{x_0} : \mathcal{S} \to \mathbb{C}^n , \qquad E_{x_0}(\varphi) = \bigl(\varphi(x_0), \varphi'(x_0), \dots, \varphi^{(n-1)}(x_0)\bigr) .

The existence theorem says $E_{x_0}$ is onto; the uniqueness theorem says its kernel is $\{0\}$ . So $E_{x_0}$ is an isomorphism and $\dim\mathcal{S} = n$ . The Wronskian at $x_0$ is precisely the determinant of the $n$ vectors $E_{x_0}(\varphi_1), \dots, E_{x_0}(\varphi_n)$ , so the Wronskian criterion is nothing but the statement that $n$ vectors in $\mathbb{C}^n$ are independent exactly when their determinant is nonzero — transported through the isomorphism. Abel's formula is separate: differentiating a determinant gives a sum of $n$ determinants, all but the last have a repeated row, and in the last the equation lets us replace the bottom row $\varphi_k^{(n)}$ by $-a_1\varphi_k^{(n-1)} - \cdots - a_n \varphi_k$ , of which only the $-a_1$ term survives row operations. So $W' = -a_1 W$ , a first-order equation.

Proof(Superposition and the dimension of the solution space)AdaptedAdapted from Coddington (1989), Ch. 3, Thms 4, 5 and 7, pp. 106-113

Step 1 (part 1). For functions with $n$ derivatives, $L(c_1\varphi_1 + \cdots + c_m\varphi_m) = c_1 L(\varphi_1) + \cdots + c_m L(\varphi_m)$ , because each derivative is linear and the coefficients multiply pointwise. If every $L(\varphi_j) = 0$ , the right side is $0$ . Hence $\mathcal{S}$ is closed under linear combinations, and it contains the zero function, so it is a vector space.

Step 2 (part 2). Existence of a solution on all of $I$ with prescribed data at $x_0$ , and its uniqueness, is the existence and uniqueness theorem applied to the first-order system for $(y, y', \dots, y^{(n-1)})$ ; it uses the continuity of $a_1, \dots, a_n$ on $I$ , which is available because we divided by $a_0 \neq 0$ . (Uniqueness alone is elementary: if $\varphi, \psi$ are two solutions with the same data, $\chi = \varphi - \psi$ solves $L(y) = 0$ with zero data at $x_0$ , and a Gronwall-type estimate on each closed bounded subinterval forces $\chi \equiv 0$ .)

Step 3 ( $E_{x_0}$ is a linear isomorphism). The map $E_{x_0}$ of the idea is linear, since differentiation and evaluation are. Part 2 says exactly that it is a bijection of $\mathcal{S}$ onto $\mathbb{C}^n$ . Hence $\dim\mathcal{S} = \dim\mathbb{C}^n = n$ .

Step 4 (existence of $n$ independent solutions). Let $e_1, \dots, e_n$ be the standard basis of $\mathbb{C}^n$ and let $\psi_k = E_{x_0}^{-1}(e_k)$ , that is, the solution with $\psi_k^{(k-1)}(x_0) = 1$ and all other derivatives of order $< n$ vanishing at $x_0$ . If $c_1\psi_1 + \cdots + c_n\psi_n = 0$ on $I$ , applying the linear map $E_{x_0}$ gives $c_1e_1 + \cdots + c_ne_n = 0$ in $\mathbb{C}^n$ , so every $c_k = 0$ . These $n$ solutions are therefore linearly independent.

Step 5 (any $n$ independent solutions are a basis). Let $\varphi_1, \dots, \varphi_n$ be solutions, linearly independent on $I$ . Linear independence in $\mathcal{S}$ means precisely linear independence as functions on $I$ — the defining relation $\sum_k c_k\varphi_k = 0$ is the same in both readings. By Step 3, $\mathcal{S}$ has dimension $n$ , and $n$ independent vectors in an $n$ -dimensional space span it. So every solution is a linear combination $\sum_k c_k\varphi_k$ , and the coefficients are unique because the $\varphi_k$ are independent. $\square$

Proof(The Wronskian criterion)AdaptedAdapted from Coddington (1989), Ch. 3, Thm 6, pp. 111-112; Stakgold & Holst (2011), Thm 3.1.2, pp. 187-188

Fix $x \in I$ and let $M(x)$ be the matrix in the definition of the Wronskian, built from $\varphi_1, \dots, \varphi_n$ ; its $k$ -th column is $E_x(\varphi_k)$ , so $W(x) = \det M(x)$ .

Step 1 ( $3 \Rightarrow 1$ ). Suppose $W(x_0) \ne 0$ and $c_1\varphi_1 + \cdots + c_n\varphi_n = 0$ on $I$ . Differentiating this identity $n-1$ times gives

c_1\varphi_1^{(j)}(x) + \cdots + c_n\varphi_n^{(j)}(x) = 0, \qquad j = 0, 1, \dots, n-1,

for every $x \in I$ . At $x = x_0$ this is the homogeneous linear system $M(x_0)\,c = 0$ for $c = (c_1, \dots, c_n)$ . Since $\det M(x_0) = W(x_0) \ne 0$ , the determinant criterion for invertibility gives $c = 0$ . So the solutions are linearly independent on $I$ .

Step 2 ( $1 \Rightarrow 2$ ), by contraposition. Suppose $W(x_1) = 0$ for some $x_1 \in I$ . Then $M(x_1)$ is singular, so there is a nonzero $c \in \mathbb{C}^n$ with $M(x_1)c = 0$ . Put $\psi = c_1\varphi_1 + \cdots + c_n\varphi_n$ . By superposition $\psi$ solves $L(y) = 0$ , and $M(x_1)c = 0$ says that $\psi$ and its first $n-1$ derivatives vanish at $x_1$ . The zero function is also such a solution, so by the uniqueness in part 2 of the structure theorem, $\psi \equiv 0$ on $I$ . Since not all $c_k$ vanish, $\varphi_1, \dots, \varphi_n$ are linearly dependent on $I$ . Hence independence forces $W$ to be nonvanishing at every point of $I$ .

Step 3 (the three statements are equivalent). Step 1 is $3 \Rightarrow 1$ , Step 2 is $1 \Rightarrow 2$ , and $2 \Rightarrow 3$ is immediate because $x_0 \in I$ . $\square$

The point of the equivalence of 2 and 3 is practical: the Wronskian of solutions cannot vanish at some points of $I$ and not at others, so one evaluation suffices. Abel's formula below shows the same thing directly, since an exponential never vanishes.

Proof(Abel's formula)AdaptedAdapted from Coddington (1989), Ch. 3, Thm 8, pp. 113-115

Step 1 (differentiate the determinant). A determinant is a sum of products of one entry from each row, so by the product rule its derivative is the sum of the $n$ determinants obtained by differentiating one row at a time:

W' \;=\; V_1 + V_2 + \cdots + V_n ,

where $V_k$ agrees with $W$ except that its $k$ -th row is differentiated.

Step 2 (all but the last term vanish). For $k \le n-1$ , differentiating row $k$ of $W$ (whose entries are $\varphi_j^{(k-1)}$ ) produces the row with entries $\varphi_j^{(k)}$ , which is already row $k+1$ of $W$ . So $V_k$ has two equal rows and $V_k = 0$ . Only $V_n$ survives, whose bottom row has entries $\varphi_j^{(n)}$ .

Step 3 (use the equation). Here we use that every $\varphi_j$ solves $L(y) = 0$ :

\varphi_j^{(n)} \;=\; -a_1\varphi_j^{(n-1)} - a_2\varphi_j^{(n-2)} - \cdots - a_n\varphi_j .

Substituting this bottom row into $V_n$ , and then adding to the bottom row $a_n$ times the first row, $a_{n-1}$ times the second row, ..., $a_2$ times the $(n-1)$ -st row — operations that do not change a determinant — cancels every term except the first, leaving the bottom row $(-a_1\varphi_1^{(n-1)}, \dots, -a_1\varphi_n^{(n-1)})$ . That row is $-a_1$ times row $n$ of $W$ , so

W' \;=\; V_n \;=\; -a_1(x)\,W .

Step 4 (solve the first-order equation). $a_1$ is continuous on $I$ , so $A(x) = \int_{x_0}^{x} a_1(t)\,dt$ is defined and differentiable on $I$ . Then $\bigl(e^{A}W\bigr)' = e^{A}(W' + a_1W) = 0$ , so $e^{A}W$ is constant on $I$ ; evaluating at $x = x_0$ , where $A(x_0) = 0$ , that constant is $W(x_0)$ . Hence $W(x) = e^{-A(x)}W(x_0)$ , which is the formula. $\square$

For $n = 2$ Step 3 is a two-line computation: $W = \varphi_1\varphi_2' - \varphi_1'\varphi_2$ gives $W' = \varphi_1\varphi_2'' - \varphi_1''\varphi_2$ , and substituting $\varphi_j'' = -a_1\varphi_j' - a_2\varphi_j$ makes the $a_2$ terms cancel, leaving $W' = -a_1W$ .

Theorem(Inhomogeneous equations)AdaptedAdapted from Coddington (1989), Ch. 3, §6, p. 122 and Thm 11, p. 123; Riley, Hobson & Bence (2006), §15.1, Eq. (15.7), p. 492

Let $b$ be continuous on $I$ , let $\varphi_p$ be one solution of $L(y) = b$ on $I$ , and let $\varphi_1, \dots, \varphi_n$ be a fundamental set for $L(y) = 0$ on $I$ . Then the solutions of $L(y) = b$ on $I$ are exactly the functions

\varphi \;=\; \varphi_p + c_1\varphi_1 + \cdots + c_n\varphi_n ,

$c_1, \dots, c_n$ constants, and the constants are uniquely determined by $\varphi$ . In the physicists' names: the general solution is a particular integral plus the complementary function, the latter carrying the $n$ arbitrary constants that the initial or boundary conditions fix.

ProofAdaptedAdapted from Coddington (1989), Ch. 3, §6, p. 122

If $\varphi$ solves $L(y) = b$ , then by linearity $L(\varphi - \varphi_p) = b - b = 0$ , so $\varphi - \varphi_p$ is a solution of the homogeneous equation and the structure theorem writes it as $\sum_k c_k\varphi_k$ with unique constants. Conversely, if $\varphi = \varphi_p + \sum_k c_k\varphi_k$ then $L(\varphi) = L(\varphi_p) + \sum_k c_k L(\varphi_k) = b + 0 = b$ . $\square$

Remark(Where the hypotheses are used)Original

$a_0 \ne 0$ on $I$ is essential, not cosmetic. Take $x^2y'' - 2xy' + 2y = 0$ on $I = \mathbb{R}$ , whose leading coefficient vanishes at $x = 0$ . Both $\varphi_1(x) = x$ and $\varphi_2(x) = x^2$ are solutions on $\mathbb{R}$ , and they are linearly independent there, yet
$W(\varphi_1, \varphi_2)(x) = x\cdot 2x - 1 \cdot x^2 = x^2$
vanishes at $x = 0$ : statements 2 and 3 of the Wronskian criterion fail. Part 2 of the structure theorem fails in both directions at the same point. Existence fails because setting $x = 0$ in the equation leaves $2y(0) = 0$ , so every solution on an interval containing $0$ — not merely those of the form $c_1x + c_2x^2$ — has $\varphi(0) = 0$ , and the initial-value problem with $\varphi(0) = 1$ has no solution. Uniqueness fails because the data $\varphi(0) = \varphi'(0) = 0$ are satisfied by $c_2x^2$ for every $c_2$ . Away from the singular point, say on $I = (0, \infty)$ , the theory applies in full: there $W = x^2 \neq 0$ and $\{x, x^2\}$ is a fundamental set (see Euler-Cauchy equations).
Continuity of the $a_k$ on the whole of $I$ is what makes Step 2 give a solution on all of $I$ rather than only near $x_0$ . This is special to linear equations: nonlinear equations can lose their solutions in finite time even with smooth data.
The criterion is about solutions of one equation. For functions that are not solutions of a common linear equation, a Wronskian vanishing identically does not imply dependence. On $I = (-1, 1)$ put $f_1(x) = x^2$ and $f_2(x) = x|x|$ . Both are differentiable, $f_2' = 2|x|$ , so
$W(f_1, f_2) = x^2 \cdot 2|x| - 2x \cdot x|x| = 0$
identically, yet $f_1, f_2$ are independent on $I$ : $c_1f_1 + c_2f_2 = 0$ at $x = \tfrac12$ and at $x = -\tfrac12$ gives $c_1 + c_2 = 0$ and $c_1 - c_2 = 0$ . The two functions cannot both be solutions of a normalised second-order equation on $I$ , and indeed the equation they do satisfy is the singular one above, $x^2y'' - 2xy' + 2y = 0$ — and only for $x \neq 0$ , since $f_2 = x|x|$ has no second derivative at the origin. This example is Stakgold and Holst's (Stakgold & Holst (2011), p. 187).
How the books phrase that warning. We follow (Coddington (1989), Ch. 3, Thm 6, pp. 111-112) and (Stakgold & Holst (2011), Thm 3.1.2, pp. 187-188) in stating the equivalence only for $n$ solutions of $L(y) = 0$ . Riley, Hobson and Bence (Riley, Hobson & Bence (2006), p. 491) give the same direction we use (a nonvanishing Wronskian implies independence) and warn that a Wronskian which vanishes is not by itself evidence of dependence. Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §7.6, pp. 359-360) do assert that functions whose Wronskian is zero over the entire range are linearly dependent, but their footnote to that assertion assumes continuous derivatives and that at least one minor of the bottom row does not vanish on the interval; the $x^2, x|x|$ pair fails that extra condition. Hassani (Hassani (2013), Thm 14.4.3, p. 425) states the equivalence for two differentiable functions that are nonzero in the interval, a hypothesis the pair also fails.
Abel's formula needs no independence. It holds for any $n$ solutions. If they are dependent, $W(x_0) = 0$ and the formula correctly returns $W \equiv 0$ .

Example(A fundamental set for Legendre's equation with

\ell = 1

)Original

Take

(1 - x^2)\,y'' - 2x\,y' + 2y = 0 \qquad \text{on } I = (-1, 1) .

The leading coefficient $1 - x^2$ is nonzero on $I$ , and dividing by it gives the normalised coefficients

a_1(x) = \frac{-2x}{1 - x^2}, \qquad a_2(x) = \frac{2}{1 - x^2} ,

both continuous on $I$ . So the solution space on $I$ has dimension $2$ .

Two solutions. By inspection $\varphi_1(x) = x$ solves the equation: $y'' = 0$ and $-2x + 2x = 0$ . A second solution, obtained from the first by reduction of order, is
$\varphi_2(x) = \frac{x}{2}\,\ln\frac{1 + x}{1 - x} - 1 .$
Check it directly: with $L = \ln\frac{1+x}{1-x}$ one has $L' = \frac{2}{1 - x^2}$ , so $\varphi_2' = \tfrac12 L + \frac{x}{1 - x^2}$ and $\varphi_2'' = \frac{1}{1-x^2} + \frac{1 + x^2}{(1 - x^2)^2}$ , and substituting gives $0$ .
Wronskian, computed directly.
$W = \varphi_1\varphi_2' - \varphi_1'\varphi_2 = x\left(\tfrac12 L + \frac{x}{1 - x^2}\right) - \left(\tfrac{x}{2}L - 1\right) = \frac{x^2}{1 - x^2} + 1 = \frac{1}{1 - x^2} .$
The $L$ terms cancel, as they must: the Wronskian of two solutions is far simpler than either solution.
Wronskian, from Abel's formula. With $x_0 = 0$ ,
$-\int_0^x a_1(t)\,dt = \int_0^x \frac{2t}{1 - t^2}\,dt = -\ln(1 - x^2) ,$
so $W(x) = W(0)/(1 - x^2)$ . Since $\varphi_1(0) = 0$ , $\varphi_1'(0) = 1$ , $\varphi_2(0) = -1$ , $\varphi_2'(0) = 0$ , we get $W(0) = 0\cdot 0 - 1 \cdot(-1) = 1$ and $W(x) = 1/(1 - x^2)$ , agreeing with step 2. This is also the self-adjoint case of Abel's formula: the equation is $\bigl((1 - x^2)y'\bigr)' + 2y = 0$ , so $W = C/p$ with $p = 1 - x^2$ and $C = 1$ .
Fundamental set. $W(0) = 1 \neq 0$ , so by the Wronskian criterion $\{\varphi_1, \varphi_2\}$ is a fundamental set on $(-1,1)$ , and $W$ is nonvanishing on all of $(-1,1)$ — as the closed form confirms.
Solving an initial-value problem. For $\varphi(0) = 2$ , $\varphi'(0) = 3$ , write $\varphi = c_1\varphi_1 + c_2\varphi_2$ . The data give $-c_2 = 2$ and $c_1 = 3$ , so
$\varphi(x) = 3x - 2\left(\frac{x}{2}\ln\frac{1+x}{1-x} - 1\right) = 3x + 2 - x\,\ln\frac{1+x}{1-x} .$
At $x = \tfrac12$ this is $3.5 - \tfrac12\ln 3 = 2.950\,69$ to six significant figures ( $2.950\,693\,856$ to ten).
What breaks at $x = \pm 1$ . The normalised coefficients blow up there, the solution space on an interval containing $x = 1$ is not described by the theorem, and indeed $\varphi_2$ is unbounded as $x \to 1^-$ while $\varphi_1$ is not. The Wronskian $1/(1 - x^2)$ diverges at the same places, which is Abel's formula telling us where the trouble is without our having to find the solutions.

ApplicationOriginal

Three uses, in increasing order of depth.

Counting the constants. A physicist writing down a linear model of order $n$ — a beam equation, an RLC circuit, a linearised normal-mode problem — knows before doing any work that the general solution carries exactly $n$ arbitrary constants, and therefore that exactly $n$ conditions (initial values, or boundary values at two ends) are needed and that more than $n$ will generally be inconsistent. The dimension theorem is what licenses the superposition of normal modes: because the solution space is finite-dimensional with a basis, decomposing an arbitrary motion into modes is a change of basis, not an approximation.

The Wronskian in the denominator. The Wronskian of a fundamental set is the determinant that must be inverted when a particular solution is built from the homogeneous ones, so it appears in the denominator of variation of parameters and in the normalisation of the one-dimensional Green's function. Abel's formula makes that denominator explicit before the solutions are known: for the self-adjoint equation $(p\,y')' + q\,y = 0$ , $W = C/p$ , which is why Green's functions for Sturm-Liouville operators carry a factor $1/(p\,W)$ that is simply a constant.

Constancy of the Wronskian in quantum mechanics. The time-independent Schrodinger equation $-\frac{\hbar^2}{2m}\psi'' + V(x)\psi = E\psi$ has no first-derivative term, so $a_1 = 0$ and the Wronskian of any two solutions is constant in $x$ . Several standard facts are this statement in disguise: the probability current of a stationary state is $x$ -independent; the transmission and reflection amplitudes of a one-dimensional barrier satisfy $|t|^2 + |r|^2 = 1$ , obtained by equating the constant Wronskian of the incoming and outgoing asymptotic forms on the two sides of the barrier; and a bound state, which decays at both ends, forces the Wronskian to vanish, so the decaying solutions at the two ends must be proportional — which is exactly the quantisation condition that selects the allowed energies.

Method

Constant-coefficient equations

Adapted

Add to PDF

One identity, $L(e^{\lambda x}q) = e^{\lambda x}\sum_j p^{(j)}(\lambda)q^{(j)}/j!$ , delivers the whole method — a fundamental set $x^{r}e^{\lambda_j x}$ read off the roots of the characteristic polynomial, the real form for a complex pair, and a polynomial-times-exponential particular solution whose degree rises by the multiplicity of the forcing exponent, which is resonance.

Prerequisites: Structure of linear equations and the Wronskian

Definition(The characteristic polynomial)AdaptedAdapted from Coddington (1989), Ch. 2, §7, p. 71; Hassani (2013), §14.7, Eqs. (14.37)-(14.39), p. 446; Riley, Hobson & Bence (2006), §15.1, Eqs. (15.8)-(15.10), p. 492

Let $a_1, \dots, a_n$ be constants (real or complex) and

L(y) \;=\; y^{(n)} + a_1 y^{(n-1)} + a_2 y^{(n-2)} + \cdots + a_n y .

Its characteristic polynomial is

p(\lambda) \;=\; \lambda^{n} + a_1\lambda^{n-1} + a_2\lambda^{n-2} + \cdots + a_n ,

so that $L$ is obtained from $p$ by replacing $\lambda^{k}$ with $d^{k}/dx^{k}$ . The equation $p(\lambda) = 0$ is called the characteristic equation, or (in much of the physics literature) the auxiliary equation. By the fundamental theorem of algebra it factorises as

p(\lambda) \;=\; (\lambda - \lambda_1)^{k_1}(\lambda - \lambda_2)^{k_2}\cdots(\lambda - \lambda_m)^{k_m} , \qquad k_1 + k_2 + \cdots + k_m = n ,

with $\lambda_1, \dots, \lambda_m$ the distinct roots and $k_j$ the multiplicity of $\lambda_j$ . Since the coefficients are constants they are continuous on every interval, so the whole of Structure of linear equations and the Wronskian applies on $I = \mathbb{R}$ : the solution space of $L(y) = 0$ has dimension $n$ there.

Method(The exponential-shift identity)AdaptedAdapted from Coddington (1989), Ch. 2, §7, Eq. (7.1) and the computation of $\\partial^k L(e^{rx})/\\partial r^k$, pp. 71-72

For every $\lambda \in \mathbb{C}$ and every function $q$ with $n$ derivatives,

L\bigl(e^{\lambda x}q(x)\bigr) \;=\; e^{\lambda x}\sum_{j=0}^{n}\frac{p^{(j)}(\lambda)}{j!}\,q^{(j)}(x) . \tag{$\ast$}

In particular $L(e^{\lambda x}) = p(\lambda)e^{\lambda x}$ , and for an integer $k \ge 0$

L\bigl(x^{k}e^{\lambda x}\bigr) \;=\; e^{\lambda x}\sum_{j=0}^{k}\binom{k}{j}\,p^{(j)}(\lambda)\,x^{k-j} .

Method(A fundamental set from the roots)AdaptedAdapted from Coddington (1989), Ch. 2, §7, Thms 11 and 12, pp. 71-73; Hassani (2013), §14.7, Thm 14.7.1, p. 447

Let $\lambda_1, \dots, \lambda_m$ be the distinct roots of $p$ , with multiplicities $k_1, \dots, k_m$ . Then the $n$ functions

x^{r}e^{\lambda_j x}, \qquad r = 0, 1, \dots, k_j - 1, \quad j = 1, \dots, m ,

are solutions of $L(y) = 0$ , they are linearly independent on every interval $I$ , and they are therefore a fundamental set: every solution of $L(y) = 0$ is

y(x) \;=\; \sum_{j=1}^{m}\bigl(c_{j,0} + c_{j,1}x + \cdots + c_{j,k_j-1}x^{k_j-1}\bigr)e^{\lambda_j x}

with uniquely determined constants $c_{j,r}$ .

Method(Undetermined coefficients for polynomial-times-exponential forcing)AdaptedAdapted from Hassani (2013), §14.7, Thm 14.7.3, p. 448

Let $S$ be a polynomial of degree $d$ . Then $L(y) = e^{\lambda x}S(x)$ always has a particular solution $y_p = e^{\lambda x}q(x)$ with $q$ again a polynomial, and the degree of $q$ is

\deg q \;=\; \begin{cases} d , & p(\lambda) \ne 0 , \\[2pt] d + k_j , & \lambda = \lambda_j \text{ a root of } p \text{ of multiplicity } k_j . \end{cases}

So resonance costs exactly the multiplicity in extra powers of $x$ . Two riders, both visible in the proof: for $p(\lambda) \ne 0$ the polynomial $q$ is unique, while in the resonant case the $k_j$ lowest coefficients of $q$ are left free by the equation and may be set to zero, since $x^{r}e^{\lambda_j x}$ with $r < k_j$ is already part of the complementary function.

Idea of the proofOriginal

Differentiating $e^{\lambda x}$ reproduces it, so $L$ acting on $e^{\lambda x}$ can only multiply it by a number, and that number is $p(\lambda)$ . Everything else follows from watching how this degenerates. Reading $L(e^{\lambda x}) = p(\lambda)e^{\lambda x}$ as an identity in $\lambda$ and differentiating it $k$ times in $\lambda$ — legitimate because $L$ differentiates in $x$ , and the two derivatives commute — produces $x^{k}e^{\lambda x}$ on the left and the derivatives $p^{(j)}(\lambda)$ on the right. A root of multiplicity $k$ is exactly a $\lambda$ at which $p$ and its first $k-1$ derivatives vanish, so the right-hand side vanishes for $x^{0}e^{\lambda x}, \dots, x^{k-1}e^{\lambda x}$ : the multiplicity of the root is the number of extra powers of $x$ it supplies. Collecting the identities for all $k$ gives one formula, $(\ast)$ , which says that $L$ acts on $e^{\lambda x}\times\{\text{polynomials}\}$ as the constant-coefficient operator $\sum_j p^{(j)}(\lambda)D^{j}/j!$ on the polynomial alone. That operator is triangular in the degree, with $p(\lambda)$ on the diagonal; it is invertible when $p(\lambda) \ne 0$ , and when $\lambda$ is a root of multiplicity $k$ it kills polynomials of degree $< k$ and is " $k$ integrations" on the rest. Hence the degree of the answer goes up by exactly $k$ . That degree increase is resonance.

Proof(The exponential-shift identity)AdaptedAdapted from Coddington (1989), Ch. 2, §7, pp. 71-72

Step 1 (the bare exponential). $\frac{d^{k}}{dx^{k}}e^{\lambda x} = \lambda^{k}e^{\lambda x}$ , so

L(e^{\lambda x}) = \bigl(\lambda^{n} + a_1\lambda^{n-1} + \cdots + a_n\bigr)e^{\lambda x} = p(\lambda)e^{\lambda x} .

Step 2 (differentiate in $\lambda$ ). Both sides are entire in $\lambda$ and polynomial in $x$ times $e^{\lambda x}$ , so we may differentiate $k$ times with respect to $\lambda$ and interchange $\partial_\lambda^{k}$ with the $x$ -derivatives in $L$ . On the left, $\partial_\lambda^{k}e^{\lambda x} = x^{k}e^{\lambda x}$ . On the right, Leibniz's rule for the product $p(\lambda)\cdot e^{\lambda x}$ gives

L\bigl(x^{k}e^{\lambda x}\bigr) \;=\; \partial_\lambda^{k}\bigl(p(\lambda)e^{\lambda x}\bigr) \;=\; \sum_{j=0}^{k}\binom{k}{j}p^{(j)}(\lambda)\,x^{k-j}e^{\lambda x} ,

the sum stopping at $j = k$ because $\partial_\lambda^{k-j}e^{\lambda x} = x^{k-j}e^{\lambda x}$ .

Step 3 (rewrite as an operator on $q$ ). For $q(x) = x^{k}$ we have $q^{(j)}(x) = \frac{k!}{(k-j)!}x^{k-j}$ for $j \le k$ and $q^{(j)} = 0$ for $j > k$ , so

\frac{p^{(j)}(\lambda)}{j!}\,q^{(j)}(x) = \frac{k!}{j!\,(k-j)!}\,p^{(j)}(\lambda)\,x^{k-j} = \binom{k}{j}p^{(j)}(\lambda)\,x^{k-j} .

Thus $(\ast)$ holds for every monomial, hence by linearity of both sides for every polynomial $q$ . Finally, $p^{(j)} \equiv 0$ for $j > n$ , so for a general $n$ -times differentiable $q$ the sum in $(\ast)$ is the finite sum written, and the identity extends to it: both sides are linear in $q$ and, expanded by the Leibniz rule, the coefficient of $q^{(j)}$ on the left is $\sum_{i \ge j}\binom{i}{j}a_{n-i}\lambda^{i-j} = p^{(j)}(\lambda)/j!$ , where $a_0 := 1$ . $\square$

Proof(A fundamental set from the roots)AdaptedAdapted from Coddington (1989), Ch. 2, §7, Thms 11 and 12, pp. 71-73

Step 1 (they are solutions). Let $\lambda_j$ have multiplicity $k_j$ . Writing $p(\lambda) = (\lambda - \lambda_j)^{k_j}g(\lambda)$ and differentiating, every derivative $p^{(i)}$ with $i < k_j$ still carries a factor $(\lambda - \lambda_j)$ , so

p(\lambda_j) = p'(\lambda_j) = \cdots = p^{(k_j-1)}(\lambda_j) = 0 .

For $0 \le r \le k_j - 1$ , the identity of Step 2 of the previous proof gives $L(x^{r}e^{\lambda_j x}) = \sum_{i=0}^{r}\binom{r}{i}p^{(i)}(\lambda_j)x^{r-i}e^{\lambda_j x}$ , and every index $i$ occurring satisfies $i \le r \le k_j - 1$ , so every term vanishes. Hence $x^{r}e^{\lambda_j x}$ solves $L(y) = 0$ .

Step 2 (there are $n$ of them). Summing the multiplicities, $k_1 + \cdots + k_m = \deg p = n$ .

Step 3 (independence: reduce to a statement about polynomials). Suppose constants $c_{j,r}$ satisfy

P_1(x)e^{\lambda_1 x} + P_2(x)e^{\lambda_2 x} + \cdots + P_m(x)e^{\lambda_m x} = 0 \qquad\text{on an interval } I,

where $P_j(x) = \sum_{r=0}^{k_j-1}c_{j,r}x^{r}$ . We must show every $P_j$ is the zero polynomial. Suppose not, and (relabelling) let $P_1 \not\equiv 0$ . Multiply the identity by $e^{-\lambda_1 x}$ , which never vanishes:

P_1(x) + P_2(x)e^{(\lambda_2-\lambda_1)x} + \cdots + P_m(x)e^{(\lambda_m-\lambda_1)x} = 0 .

Step 4 (differentiate away one exponential at a time). Differentiate this identity $\deg P_1 + 1$ times. The first term dies. In every other term, $\frac{d}{dx}\bigl(Q(x)e^{\mu x}\bigr) = \bigl(\mu Q + Q'\bigr)e^{\mu x}$ with $\mu = \lambda_i - \lambda_1 \ne 0$ , so each differentiation replaces the polynomial $Q$ by $\mu Q + Q'$ , which has the same degree and the same leading coefficient scaled by the nonzero $\mu$ ; in particular a nonzero polynomial stays nonzero. We are left with

\widetilde{P}_2(x)e^{(\lambda_2-\lambda_1)x} + \cdots + \widetilde{P}_m(x)e^{(\lambda_m-\lambda_1)x} = 0

with $\deg \widetilde{P}_i = \deg P_i$ and $\widetilde{P}_i \not\equiv 0$ whenever $P_i \not\equiv 0$ . One exponential has been eliminated at no cost.

Step 5 (iterate and conclude). Repeating Steps 3 and 4 at most $m-1$ more times leaves a single term $R(x)e^{\mu x} = 0$ on $I$ with $\deg R = \deg P_{i_0}$ for some index $i_0$ with $P_{i_0} \not\equiv 0$ . Since $e^{\mu x} \ne 0$ , $R$ vanishes at every point of the interval $I$ , which has infinitely many points, so $R \equiv 0$ — contradicting that $R$ is nonzero of the same degree. Hence all $P_j \equiv 0$ , i.e. all $c_{j,r} = 0$ : the $n$ functions are linearly independent on $I$ .

Step 6 (they span). They are $n$ linearly independent solutions, so by the dimension theorem they form a basis of the $n$ -dimensional solution space, and the coefficients in the expansion of a given solution are unique. $\square$

Proof(Undetermined coefficients)Original

By $(\ast)$ , for a polynomial $q$ ,

L\bigl(e^{\lambda x}q\bigr) = e^{\lambda x}\,M_\lambda[q], \qquad M_\lambda \;=\; \sum_{j\ge 0}\frac{p^{(j)}(\lambda)}{j!}\,D^{j} , \quad D = \frac{d}{dx} .

So we must solve $M_\lambda[q] = S$ for a polynomial $q$ . Let $d = \deg S$ .

Step 1 (non-resonant case $p(\lambda) \ne 0$ ). $M_\lambda$ maps the space $\mathcal{P}_d$ of polynomials of degree $\le d$ into itself, since $D$ lowers degree. In the basis $1, x, \dots, x^{d}$ its matrix is triangular with every diagonal entry equal to $p(\lambda)$ , because the $j = 0$ term is multiplication by $p(\lambda)$ and all other terms strictly lower the degree. Its determinant is $p(\lambda)^{d+1} \neq 0$ , so $M_\lambda$ is a bijection of $\mathcal{P}_d$ : there is exactly one $q \in \mathcal{P}_d$ with $M_\lambda[q] = S$ . Its degree is exactly $d$ , since the coefficient of $x^{d}$ in $M_\lambda[q]$ is $p(\lambda)$ times that of $q$ .

Step 2 (resonant case: $\lambda$ a root of multiplicity $k$ ). Then $p(\lambda) = \cdots = p^{(k-1)}(\lambda) = 0$ and $p^{(k)}(\lambda) \ne 0$ , so the sum defining $M_\lambda$ starts at $j = k$ :

M_\lambda[q] \;=\; \sum_{j \ge k}\frac{p^{(j)}(\lambda)}{j!}\,q^{(j)} \;=\; N_\lambda\bigl[q^{(k)}\bigr], \qquad N_\lambda = \sum_{i\ge 0}\frac{p^{(k+i)}(\lambda)}{(k+i)!}\,D^{i} .

$N_\lambda$ has the form treated in Step 1 with $p(\lambda)$ replaced by the nonzero number $p^{(k)}(\lambda)/k!$ , so there is a unique polynomial $Q$ of degree $d$ with $N_\lambda[Q] = S$ .

Step 3 (undo the $k$ derivatives). Choose any polynomial $q$ with $q^{(k)} = Q$ , for instance by integrating $Q$ $k$ times with zero constants. Then $\deg q = d + k$ and $M_\lambda[q] = N_\lambda[Q] = S$ , so $y_p = e^{\lambda x}q$ solves $L(y) = e^{\lambda x}S$ . The $k$ constants of integration are free: they change $q$ by a polynomial of degree $< k$ , that is, they change $y_p$ by a member of the complementary function $x^{r}e^{\lambda x}$ , $r < k$ . Choosing them zero gives $\deg q = d + k$ exactly, and no $q$ of lower degree can work, since $q^{(k)}$ must equal the degree- $d$ polynomial $Q$ . $\square$

Remark(Real equations and complex roots)AdaptedAdapted from Riley, Hobson & Bence (2006), §15.1.1, Eq. (15.12), p. 493; Arfken, Weber & Harris (2013), §7.3, p. 343; Hassani (2013), §14.7, p. 447

Suppose $a_1, \dots, a_n$ are real, so that $p$ has real coefficients. Then non-real roots come in conjugate pairs of equal multiplicity: if $\lambda = \alpha + i\beta$ with $\beta \ne 0$ is a root of multiplicity $k$ , so is $\bar\lambda = \alpha - i\beta$ . For each $r < k$ the two complex solutions $x^{r}e^{(\alpha\pm i\beta)x}$ may be replaced by the two real ones

x^{r}e^{\alpha x}\cos\beta x , \qquad x^{r}e^{\alpha x}\sin\beta x ,

their half-sum and half-difference-over- $i$ ; the replacement is invertible, so it carries a fundamental set to a fundamental set. For $r = 0$ , Riley, Hobson and Bence (Riley, Hobson & Bence (2006), Eq. (15.12), p. 493) write the same pair in amplitude-phase form,

c_1e^{(\alpha+i\beta)x} + c_2e^{(\alpha-i\beta)x} = e^{\alpha x}\bigl(d_1\cos\beta x + d_2\sin\beta x\bigr) = Ae^{\alpha x}\sin(\beta x + \phi) ,

with $A$ and $\phi$ arbitrary constants; this is the form a physicist wants, since $A$ is the amplitude and $\phi$ the phase. A real-valued solution is obtained from real $d_1, d_2$ , equivalently from $c_2 = \bar c_1$ .

Method(Solving

L(y) = f

with constant coefficients)AdaptedAdapted from Riley, Hobson & Bence (2006), §15.1.1-15.1.3, pp. 492-495; Hassani (2013), §14.7, pp. 446-449

Write down $p(\lambda)$ by replacing $d^{k}y/dx^{k}$ with $\lambda^{k}$ , and find its roots with their multiplicities. (For $n \le 4$ this is algebra; beyond that, numerically.)
Write the complementary function. Each root $\lambda_j$ of multiplicity $k_j$ contributes $(c_{j,0} + c_{j,1}x + \cdots + c_{j,k_j-1}x^{k_j-1})e^{\lambda_j x}$ . If the equation is real and you want a real answer, convert each conjugate pair to $e^{\alpha x}(\cdots\cos\beta x + \cdots\sin\beta x)$ as in the remark.
Split the forcing $f$ into terms of the form $e^{\lambda x}S(x)$ with $S$ a polynomial, using superposition: a particular solution of $L(y) = f_1 + f_2$ is the sum of particular solutions for $f_1$ and $f_2$ . A pure polynomial is the case $\lambda = 0$ ; $\cos\omega x$ and $\sin\omega x$ are handled either by taking real and imaginary parts of the case $\lambda = i\omega$ , or by the trial pair $b_1\cos\omega x + b_2\sin\omega x$ .
Choose the trial form for each piece: $y_p = e^{\lambda x}q(x)$ with $\deg q = \deg S$ if $p(\lambda) \ne 0$ , and $\deg q = \deg S + k$ if $\lambda$ is a root of multiplicity $k$ . In the resonant case drop the $k$ lowest powers of $x$ from $q$ (they duplicate the complementary function), i.e. take $q(x) = x^{k}\bigl(b_0 + b_1x + \cdots + b_{\deg S}x^{\deg S}\bigr)$ . Equivalently, in the recipe of (Riley, Hobson & Bence (2006), pp. 494-495): take the trial function suggested by $f$ and multiply it by the smallest power of $x$ that leaves no term already present in the complementary function.
Fix the coefficients by substituting the trial form into the equation and matching powers of $x$ in each exponential. The computation is $\sum_{j}p^{(j)}(\lambda)q^{(j)}/j! = S$ , which is quickest to use directly: the derivatives $p^{(j)}(\lambda)$ are a few numbers, computed once.
Add and impose the data. The general solution is complementary function plus particular integral (Structure of linear equations and the Wronskian); then solve the linear system that the $n$ initial or boundary conditions impose on $c_{j,r}$ .

Example(A third-order equation with a resonant forcing)Original

Solve

y''' - 3y'' + 4y \;=\; e^{2x} + 12x , \qquad y(0) = y'(0) = y''(0) = 0 .

Characteristic polynomial. $p(\lambda) = \lambda^{3} - 3\lambda^{2} + 4 = (\lambda - 2)^{2}(\lambda + 1)$ . (Check: $(\lambda^2 - 4\lambda + 4)(\lambda + 1) = \lambda^3 - 3\lambda^2 + 4$ .) Roots: $\lambda = 2$ with multiplicity $2$ , and $\lambda = -1$ , simple.
Complementary function. $y_c = (c_1 + c_2x)e^{2x} + c_3e^{-x}$ .
Derivatives of $p$ , needed in step 5 of the method: $p'(\lambda) = 3\lambda^{2} - 6\lambda$ , $p''(\lambda) = 6\lambda - 6$ , $p'''(\lambda) = 6$ . Hence
$p(2) = 0, \quad p'(2) = 0, \quad p''(2) = 6, \quad p'''(2) = 6; \qquad p(0) = 4, \quad p'(0) = 0 .$
The resonant piece $e^{2x}$ . Here $\lambda = 2$ is a double root, $S = 1$ has degree $0$ , so $q$ has degree $0 + 2 = 2$ and we may take $q = Ax^{2}$ . The equation $\sum_j p^{(j)}(2)q^{(j)}/j! = 1$ reduces, because $p(2) = p'(2) = 0$ , to
$\frac{p''(2)}{2!}q'' + \frac{p'''(2)}{3!}q''' \;=\; \frac{6}{2}(2A) + 0 \;=\; 6A \;=\; 1 ,$
so $A = \tfrac16$ and $y_{p,1} = \tfrac16 x^{2}e^{2x}$ .
The polynomial piece $12x$ . This is $\lambda = 0$ , and $p(0) = 4 \ne 0$ , so $q = ax + b$ has the same degree as $S = 12x$ . The equation $p(0)q + p'(0)q' = 12x$ reads $4(ax + b) = 12x$ , giving $a = 3$ , $b = 0$ and $y_{p,2} = 3x$ .
General solution.
$y(x) \;=\; (c_1 + c_2x)e^{2x} + c_3e^{-x} + \tfrac16x^{2}e^{2x} + 3x .$
Direct substitution confirms that $\tfrac16x^{2}e^{2x} + 3x$ solves the equation.
Impose the data. With $y_p = \tfrac16x^2e^{2x} + 3x$ one finds $y_p(0) = 0$ , $y_p'(0) = 3$ , $y_p''(0) = \tfrac13$ . So the conditions become
$c_1 + c_3 = 0, \qquad 2c_1 + c_2 - c_3 = -3, \qquad 4c_1 + 4c_2 + c_3 = -\tfrac13 ,$
with solution $c_1 = -\tfrac{35}{27}$ , $c_2 = \tfrac89$ , $c_3 = \tfrac{35}{27}$ . Hence
$y(x) \;=\; \left(-\frac{35}{27} + \frac{8}{9}x + \frac{1}{6}x^{2}\right)e^{2x} + \frac{35}{27}e^{-x} + 3x .$
A check. At $x = 1$ this gives $y(1) = -\tfrac{13}{54}e^{2} + \tfrac{35}{27}e^{-1} + 3 = 1.698\,034$ to six decimal places. Differentiating the closed form three times and substituting returns $e^{2x} + 12x$ exactly, and $y, y', y''$ all vanish at $x = 0$ .
Where the resonance shows. The forcing $e^{2x}$ is itself a solution of the homogeneous equation, and so is $xe^{2x}$ ; the response is therefore $x^{2}e^{2x}$ , larger by two powers of $x$ — the multiplicity of the root. Had the forcing been $e^{3x}$ , with $p(3) = 4 \neq 0$ , the response would have been the plain $\tfrac14e^{3x}$ .

ApplicationOriginal

Circuits and mechanics. The series $RLC$ circuit obeys $L\ddot q + R\dot q + q/C = V(t)$ and the mass-spring-damper obeys $m\ddot x + c\dot x + kx = F(t)$ ; these are the same equation, so the characteristic polynomial $m\lambda^{2} + c\lambda + k$ decides the behaviour of both. Its two roots being real and distinct, equal, or complex is exactly the over-, critically- and under-damped case, worked out in The damped and driven harmonic oscillator. The algebraic condition $c^{2} = 4mk$ for critical damping — a double root — is the reason instrument movements and door closers are specified by a dimensionless damping ratio: it is the quantity that says where you are relative to the double root.

Stability from the roots. Since each solution is a sum of terms $x^{r}e^{\lambda_j x}$ , every solution of $L(y) = 0$ tends to $0$ as $x \to +\infty$ precisely when every root has negative real part, and solutions stay bounded when no root has positive real part and the purely imaginary roots are simple. This is why control engineering reads stability off the locations of the roots of a polynomial, and why criteria (Routh-Hurwitz, Nyquist) that decide where the roots lie without computing them are useful: the linear algebra is already done by the theorem above.

Resonance and secular growth. The degree increase in the resonant case is the mathematical content of resonance. Driving an undamped oscillator $\ddot x + \omega_0^{2}x = F\cos\omega_0 t$ at its own frequency means forcing with $e^{\pm i\omega_0 t}$ , which is a simple root of $\lambda^{2} + \omega_0^{2}$ ; the response therefore gains one power of $t$ , $x_p = (F/2\omega_0)\,t\sin\omega_0 t$ , and grows without bound. The same count appears in perturbation theory as the secular terms that must be removed by renormalising the frequency, and in quantum mechanics as the linear-in-time growth of a transition amplitude on resonance that underlies Fermi's golden rule.

Discrete analogues. Linear recurrence relations with constant coefficients, which arise in numerical schemes and in statistical physics transfer matrices, have the same structure with $e^{\lambda x}$ replaced by $\mu^{n}$ : a root of multiplicity $k$ of the characteristic polynomial contributes $\mu^{n}, n\mu^{n}, \dots, n^{k-1}\mu^{n}$ , and a forcing term that duplicates a root gains a factor $n$ . The algebra below the surface is the same shift identity.

Method

Variation of parameters

Adapted

Add to PDF

Given a fundamental set of the homogeneous equation, replacing its constants by functions and imposing $n-1$ convenient constraints turns the inhomogeneous equation into a linear system whose determinant is the Wronskian; it yields a particular solution as an integral, for any continuous forcing and any order.

Prerequisites: Structure of linear equations and the Wronskian

Note(The problem)Original

Suppose we have solved the homogeneous equation and want a particular solution of the inhomogeneous one. The method of undetermined coefficients works only for constant coefficients and for special forcings. Variation of parameters needs neither: it works for variable coefficients and for any continuous forcing, at the price of two integrals (for $n = 2$ ) that may or may not be elementary. It is the one general-purpose method, and it is the ancestor of the Green's function.

Definition(The set-up)AdaptedAdapted from Coddington (1989), Ch. 3, §6, p. 122; Riley, Hobson & Bence (2006), §15.2.4, Eq. (15.53), p. 508

Let $I$ be an interval, let $a_1, \dots, a_n$ and $b$ be continuous on $I$ , let

L(y) \;=\; y^{(n)} + a_1(x)y^{(n-1)} + \cdots + a_n(x)y ,

and let $\varphi_1, \dots, \varphi_n$ be a fundamental set for $L(y) = 0$ on $I$ , with Wronskian $W = W(\varphi_1, \dots, \varphi_n)$ , which is nowhere zero on $I$ (Structure of linear equations and the Wronskian). We look for a particular solution of $L(y) = b$ in the form

y_p \;=\; u_1\varphi_1 + u_2\varphi_2 + \cdots + u_n\varphi_n ,

where the constants of the complementary function have been replaced by unknown functions $u_1, \dots, u_n$ on $I$ — hence the names variation of parameters and variation of constants.

If the equation is presented unnormalised as $A_n(x)y^{(n)} + \cdots + A_0(x)y = f(x)$ with $A_n \ne 0$ on $I$ , then $b = f/A_n$ ; this factor is easy to forget.

Method(Second order)AdaptedAdapted from Arfken, Weber & Harris (2013), §7.7, Eqs. (7.95)-(7.98), pp. 375-376 and Exercise 7.7.1, p. 377; Coddington (1989), Ch. 2, §6, Eq. (6.3) and Thm 10, pp. 68

For $y'' + p(x)y' + q(x)y = b(x)$ with $p, q, b$ continuous on $I$ and $\{y_1, y_2\}$ a fundamental set of the homogeneous equation, put $W = y_1y_2' - y_1'y_2$ . Then $u_1' = -y_2b/W$ and $u_2' = y_1b/W$ , so for any $x_0 \in I$

y_p(x) \;=\; -\,y_1(x)\int_{x_0}^{x}\frac{y_2(t)\,b(t)}{W(t)}\,dt \;+\; y_2(x)\int_{x_0}^{x}\frac{y_1(t)\,b(t)}{W(t)}\,dt

is a solution of $y'' + py' + qy = b$ on $I$ . Equivalently, with the two integrals combined under one integral sign,

y_p(x) \;=\; \int_{x_0}^{x} K(x,t)\,b(t)\,dt , \qquad K(x,t) \;=\; \frac{y_1(t)\,y_2(x) - y_1(x)\,y_2(t)}{W(t)} .

This $y_p$ satisfies $y_p(x_0) = y_p'(x_0) = 0$ . The general solution of the equation is $y_p + c_1y_1 + c_2y_2$ .

Method(General order)AdaptedAdapted from Coddington (1989), Ch. 3, §6, Thm 11 and Eq. (6.2), pp. 122-123; Riley, Hobson & Bence (2006), §15.2.4, Eq. (15.55), pp. 508-509

With the set-up above, the functions $u_1', \dots, u_n'$ are determined at each $x \in I$ by the $n$ linear equations

\begin{aligned} u_1'\varphi_1 + u_2'\varphi_2 + \cdots + u_n'\varphi_n &= 0, \\ u_1'\varphi_1' + u_2'\varphi_2' + \cdots + u_n'\varphi_n' &= 0, \\ &\ \ \vdots \\ u_1'\varphi_1^{(n-2)} + \cdots + u_n'\varphi_n^{(n-2)} &= 0, \\ u_1'\varphi_1^{(n-1)} + \cdots + u_n'\varphi_n^{(n-1)} &= b . \end{aligned}

The determinant of this system is $W \ne 0$ , so it has a unique solution, and by Cramer's rule

u_k'(x) \;=\; \frac{W_k(x)\,b(x)}{W(x)} ,

and here $W_k$ denotes what $W$ becomes when its $k$ -th column $\bigl(\varphi_k, \varphi_k', \dots, \varphi_k^{(n-1)}\bigr)$ is overwritten by $(0, 0, \dots, 0, 1)$ . Hence, for any $x_0 \in I$ ,

y_p(x) \;=\; \sum_{k=1}^{n}\varphi_k(x)\int_{x_0}^{x}\frac{W_k(t)\,b(t)}{W(t)}\,dt

is a solution of $L(y) = b$ on $I$ , and every solution of $L(y) = b$ is $y_p + c_1\varphi_1 + \cdots + c_n\varphi_n$ for constants $c_k$ .

Idea of the proofOriginal

We have $n$ unknown functions and only one equation to satisfy, so there is a large amount of freedom; the trick is to spend it on making the computation collapse. Differentiate $y_p = \sum_k u_k\varphi_k$ once: the answer is what a genuine solution would give, $\sum_k u_k\varphi_k'$ , plus the unwanted $\sum_k u_k'\varphi_k$ . Declare that unwanted sum to be zero. Differentiate again and do the same. After $n-1$ such declarations, $y_p$ and its first $n-1$ derivatives look exactly as if the $u_k$ were constants, so when we finally substitute into $L$ , every term groups into $\sum_k u_k L(\varphi_k) = 0$ except the single leftover $\sum_k u_k'\varphi_k^{(n-1)}$ , which must therefore equal $b$ . The $n$ conditions are a linear system for $u_1', \dots, u_n'$ whose coefficient matrix is the Wronskian matrix — invertible precisely because the $\varphi_k$ are a fundamental set. So the freedom was spent legitimately: nothing was over-determined, and the remaining work is $n$ integrations.

ProofAdaptedAdapted from Riley, Hobson & Bence (2006), §15.2.4, pp. 508-509; Arfken, Weber & Harris (2013), §7.7, pp. 375-376

Step 1 (the first $n-1$ constraints). Let $u_1, \dots, u_n$ be any differentiable functions on $I$ and set $y_p = \sum_k u_k\varphi_k$ . Differentiating,

y_p' \;=\; \sum_k u_k\varphi_k' \;+\; \sum_k u_k'\varphi_k .

Impose $\sum_k u_k'\varphi_k = 0$ . Then $y_p' = \sum_k u_k\varphi_k'$ , which has the same shape as $y_p$ with $\varphi_k$ replaced by $\varphi_k'$ . Differentiating that and imposing $\sum_k u_k'\varphi_k' = 0$ gives $y_p'' = \sum_k u_k\varphi_k''$ . Repeating, after imposing the $m$ -th constraint $\sum_k u_k'\varphi_k^{(m-1)} = 0$ for $m = 1, \dots, n-1$ we have

y_p^{(m)} \;=\; \sum_k u_k\varphi_k^{(m)} , \qquad m = 0, 1, \dots, n-1 .

Step 2 (the last derivative). One more differentiation, with no constraint left to impose, gives

y_p^{(n)} \;=\; \sum_k u_k\varphi_k^{(n)} \;+\; \sum_k u_k'\varphi_k^{(n-1)} .

Step 3 (substitute). Insert these into $L(y_p) = y_p^{(n)} + a_1y_p^{(n-1)} + \cdots + a_ny_p$ and collect the terms containing $u_k$ rather than $u_k'$ :

L(y_p) \;=\; \sum_{k=1}^{n}u_k\Bigl[\varphi_k^{(n)} + a_1\varphi_k^{(n-1)} + \cdots + a_n\varphi_k\Bigr] \;+\; \sum_{k=1}^{n}u_k'\varphi_k^{(n-1)} .

Every bracket is $L(\varphi_k) = 0$ , because each $\varphi_k$ solves the homogeneous equation — this is where that hypothesis enters. Hence $L(y_p) = b$ if and only if

\sum_{k=1}^{n}u_k'\varphi_k^{(n-1)} \;=\; b ,

which is the last equation of the system. So the $n$ displayed equations are together equivalent to " $y_p$ solves $L(y) = b$ and the $n-1$ constraints hold".

Step 4 (the system is solvable). For each fixed $x$ , the $n$ equations are linear in the numbers $u_1'(x), \dots, u_n'(x)$ , with coefficient matrix the Wronskian matrix of $\varphi_1, \dots, \varphi_n$ at $x$ . Since the $\varphi_k$ are a fundamental set, $W(x) \ne 0$ for every $x \in I$ (the Wronskian criterion), so by Cramer's rule there is exactly one solution, namely $u_k'(x) = W_k(x)b(x)/W(x)$ with $W_k$ as described: the right-hand side column is $(0, \dots, 0, b)$ , and factoring $b$ out of the determinant in the numerator leaves $W_k$ .

Step 5 (integrate). $W_k$ is continuous (a determinant of continuous functions), $W$ is continuous and nonvanishing, and $b$ is continuous, so $u_k' = W_kb/W$ is continuous on $I$ and $u_k(x) = \int_{x_0}^{x}W_kb/W\,dt$ is a well-defined antiderivative on $I$ . With these $u_k$ , Steps 1-3 run forwards and give $L(y_p) = b$ . Adding any constants $c_k$ to the $u_k$ adds $\sum_k c_k\varphi_k$ to $y_p$ , which is the complementary function; by the structure theorem these are all the solutions.

Step 6 (the case $n = 2$ ). Here the system is $u_1'y_1 + u_2'y_2 = 0$ , $u_1'y_1' + u_2'y_2' = b$ , and Cramer's rule gives

u_1' = \frac{1}{W}\begin{vmatrix}0 & y_2\\ b & y_2'\end{vmatrix} = -\frac{y_2b}{W}, \qquad u_2' = \frac{1}{W}\begin{vmatrix}y_1 & 0\\ y_1' & b\end{vmatrix} = \frac{y_1b}{W} .

Integrating from $x_0$ and substituting into $y_p = u_1y_1 + u_2y_2$ gives both displayed forms of the second-order formula; the kernel form follows by putting the two integrals over a common integrand. Finally $y_p(x_0) = 0$ because both integrals vanish at $x_0$ , and $y_p'(x_0) = u_1(x_0)y_1'(x_0) + u_2(x_0)y_2'(x_0) = 0$ by Step 1. $\square$

Remark(Where the hypotheses are used, and practical points)Original

A full fundamental set is needed, not just one solution: $W \ne 0$ in Step 4 is what makes the system solvable. If only one solution is known, use reduction of order first, or the substitution $y = u(x)v(x)$ of (Riley, Hobson & Bence (2006), §15.2.3, p. 507), which lowers the order by one without needing the rest of the set.
Normalise first. The last equation of the system has $b$ on the right, where $b$ is the forcing after dividing by the coefficient of $y^{(n)}$ . In (Riley, Hobson & Bence (2006), Eq. (15.55), p. 509) this appears as $f(x)/a_n(x)$ .
Continuity of $b$ is used only in Step 5, to integrate. Piecewise-continuous forcing (a switched voltage, an impulse) is handled by applying the formula on each interval of continuity, or directly: the integral formula still makes sense and produces a solution with $n-1$ continuous derivatives.
The constants of integration are free. Choosing them zero gives the $y_p$ above, which vanishes together with its first $n-1$ derivatives at $x_0$ — convenient for an initial-value problem at $x_0$ . Keeping them simply redefines the arbitrary constants of the complementary function, as (Riley, Hobson & Bence (2006), p. 510) notes.
A particular solution may be simplified. Any multiple of a $\varphi_k$ may be dropped from $y_p$ , since it belongs to the complementary function. This is worth doing: the raw formula often produces terms that cancel or duplicate the homogeneous solutions.
The kernel is a Green's function. For $n = 2$ , $K(x,t)$ depends on the equation alone, not on $b$ : once computed, it solves the equation for every forcing. Extended by $K = 0$ for $t > x$ it is the causal Green's function of Causal Green's functions for initial-value problems, and the method is the one-dimensional ancestor of Duhamel's principle.

Method(Variation of parameters)AdaptedAdapted from Riley, Hobson & Bence (2006), §15.2.4, Solution method, p. 510; Arfken, Weber & Harris (2013), §7.7, pp. 375-376

Normalise: divide by the coefficient of $y^{(n)}$ so that the equation reads $L(y) = b$ with leading coefficient $1$ .
Get a fundamental set $\varphi_1, \dots, \varphi_n$ of $L(y) = 0$ and compute $W = W(\varphi_1, \dots, \varphi_n)$ . Abel's formula is often the quickest route to $W$ (Structure of linear equations and the Wronskian).
Write the system: all the derivative orders $0$ to $n-2$ give $\sum_k u_k'\varphi_k^{(m)} = 0$ , and order $n-1$ gives $\sum_k u_k'\varphi_k^{(n-1)} = b$ .
Solve it for $u_1', \dots, u_n'$ . For $n = 2$ use $u_1' = -y_2b/W$ , $u_2' = y_1b/W$ directly; for larger $n$ use $u_k' = W_kb/W$ or just eliminate by hand.
Integrate each $u_k'$ , dropping the constants of integration.
Assemble $y_p = \sum_k u_k\varphi_k$ , discard any part of it that is a multiple of some $\varphi_k$ , add the complementary function, and impose the data.

Example(A second-order equation with variable coefficients)Original

Solve

x^{2}y'' - 2x\,y' + 2y \;=\; x^{3}\ln x \qquad \text{on } x > 0, \qquad y(1) = y'(1) = 0 .

Normalise. Dividing by $x^{2}$ : $y'' - \frac{2}{x}y' + \frac{2}{x^{2}}y = x\ln x$ , so $b(x) = x\ln x$ — not $x^{3}\ln x$ .
Fundamental set. The homogeneous equation is an Euler equation (Euler-Cauchy equations); trying $y = x^{s}$ gives $s(s-1) - 2s + 2 = (s-1)(s-2) = 0$ , so $y_1 = x$ and $y_2 = x^{2}$ . Their Wronskian is
$W = y_1y_2' - y_1'y_2 = x\cdot 2x - 1\cdot x^{2} = x^{2} ,$
nonzero for $x > 0$ , so this is a fundamental set there. (Abel's formula agrees: $a_1 = -2/x$ gives $W \propto \exp\int(2/x)\,dx = x^{2}$ .)
Solve for $u_1', u_2'$ .
$u_1' = -\frac{y_2b}{W} = -\frac{x^{2}\cdot x\ln x}{x^{2}} = -x\ln x , \qquad u_2' = \frac{y_1b}{W} = \frac{x\cdot x\ln x}{x^{2}} = \ln x .$
Integrate (by parts, constants dropped):
$u_1 = -\int x\ln x\,dx = -\frac{x^{2}}{2}\ln x + \frac{x^{2}}{4} , \qquad u_2 = \int \ln x\,dx = x\ln x - x .$
Assemble.
$y_p = u_1y_1 + u_2y_2 = \left(-\frac{x^{2}\ln x}{2} + \frac{x^{2}}{4}\right)x + \bigl(x\ln x - x\bigr)x^{2} = \frac{x^{3}\ln x}{2} - \frac{3x^{3}}{4} .$
Substituting back: $y_p' = \frac{3x^{2}\ln x}{2} + \frac{x^{2}}{2} - \frac{9x^{2}}{4}$ and $y_p'' = 3x\ln x + \frac{3x}{2} + x - \frac{9x}{2} = 3x\ln x - 2x$ , and then $x^{2}y_p'' - 2xy_p' + 2y_p = x^{3}\ln x$ exactly.
General solution and data.
$y = c_1x + c_2x^{2} + \frac{x^{3}}{4}\bigl(2\ln x - 3\bigr) .$
At $x = 1$ : $y_p(1) = -\tfrac34$ and $y_p'(1) = \tfrac12 - \tfrac94 = -\tfrac74$ , so $c_1 + c_2 = \tfrac34$ and $c_1 + 2c_2 = \tfrac74$ , giving $c_2 = 1$ , $c_1 = -\tfrac14$ :
$y(x) = -\frac{x}{4} + x^{2} + \frac{x^{3}}{4}\bigl(2\ln x - 3\bigr) .$
A number to check against. $y(2) = -\tfrac12 + 4 + 2(2\ln 2 - 3) = 3.5 + 2(2\ln 2 - 3) = 0.272\,589$ to six decimal places.

Example(Third order, with the

W_k

formula)Original

Solve $y''' - y' = e^{2x}$ on $\mathbb{R}$ .

Fundamental set. $p(\lambda) = \lambda^{3} - \lambda = \lambda(\lambda-1)(\lambda+1)$ , so $\varphi_1 = 1$ , $\varphi_2 = e^{x}$ , $\varphi_3 = e^{-x}$ (Constant-coefficient equations).
Wronskian.
$W = \begin{vmatrix} 1 & e^{x} & e^{-x} \\ 0 & e^{x} & -e^{-x} \\ 0 & e^{x} & e^{-x}\end{vmatrix} = 1\cdot\bigl(e^{x}e^{-x} + e^{-x}e^{x}\bigr) = 2 ,$
constant, as Abel's formula predicts since the coefficient of $y''$ is $0$ .
The modified determinants. Replacing the $k$ -th column by $(0,0,1)^{\mathsf T}$ :
$W_1 = \begin{vmatrix} 0 & e^{x} & e^{-x} \\ 0 & e^{x} & -e^{-x} \\ 1 & e^{x} & e^{-x}\end{vmatrix} = -2 , \qquad W_2 = \begin{vmatrix} 1 & 0 & e^{-x} \\ 0 & 0 & -e^{-x} \\ 0 & 1 & e^{-x}\end{vmatrix} = e^{-x} , \qquad W_3 = \begin{vmatrix} 1 & e^{x} & 0 \\ 0 & e^{x} & 0 \\ 0 & e^{x} & 1\end{vmatrix} = e^{x} .$
Integrate $u_k' = W_kb/W$ with $b = e^{2x}$ .
$u_1' = -e^{2x} \Rightarrow u_1 = -\tfrac12e^{2x}, \qquad u_2' = \tfrac12e^{x} \Rightarrow u_2 = \tfrac12e^{x}, \qquad u_3' = \tfrac12e^{3x} \Rightarrow u_3 = \tfrac16e^{3x} .$
Assemble.
$y_p = -\tfrac12e^{2x} + \tfrac12e^{x}\cdot e^{x} + \tfrac16e^{3x}\cdot e^{-x} = \tfrac16e^{2x} .$
The first two terms cancel, which is typical: the raw formula produces more structure than the answer has. Check: $y_p''' - y_p' = \tfrac{8}{6}e^{2x} - \tfrac{2}{6}e^{2x} = e^{2x}$ .
Cross-check with undetermined coefficients. Since $p(2) = 8 - 2 = 6 \ne 0$ , that method predicts $y_p = e^{2x}/p(2) = e^{2x}/6$ — the same answer, obtained in one line. Variation of parameters is the method to reach for when the forcing is not of that special form: replacing $e^{2x}$ by, say, $1/(1+x^{2})$ leaves steps 1-3 untouched and changes only the integrals.

ApplicationOriginal

Impulse response and causality. Read the kernel form backwards: $y_p(x) = \int_{x_0}^{x}K(x,t)b(t)\,dt$ says that the response at $x$ is a superposition of the responses to the forcing applied at each earlier instant $t$ , with weight $K(x,t)$ . For fixed $t$ , the function $x \mapsto K(x,t)$ solves the homogeneous equation and satisfies $K(t,t) = 0$ and $\partial_xK(x,t)\big|_{x=t} = \bigl(y_1(t)y_2'(t) - y_1'(t)y_2(t)\bigr)/W(t) = 1$ : it is the free motion started from rest with unit velocity at time $t$ , that is, the response to a unit impulse delivered at $t$ . In circuit and control engineering that function is the impulse response; its Fourier or Laplace transform is the transfer function. The construction here is where it comes from, and it shows that a linear system is characterised by a single two-variable kernel.

Forced oscillators with awkward forcing. For $\ddot x + \omega_0^{2}x = F(t)/m$ the fundamental set is $\cos\omega_0t$ , $\sin\omega_0t$ with $W = \omega_0$ , and the formula gives

x_p(t) = \frac{1}{m\omega_0}\int_{0}^{t}\sin\bigl(\omega_0(t - s)\bigr)F(s)\,ds ,

valid for any measured or switched force — a seismometer trace, a shock pulse — where no trial function would help. The same integral, with $\sin$ replaced by the appropriate decaying solution, gives the response of a damped oscillator and is the standard starting point for response spectra in earthquake engineering.

Perturbation theory. In first-order perturbation theory one knows two solutions of an unperturbed equation and treats the perturbation as forcing. For the radial Schrodinger equation $u'' + (k^{2} - V)u = 0$ , writing the perturbing potential as the forcing $b = \delta V\,u$ and applying this formula produces the integral equation whose iteration is the Born series, and the first term is the Born approximation (Scattering and the Born approximation). The appearance of $1/W$ there is the normalisation of the regular and irregular radial solutions.

Systems and Duhamel. The same calculation for a first-order system $X' = AX + G(t)$ reads $X(t) = e^{tA}\bigl(X_0 + \int_0^{t}e^{-sA}G(s)\,ds\bigr)$ , as Hirsch, Smale and Devaney (Hirsch, Smale & Devaney (2013), §6.5, pp. 131-132) give it: the fundamental matrix $e^{tA}$ plays the role of the $\varphi_k$ and its inverse the role of $1/W$ . That form generalises directly to Duhamel's principle and to forced evolution equations in partial differential equations.

Method

Reduction of order

Adapted

Add to PDF

One known solution $\varphi_1$ of a linear homogeneous equation, nonvanishing on an interval, lowers the order by one through the substitution $y = u\varphi_1$ ; for second-order equations this gives the second solution in closed form, $\varphi_2 = \varphi_1\int \varphi_1^{-2}e^{-\int a_1}$ , which is Abel's formula integrated once.

Prerequisites: Structure of linear equations and the Wronskian

Note(The problem)Original

It is common to find one solution of a linear homogeneous equation by luck, inspection or a series. Reduction of order converts that one solution into a lowering of the order by one — and for a second-order equation, directly into the missing second solution, in closed form. The idea is the same as variation of parameters: let a constant become a function. The difference is that here the equation is homogeneous, so instead of a particular solution we gain a smaller homogeneous equation.

Method(Reduction of order, general order)AdaptedAdapted from Coddington (1989), Ch. 3, §5, Eq. (5.1) and Thm 9, pp. 118-119

Let $a_1, \dots, a_n$ be continuous on an interval $I$ , let

L(y) = y^{(n)} + a_1(x)y^{(n-1)} + \cdots + a_n(x)y ,

and let $\varphi_1$ be a solution of $L(y) = 0$ on $I$ with

\varphi_1(x) \ne 0 \quad\text{for every } x \in I .

Substituting $y = u\varphi_1$ into $L(y) = 0$ and writing $v = u'$ gives a linear homogeneous equation of order $n-1$ for $v$ ,

\varphi_1\,v^{(n-1)} + \cdots + \bigl[n\varphi_1^{(n-1)} + (n-1)a_1\varphi_1^{(n-2)} + \cdots + a_{n-1}\varphi_1\bigr]v \;=\; 0 , \tag{R}

whose leading coefficient $\varphi_1$ does not vanish on $I$ , so that $(\mathrm{R})$ has $n-1$ linearly independent solutions $v_2, \dots, v_n$ on $I$ . If $x_0 \in I$ and

u_k(x) = \int_{x_0}^{x} v_k(t)\,dt , \qquad k = 2, \dots, n ,

then

\varphi_1, \quad u_2\varphi_1, \quad \dots, \quad u_n\varphi_1

is a basis for the solutions of $L(y) = 0$ on $I$ .

Method(Second order: the second solution in closed form)AdaptedAdapted from Coddington (1989), Ch. 3, §5, Eq. (5.7) and Thm 10, p. 120; Hassani (2013), §14.4.1, Thm 14.4.4, p. 426; Arfken, Weber & Harris (2013), §7.6, Eq. (7.67), p. 363; Riley, Hobson & Bence (2006), §16.4.1, Eq. (16.25), p. 544

Let $a_1, a_2$ be continuous on $I$ , consider

y'' + a_1(x)y' + a_2(x)y = 0 ,

and let $\varphi_1$ be a solution with $\varphi_1(x) \ne 0$ on $I$ . Fix $x_0 \in I$ . Then

\varphi_2(x) \;=\; \varphi_1(x)\int_{x_0}^{x}\frac{1}{\varphi_1(s)^{2}}\exp\left[-\int_{x_0}^{s}a_1(t)\,dt\right]ds

is a solution on $I$ , and $\{\varphi_1, \varphi_2\}$ is a basis for the solutions on $I$ . Its Wronskian is

W(\varphi_1, \varphi_2)(x) \;=\; \exp\left[-\int_{x_0}^{x}a_1(t)\,dt\right] ,

in particular $W(x_0) = 1$ : the formula is normalised so that independence is automatic.

Two special cases are worth remembering. If $a_1 \equiv 0$ (as in $y'' + a_2y = 0$ , the form every second-order equation can be put into), then $\varphi_2 = \varphi_1\int \varphi_1^{-2}$ . And changing $x_0$ , or adding a constant to the outer integral, changes $\varphi_2$ only by a multiple of $\varphi_1$ , so any convenient antiderivative may be used and any resulting multiple of $\varphi_1$ discarded.

Idea of the proofOriginal

Two ideas, and they meet.

The substitution. If we look for further solutions in the form $y = u\varphi_1$ , then the terms in which $u$ appears undifferentiated assemble into $u\,L(\varphi_1) = 0$ and disappear. What survives involves only $u', u'', \dots, u^{(n)}$ , which is an equation of order $n-1$ in $v = u'$ . We have not solved anything, but we have one order less to solve, and we needed $\varphi_1 \neq 0$ only to keep the new leading coefficient nonvanishing.

Abel's formula. For $n = 2$ the reduced equation is first-order and so can be integrated explicitly; but there is a shortcut that bypasses it. Abel's formula already tells us the Wronskian of the pair $\{\varphi_1, \varphi_2\}$ before we know $\varphi_2$ . And the Wronskian is exactly $\varphi_1^{2}(\varphi_2/\varphi_1)'$ . So $(\varphi_2/\varphi_1)'$ is known; one integration gives $\varphi_2/\varphi_1$ , and the formula drops out. This is why the second solution of a second-order equation costs only an integral: the Wronskian does the work.

Proof(Reduction of order, general order)AdaptedAdapted from Coddington (1989), Ch. 3, §5, pp. 118-119

Step 1 (substitute). Let $u$ be any $n$ -times differentiable function on $I$ and expand $L(u\varphi_1)$ by the Leibniz rule. For $0 \le m \le n$ ,

(u\varphi_1)^{(m)} = \sum_{i=0}^{m}\binom{m}{i}u^{(i)}\varphi_1^{(m-i)} ,

so, with $a_0 := 1$ ,

L(u\varphi_1) = \sum_{m=0}^{n}a_{n-m}\,(u\varphi_1)^{(m)} = \sum_{i=0}^{n}u^{(i)}\underbrace{\left[\sum_{m\ge i}\binom{m}{i}a_{n-m}\varphi_1^{(m-i)}\right]}_{=\,b_i} .

Step 2 (the $u$ term vanishes). The coefficient of $u^{(0)} = u$ is $b_0 = \sum_m a_{n-m}\varphi_1^{(m)} = L(\varphi_1) = 0$ , because $\varphi_1$ solves the equation — this is the only place that hypothesis is used. Hence

L(u\varphi_1) = \sum_{i=1}^{n}b_i\,u^{(i)} = \sum_{j=0}^{n-1}b_{j+1}\,v^{(j)}, \qquad v = u' .

So $u\varphi_1$ solves $L(y) = 0$ exactly when $v = u'$ solves the order- $(n-1)$ linear homogeneous equation $\sum_{j}b_{j+1}v^{(j)} = 0$ , which is $(\mathrm{R})$ .

Step 3 (the reduced equation is non-singular). Its leading coefficient is $b_n = \binom{n}{n}a_0\varphi_1 = \varphi_1$ , nonzero on $I$ by hypothesis. Dividing by it gives a normalised equation of order $n-1$ with coefficients continuous on $I$ , so by the dimension theorem it has $n-1$ linearly independent solutions $v_2, \dots, v_n$ on $I$ — and no more.

Step 4 (the $n$ functions solve the equation). $\varphi_1 = 1\cdot\varphi_1$ corresponds to $u \equiv 1$ , $v \equiv 0$ . Each $u_k = \int_{x_0}^{x}v_k$ is differentiable with $u_k' = v_k$ , so $u_k\varphi_1$ solves $L(y) = 0$ by Step 2.

Step 5 (independence). Suppose $c_1\varphi_1 + c_2u_2\varphi_1 + \cdots + c_nu_n\varphi_1 = 0$ on $I$ . Dividing by $\varphi_1$ , which is nonvanishing,

c_1 + c_2u_2(x) + \cdots + c_nu_n(x) = 0 \qquad (x \in I) .

Differentiating, $c_2v_2 + \cdots + c_nv_n = 0$ on $I$ , and the independence of $v_2, \dots, v_n$ forces $c_2 = \cdots = c_n = 0$ . The displayed identity then gives $c_1 = 0$ .

Step 6 (basis). These are $n$ linearly independent solutions of an equation of order $n$ , hence a basis (Structure of linear equations and the Wronskian). $\square$

Proof(Second order, by Abel's formula)AdaptedAdapted from Arfken, Weber & Harris (2013), §7.6, Eqs. (7.63)-(7.67), pp. 362-363; Riley, Hobson & Bence (2006), §16.4.1, p. 544

Step 1 (the quotient derivative). Wherever $\varphi_1 \ne 0$ , for any differentiable $\varphi_2$ ,

\left(\frac{\varphi_2}{\varphi_1}\right)' = \frac{\varphi_1\varphi_2' - \varphi_1'\varphi_2}{\varphi_1^{2}} = \frac{W(\varphi_1,\varphi_2)}{\varphi_1^{2}} .

Step 2 (what the Wronskian must be). Suppose $\varphi_2$ is a solution. By Abel's formula, $W(\varphi_1,\varphi_2)(x) = W(x_0)\,e^{-A(x)}$ where $A(x) = \int_{x_0}^{x}a_1$ . We are free to look for the solution with $W(x_0) = 1$ , since multiplying $\varphi_2$ by a constant rescales $W$ . So we require

\left(\frac{\varphi_2}{\varphi_1}\right)'(x) = \frac{e^{-A(x)}}{\varphi_1(x)^{2}} .

Step 3 (integrate). The right-hand side is continuous on $I$ ( $a_1$ is continuous, $\varphi_1$ is continuous and nonvanishing), so it has the antiderivative $\int_{x_0}^{x}\varphi_1(s)^{-2}e^{-A(s)}\,ds$ , and

\varphi_2(x) = \varphi_1(x)\int_{x_0}^{x}\frac{e^{-A(s)}}{\varphi_1(s)^{2}}\,ds

is forced, up to adding a constant multiple of $\varphi_1$ .

Step 4 (it really is a solution). So far we have only derived a necessary form; now verify it. Put $g(x) = \int_{x_0}^{x}\varphi_1^{-2}e^{-A}$ , so $\varphi_2 = \varphi_1 g$ with $g' = \varphi_1^{-2}e^{-A}$ . Then

\varphi_2' = \varphi_1'g + \varphi_1g', \qquad \varphi_2'' = \varphi_1''g + 2\varphi_1'g' + \varphi_1g'' ,

and therefore

\varphi_2'' + a_1\varphi_2' + a_2\varphi_2 = \underbrace{\bigl(\varphi_1'' + a_1\varphi_1' + a_2\varphi_1\bigr)}_{=\,0}g \;+\; 2\varphi_1'g' + \varphi_1g'' + a_1\varphi_1g' .

For the remainder, differentiate $g' = \varphi_1^{-2}e^{-A}$ : since $A' = a_1$ ,

g'' = -2\frac{\varphi_1'}{\varphi_1^{3}}e^{-A} - a_1\frac{e^{-A}}{\varphi_1^{2}} = -\left(\frac{2\varphi_1'}{\varphi_1} + a_1\right)g' ,

so $\varphi_1g'' = -2\varphi_1'g' - a_1\varphi_1g'$ , which cancels the remaining three terms exactly. Hence $\varphi_2$ solves the equation.

Step 5 (basis). By Step 1, $W(\varphi_1,\varphi_2) = \varphi_1^{2}g' = e^{-A}$ , which is nonzero everywhere on $I$ ; by the Wronskian criterion, $\{\varphi_1,\varphi_2\}$ is a basis. $\square$

The same formula comes out of the substitution route of the previous proof: for $n = 2$ the reduced equation $(\mathrm{R})$ is $\varphi_1v' + (2\varphi_1' + a_1\varphi_1)v = 0$ , whose integrating factor is $\varphi_1e^{A}$ , giving $\varphi_1^{2}e^{A}v = \text{const}$ and hence $v = u' = \varphi_1^{-2}e^{-A}$ as above.

Remark(Where the hypotheses are used)Original

$\varphi_1 \ne 0$ on $I$ enters twice: it keeps the leading coefficient of the reduced equation away from zero (Step 3 of the first proof), and it is needed to divide by $\varphi_1^{2}$ . The hypothesis is genuinely restrictive. For $y'' + y = 0$ with $\varphi_1 = \sin x$ , the formula $\varphi_2 = \sin x\int\csc^{2} = -\cos x$ is valid on $(0,\pi)$ and the answer happens to extend to $\mathbb{R}$ , but the derivation does not: it must be run on each interval between zeros of $\sin$ and the pieces matched. In general the second solution exists on all of $I$ (by the existence theorem) even though this formula for it does not.
A zero of $\varphi_1$ is where the second solution misbehaves. Near a simple zero at $x_1$ , $\varphi_1(s)^{2} \sim c(s - x_1)^{2}$ , so the inner integral diverges like $1/(x - x_1)$ and $\varphi_2 \to$ a finite nonzero limit; near a point where $\varphi_1$ vanishes to higher order (which happens only at a singular point of the equation) the divergence is stronger. This is the mechanism behind the logarithms and inverse powers in the irregular solutions of Legendre's and Bessel's equations.
Continuity of $a_1$ is what makes $e^{-A}$ well defined and continuous; $a_2$ is never used explicitly in the second-order formula, which is striking — all of its influence is already encoded in $\varphi_1$ .
Constants are free, as the books stress. Hassani (Hassani (2013), Thm 14.4.4 and the preceding paragraph, p. 426) writes the general second solution as $\varphi_1\bigl(C + K\int\cdots\bigr)$ and sets $C = 0$ because that term is proportional to $\varphi_1$ ; Arfken, Weber and Harris (Arfken, Weber & Harris (2013), pp. 362-363) drop the lower limits for the same reason and normalise $W = 1$ . We fix $W(x_0) = 1$ by taking the lower limit $x_0$ in both integrals, which is the same choice.
The two routes are the same calculation. (Riley, Hobson & Bence (2006), §16.4.1, p. 545) notes that assuming $y_2 = u(z)y_1(z)$ and the Wronskian route give identical results; the last paragraph of the proof above is that identity.

Method(Getting a second solution)AdaptedAdapted from Coddington (1989), Ch. 3, §5, pp. 119-121; Arfken, Weber & Harris (2013), §7.6, p. 363

Normalise to $y'' + a_1y' + a_2y = 0$ and fix an interval $I$ on which $\varphi_1$ has no zero.
Compute the Wronskian factor $e^{-A}$ , $A(x) = \int_{x_0}^{x}a_1$ . For $a_1 = c/x$ this is a power of $x$ ; for $a_1$ constant it is an exponential; for $a_1 = 0$ it is $1$ .
Integrate $\varphi_1^{-2}e^{-A}$ once. Partial fractions usually does it when $\varphi_1$ is a polynomial, and this step is where logarithms are born.
Multiply by $\varphi_1$ and simplify, discarding any multiple of $\varphi_1$ that the constants of integration produced.
Check by substitution, or by confirming $W(\varphi_1,\varphi_2) = e^{-A}$ , which is a weaker but much cheaper test.
For order $n > 2$ , substitute $y = u\varphi_1$ directly, set $v = u'$ , and solve the resulting equation of order $n-1$ by whatever means; then integrate each of its solutions once and multiply by $\varphi_1$ .

Example(Legendre's equation with

\ell = 1

: the logarithmic second solution)Original

Take $(1 - x^{2})y'' - 2xy' + 2y = 0$ with the known solution $\varphi_1(x) = x$ , and work on $I = (0,1)$ , where $\varphi_1 \ne 0$ .

Normalise: $a_1(x) = -2x/(1-x^{2})$ , $a_2(x) = 2/(1-x^{2})$ .
Wronskian factor. With $x_0 = \tfrac12$ ,
$A(s) = \int_{1/2}^{s}\frac{-2t}{1 - t^{2}}\,dt = \Bigl[\ln(1 - t^{2})\Bigr]_{1/2}^{s} = \ln\frac{1 - s^{2}}{3/4}, \qquad e^{-A(s)} = \frac{3}{4}\cdot\frac{1}{1 - s^{2}} .$
Integrate. By partial fractions, $\dfrac{1}{s^{2}(1-s^{2})} = \dfrac{1}{s^{2}} + \dfrac{1}{1-s^{2}}$ , so
$\int\frac{ds}{s^{2}(1-s^{2})} = -\frac{1}{s} + \frac{1}{2}\ln\frac{1+s}{1-s} .$
Assemble. Up to the constant factor $\tfrac34$ and a multiple of $\varphi_1 = x$ coming from the lower limit,
$\varphi_2(x) \;=\; x\left(-\frac{1}{x} + \frac{1}{2}\ln\frac{1+x}{1-x}\right) \;=\; \frac{x}{2}\ln\frac{1+x}{1-x} - 1 .$
This is the second solution used in Structure of linear equations and the Wronskian, where it is checked by substitution and its Wronskian with $x$ is computed to be $1/(1-x^{2})$ .
The formula's own prediction. $W(\varphi_1,\varphi_2)$ should be $e^{-A(x)} = \tfrac34(1-x^{2})^{-1}$ for the normalised $\varphi_2$ ; our $\varphi_2$ has had the factor $\tfrac34$ removed, so $W = (1-x^{2})^{-1}$ , as found independently there. The two computations agree.
Beyond $(0,1)$ . The resulting $\varphi_2$ in fact solves the equation on all of $(-1,1)$ , as direct substitution shows, even though the derivation needed $x \ne 0$ ; at $x = 0$ the inner integral diverges but the product with $\varphi_1 = x$ stays finite, with $\varphi_2(0) = -1$ . At $x = \pm1$ the equation is singular and $\varphi_2$ blows up logarithmically — the generic fate of the second solution at a singular point.

Example(A polynomial second solution)Original

Take

x\,y'' - (x+2)\,y' + 2y = 0 \qquad \text{on } I = (0,\infty) .

A solution by inspection. Try $\varphi_1 = e^{x}$ : then $xe^{x} - (x+2)e^{x} + 2e^{x} = 0$ . It never vanishes, so the hypothesis holds on all of $I$ — indeed on $\mathbb{R}$ , but the equation is singular at $x = 0$ , so we stay on $x > 0$ .
Normalise: $y'' - \dfrac{x+2}{x}y' + \dfrac{2}{x}y = 0$ , so $a_1(x) = -1 - \dfrac{2}{x}$ .
Wronskian factor. With $x_0 = 1$ ,
$A(s) = \int_{1}^{s}\left(-1 - \frac{2}{t}\right)dt = -(s - 1) - 2\ln s , \qquad e^{-A(s)} = e^{\,s-1}s^{2} .$
Integrate.
$\frac{e^{-A(s)}}{\varphi_1(s)^{2}} = \frac{e^{\,s-1}s^{2}}{e^{2s}} = e^{-1}s^{2}e^{-s} , \qquad \int_{1}^{x} s^{2}e^{-s}\,ds = \Bigl[-(s^{2}+2s+2)e^{-s}\Bigr]_{1}^{x} = 5e^{-1} - (x^{2}+2x+2)e^{-x} .$
Assemble.
$\varphi_2(x) = e^{x}\cdot e^{-1}\Bigl[5e^{-1} - (x^{2}+2x+2)e^{-x}\Bigr] = 5e^{\,x-2} - e^{-1}\bigl(x^{2}+2x+2\bigr) .$
The first term is a multiple of $\varphi_1$ , so discard it, and drop the constant $-e^{-1}$ : a second solution is
$\varphi_2(x) = x^{2} + 2x + 2 .$
Checks. Substituting: $x\cdot 2 - (x+2)(2x+2) + 2(x^{2}+2x+2) = 2x - 2x^{2} - 6x - 4 + 2x^{2} + 4x + 4 = 0$ . And $W(e^{x}, x^{2}+2x+2) = e^{x}(2x+2) - e^{x}(x^{2}+2x+2) = -x^{2}e^{x}$ , which is a constant multiple of $e^{-A(x)} = x^{2}e^{\,x-1}$ as Abel's formula requires. The general solution on $x > 0$ is $y = c_1e^{x} + c_2(x^{2}+2x+2)$ .
A remark worth the attention. An equation with a non-polynomial solution can have a polynomial one as well, and here the polynomial is the second solution, produced by an integration rather than by a series. Note also that neither solution is a polynomial times the other, and that the Wronskian $-x^{2}e^{x}$ vanishes at the singular point $x = 0$ , exactly as Structure of linear equations and the Wronskian predicts.

ApplicationOriginal

Second solutions of the equations of physics. The standard special functions come in pairs, and the second member of each pair is produced this way. Legendre's equation gives $P_\ell$ by a series and then $Q_\ell$ by exactly the computation of the first example above, with its $\ln\frac{1+x}{1-x}$ ; Bessel's equation gives $J_n$ and then $Y_n$ , whose $\ln x$ has the same origin. In a physical problem the irregular solution is usually the one discarded by a boundary condition — but one must know it to know that it has been discarded, and in exterior or scattering problems it is the one that is kept.

The repeated-root case in the Frobenius method. If the two exponents at a regular singular point coincide, or if their difference is a whole number, the series method yields only one solution (The Frobenius method). Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §16.4.1, pp. 544-545) use precisely this formula to manufacture the missing one; it is also the cleanest explanation of why a logarithm must appear, since $\int \varphi_1^{-2}e^{-A}$ produces $\int ds/s$ exactly when the exponents are equal or differ by an integer.

Critical damping. For $\ddot x + 2\gamma\dot x + \gamma^{2}x = 0$ the characteristic polynomial has the double root $-\gamma$ , so inspection gives only $\varphi_1 = e^{-\gamma t}$ . Here $a_1 = 2\gamma$ , $e^{-A} = e^{-2\gamma t}$ and $\varphi_1^{-2} = e^{2\gamma t}$ , so the integrand is $1$ and $\varphi_2 = t\,e^{-\gamma t}$ . The factor $t$ that appears in every repeated-root formula is this trivial integral; the critically damped oscillator of The damped and driven harmonic oscillator is its most familiar instance.

Numerics. When a solution has been computed numerically on a grid, this formula gives the second one by quadrature, which is far more stable than integrating the equation a second time: a numerical solution of a second-order equation started near the "wrong" solution is contaminated by the dominant one, whereas the quadrature has no such instability as long as $\varphi_1$ stays away from zero.

Method

Euler-Cauchy equations

Adapted

Add to PDF

An equation in which each term carries $x^{k}$ with the $k$ -th derivative is scale invariant, so powers $|x|^{r}$ solve it; the exponents are the roots of an indicial polynomial, a root of multiplicity $m$ contributes $|x|^{r}$ times $1, \log|x|, \dots, (\log|x|)^{m-1}$ , and the substitution $x = e^{t}$ turns the equation into one with constant coefficients.

Prerequisites: Constant-coefficient equations

Definition(Euler-Cauchy equations)AdaptedAdapted from Coddington (1989), Ch. 4, §2, Eqs. (2.1) and (2.5), pp. 145, 148; Riley, Hobson & Bence (2006), §15.2.1, Eq. (15.38), p. 504

An Euler equation (also Euler-Cauchy, or equidimensional) of order $n$ is

L(y) \;=\; x^{n}y^{(n)} + a_1x^{n-1}y^{(n-1)} + \cdots + a_{n-1}xy' + a_ny \;=\; f(x) ,

with $a_1, \dots, a_n$ constants. The defining feature is that the $k$ -th derivative appears multiplied by $x^{k}$ . The consequence is a scaling symmetry: for any $c > 0$ , if $y(x)$ solves $L(y) = 0$ then so does $y(cx)$ , because each term $x^{k}y^{(k)}$ is unchanged in form when $x$ is rescaled. This is the reason for the name equidimensional: every term has the same dimension if $x$ and $y$ are given dimensions.

The leading coefficient $x^{n}$ vanishes at $x = 0$ , so $x = 0$ is a singular point and the theory of Structure of linear equations and the Wronskian applies only on intervals not containing it. In this unit $I$ is any interval with $0 \notin I$ ; the two cases $x > 0$ and $x < 0$ are separate problems with identical answers.

For $r \in \mathbb{C}$ and $x \neq 0$ define $|x|^{r} = e^{r\log|x|}$ , with $\log$ the real logarithm of the positive number $|x|$ . The second-order case,

x^{2}y'' + a\,xy' + b\,y = 0 ,

is the one that occurs most often; its indicial polynomial is $q(r) = r(r-1) + ar + b$ .

Method(The indicial polynomial)AdaptedAdapted from Coddington (1989), Ch. 4, §2, Eqs. (2.2) and (2.6), pp. 145, 148

For every $r \in \mathbb{C}$ and every $x \ne 0$ ,

L\bigl(|x|^{r}\bigr) \;=\; q(r)\,|x|^{r} ,

where $q$ , the indicial polynomial of $L$ , is the polynomial of degree $n$

q(r) \;=\; r(r-1)\cdots(r-n+1) \;+\; a_1\,r(r-1)\cdots(r-n+2) \;+\; \cdots \;+\; a_{n-1}r \;+\; a_n .

For $n = 2$ this is $q(r) = r(r-1) + ar + b$ , and $q(r) = 0$ is called the indicial equation. Equivalently: $q$ is obtained from $L$ by replacing $x^{k}d^{k}/dx^{k}$ with the falling factorial $r(r-1)\cdots(r-k+1)$ .

Method(A basis of solutions)AdaptedAdapted from Coddington (1989), Ch. 4, §2, Thm 1, p. 147 and Thm 2, p. 149; Riley, Hobson & Bence (2006), §15.2.1, p. 504

Write the indicial polynomial $q$ of the $n$ -th order Euler equation as $q(r) = (r-r_1)^{m_1}\cdots(r-r_s)^{m_s}$ with $r_1, \dots, r_s$ distinct, so that $m_1 + \cdots + m_s = n$ . Then the $n$ functions

|x|^{r_j},\quad |x|^{r_j}\log|x|,\quad \dots,\quad |x|^{r_j}\bigl(\log|x|\bigr)^{m_j-1}, \qquad j = 1, \dots, s ,

are a basis for the solution space of $L(y) = 0$ on every interval that avoids $x = 0$ .

In particular, for the second-order equation $x^{2}y'' + axy' + by = 0$ : if $q$ has distinct roots $r_1 \ne r_2$ a basis is $|x|^{r_1}, |x|^{r_2}$ ; if $r_1$ is a double root a basis is $|x|^{r_1}, |x|^{r_1}\log|x|$ .

Method(Real form for a complex pair)AdaptedAdapted from Coddington (1989), Ch. 4, §2, Ex. 4, pp. 149-150

Suppose $a_1, \dots, a_n$ are real, so $q$ has real coefficients, and let $r = \sigma + i\tau$ with $\tau \neq 0$ be a root of multiplicity $m$ ; then $\bar r = \sigma - i\tau$ is a root of the same multiplicity. For each $k = 0, \dots, m-1$ the two complex solutions $|x|^{\sigma\pm i\tau}(\log|x|)^{k}$ may be replaced by the real pair

|x|^{\sigma}\bigl(\log|x|\bigr)^{k}\cos\bigl(\tau\log|x|\bigr), \qquad |x|^{\sigma}\bigl(\log|x|\bigr)^{k}\sin\bigl(\tau\log|x|\bigr),

and the result is again a basis. For $n = 2$ with $r = \sigma \pm i\tau$ this gives the basis $|x|^{\sigma}\cos(\tau\log|x|)$ , $|x|^{\sigma}\sin(\tau\log|x|)$ : a power law modulated by an oscillation that is periodic in $\log|x|$ , not in $x$ .

Method(Reduction to constant coefficients)AdaptedAdapted from Riley, Hobson & Bence (2006), §15.2.1, Eqs. (15.36)-(15.38), pp. 503-504; Coddington (1989), Ch. 4, §2, Ex. 3, p. 149

On $x > 0$ substitute $x = e^{t}$ , that is, $t = \log x$ , and set $Y(t) = y(e^{t})$ . Writing $\theta = d/dt$ ,

x^{k}\frac{d^{k}y}{dx^{k}} \;=\; \theta(\theta - 1)\cdots(\theta - k + 1)\,Y , \qquad k = 0, 1, 2, \dots

Hence $L(y) = f(x)$ becomes the constant-coefficient equation

\Bigl[\theta(\theta-1)\cdots(\theta-n+1) + a_1\theta(\theta-1)\cdots(\theta-n+2) + \cdots + a_n\Bigr]Y \;=\; f(e^{t}) ,

whose characteristic polynomial is exactly the indicial polynomial $q$ . The methods of Constant-coefficient equations apply, and translating back with $t = \log x$ , $e^{rt} = x^{r}$ , $t^{k} = (\log x)^{k}$ reproduces the basis above.

More generally, Riley, Hobson and Bence (Riley, Hobson & Bence (2006), Eq. (15.36), p. 503) call

A_n(\alpha x + \beta)^{n}y^{(n)} + \cdots + A_1(\alpha x + \beta)y' + A_0y = f(x)

Legendre's linear equation and reduce it the same way, by $\alpha x + \beta = e^{t}$ .

Idea of the proofOriginal

A scale-invariant operator should have power-law eigenfunctions, and it does: $x^{k}d^{k}/dx^{k}$ applied to $x^{r}$ returns a number times $x^{r}$ , so $L(x^{r}) = q(r)x^{r}$ and the exponents we want are the roots of $q$ . The resemblance to the constant-coefficient case — where $L(e^{\lambda x}) = p(\lambda)e^{\lambda x}$ — is not an accident: $\log x$ converts scalings of $x$ into translations of $t$ , and a translation-invariant operator is one with constant coefficients. That one change of variable explains everything, including where the logarithms come from: a repeated root of a constant-coefficient equation contributes $t^{k}e^{\lambda t}$ , and $t = \log x$ . One can also stay in $x$ and differentiate the identity $L(|x|^{r}) = q(r)|x|^{r}$ with respect to $r$ , which produces $\log|x|$ on the left and the derivatives $q^{(j)}(r)$ on the right — the same bookkeeping as for repeated roots in the constant-coefficient case.

Proof(The indicial polynomial)AdaptedAdapted from Coddington (1989), Ch. 4, §2, pp. 145, 148

Step 1 ( $x > 0$ ). For $x > 0$ , $|x|^{r} = x^{r} = e^{r\log x}$ and $\frac{d}{dx}x^{r} = rx^{r-1}$ , which iterates to

\frac{d^{k}}{dx^{k}}x^{r} = r(r-1)\cdots(r-k+1)\,x^{r-k} , \qquad\text{so}\qquad x^{k}\frac{d^{k}}{dx^{k}}x^{r} = r(r-1)\cdots(r-k+1)\,x^{r} .

Step 2 (sum the terms). The $j$ -th term of $L$ is $a_j x^{n-j}d^{n-j}/dx^{n-j}$ (with $a_0 = 1$ ), so applying it to $x^{r}$ gives $a_j\,r(r-1)\cdots(r-n+j+1)x^{r}$ . Adding over $j = 0, \dots, n$ gives $L(x^{r}) = q(r)x^{r}$ with $q$ as stated.

Step 3 ( $x < 0$ ). For $x < 0$ , $|x| = -x$ , so $\frac{d}{dx}|x|^{r} = \frac{d}{dx}(-x)^{r} = -r(-x)^{r-1}$ and more generally $\frac{d^{k}}{dx^{k}}|x|^{r} = (-1)^{k}r(r-1)\cdots(r-k+1)(-x)^{r-k}$ . Multiplying by $x^{k} = (-1)^{k}(-x)^{k}$ the two sign factors cancel, and again

x^{k}\frac{d^{k}}{dx^{k}}|x|^{r} = r(r-1)\cdots(r-k+1)\,|x|^{r} .

So the identity, and hence everything derived from it, holds verbatim on either side of the origin. $\square$

Proof(A basis of solutions)AdaptedAdapted from Coddington (1989), Ch. 4, §2, pp. 145-149; Riley, Hobson & Bence (2006), §15.2.1, pp. 503-504

Step 1 (change variables). Work on $x > 0$ ; the case $x < 0$ follows by replacing $x$ with $|x|$ throughout, by Step 3 of the previous proof. Put $t = \log x$ , $Y(t) = y(e^{t})$ . From $dt/dx = 1/x$ we get $x\,dy/dx = dY/dt = \theta Y$ . Assume inductively that $x^{k}y^{(k)} = \theta(\theta-1)\cdots(\theta-k+1)Y$ and write $G = \theta(\theta-1)\cdots(\theta-k+1)Y$ for brevity. Differentiating $x^{k}y^{(k)} = G$ in $x$ ,

kx^{k-1}y^{(k)} + x^{k}y^{(k+1)} = \frac{dG}{dt}\cdot\frac{1}{x} ,

and multiplying by $x$ ,

x^{k+1}y^{(k+1)} = \theta G - k\,x^{k}y^{(k)} = (\theta - k)\,\theta(\theta-1)\cdots(\theta-k+1)Y ,

which is the statement for $k+1$ . The base case $k = 0$ is trivial. This proves the operator identities, hence that $L(y) = f$ becomes the constant-coefficient equation displayed above, with characteristic polynomial

\theta(\theta-1)\cdots(\theta-n+1) + a_1\theta(\theta-1)\cdots(\theta-n+2) + \cdots + a_n \Big|_{\theta \to r} \;=\; q(r) .

Step 2 (solve it). By Constant-coefficient equations, a root $r_j$ of $q$ of multiplicity $m_j$ contributes the $m_j$ solutions $t^{k}e^{r_jt}$ , $k = 0, \dots, m_j - 1$ , and together these $n$ functions are a basis of the solution space in $t$ .

Step 3 (translate back). $t = \log x$ turns $e^{r_jt}$ into $x^{r_j}$ and $t^{k}$ into $(\log x)^{k}$ , so the $n$ functions $x^{r_j}(\log x)^{k}$ solve $L(y) = 0$ on $x > 0$ . The map $y \mapsto Y$ , $Y(t) = y(e^{t})$ , is a linear bijection between functions on $(0,\infty)$ and functions on $\mathbb{R}$ with inverse $y(x) = Y(\log x)$ ; a linear bijection carries a basis to a basis and independent functions to independent functions. Hence these $n$ functions are a basis of the solutions of $L(y) = 0$ on any interval of $(0,\infty)$ .

Step 4 (a check that avoids the substitution). The same conclusion follows from the indicial identity directly. Differentiating $L(|x|^{r}) = q(r)|x|^{r}$ $k$ times with respect to $r$ , and using $\partial_r|x|^{r} = |x|^{r}\log|x|$ together with the commuting of $\partial_r$ and $d/dx$ ,

L\Bigl(|x|^{r}\bigl(\log|x|\bigr)^{k}\Bigr) \;=\; \sum_{i=0}^{k}\binom{k}{i}q^{(i)}(r)\,|x|^{r}\bigl(\log|x|\bigr)^{k-i} .

If $r_j$ has multiplicity $m_j$ then $q(r_j) = \cdots = q^{(m_j-1)}(r_j) = 0$ , so the right-hand side vanishes for every $k \le m_j - 1$ . $\square$

Proof(Real form for a complex pair)Original

Since $q$ has real coefficients, $q(\bar r) = \overline{q(r)}$ , so $\bar r$ is a root whenever $r$ is, with the same multiplicity (apply the same conjugation to $q^{(i)}$ ). With $r = \sigma + i\tau$ and $u = \log|x|$ ,

|x|^{\sigma\pm i\tau} = e^{(\sigma\pm i\tau)u} = |x|^{\sigma}\bigl(\cos\tau u \pm i\sin\tau u\bigr) ,

|x|^{\sigma}\cos(\tau u) = \tfrac12\bigl(|x|^{r} + |x|^{\bar r}\bigr), \qquad |x|^{\sigma}\sin(\tau u) = \tfrac{1}{2i}\bigl(|x|^{r} - |x|^{\bar r}\bigr) ,

and multiplying both identities by $u^{k}$ gives the claimed pair as combinations of the two complex solutions. The $2\times 2$ matrix $\begin{pmatrix}1/2 & 1/2\\ 1/(2i) & -1/(2i)\end{pmatrix}$ has determinant $-1/(2i) \neq 0$ , so the replacement is invertible and the new $n$ functions span the same space; being $n$ spanning vectors of an $n$ -dimensional space they are a basis. $\square$

Remark(Where the hypotheses are used)Original

$x \ne 0$ throughout. The indicial identity $x^{k}(|x|^{r})^{(k)} = r(r-1)\cdots(r-k+1)|x|^{r}$ is derived in Steps 1 and 3 for $x \neq 0$ only, and that is all one can ask: the leading coefficient $x^{n}$ vanishes at the origin, so the equation is singular there, and the functions themselves misbehave — $|x|^{r}$ with $r$ non-integer is not differentiable enough times at $0$ (for $0 < r < 1$ not even once), and $\log|x|$ is not defined there at all. Coddington's theorems are stated, carefully, on any interval not containing the origin; a solution on $(0,\infty)$ and one on $(-\infty,0)$ have no reason to join up. For $x^{2}y'' - 2y = 0$ ( $q(r) = r^{2} - r - 2$ , roots $2$ and $-1$ ) the solution $|x|^{-1}$ exists on each side and blows up at the origin, and $x^{2}$ extends smoothly but is not part of a two-dimensional solution space at $0$ : compare the warning in Structure of linear equations and the Wronskian.
$|x|^{r}$ rather than $x^{r}$ . On $x < 0$ and for complex or non-integer $r$ , $x^{r}$ needs a branch of $\log$ and becomes complex-valued ((Coddington (1989), Ch. 4, §2, Ex. 5, p. 150) works this out), whereas $|x|^{r}$ is real for real $r$ and is what the identity of Step 3 delivers. On $x > 0$ the two agree and we write $x^{r}$ and $\log x$ .
Constant $a_k$ are essential. If one $a_k$ is replaced by a function the scaling symmetry is lost, $L(x^{r})$ is no longer a multiple of $x^{r}$ , and only the leading behaviour near $0$ survives — which is precisely the content of the Frobenius method (The Frobenius method): near a regular singular point the solutions are the Euler solutions multiplied by convergent power series. Strauss (Strauss (2008), Appendix A.4, p. 424) states this as the main theorem about regular singular points, with the same indicial equation $x(x-1) + \beta x + \gamma = 0$ .
Multiplicity, not root separation, produces the logarithm. For $n = 2$ the double root occurs when $(a-1)^{2} = 4b$ ; there $|x|^{r_1}$ and $|x|^{r_1}\log|x|$ are the basis. The logarithm can also be obtained from reduction of order: with $\varphi_1 = x^{r_1}$ and $a_1 = a/x$ one gets $e^{-A} = x^{-a}$ and $\varphi_1^{-2}e^{-A} = x^{-2r_1-a} = x^{-1}$ exactly when $2r_1 + a = 1$ , which is the condition $r_1 = (1-a)/2$ for a double root. The integral $\int dx/x$ is the logarithm.
Forcing by a power is easy. Since $L(cx^{k}) = c\,q(k)x^{k}$ , the equation $L(y) = x^{k}$ has the particular solution $x^{k}/q(k)$ whenever $q(k) \ne 0$ ; if $k$ is a root of $q$ of multiplicity $m$ , the trial $c\,x^{k}(\log x)^{m}$ works, with $c = 1/q^{(m)}(k)$ . The reason is the differentiated identity of Step 4 above with the exponent $m$ : in $L\Bigl(x^{k}(\log x)^{m}\Bigr) = \sum_{i=0}^{m}\binom{m}{i}q^{(i)}(k)\,x^{k}(\log x)^{m-i}$ every term with $i \le m-1$ carries a factor $q^{(i)}(k) = 0$ , and the surviving $i = m$ term is $\binom{m}{m}q^{(m)}(k)x^{k} = q^{(m)}(k)\,x^{k}$ . So $L\bigl(x^{k}(\log x)^{m}\bigr) = q^{(m)}(k)x^{k}$ exactly, and $q^{(m)}(k) \ne 0$ because the multiplicity is exactly $m$ . (For $m = 1$ this reads $c = 1/q'(k)$ ; for the triple root of the first example, $q = (r-1)^{3}$ and $q'''(1) = 6$ , so forcing $x$ would need $y_p = \tfrac16 x(\log x)^{3}$ , and indeed $x^{3}y''' + xy' - y$ applied to $x(\log x)^{3}$ gives $6x$ .) This is the exact analogue of resonance in Constant-coefficient equations, and (Coddington (1989), Ch. 4, §2, Ex. 6, p. 150) sets it as an exercise for $n=2$ .

Method(Solving an Euler equation)AdaptedAdapted from Riley, Hobson & Bence (2006), §15.2.1, Solution method, pp. 504-505; Coddington (1989), Ch. 4, §2, pp. 145-149

Recognise it: each $y^{(k)}$ multiplied by $x^{k}$ (times a constant). If the variable appears as $\alpha x + \beta$ throughout, substitute $u = \alpha x + \beta$ first.
Write the indicial polynomial by replacing $x^{k}y^{(k)}$ with $r(r-1)\cdots(r-k+1)$ , and find its roots with multiplicities. For $n = 2$ : $q(r) = r(r-1) + ar + b$ .
Write the basis: each root $r_j$ of multiplicity $m_j$ contributes $|x|^{r_j}\bigl(\log|x|\bigr)^{k}$ , $k = 0, \dots, m_j-1$ . Convert a complex pair to the real $|x|^{\sigma}\cos(\tau\log|x|)$ , $|x|^{\sigma}\sin(\tau\log|x|)$ if a real answer is wanted.
For a nonzero right-hand side, either (a) try $y_p = c\,x^{k}$ for forcing $x^{k}$ , using $L(cx^{k}) = c\,q(k)x^{k}$ , with a factor $(\log x)^{m}$ if $k$ is a root of multiplicity $m$ ; or (b) substitute $x = e^{t}$ and use undetermined coefficients or variation of parameters on the resulting constant-coefficient equation; or (c) use variation of parameters directly. For $n = 2$ the Wronskian needed there is, for distinct roots $r_1 \ne r_2$ with basis $|x|^{r_1}, |x|^{r_2}$ , $W = (r_2 - r_1)\,|x|^{r_1 + r_2 - 1}\operatorname{sgn} x ,$ which is nonvanishing exactly because $r_1 \ne r_2$ ; and for a double root $r_1$ , with basis $|x|^{r_1}, |x|^{r_1}\log|x|$ , $W = |x|^{2r_1 - 1}\operatorname{sgn} x .$
State the interval. The answer is valid on an interval that excludes $x = 0$ ; say which one.

Example(A triple root, with forcing)Original

Solve

x^{3}y''' + x\,y' - y \;=\; x^{2} \quad\text{on } x > 0 , \qquad y(1) = y'(1) = y''(1) = 0 .

Indicial polynomial. Replace $x^{3}y'''$ by $r(r-1)(r-2)$ , $xy'$ by $r$ , and $y$ by $1$ :
$q(r) = r(r-1)(r-2) + r - 1 = r^{3} - 3r^{2} + 2r + r - 1 = r^{3} - 3r^{2} + 3r - 1 = (r-1)^{3} .$
So $r = 1$ is a root of multiplicity $3$ .
Complementary function. $y_c = \bigl(c_1 + c_2\log x + c_3(\log x)^{2}\bigr)x$ . Each piece checks out directly; for instance $y = x(\log x)^{2}$ has $y' = (\log x)^{2} + 2\log x$ , $y'' = (2\log x + 2)/x$ , $y''' = -2\log x/x^{2}$ , and
$x^{3}\left(\frac{-2\log x}{x^{2}}\right) + x\bigl((\log x)^{2} + 2\log x\bigr) - x(\log x)^{2} = -2x\log x + 2x\log x = 0 .$
Their Wronskian at $x = 1$ is the determinant of the matrix with rows $(1,0,0)$ , $(1,1,0)$ , $(0,1,2)$ , namely $2 \neq 0$ ; and since the normalised equation has no $y''$ term, Abel's formula makes $W$ constant, so $W \equiv 2$ on $x>0$ .
Particular integral. The forcing is $x^{2}$ and $q(2) = (2-1)^{3} = 1 \ne 0$ , so $y_p = x^{2}/q(2) = x^{2}$ . Check: $x^{3}\cdot 0 + x\cdot 2x - x^{2} = x^{2}$ .
General solution. $y = \bigl(c_1 + c_2\log x + c_3(\log x)^{2}\bigr)x + x^{2}$ .
Impose the data at $x = 1$ , where $\log x = 0$ . With $y_h = x(c_1 + c_2L + c_3L^{2})$ , $L = \log x$ :
$y_h' = c_1 + c_2 + (c_2 + 2c_3)L + c_3L^{2} , \qquad y_h'' = \frac{c_2 + 2c_3 + 2c_3L}{x} ,$
so $y_h(1) = c_1$ , $y_h'(1) = c_1 + c_2$ , $y_h''(1) = c_2 + 2c_3$ . Since $y_p(1) = 1$ , $y_p'(1) = 2$ , $y_p''(1) = 2$ , the conditions give $c_1 = -1$ , then $c_2 = -1$ , then $c_3 = -\tfrac12$ :
$y(x) \;=\; x^{2} - x\left(1 + \log x + \tfrac12(\log x)^{2}\right) .$
Checks. Substituting this into the equation returns $x^{2}$ , and $y, y', y''$ vanish at $x = 1$ . At $x = 2$ : $y(2) = 4 - 2\bigl(1 + \ln 2 + \tfrac12\ln^{2}2\bigr) = 0.133\,253$ to six decimal places.

Example(A complex pair: oscillation in

\log x

)Original

Solve $x^{2}y'' + 3xy' + 5y = 0$ for $x \ne 0$ .

Indicial equation. $q(r) = r(r-1) + 3r + 5 = r^{2} + 2r + 5$ , with roots $r = -1 \pm 2i$ (so $\sigma = -1$ , $\tau = 2$ ).
Real basis.
$\varphi_1(x) = \frac{\cos\bigl(2\log|x|\bigr)}{|x|} , \qquad \varphi_2(x) = \frac{\sin\bigl(2\log|x|\bigr)}{|x|} .$
Check $\varphi_1$ on $x > 0$ . With $u = \log x$ , so that $du/dx = 1/x$ ,
$\varphi_1 = \frac{\cos 2u}{x}, \qquad \varphi_1' = -\frac{\cos 2u + 2\sin 2u}{x^{2}}, \qquad \varphi_1'' = \frac{2\bigl(3\sin 2u - \cos 2u\bigr)}{x^{3}} ,$
so that
$x^{2}\varphi_1'' + 3x\varphi_1' + 5\varphi_1 = \frac{6\sin 2u - 2\cos 2u - 3\cos 2u - 6\sin 2u + 5\cos 2u}{x} = 0 .$
The computation for $\varphi_2$ is the same with $\cos$ and $\sin$ exchanged and a sign flipped.
What the solution looks like. As $x \to 0^{+}$ the amplitude $1/x$ grows while the phase $2\log x \to -\infty$ : the solution oscillates infinitely often in any neighbourhood of the origin, with the zeros of $\varphi_1$ at $x = e^{\pi/4 + k\pi/2}$ , $k \in \mathbb{Z}$ , accumulating geometrically at the origin. The oscillation is periodic in $\log x$ with period $\pi$ , so consecutive zeros are in the fixed ratio $e^{-\pi/2} = 0.2079$ as one moves towards the origin. This log-periodic behaviour is the signature of a complex pair of indicial roots.
A definite solution. The solution with $y(1) = 1$ , $y'(1) = 0$ is $y = x^{-1}\bigl(\cos(2\log x) + \tfrac12\sin(2\log x)\bigr)$ , since $\varphi_1(1) = 1$ , $\varphi_1'(1) = -1$ , $\varphi_2(1) = 0$ , $\varphi_2'(1) = 2$ . At $x = 2$ : $\tfrac12\bigl(\cos(2\ln 2) + \tfrac12\sin(2\ln 2)\bigr) = 0.337\,485$ to six decimal places.

ApplicationOriginal

Laplace's equation in polar and spherical coordinates. Separating $\nabla^{2}u = 0$ in the plane as $u = R(r)\Theta(\theta)$ gives $r^{2}R'' + rR' - m^{2}R = 0$ — an Euler equation with $q(\rho) = \rho^{2} - m^{2}$ , hence $R = r^{m}, r^{-m}$ for $m \ne 0$ and the double root $m = 0$ giving $R = 1, \log r$ . The logarithmic potential of a line charge is the double-root case of this unit. In three dimensions the radial equation $r^{2}R'' + 2rR' - \ell(\ell+1)R = 0$ has $q(\rho) = \rho(\rho+1) - \ell(\ell+1)$ , with roots $\ell$ and $-\ell-1$ : the interior and exterior multipole solutions $r^{\ell}$ and $r^{-\ell-1}$ of electrostatics are nothing but the exponents of an Euler equation.

Scale-free problems and power laws. Whenever a linear problem has no intrinsic length scale, the governing equation is equidimensional and the solutions are power laws, possibly with logarithms. Examples: the stress field near a crack tip or a wedge corner in elasticity (the Williams expansion, whose exponents solve a transcendental version of an indicial equation), the velocity profile in a corner flow, and the self-similar solutions of diffusion problems. The appearance of a logarithm always signals a degenerate (repeated) exponent, and it is the mathematical origin of such anomalies as the logarithmic corrections in two-dimensional problems.

The inverse-square potential. The radial Schrodinger equation with $V = -g/r^{2}$ reduces, for the leading behaviour, to $r^{2}u'' + \bigl(g - \ell(\ell+1)\bigr)u = 0$ , an Euler equation with $q(\rho) = \rho(\rho-1) + g - \ell(\ell+1)$ . When $g$ exceeds a critical value the indicial roots become a complex pair and the solution becomes log-periodic, $r^{1/2}\cos(\tau\log r + \delta)$ : this is the celebrated "fall to the centre", where no ground state exists and the spectrum must be regularised. The same mathematics governs the Efimov effect in three-body physics, where the log-periodicity appears as a geometric sequence of bound states with ratio $e^{\pi/\tau}$ — the ratio seen in step 3 of the second example.

The model for regular singular points. Finally, the Euler equation is the reference case for the whole of the series-solution theory: near a regular singular point the coefficients $xp(x)$ and $x^{2}q(x)$ tend to constants, and the equation is asymptotically an Euler equation whose indicial roots give the leading exponents of the two solutions (The Frobenius method, Classifying singular points, including infinity). Every indicial equation in the theory of Bessel, Legendre and hypergeometric functions is the indicial equation of an Euler equation.

Application

The damped and driven harmonic oscillator

Adapted

Add to PDF

The equation $\ddot x + 2\gamma\dot x + \omega_0^2 x = f(t)$ worked out in full: the three free regimes fixed by the sign of $\gamma - \omega_0$ , the decaying transient plus periodic steady state under harmonic driving, the amplitude and phase as functions of frequency, and the quality factor $Q = \omega_0/2\gamma$ read three ways.

Prerequisites: Constant-coefficient equations, Structure of linear equations and the Wronskian

Definition(The oscillator and its parameters)AdaptedAdapted from Hirsch, Smale & Devaney (2013), §2.2, pp. 25-26; Hassani (2013), §14.7, Example 14.7.2, pp. 447-448

A mass $m$ on a spring of stiffness $k > 0$ , with a frictional force proportional to the velocity through a damping constant $b \ge 0$ and an external force $F(t)$ , obeys Newton's law in the form

m\ddot x + b\dot x + kx \;=\; F(t) .

Dividing by $m$ and naming

\omega_0 = \sqrt{\frac{k}{m}} \;>\; 0 , \qquad \gamma = \frac{b}{2m} \;\ge\; 0 , \qquad f(t) = \frac{F(t)}{m} ,

gives the standard form used throughout this unit,

\ddot x + 2\gamma\dot x + \omega_0^{2}x \;=\; f(t) .

Here $\omega_0$ is the undamped natural frequency and $\gamma$ the damping rate (of dimension $1/\text{time}$ ). The characteristic polynomial is $p(\lambda) = \lambda^{2} + 2\gamma\lambda + \omega_0^{2}$ , with roots

\lambda_{\pm} \;=\; -\gamma \pm \sqrt{\gamma^{2} - \omega_0^{2}} .

The series $RLC$ circuit $L\ddot q + R\dot q + q/C = V(t)$ is the same equation: Hassani (Hassani (2013), p. 448) records the correspondence $2\gamma = R/L$ and $\omega_0^{2} = 1/(LC)$ . Everything below therefore applies verbatim to circuits, with charge for displacement and voltage for force.

Application(Free motion: three regimes)AdaptedAdapted from Hassani (2013), §14.7, Example 14.7.2, pp. 447-448; Arfken, Weber & Harris (2013), §20.8, Example 20.8.5, Eq. (20.161), p. 1021 and Ex. 20.8.6, p. 1030

Let $f \equiv 0$ and $\gamma, \omega_0 > 0$ . The general real solution of $\ddot x + 2\gamma\dot x + \omega_0^{2}x = 0$ is determined by the sign of $\gamma^{2} - \omega_0^{2}$ .

Underdamped, $\gamma < \omega_0$ . The roots are the complex pair $-\gamma \pm i\omega_d$ with
$\omega_d = \sqrt{\omega_0^{2} - \gamma^{2}} ,$
and
$x(t) = Ae^{-\gamma t}\cos(\omega_dt - \phi) ,$
with $A \ge 0$ and $\phi$ fixed by the data. The motion oscillates at the damped frequency $\omega_d < \omega_0$ inside a decaying envelope $Ae^{-\gamma t}$ , and has infinitely many zeros. For the data $x(0) = x_0$ , $\dot x(0) = 0$ ,
$x(t) = x_0\,\frac{\omega_0}{\omega_d}\,e^{-\gamma t}\cos(\omega_dt - \phi), \qquad \tan\phi = \frac{\gamma}{\omega_d} .$
Critically damped, $\gamma = \omega_0$ . The root $-\gamma$ is double and
$x(t) = (c_1 + c_2t)\,e^{-\gamma t} .$
For $x(0) = x_0$ , $\dot x(0) = 0$ this is $x(t) = x_0(1 + \gamma t)e^{-\gamma t}$ . The motion crosses zero at most once.
Overdamped, $\gamma > \omega_0$ . Both roots
$\lambda_{\pm} = -\gamma \pm \sqrt{\gamma^{2} - \omega_0^{2}}$
are real and negative, with $\lambda_- < -\gamma < \lambda_+ < 0$ , and
$x(t) = c_+e^{\lambda_+t} + c_-e^{\lambda_-t} .$
The motion crosses zero at most once, and for large $t$ decays like $e^{\lambda_+t}$ , that is, more slowly than in the critical case.

In all three cases every solution tends to $0$ as $t \to \infty$ , since both roots have negative real part. In circuit language, Hassani (Hassani (2013), p. 448) states case 3 as: for $R \ge 2\sqrt{L/C}$ the circuit does not oscillate.

Application(Harmonic driving: transient and steady state)AdaptedAdapted from Hirsch, Smale & Devaney (2013), §6.5, pp. 132-133

Let $\gamma > 0$ and $\omega_0 > 0$ , and drive the oscillator harmonically:

\ddot x + 2\gamma\dot x + \omega_0^{2}x \;=\; f_0\cos\omega t , \qquad \omega > 0 .

Steady state. There is exactly one solution that is periodic with the period $2\pi/\omega$ of the driving, namely
$x_{\mathrm{ss}}(t) = A(\omega)\cos\bigl(\omega t - \delta(\omega)\bigr) ,$ $A(\omega) = \frac{f_0}{\sqrt{(\omega_0^{2} - \omega^{2})^{2} + 4\gamma^{2}\omega^{2}}} , \qquad \cos\delta = \frac{\omega_0^{2}-\omega^{2}}{\sqrt{(\omega_0^{2}-\omega^{2})^{2}+4\gamma^{2}\omega^{2}}}, \quad \sin\delta = \frac{2\gamma\omega}{\sqrt{(\omega_0^{2}-\omega^{2})^{2}+4\gamma^{2}\omega^{2}}} ,$
so that $\delta \in (0,\pi)$ is the phase lag of the response behind the force.
Every solution approaches it. Any solution is $x = x_{\mathrm{ss}} + x_{\mathrm{tr}}$ where the transient $x_{\mathrm{tr}}$ is a solution of the free equation, hence $x(t) - x_{\mathrm{ss}}(t) \to 0$ as $t \to \infty$ . The data affect only the transient.
Limits. As $\omega \to 0^{+}$ , $A(\omega) \to f_0/\omega_0^{2}$ and $\delta(\omega) \to 0$ ; we write $A(0) = f_0/\omega_0^{2}$ for this limit, the static deflection, which is also what the displayed formula for $A$ returns at $\omega = 0$ . (The phase lag itself is defined only for $\omega > 0$ , where $\sin\delta > 0$ puts it strictly inside $(0,\pi)$ ; the values $0$ and $\pi$ are approached but never attained.) As $\omega \to \infty$ , $A(\omega) \sim f_0/\omega^{2}$ and $\delta(\omega) \to \pi$ . At $\omega = \omega_0$ , $A = f_0/(2\gamma\omega_0)$ with $\delta = \pi/2$ exactly.

Definition(Quality factor)Original

For $\gamma > 0$ the quality factor of the oscillator is the dimensionless number

Q \;=\; \frac{\omega_0}{2\gamma} \;=\; \frac{\sqrt{mk}}{b} \;=\; \frac{1}{R}\sqrt{\frac{L}{C}} .

The three regimes are $Q > \tfrac12$ (underdamped), $Q = \tfrac12$ (critical) and $Q < \tfrac12$ (overdamped). Large $Q$ means light damping.

Application(Resonance and three readings of

Q

)Original

With $\gamma > 0$ and $\omega_0 > 0$ :

Amplitude resonance. If $\gamma < \omega_0/\sqrt2$ (that is, $Q > 1/\sqrt2$ ), $A(\omega)$ has a unique maximum at
$\omega_r = \sqrt{\omega_0^{2} - 2\gamma^{2}} , \qquad A(\omega_r) = \frac{f_0}{2\gamma\sqrt{\omega_0^{2}-\gamma^{2}}} = \frac{f_0}{2\gamma\omega_d} .$
If $\gamma \ge \omega_0/\sqrt2$ , $A$ is strictly decreasing on $(0,\infty)$ : there is no resonant peak at all.
Velocity resonance. The velocity amplitude $\omega A(\omega)$ , and hence the time-averaged power $\overline{P} = \tfrac12 m(2\gamma)\,\omega^{2}A^{2}$ absorbed from the force, is maximal exactly at $\omega = \omega_0$ for every $\gamma > 0$ , where $\delta = \pi/2$ .
Amplification. $A(\omega_0) = Q\cdot A(0)$ : at the natural frequency the response is $Q$ times the static deflection.
Bandwidth. The two frequencies at which $\overline{P}$ is half its maximum are $\omega_\pm = \sqrt{\omega_0^{2}+\gamma^{2}} \pm \gamma$ , so the full width at half maximum is exactly
$\Delta\omega = \omega_+ - \omega_- = 2\gamma , \qquad\text{hence}\qquad Q = \frac{\omega_0}{\Delta\omega} .$
Energy decay. In free underdamped motion the energy $E = \tfrac12m\dot x^{2} + \tfrac12kx^{2}$ has envelope proportional to $e^{-2\gamma t}$ , so the fraction lost per cycle is $1 - e^{-4\pi\gamma/\omega_d}$ , and
$2\pi\,\frac{E}{|\Delta E|_{\text{per cycle}}} \;=\; \frac{\omega_d}{2\gamma}\Bigl(1 + O\bigl(\gamma/\omega_d\bigr)\Bigr) \;\longrightarrow\; Q \qquad\text{as } \gamma/\omega_0 \to 0 .$

Idea of the proofOriginal

Two ingredients do all the work. First, the free motion is read straight off the characteristic polynomial $\lambda^{2}+2\gamma\lambda+\omega_0^{2}$ : the discriminant $4(\gamma^{2}-\omega_0^{2})$ changes sign at $\gamma = \omega_0$ , and the three regimes are the three cases of Constant-coefficient equations — complex pair, double root, two real roots. Second, for harmonic driving it is easier to drive with $e^{i\omega t}$ than with $\cos\omega t$ : since $p(i\omega) \ne 0$ whenever $\gamma > 0$ , the complex response is simply $f_0e^{i\omega t}/p(i\omega)$ , and taking real parts gives amplitude $f_0/|p(i\omega)|$ and phase lag $\arg p(i\omega)$ . Every statement about resonance is then a statement about where the modulus of the quadratic $p(i\omega) = \omega_0^{2}-\omega^{2} + 2i\gamma\omega$ is small — and $|p(i\omega)|^{2}$ is a quadratic in $\omega^{2}$ , so its minimum is found by completing a square, with no calculus at all.

Proof(Free motion)AdaptedAdapted from Hassani (2013), §14.7, Example 14.7.2, p. 447

Step 1 (roots). $p(\lambda) = \lambda^{2} + 2\gamma\lambda + \omega_0^{2} = (\lambda+\gamma)^{2} - (\gamma^{2}-\omega_0^{2})$ , giving $\lambda_\pm = -\gamma \pm \sqrt{\gamma^{2}-\omega_0^{2}}$ .

Step 2 (the three cases). By Constant-coefficient equations: if $\gamma < \omega_0$ the square root is $i\omega_d$ with $\omega_d = \sqrt{\omega_0^{2}-\gamma^{2}} > 0$ , and the real basis is $e^{-\gamma t}\cos\omega_dt$ , $e^{-\gamma t}\sin\omega_dt$ , whose general real combination is $Ae^{-\gamma t}\cos(\omega_dt-\phi)$ . If $\gamma = \omega_0$ the root $-\gamma$ is double, giving the basis $e^{-\gamma t}$ , $te^{-\gamma t}$ . If $\gamma > \omega_0$ the two roots are real and distinct.

Step 3 (signs in the overdamped case). Let $s = \sqrt{\gamma^{2}-\omega_0^{2}}$ , so $0 < s < \gamma$ because $\omega_0 > 0$ . Then $\lambda_+ = -\gamma + s < 0$ and $\lambda_- = -\gamma - s < -\gamma < \lambda_+$ . Both roots are negative, so both exponentials decay; the slower one is $e^{\lambda_+t}$ . Since $\lambda_+ = -\omega_0^{2}/(\gamma + s) \to 0$ as $\gamma \to \infty$ , heavy damping makes the return to equilibrium arbitrarily slow.

Step 4 (the stated initial-value solutions). Underdamped with $x(0)=x_0$ , $\dot x(0) = 0$ : writing $x = e^{-\gamma t}(a\cos\omega_dt + c\sin\omega_dt)$ gives $a = x_0$ and $\dot x(0) = -\gamma a + \omega_dc = 0$ , so $c = \gamma x_0/\omega_d$ . Then

a\cos\omega_dt + c\sin\omega_dt = \sqrt{a^{2}+c^{2}}\,\cos(\omega_dt - \phi), \qquad \tan\phi = c/a = \gamma/\omega_d ,

and $\sqrt{a^{2}+c^{2}} = x_0\sqrt{1 + \gamma^{2}/\omega_d^{2}} = x_0\omega_0/\omega_d$ , since $\omega_d^{2}+\gamma^{2} = \omega_0^{2}$ . This is the form Arfken, Weber and Harris (Arfken, Weber & Harris (2013), Eq. (20.161), p. 1021) obtain by Laplace transform. Critically damped: $x = (c_1+c_2t)e^{-\gamma t}$ gives $c_1 = x_0$ and $\dot x(0) = c_2 - \gamma c_1 = 0$ , so $x = x_0(1+\gamma t)e^{-\gamma t}$ .

Step 5 (counting zeros). In case 2, $x(t) = 0$ forces $c_1 + c_2t = 0$ : at most one zero. In case 3, $x(t) = 0$ forces $e^{(\lambda_+-\lambda_-)t} = -c_-/c_+$ , and $t \mapsto e^{(\lambda_+-\lambda_-)t}$ is strictly increasing, so again at most one zero. In case 1 the zeros are where $\omega_dt - \phi$ is an odd multiple of $\pi/2$ : infinitely many, spaced $\pi/\omega_d$ apart. $\square$

Proof(Steady state, amplitude and phase)AdaptedAdapted from Hirsch, Smale & Devaney (2013), §6.5, pp. 132-133

Step 1 (complexify). Consider $\ddot z + 2\gamma\dot z + \omega_0^{2}z = f_0e^{i\omega t}$ ; if $z$ solves it then $x = \operatorname{Re}z$ solves the real equation, because the coefficients are real. Since

p(i\omega) = (i\omega)^{2} + 2\gamma(i\omega) + \omega_0^{2} = \bigl(\omega_0^{2}-\omega^{2}\bigr) + 2i\gamma\omega

has $\operatorname{Im}p(i\omega) = 2\gamma\omega \ne 0$ for $\gamma > 0$ and $\omega > 0$ , we have $p(i\omega) \neq 0$ : the driving is never resonant in the sense of Constant-coefficient equations. Hence, by the exponential-shift identity $L(e^{i\omega t}) = p(i\omega)e^{i\omega t}$ ,

z(t) = \frac{f_0}{p(i\omega)}e^{i\omega t}

is a solution.

Step 2 (read off amplitude and phase). Write $p(i\omega) = |p(i\omega)|e^{i\delta}$ with

|p(i\omega)| = \sqrt{(\omega_0^{2}-\omega^{2})^{2} + 4\gamma^{2}\omega^{2}} , \qquad \cos\delta = \frac{\omega_0^{2}-\omega^{2}}{|p(i\omega)|}, \quad \sin\delta = \frac{2\gamma\omega}{|p(i\omega)|} .

Since $\sin\delta > 0$ , $\delta$ lies in $(0,\pi)$ and is uniquely determined. Then $z = (f_0/|p|)e^{i(\omega t - \delta)}$ and

x_{\mathrm{ss}} = \operatorname{Re}z = A(\omega)\cos(\omega t - \delta), \qquad A(\omega) = \frac{f_0}{|p(i\omega)|} ,

which is the stated formula. It is periodic with period $2\pi/\omega$ .

Step 3 (uniqueness of the periodic solution). If $x_1, x_2$ are both solutions of period $2\pi/\omega$ , then $u = x_1 - x_2$ solves the free equation and is periodic. But every solution of the free equation tends to $0$ (Step 1 of the previous proof: both roots have real part $\le -\min(\gamma, \omega_0^{2}/2\gamma) < 0$ ), and a periodic function tending to $0$ is identically $0$ . So $x_1 = x_2$ . This is the conclusion Hirsch, Smale and Devaney (Hirsch, Smale & Devaney (2013), §6.5, pp. 132-133) reach for the system form of the equation, by showing that $e^{2\pi A}-I$ is invertible.

Step 4 (transients die). Any solution $x$ has $x - x_{\mathrm{ss}}$ solving the free equation, hence decaying to $0$ ; so $x$ approaches $x_{\mathrm{ss}}$ , and the initial data influence only the transient part.

Step 5 (limits). At $\omega = 0$ : $|p| = \omega_0^{2}$ , $\delta = 0$ . As $\omega \to \infty$ : $|p| \sim \omega^{2}$ and $\cos\delta \to -1$ , so $\delta \to \pi$ — the response is in antiphase with the force. At $\omega = \omega_0$ : $\cos\delta = 0$ and $\sin\delta = 1$ , so $\delta = \pi/2$ and $|p| = 2\gamma\omega_0$ . $\square$

Proof(Resonance and the readings of

Q

)Original

Step 1 (amplitude peak). $A$ is largest where $|p(i\omega)|^{2}$ is smallest. Put $v = \omega^{2} \ge 0$ :

|p|^{2} = (\omega_0^{2}-v)^{2} + 4\gamma^{2}v = v^{2} - 2(\omega_0^{2}-2\gamma^{2})v + \omega_0^{4} = \bigl(v - (\omega_0^{2}-2\gamma^{2})\bigr)^{2} + 4\gamma^{2}\bigl(\omega_0^{2}-\gamma^{2}\bigr) .

This is a parabola in $v$ with vertex at $v = \omega_0^{2}-2\gamma^{2}$ . If $\gamma < \omega_0/\sqrt2$ the vertex lies in $v > 0$ and is attained at $\omega_r = \sqrt{\omega_0^{2}-2\gamma^{2}}$ , where $|p|^{2} = 4\gamma^{2}(\omega_0^{2}-\gamma^{2}) = 4\gamma^{2}\omega_d^{2}$ , giving $A(\omega_r) = f_0/(2\gamma\omega_d)$ . If $\gamma \ge \omega_0/\sqrt2$ the vertex is at $v \le 0$ , so $|p|^{2}$ is strictly increasing on $v > 0$ and $A$ strictly decreasing: no peak.

Step 2 (velocity peak). The velocity amplitude is $\omega A$ , and

\frac{1}{\omega^{2}A^{2}} = \frac{|p|^{2}}{f_0^{2}\omega^{2}} = \frac{1}{f_0^{2}}\left[\left(\frac{\omega_0^{2}-\omega^{2}}{\omega}\right)^{2} + 4\gamma^{2}\right] .

The bracket is minimal exactly when $\omega_0^{2}-\omega^{2} = 0$ , for any $\gamma > 0$ , giving $\omega = \omega_0$ and $\omega A = f_0/2\gamma$ . The average power delivered by the force $mf_0\cos\omega t$ to the velocity $-\omega A\sin(\omega t-\delta)$ over one period is

\overline{P} = \tfrac12\,m f_0\,\omega A\sin\delta = \tfrac12\,m f_0\,\omega A\cdot\frac{2\gamma\omega}{|p|} = m\gamma\,\omega^{2}A^{2} ,

which is proportional to $(\omega A)^{2}$ , so it peaks at $\omega_0$ as well.

Step 3 (amplification). $A(\omega_0)/A(0) = \omega_0^{2}/(2\gamma\omega_0) = \omega_0/2\gamma = Q$ .

Step 4 (bandwidth). $\overline{P}$ is half its maximum where the bracket in Step 2 doubles, that is where $\bigl|(\omega_0^{2}-\omega^{2})/\omega\bigr| = 2\gamma$ . This gives two quadratics, $\omega^{2} + 2\gamma\omega - \omega_0^{2} = 0$ and $\omega^{2} - 2\gamma\omega - \omega_0^{2} = 0$ , with the positive roots

\omega_- = \sqrt{\omega_0^{2}+\gamma^{2}} - \gamma , \qquad \omega_+ = \sqrt{\omega_0^{2}+\gamma^{2}} + \gamma ,

whose difference is exactly $2\gamma$ : no approximation is involved. Hence $Q = \omega_0/\Delta\omega$ .

Step 5 (energy). For $x = Ae^{-\gamma t}\cos(\omega_dt-\phi)$ ,

E = \tfrac12 m\dot x^{2} + \tfrac12 m\omega_0^{2}x^{2} = \tfrac12 mA^{2}e^{-2\gamma t}\Bigl[\bigl(\gamma\cos\theta + \omega_d\sin\theta\bigr)^{2} + \omega_0^{2}\cos^{2}\theta\Bigr],\quad \theta = \omega_dt-\phi ,

The bracket is a bounded function of $\theta$ alone, so $E(t)$ is the decaying exponential $e^{-2\gamma t}$ times a bounded oscillation. Averaging the bracket over one cycle with $\overline{\cos^{2}\theta} = \overline{\sin^{2}\theta} = \tfrac12$ and $\overline{\sin\theta\cos\theta} = 0$ gives $\tfrac12\gamma^{2}+\tfrac12\omega_d^{2}+\tfrac12\omega_0^{2} = \omega_0^{2}$ , because $\gamma^{2}+\omega_d^{2} = \omega_0^{2}$ . Treating $e^{-2\gamma t}$ as constant over a cycle — legitimate to leading order when $\gamma \ll \omega_d$ — the cycle-averaged energy is $\tfrac12 mA^{2}\omega_0^{2}e^{-2\gamma t}$ , and its fractional loss over one period $T = 2\pi/\omega_d$ is

\frac{|\Delta E|}{E} \;=\; 1 - e^{-2\gamma T} \;=\; 1 - e^{-4\pi\gamma/\omega_d} \;=\; \frac{4\pi\gamma}{\omega_d}\Bigl(1 + O\bigl(\gamma/\omega_d\bigr)\Bigr) .

Therefore $2\pi E/|\Delta E| = \bigl(\omega_d/2\gamma\bigr)\bigl(1+O(\gamma/\omega_d)\bigr)$ , and since $\omega_d \to \omega_0$ as $\gamma/\omega_0 \to 0$ , this tends to $Q$ . $\square$

Remark(Reading the regimes, and what the hypotheses buy)Original

$k > 0$ is needed for oscillation at all. If $k < 0$ (an inverted pendulum) the roots are real of opposite signs and the equilibrium is a saddle: solutions grow. All of the above assumes $\omega_0^{2} > 0$ .
$\gamma > 0$ is what makes the steady state unique and attracting. For $\gamma = 0$ the free solutions do not decay, the "transient" never dies, and driving at $\omega = \omega_0$ is resonant in the strict sense of Constant-coefficient equations: $p(i\omega_0) = 0$ and the response is $x_p = (f_0/2\omega_0)\,t\sin\omega_0t$ , growing linearly for ever. The physical resonance of a damped oscillator is the finite remnant of that: $A(\omega_0) = f_0/(2\gamma\omega_0)$ diverges as $\gamma \to 0$ .
Critical damping is the fastest return, but only just. Overdamping is slower (Step 3 of the first proof: $\lambda_+ \to 0$ as $\gamma \to \infty$ ), and underdamping overshoots. This is why galvanometers, door closers and servo loops are specified near $Q = \tfrac12$ . The optimum is shallow, however: the decay rate is $\gamma$ for $\gamma \le \omega_0$ and $\gamma - \sqrt{\gamma^{2}-\omega_0^{2}}$ above it, which is continuous and has a corner, not a sharp peak, at $\gamma = \omega_0$ .
In the phase plane. As a planar system $\dot x = y$ , $\dot y = -\omega_0^{2}x - 2\gamma y$ , the matrix has trace $\tau = -2\gamma$ and determinant $\Delta = \omega_0^{2}$ , so $\tau^{2}-4\Delta = 4(\gamma^{2}-\omega_0^{2})$ . In the classification of Strogatz (Strogatz (1994), §5.2, Fig. 5.2.8, p. 137), the underdamped case is a stable spiral, the overdamped case a stable node, and the critical case the degenerate node that lies on the parabola $\tau^{2} = 4\Delta$ separating them — the curve he calls "the borderline between nodes and spirals". The three regimes of this unit are the three regions of that diagram with $\Delta > 0$ , $\tau < 0$ .
Two resonant frequencies, not one. The amplitude peaks at $\omega_r = \sqrt{\omega_0^{2}-2\gamma^{2}}$ , the velocity and the absorbed power at $\omega_0$ , and the free motion rings at $\omega_d = \sqrt{\omega_0^{2}-\gamma^{2}}$ . The three are distinct, ordered $\omega_r < \omega_d < \omega_0$ , and agree to $O(\gamma^{2}/\omega_0^{2}) = O(1/4Q^{2})$ . For $Q = 10$ , where $\gamma/\omega_0 = 1/(2Q) = 0.05$ , the spread is $\omega_r/\omega_0 = \sqrt{0.995} = 0.997\,497$ against $\omega_d/\omega_0 = \sqrt{0.9975} = 0.998\,749$ and $1$ : about $0.25\%$ from end to end, as $1/4Q^{2}$ predicts. At the $Q = 5$ of the example below it is about $1\%$ . Loose talk of "the resonant frequency" means $\omega_0$ , and the quantity actually measured in a swept-frequency experiment depends on whether amplitude or power is recorded.
$Q$ measured three ways. Items 3, 4 and 5 of the statement are three experiments: compare the peak response with the static response; measure the half-power bandwidth; or ring the oscillator and count cycles until the energy falls by $e^{-2\pi}$ . Only item 4 is exact for all $\gamma$ ; items 3 and 5 are exact only in the stated senses ( $A(\omega_0)$ exactly, and the energy ratio only as $\gamma/\omega_0 \to 0$ ).

Example(A complete worked oscillator)Original

Take $m = 0.5\ \mathrm{kg}$ , $k = 50\ \mathrm{N\,m^{-1}}$ , $b = 1\ \mathrm{kg\,s^{-1}}$ , driven by $F(t) = 2\cos\omega t$ newtons.

Parameters. $\omega_0 = \sqrt{k/m} = \sqrt{100} = 10\ \mathrm{s^{-1}}$ , $\gamma = b/2m = 1\ \mathrm{s^{-1}}$ , $f_0 = F_0/m = 4\ \mathrm{m\,s^{-2}}$ . So $Q = \omega_0/2\gamma = 5$ , and the oscillator is underdamped.
Free motion. $\omega_d = \sqrt{100-1} = \sqrt{99} = 9.949\,874\ \mathrm{s^{-1}}$ , only $0.5\%$ below $\omega_0$ . Released from $x_0$ at rest, $x = x_0(10/\sqrt{99})e^{-t}\cos(\sqrt{99}\,t-\phi)$ with $\tan\phi = 1/\sqrt{99}$ , i.e. $\phi = 0.1002\ \mathrm{rad}$ . The envelope falls by $e^{-1}$ in $1\ \mathrm{s}$ , i.e. in $\omega_d/2\pi \approx 1.58$ cycles, and by $e^{-5} = 0.0067$ in $5\ \mathrm{s}$ .
Static and resonant response. $A(0) = f_0/\omega_0^{2} = 0.04\ \mathrm{m}$ . At $\omega = \omega_0$ : $A = f_0/(2\gamma\omega_0) = 4/20 = 0.2\ \mathrm{m} = Q\,A(0)$ , with phase lag exactly $\pi/2$ . The amplitude peak is at $\omega_r = \sqrt{100-2} = 9.899\,495$ , where $A = f_0/(2\gamma\omega_d) = 4/(2\sqrt{99}) = 0.201\,008\ \mathrm{m}$ — larger than $A(\omega_0)$ by only $0.5\%$ .
Bandwidth. $\omega_\pm = \sqrt{101}\pm1 = 11.049\,876$ and $9.049\,876$ , so $\Delta\omega = 2.000\,000 = 2\gamma$ and $\omega_0/\Delta\omega = 5 = Q$ , as the statement requires exactly.
A phase. At $\omega = 9$ : $\tan\delta = 2\gamma\omega/(\omega_0^{2}-\omega^{2}) = 18/19$ , so $\delta = 0.758\,378\ \mathrm{rad} = 43.45^{\circ}$ ; the response still leads the $\pi/2$ of resonance.
A full solution, transient included. Drive at $\omega = \omega_0 = 10$ from rest at the origin: $\ddot x + 2\dot x + 100x = 4\cos 10t$ , $x(0) = \dot x(0) = 0$ . The steady state is $x_{\mathrm{ss}} = 0.2\cos(10t - \pi/2) = 0.2\sin 10t$ . Adding the transient $e^{-t}\bigl(c_1\cos\sqrt{99}t + c_2\sin\sqrt{99}t\bigr)$ and imposing the data: $x(0) = c_1 = 0$ , and $\dot x(0) = 2 + c_2\sqrt{99} = 0$ gives $c_2 = -2/\sqrt{99}$ . Hence
$x(t) \;=\; \frac{1}{5}\sin 10t \;-\; \frac{2}{\sqrt{99}}\,e^{-t}\sin\bigl(\sqrt{99}\,t\bigr) .$
Substituting this back into the equation returns $4\cos 10t$ identically, and $x(0) = \dot x(0) = 0$ .
Numbers to check against. $x(1) = -0.071\,735\ \mathrm{m}$ and, once the transient has gone, the motion is $0.2\sin 10t$ : the amplitude builds up to $0.2\ \mathrm{m}$ on the timescale $1/\gamma = 1\ \mathrm{s}$ , which is $Q/\pi \approx 1.6$ cycles. The general rule visible here is that a high- $Q$ oscillator takes about $Q/\pi$ cycles to reach its resonant amplitude — the same $Q$ that sets the bandwidth, as it must, since a narrow resonance is a slow one.

ApplicationOriginal

Measurement instruments. A moving-coil galvanometer, a seismometer and an accelerometer are all this equation, and the designer's choice is a point in the $(\gamma, \omega_0)$ plane. An instrument meant to follow a signal faithfully is damped near critical, so that it neither rings nor lags; one meant to select a frequency is made high- $Q$ , so that the bandwidth $\Delta\omega = 2\gamma$ is narrow. The same trade-off, exactly, sets the ringing and the settling time of a servo loop.

Spectroscopy and linewidths. A classical model of a radiating atom is a damped oscillator, and item 4 above says that its absorption line has full width at half maximum $\Delta\omega = 2\gamma$ in angular frequency. Measuring a linewidth therefore measures a damping rate, and the Lorentzian line shape $\overline{P}(\omega) \propto \omega^{2}/\bigl[(\omega_0^{2}-\omega^{2})^{2}+4\gamma^{2}\omega^{2}\bigr]$ , which near $\omega_0$ becomes $\propto 1/\bigl[(\omega-\omega_0)^{2}+\gamma^{2}\bigr]$ , is the standard fitting function for a resonance in optics, nuclear physics (the Breit-Wigner shape) and nuclear magnetic resonance.

Tuned circuits. For the series $RLC$ circuit, $Q = R^{-1}\sqrt{L/C}$ is the figure of merit of a tuned stage: the selectivity of a radio receiver is $\omega_0/\Delta\omega$ , so a $Q$ of a few hundred separates stations a fraction of a percent apart in frequency. The statement that the circuit does not oscillate for $R \ge 2\sqrt{L/C}$ — that is, for $Q \le \tfrac12$ — is Hassani's (Hassani (2013), p. 448) version of the overdamping condition.

Structures and vehicles. The collapse of structures at resonance is the $A(\omega_r) = f_0/(2\gamma\omega_d)$ of item 1: lightly damped structures amplify periodic loads by a factor $Q$ . Engineering practice is to move $\omega_0$ away from the expected forcing, or to add damping, and the design quantity is the dimensionless damping ratio $\zeta = \gamma/\omega_0 = 1/(2Q)$ . Vehicle suspensions are tuned to $\zeta \approx 0.2$ to $0.3$ , well underdamped, because comfort penalises a slow return more than a small overshoot.

Normal modes of larger systems. A system of many coupled linear oscillators diagonalises into independent modes (see Spectral theorem for Hermitian matrices), each obeying this very equation with its own $\omega_0$ and $\gamma$ . The whole response of a molecule, a crystal lattice, an optical cavity or a bridge is then a sum of terms of the form $f_0/p_j(i\omega)$ , and the measured spectrum is a sum of Lorentzians. This is why the single-oscillator analysis above is worth doing in such detail: it is the building block of linear response theory, where $1/p(i\omega)$ is the susceptibility and the relation between its real and imaginary parts is a Kramers-Kronig relation.

Series solutions

Method

Power-series solutions at ordinary points

Adapted

Add to PDF

At a point where the normalised coefficients are analytic, every solution is analytic: substituting a power series turns the equation into a recurrence for the coefficients, the first $n$ of which are the initial data, and the series converges at least as far as the coefficient series do.

Prerequisites: Structure of linear equations and the Wronskian, Taylor series of analytic functions

Definition(Ordinary and singular points)AdaptedAdapted from Riley, Hobson & Bence (2006), §16.1.1, p. 533; Coddington (1989), Ch. 3, §7, p. 126; Arfken, Weber & Harris (2013), §7.4, p. 344

A function $g$ defined near $x_0$ is analytic at $x_0$ if it can be written as a power series

g(x) = \sum_{k=0}^{\infty} g_k\,(x - x_0)^{k}

convergent on $|x - x_0| < r_0$ for some $r_0 > 0$ . Such a $g$ has derivatives of every order there, computed by differentiating term by term, and the differentiated series have the same radius (Taylor series of analytic functions).

Write the equation in normalised form,

y^{(n)} + a_1(x)\,y^{(n-1)} + \cdots + a_n(x)\,y = 0 ,

and for $n = 2$ ,

y'' + p(x)\,y' + q(x)\,y = 0 .

The point $x_0$ is an ordinary point of the equation if every coefficient $a_1, \dots, a_n$ (for $n = 2$ : both $p$ and $q$ ) is analytic at $x_0$ , and a singular point otherwise. A point where the original, unnormalised leading coefficient vanishes is a singular point unless the division cancels it.

By translating, $x_0 = 0$ throughout what follows; for a general $x_0$ substitute $X = x - x_0$ .

Method(Analytic solutions at an ordinary point)AdaptedAdapted from Coddington (1989), Ch. 3, §7, Thm 12, p. 129; Hassani (2013), §14.6.1, Thms 14.6.2 and 14.6.4, pp. 440-441

Suppose $a_1, \dots, a_n$ have power-series expansions in powers of $x - x_0$ convergent on $|x - x_0| < r_0$ with $r_0 > 0$ , and let $\alpha_1, \dots, \alpha_n$ be any constants. Then the initial-value problem

L(y) = 0 , \qquad y(x_0) = \alpha_1, \ \dots, \ y^{(n-1)}(x_0) = \alpha_n

has a solution given by a power series

\varphi(x) = \sum_{k=0}^{\infty} c_k\,(x-x_0)^{k}

convergent on $|x - x_0| < r_0$ . Its first $n$ coefficients are the data,

k!\,c_k = \alpha_{k+1} , \qquad k = 0, 1, \dots, n-1 ,

and $c_k$ for $k \ge n$ is computed in terms of $c_0, \dots, c_{k-1}$ by substituting the series into $L(y) = 0$ . Since the initial-value problem has only one solution (Structure of linear equations and the Wronskian), every solution of $L(y) = 0$ on $|x-x_0| < r_0$ is given by such a convergent series: at an ordinary point all solutions are analytic.

For $n = 2$ , with $p = \sum_m a_mx^{m}$ and $q = \sum_m b_mx^{m}$ about $x_0 = 0$ , substituting $y = \sum_k c_kx^{k}$ gives, for every $n \ge 0$ , the recurrence

(n+2)(n+1)\,c_{n+2} \;=\; -\sum_{m=0}^{n}\Bigl[(n-m+1)\,a_m\,c_{n-m+1} \;+\; b_m\,c_{n-m}\Bigr] ,

which determines $c_2, c_3, \dots$ uniquely from $c_0$ and $c_1$ , and the resulting series has radius of convergence at least the smaller of the radii of the series for $p$ and $q$ .

Note(Fuchs' theorem and the two exponents)AdaptedAdapted from Arfken, Weber & Harris (2013), §7.5, Fuchs' theorem and Summary, p. 355; Arfken, Weber & Harris (2013), §7.5, Exercise 7.5.2, p. 356

Arfken, Weber and Harris (Arfken, Weber & Harris (2013), p. 355) record the scope of the method as Fuchs' theorem: at least one power-series solution is always available, so long as the expansion point is either ordinary or a regular singular point, the worst case the theorem allows. About an irregular singular point the method may fail outright.

The more general substitution $y = x^{s}\sum_{j\ge0}a_jx^{j}$ with $a_0 \ne 0$ and $s$ undetermined — the Frobenius substitution of The Frobenius method — needs no modification at an ordinary point: there the indicial equation has roots $s = 0$ and $s = 1$ ((Arfken, Weber & Harris (2013), Exercise 7.5.2, p. 356)), so one solution is a plain power series and the other is $x$ times a power series, which is again a plain power series. Nothing is gained by carrying $s$ , and the plain substitution $y = \sum_k c_kx^{k}$ suffices. This is the practical difference between this unit and the next.

Idea of the proofOriginal

The substitution is bookkeeping, and the only real question is whether the bookkeeping converges. Substituting $y = \sum c_kx^{k}$ into the equation produces, at order $x^{n}$ , one term containing $c_{n+2}$ with coefficient $(n+2)(n+1)$ — which is never zero — and finitely many terms containing lower coefficients. So every coefficient after the first two is forced, and the equation becomes an explicit recurrence: the differential equation has been traded for an algebraic one. The two free constants $c_0, c_1$ are exactly $y(0)$ and $y'(0)$ , which is why a second-order equation has a two-dimensional solution space, visible here as the two free coefficients.

Convergence is the part that needs an estimate. The recurrence expresses $c_{n+2}$ as a sum of up to $n$ terms, so a naive bound on $|c_{n+2}|$ grows; but the divisor $(n+2)(n+1)$ grows too, and the coefficients $a_m, b_m$ decay geometrically because $p$ and $q$ have positive radii. Comparing the recurrence with a geometric majorant turns this into an elementary induction, and gives convergence in some disc. That the disc may be taken as large as the coefficients' own disc is a sharper fact, proved in the sources by a complex-analytic argument.

Proof(The recurrence, and convergence in some disc)AdaptedAdapted from Hassani (2013), §14.6.1, Eq. (14.28) and Thm 14.6.2, pp. 439-440; Riley, Hobson & Bence (2006), §16.2, Eqs. (16.9)-(16.11), p. 536

Take $n = 2$ and $x_0 = 0$ ; the general order is identical with longer sums.

Step 1 (differentiate the series). If $y = \sum_{k\ge0}c_kx^{k}$ converges on $|x| < R$ , term-by-term differentiation (Taylor series of analytic functions) gives, after shifting the summation index so that each series is written as a coefficient of $x^{n}$ ,

y' = \sum_{n\ge0}(n+1)c_{n+1}x^{n} , \qquad y'' = \sum_{n\ge0}(n+2)(n+1)c_{n+2}x^{n} .

Step 2 (multiply the series). On the common disc of convergence, products of power series may be multiplied term by term:

p\,y' = \sum_{n\ge0}\Bigl[\sum_{m=0}^{n}a_m(n-m+1)c_{n-m+1}\Bigr]x^{n} , \qquad q\,y = \sum_{n\ge0}\Bigl[\sum_{m=0}^{n}b_mc_{n-m}\Bigr]x^{n} .

Step 3 (a power series vanishes only if all its coefficients do). Adding the three series, $y''+py'+qy = 0$ holds on a neighbourhood of $0$ if and only if every coefficient of $x^{n}$ vanishes, which is the displayed recurrence. The coefficient of $c_{n+2}$ is $(n+2)(n+1) \ne 0$ , so the recurrence determines $c_{n+2}$ from $c_0, \dots, c_{n+1}$ , and $c_0 = y(0)$ , $c_1 = y'(0)$ are free. This proves the formal part: there is exactly one formal power series with given $c_0, c_1$ satisfying the equation.

Step 4 (geometric bounds on the coefficients of $p$ and $q$ ). Let $p$ and $q$ converge on $|x| < r$ , and fix $0 < \rho < r$ . Then $\sum_m|a_m|\rho^{m}$ and $\sum_m|b_m|\rho^{m}$ converge, so their terms are bounded: there is $M \ge 1$ with

|a_m| \le M\rho^{-m} , \qquad |b_m| \le M\rho^{-m} \qquad\text{for all } m \ge 0 .

Step 5 (rescale). Put $d_k = \rho^{k}|c_k|$ . Taking absolute values in the recurrence and inserting the bounds of Step 4,

(n+2)(n+1)\,\rho^{-(n+2)}d_{n+2} \;\le\; M\sum_{m=0}^{n}\Bigl[(n-m+1)\rho^{-(n+1)}d_{n-m+1} + \rho^{-n}d_{n-m}\Bigr] ,

because $\rho^{-m}|c_{n-m+1}| = \rho^{-(n+1)}d_{n-m+1}$ and $\rho^{-m}|c_{n-m}| = \rho^{-n}d_{n-m}$ . Multiplying by $\rho^{n+2}$ and renaming the summation index ( $j = n-m+1$ , then $j = n-m$ ),

(n+2)(n+1)\,d_{n+2} \;\le\; M\Bigl[\rho\sum_{j=1}^{n+1}j\,d_j \;+\; \rho^{2}\sum_{j=0}^{n}d_j\Bigr] \;\le\; B\sum_{j=0}^{n+1}(j+1)\,d_j , \qquad B = M\rho(1+\rho) .

Step 6 (induct on partial sums). Set $e_j = (j+1)d_j$ and $S_n = \sum_{j=0}^{n}e_j$ . Since $(n+3) \le \tfrac32(n+2)(n+1)$ for every $n \ge 0$ ,

e_{n+2} = (n+3)d_{n+2} \le \frac{n+3}{(n+2)(n+1)}\,B\,S_{n+1} \le \tfrac32B\,S_{n+1} ,

so $S_{n+2} = S_{n+1} + e_{n+2} \le \lambda S_{n+1}$ with $\lambda = 1 + \tfrac32B \ge 1$ . By induction $S_n \le S_1\lambda^{\,n-1}$ for $n \ge 1$ , and hence

d_n \le e_n \le S_n \le S_1\lambda^{\,n-1} , \qquad\text{so}\qquad |c_n| = \rho^{-n}d_n \le \frac{S_1}{\lambda}\left(\frac{\lambda}{\rho}\right)^{n} .

Step 7 (conclude). Therefore $\limsup_n|c_n|^{1/n} \le \lambda/\rho$ , so the formal series converges on $|x| < \rho/\lambda$ , a disc of positive radius. On that disc all the manipulations of Steps 1-3 are valid, so the sum is a genuine solution, with $y(0) = c_0$ and $y'(0) = c_1$ . $\square$

Note(The sharp radius)Original

Step 7 gives convergence on a possibly small disc, which is all that is needed to know that the formal series is a solution. The sharp statement — that the radius is at least $r$ itself, the radius of the coefficient series — is the content of Coddington's Theorem 12 (Coddington (1989), Ch. 3, Thm 12, p. 129) and Hassani's Theorem 14.6.4 (Hassani (2013), Thm 14.6.4, p. 441), and we use it without proof: Coddington devotes his Chapter 3, §9 to it, and Hassani refers the reader to Birkhoff and Rota. Both proofs work in the complex domain, where the natural statement is that a linear equation with coefficients holomorphic on a disc has solutions holomorphic on that whole disc.

Remark(Where the hypotheses are used, and a disagreement between sources)Original

Normalisation matters. The theorem is about the equation after division by the leading coefficient. Hassani (Hassani (2013), Example 14.6.3, pp. 440-441) gives $x^{2}y' - y + x = 0$ , whose formal power series is $\sum_{n\ge1}(n-1)!\,x^{n}$ — a perfectly determinate formal solution that converges for no $x \ne 0$ . The leading coefficient vanishes at $0$ , so $0$ is not an ordinary point and the theorem does not apply. Formal solvability and convergence are different questions, and Step 3 only settles the first.
$(n+2)(n+1) \ne 0$ is what makes the recurrence work. At a singular point the coefficient multiplying the highest unknown can vanish, the recurrence stalls, and one must carry an undetermined exponent instead (The Frobenius method).
Radii: what the books say. Hassani (Hassani (2013), Thm 14.6.4, p. 441) and Coddington (Coddington (1989), Thm 12, p. 129) state a lower bound: the radius is at least that of the coefficient series. Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §16.2, p. 536) state instead that the radius equals the distance from the expansion point to the nearest singular point of the equation, adding that at that radius the series may or may not converge. We follow the lower-bound version, because particular solutions can do better: Legendre's equation has singular points at $x = \pm1$ , yet for integer $\ell$ one of its two solutions is a polynomial, whose series converges everywhere. The example below shows the same phenomenon in a two-line computation. The sharp distance-to-the-singularity statement is correct as a statement about the generic solution, and it is the right rule of thumb in practice.
Why complex singularities matter for a real problem. The radius is limited by the nearest singularity in the complex plane, even when the equation, the interval and the solution are entirely real. In the example below, the series stops converging at $|x| = 1$ because the coefficients are singular at $x = \pm i$ , which no amount of looking at the real axis would reveal. This is the practical reason the theory of series solutions is complex-analytic.

Method(Series solution at an ordinary point)AdaptedAdapted from Riley, Hobson & Bence (2006), §16.2, pp. 535-537; Arfken, Weber & Harris (2013), §7.5, pp. 347-348

Check that $x_0$ is ordinary: normalise and confirm that the coefficients are analytic at $x_0$ . Translate so that $x_0 = 0$ .
Find the radius you are entitled to: locate the equation's singular points in the complex plane and measure how far the closest of them lies from $0$ . This tells you in advance where the answer can be used.
Substitute $y = \sum_{k\ge0}c_kx^{k}$ , differentiating term by term, and do not divide by the leading coefficient first if that introduces denominators: multiply the equation through by the polynomial leading coefficient instead, which keeps the recurrence finite.
Shift indices so that every series is a sum over the coefficient of $x^{n}$ , then set each coefficient to zero. The result is a recurrence for $c_{n+2}$ (or $c_{n+k}$ ) in terms of earlier coefficients.
Generate two solutions: take $c_0 = 1$ , $c_1 = 0$ for one and $c_0 = 0$ , $c_1 = 1$ for the other. These are automatically independent, since their Wronskian at $0$ is $\det\begin{pmatrix}1&0\\0&1\end{pmatrix} = 1$ .
Identify or estimate. Look for a closed form; look for termination (a recurrence factor that vanishes gives a polynomial solution); and check the radius by the ratio test on the recurrence.

Example(A two-term recurrence with one terminating solution)Original

Solve

(1 + x^{2})\,y'' + 2x\,y' - 2y = 0

by series about $x = 0$ .

$x = 0$ is ordinary. Normalised, $p(x) = 2x/(1+x^{2})$ and $q(x) = -2/(1+x^{2})$ , both analytic at $0$ with radius exactly $1$ , since $1 + x^{2}$ vanishes at $x = \pm i$ . So the theorem guarantees solutions analytic on $|x| < 1$ , and the nearest singular points are the complex pair $\pm i$ at distance $1$ .
Substitute, without dividing. With $y = \sum_{n\ge0}c_nx^{n}$ ,
$y'' = \sum_{n\ge0}(n+2)(n+1)c_{n+2}x^{n}, \quad x^{2}y'' = \sum_{n\ge0}n(n-1)c_nx^{n}, \quad 2xy' = \sum_{n\ge0}2n\,c_nx^{n} .$
The coefficient of $x^{n}$ in the equation is therefore
$(n+2)(n+1)c_{n+2} + \bigl[n(n-1) + 2n - 2\bigr]c_n = 0 .$
Simplify the recurrence. Since $n(n-1)+2n-2 = n^{2}+n-2 = (n+2)(n-1)$ , the factor $n+2$ cancels:
$c_{n+2} \;=\; -\,\frac{n-1}{n+1}\,c_n , \qquad n \ge 0 .$
Odd and even coefficients do not mix.
The odd solution terminates. Putting $n = 1$ gives $c_3 = 0$ , and then $c_5 = c_7 = \cdots = 0$ . With $c_0 = 0$ , $c_1 = 1$ the series is the single term
$y_1(x) = x ,$
which indeed solves the equation: $0 + 2x - 2x = 0$ . Its "series" converges for every $x$ .
The even solution. With $c_0 = 1$ , $c_1 = 0$ : $n = 0$ gives $c_2 = 1$ ; $n = 2$ gives $c_4 = -\tfrac13$ ; $n = 4$ gives $c_6 = -\tfrac35 c_4 = \tfrac15$ ; $n = 6$ gives $c_8 = -\tfrac57c_6 = -\tfrac17$ . By induction $c_{2k} = (-1)^{k+1}/(2k-1)$ for $k \ge 1$ , since
$c_{2k+2} = -\frac{2k-1}{2k+1}c_{2k} = -\frac{2k-1}{2k+1}\cdot\frac{(-1)^{k+1}}{2k-1} = \frac{(-1)^{k+2}}{2k+1} .$
Hence
$y_2(x) = 1 + x^{2} - \frac{x^{4}}{3} + \frac{x^{6}}{5} - \frac{x^{8}}{7} + \cdots = 1 + \sum_{k=1}^{\infty}\frac{(-1)^{k+1}}{2k-1}\,x^{2k} .$
Closed form and radius. Since $\arctan x = \sum_{k\ge1}(-1)^{k-1}x^{2k-1}/(2k-1)$ for $|x|<1$ ,
$y_2(x) = 1 + x\arctan x .$
Direct substitution confirms it: $y_2' = \arctan x + x/(1+x^{2})$ , $y_2'' = 1/(1+x^{2}) + (1-x^{2})/(1+x^{2})^{2}$ , and
$(1+x^{2})y_2'' + 2xy_2' - 2y_2 = 1 + \frac{1-x^{2}}{1+x^{2}} + 2x\arctan x + \frac{2x^{2}}{1+x^{2}} - 2 - 2x\arctan x = 0 .$
The ratio test on $|c_{2k+2}/c_{2k}| = (2k-1)/(2k+1) \to 1$ gives radius exactly $1$ for $y_2$ — the distance to $\pm i$ , as predicted — while $y_1 = x$ has infinite radius. One solution attains the predicted radius and the other beats it, which is why the theorem is stated as a lower bound.
A definite solution and a number. $y(0) = 1$ , $y'(0) = 0$ selects $y = y_2$ . At $x = \tfrac12$ :
$y\bigl(\tfrac12\bigr) = 1 + \tfrac12\arctan\tfrac12 = 1.231\,823\,8 ,$
while the partial sum through $x^{8}$ is $1 + 0.25 - 0.020833 + 0.003125 - 0.000558 = 1.231\,734$ , already right to four decimals. The Wronskian check: from $y_1(0)=0$ , $y_1'(0)=1$ , $y_2(0)=1$ , $y_2'(0)=0$ we get $W(y_1,y_2)(0) = y_1(0)y_2'(0) - y_1'(0)y_2(0) = -1 \ne 0$ , so $\{x,\,1+x\arctan x\}$ is a fundamental set on $(-1,1)$ .

ApplicationOriginal

Where the special functions come from. Legendre's, Hermite's and Chebyshev's equations all have $x = 0$ as an ordinary point, and this method is how their solutions are constructed. The recurrences they produce share the feature seen in step 4 above: for special values of a parameter a factor vanishes and the series terminates, giving a polynomial. For Hermite's equation $y'' - 2xy' + 2\nu y = 0$ the recurrence is $c_{n+2} = 2(n-\nu)c_n/\bigl((n+2)(n+1)\bigr)$ , which terminates exactly when $\nu$ is a non-negative integer. In the quantum harmonic oscillator that condition is the quantisation of the energy: the non-terminating series grows like $e^{x^{2}/2}$ and is not normalisable, so only the terminating ones are physical, and $\nu = n$ gives $E = (n+\tfrac12)\hbar\omega$ . The same mechanism quantises angular momentum through Legendre's equation.

Airy functions and turning points. The Airy equation $y'' = xy$ has no singular point in the finite plane, so its series solutions about $0$ converge for all $x$ ; they are the standard way to compute $\mathrm{Ai}$ and $\mathrm{Bi}$ for moderate argument, and they describe the behaviour of a wavefunction at a classical turning point, where the WKB approximation fails.

Numerical work and Taylor methods. The recurrence is a recipe for generating Taylor coefficients of the solution to any order at the cost of a few operations each, which is the basis of Taylor-series and automatic-differentiation integrators for ordinary differential equations. Their step size is limited by exactly the radius this unit estimates, which is why such integrators are efficient for problems with distant singularities and poor near them.

Knowing in advance where a model can be trusted. Step 2 of the method is the practically valuable one. A model whose coefficients are singular at a nearby complex point has a series solution that is useless beyond that distance, no matter how many terms are computed, and the breakdown has nothing to do with the physical interval. Legendre's functions of the second kind diverge at $x = \pm1$ , that is at the poles of a sphere, and the series expansions of multipole fields converge only inside a sphere reaching the nearest source — a statement about radii of convergence of exactly this kind.

Method

The Frobenius method

Adapted

Add to PDF

At a regular singular point the substitution $y = x^{\sigma}\sum_n a_nx^n$ works: the lowest power gives a quadratic indicial equation for $\sigma$ , the rest a recurrence, and the larger exponent always yields a solution — while equal exponents, or exponents differing by an integer, can block the second one and force a logarithm.

Prerequisites: Power-series solutions at ordinary points, Euler-Cauchy equations, Reduction of order

Definition(Regular singular points and Frobenius series)AdaptedAdapted from Hassani (2013), §15.2, Def. 15.2.2 and the expansions on p. 465; Riley, Hobson & Bence (2006), §16.3, Eq. (16.12) and the definitions of $s$ and $t$, p. 539

Consider

y'' + p(x)\,y' + q(x)\,y = 0

on a punctured interval around $x_0$ . The point $x_0$ is a regular singular point if $p$ has at worst a simple pole and $q$ at worst a pole of order $2$ at $x_0$ — equivalently, if

s(x) = (x - x_0)\,p(x) \qquad\text{and}\qquad t(x) = (x-x_0)^{2}\,q(x)

are both analytic at $x_0$ . Take $x_0 = 0$ from now on (translate otherwise) and write

s(x) = \sum_{n\ge0}s_nx^{n} , \qquad t(x) = \sum_{n\ge0}t_nx^{n} ,

both convergent for $|x| < r$ . Multiplying the equation by $x^{2}$ puts it in the form

x^{2}y'' + x\,s(x)\,y' + t(x)\,y = 0 ,

which is an Euler equation to leading order, with $s$ and $t$ frozen at their values at $0$ .

A Frobenius series, or generalised power series, is an expression

y = x^{\sigma}\sum_{n=0}^{\infty}a_nx^{n} , \qquad a_0 \neq 0 ,

with $\sigma$ a real or complex number. The normalisation $a_0 \ne 0$ costs nothing: if $a_0$ were zero, $\sigma$ could be increased to make the leading coefficient nonzero.

Method(Indicial equation and recurrence)AdaptedAdapted from Riley, Hobson & Bence (2006), §16.3, Eqs. (16.15)-(16.16), pp. 539-540; Hassani (2013), §15.2, Eqs. (15.7)-(15.9), p. 465

Substituting a Frobenius series into $x^{2}y'' + xsy' + ty = 0$ and dividing by $x^{\sigma}$ gives

\sum_{n\ge0}\Bigl[(n+\sigma)(n+\sigma-1) + s(x)(n+\sigma) + t(x)\Bigr]a_nx^{n} \;=\; 0 .

Define the indicial polynomial

I(\sigma) \;=\; \sigma(\sigma-1) + s_0\,\sigma + t_0 , \qquad s_0 = s(0),\ t_0 = t(0) .

Then:

The coefficient of $x^{0}$ gives, since $a_0 \ne 0$ , the indicial equation
$I(\sigma) = \sigma(\sigma-1) + s_0\sigma + t_0 = 0 ,$
whose two roots $\sigma_1, \sigma_2$ are the characteristic exponents (or indices) of the regular singular point. They are the exponents of the Euler equation obtained by replacing $s, t$ by $s_0, t_0$ .
The coefficient of $x^{n}$ for $n \ge 1$ gives the recurrence
$I(n+\sigma)\,a_n \;=\; -\sum_{k=0}^{n-1}\Bigl[(k+\sigma)\,s_{n-k} + t_{n-k}\Bigr]a_k .$
Since $I(\sigma) = (\sigma-\sigma_1)(\sigma-\sigma_2)$ ,
$I(n+\sigma_1) = n\,(n + \sigma_1 - \sigma_2) , \qquad I(n+\sigma_2) = n\,(n - (\sigma_1-\sigma_2)) .$

Method(Which exponents give Frobenius solutions)AdaptedAdapted from Hassani (2013), §15.2, Thm 15.2.3, p. 466, Eq. (15.10), p. 467 and Thm 15.2.6, pp. 468-469; Riley, Hobson & Bence (2006), §16.3.1-16.3.3, pp. 540-543; Arfken, Weber & Harris (2013), §7.5, Summary, p. 355

Label the exponents so that $\operatorname{Re}\sigma_1 \ge \operatorname{Re}\sigma_2$ .

The larger exponent always works. The recurrence determines $a_1, a_2, \dots$ uniquely from $a_0$ for $\sigma = \sigma_1$ , because $I(n+\sigma_1) = n(n+\sigma_1-\sigma_2) \ne 0$ for every $n \ge 1$ . The resulting Frobenius series converges on $0 < |x| < r$ .
$\sigma_1 - \sigma_2$ not an integer. Then $I(n+\sigma_2) \ne 0$ for all $n \ge 1$ as well, and there are two Frobenius solutions,
$y_1 = x^{\sigma_1}\sum_{n\ge0}a_nx^{n} , \qquad y_2 = x^{\sigma_2}\sum_{n\ge0}b_nx^{n} ,$
which are linearly independent because $y_2/y_1 = x^{\sigma_2-\sigma_1}(\cdots)$ is not constant. This case includes a complex-conjugate pair $\sigma_2 = \bar\sigma_1$ with $\operatorname{Im}\sigma_1 \ne 0$ , since then $\sigma_1-\sigma_2 = 2i\operatorname{Im}\sigma_1$ is not a real integer.
$\sigma_1 = \sigma_2$ . Only one Frobenius series exists, and the second solution necessarily contains a logarithm.
$\sigma_1 - \sigma_2 = N$ , a positive integer. The recurrence for $\sigma_2$ runs without trouble up to $a_{N-1}$ and then meets $I(N+\sigma_2) = I(\sigma_1) = 0$ , so it demands that the right-hand side at $n = N$ vanish. Call that right-hand side $\rho_N$ (it is determined by $s, t$ and $a_0$ ). If $\rho_N = 0$ the equation $0\cdot a_N = 0$ holds, $a_N$ may be chosen freely, and a second Frobenius series exists. If $\rho_N \ne 0$ there is no second Frobenius series, and Hassani (Hassani (2013), Thm 15.2.6, p. 468) gives the basis in the form
$y_1 = x^{\sigma_1}\Bigl(1 + \sum_{k\ge1}a_kx^{k}\Bigr) , \qquad y_2 = x^{\sigma_2}\Bigl(1 + \sum_{k\ge1}b_kx^{k}\Bigr) + C\,y_1\log x ,$
with the power series convergent in a neighbourhood of $0$ .

For complex exponents, "larger" means larger real part.

Idea of the proofOriginal

The equation, written as $x^{2}y'' + x\,s(x)y' + t(x)y = 0$ , is an Euler equation whose constants have been allowed to vary analytically. Near $x = 0$ only $s(0)$ and $t(0)$ matter, so the solutions must begin like the Euler solutions $x^{\sigma}$ with $\sigma$ a root of the Euler indicial polynomial — and the higher coefficients of $s$ and $t$ then correct the power into a series. That is the whole method: guess an Euler solution and let it drag a power series behind it.

What can go wrong is visible in the recurrence. Solving for $a_n$ means dividing by $I(n+\sigma)$ , and $I$ vanishes only at $\sigma_1$ and $\sigma_2$ . Starting from the larger exponent, $n+\sigma_1$ never returns to a root, so nothing obstructs. Starting from the smaller, $n + \sigma_2$ hits $\sigma_1$ exactly when $n$ equals the integer difference, and at that step the recurrence either is consistent by accident (and gives a second series) or is contradictory (and no second series exists). In the contradictory case the missing solution is supplied by reduction of order, and the logarithm appears for a precise reason: the integrand $\varphi_1^{-2}e^{-\int p}$ behaves like $x^{-(\sigma_1-\sigma_2)-1}$ , which integrates to a power unless the exponent is $-1$ , that is unless $\sigma_1 = \sigma_2$ ; and in the integer case the $x^{-1}$ term of the expansion contributes a logarithm.

Proof(Indicial equation and recurrence)AdaptedAdapted from Riley, Hobson & Bence (2006), §16.3, pp. 539-540; Hassani (2013), §15.2, p. 465

Step 1 (differentiate). For $y = x^{\sigma}\sum_{n\ge0}a_nx^{n} = \sum_{n\ge0}a_nx^{n+\sigma}$ , term-by-term differentiation on $x > 0$ gives

y' = \sum_{n\ge0}(n+\sigma)a_nx^{n+\sigma-1} , \qquad y'' = \sum_{n\ge0}(n+\sigma)(n+\sigma-1)a_nx^{n+\sigma-2} .

Step 2 (assemble). Then $x^{2}y''$ , $x\,s(x)y'$ and $t(x)y$ all carry the common factor $x^{\sigma}$ :

x^{2}y'' + xsy' + ty = x^{\sigma}\sum_{n\ge0}\Bigl[(n+\sigma)(n+\sigma-1) + s(x)(n+\sigma) + t(x)\Bigr]a_nx^{n} .

Dividing by $x^{\sigma}$ (legitimate for $x \ne 0$ ) gives the displayed identity.

Step 3 (lowest power). Expanding $s$ and $t$ as power series and collecting the coefficient of $x^{0}$ : only $n = 0$ contributes, with $s$ and $t$ at their constant terms, giving $\bigl[\sigma(\sigma-1) + s_0\sigma + t_0\bigr]a_0 = I(\sigma)\,a_0$ . Since $a_0 \ne 0$ by the normalisation, $I(\sigma) = 0$ .

Step 4 (the general power). The coefficient of $x^{n}$ , $n \ge 1$ , receives $(n+\sigma)(n+\sigma-1)a_n$ from the first term, $\sum_{k=0}^{n}(k+\sigma)s_{n-k}a_k$ from the second and $\sum_{k=0}^{n}t_{n-k}a_k$ from the third. The $k = n$ parts of the last two sums are $(n+\sigma)s_0a_n$ and $t_0a_n$ , which combine with the first term into $I(n+\sigma)a_n$ . Setting the whole coefficient to zero and moving the $k < n$ terms to the right gives the recurrence.

Step 5 (the factorisation). $I$ is monic of degree $2$ with roots $\sigma_1, \sigma_2$ , so $I(\sigma) = (\sigma-\sigma_1)(\sigma-\sigma_2)$ and therefore $I(n+\sigma_1) = n(n+\sigma_1-\sigma_2)$ , $I(n+\sigma_2) = (n+\sigma_2-\sigma_1)n$ . $\square$

Proof(Which exponents work)AdaptedAdapted from Hassani (2013), §15.2, pp. 466-469; Riley, Hobson & Bence (2006), §16.3.1-16.3.3, pp. 540-543

Step 1 (the larger exponent). For $n \ge 1$ , $I(n+\sigma_1) = n(n+\sigma_1-\sigma_2)$ . Here $n \ne 0$ and $\operatorname{Re}(n+\sigma_1-\sigma_2) = n + \operatorname{Re}(\sigma_1-\sigma_2) \ge n \ge 1 > 0$ , so $I(n+\sigma_1) \ne 0$ . Hence the recurrence solves for $a_n$ at every step and the coefficients are uniquely determined by $a_0$ . This proves part 1 up to convergence, for which see the note below.

Step 2 (non-integer difference). $I(n+\sigma_2) = n\bigl(n - (\sigma_1-\sigma_2)\bigr)$ vanishes for some $n \ge 1$ only if $\sigma_1 - \sigma_2 = n$ is a positive integer. If $\sigma_1-\sigma_2 \notin \mathbb{Z}$ this never happens, so the recurrence for $\sigma_2$ also runs freely. Independence: if $c_1y_1 + c_2y_2 = 0$ on an interval $(0,\epsilon)$ with $c_2 \ne 0$ , dividing by $x^{\sigma_2}$ gives $c_2(a_0 + O(x)) = -c_1x^{\sigma_1-\sigma_2}(b_0 + O(x))$ ; letting $x \to 0^{+}$ , the left side tends to $c_2a_0 \ne 0$ while the right side tends to $0$ if $\operatorname{Re}(\sigma_1-\sigma_2) > 0$ , a contradiction. If $\operatorname{Re}(\sigma_1-\sigma_2) = 0$ with $\sigma_1 \neq \sigma_2$ , the exponents are a conjugate pair and $x^{\sigma_1-\sigma_2} = e^{i\tau\log x}$ oscillates without limit, so no constants can make the identity hold. Hence $y_1, y_2$ are independent.

Step 3 (equal exponents force a logarithm). Let $\sigma_1 = \sigma_2 = \sigma$ . The two roots of $I$ coincide, so $2\sigma = 1 - s_0$ by the sum of the roots. Let $y_1 = x^{\sigma}(a_0 + \cdots)$ be the Frobenius solution of Step 1 and apply reduction of order on an interval where $y_1 \ne 0$ : a second solution is $y_1\int \varphi$ with

\varphi(x) = \frac{1}{y_1(x)^{2}}\exp\left[-\int p\right] .

Now $p = s/x = s_0/x + (\text{analytic})$ , so $\exp[-\int p] = x^{-s_0}h(x)$ with $h$ analytic and $h(0) \ne 0$ ; and $y_1^{-2} = x^{-2\sigma}g(x)$ with $g$ analytic, $g(0) = a_0^{-2} \ne 0$ . Hence

\varphi(x) = x^{-2\sigma-s_0}\,g(x)h(x) = x^{-1}\,gh ,

using $2\sigma + s_0 = 1$ . Since $g(0)h(0) \ne 0$ , the Laurent expansion of $\varphi$ has a nonzero residue at $0$ , so $\int\varphi = c\log x + (\text{analytic})$ with $c \ne 0$ , and the second solution is $c\,y_1\log x + y_1\cdot(\text{analytic})$ . No Frobenius series — a function of the form $x^{\mu}\times$ (convergent power series) — can contain a $\log$ , so the second solution is not of that form, and only one Frobenius solution exists.

Step 4 (integer difference). Let $\sigma_1 - \sigma_2 = N \ge 1$ . Running the recurrence at $\sigma = \sigma_2$ , the factors $I(n+\sigma_2) = n(n-N)$ are nonzero for $1 \le n \le N-1$ , so $a_1, \dots, a_{N-1}$ are determined by $a_0$ . At $n = N$ the left-hand side is $I(\sigma_1)a_N = 0$ while the right-hand side is a definite number $\rho_N$ determined by $a_0, \dots, a_{N-1}$ . If $\rho_N \ne 0$ the recurrence is inconsistent and no Frobenius series with exponent $\sigma_2$ exists; if $\rho_N = 0$ the step imposes no condition, $a_N$ is free, and the recurrence continues (the factors $I(n+\sigma_2)$ for $n > N$ are again nonzero), producing a second Frobenius series. The same reduction-of-order computation as in Step 3, now with $2\sigma_1 + s_0 = 1 + N$ , gives $\varphi = x^{-1-N}gh$ , whose Laurent expansion may or may not have a nonzero residue at $0$ . When that residue is nonzero, integrating produces a logarithm and the second solution takes the shape Hassani (Hassani (2013), Thm 15.2.6, p. 468) gives, with $C$ the residue; when it vanishes, $\int\varphi$ is a Laurent polynomial plus a power series and the second solution is again a Frobenius series. $\square$

Note(Convergence, and what happens at an irregular point)Original

That the Frobenius series actually converge — on $0 < |x| < r$ , with $r$ the radius common to $s$ and $t$ — is a theorem we use without proof; Hassani (Hassani (2013), p. 467) refers to Birkhoff and Rota for it, and Riley, Hobson and Bence (Riley, Hobson & Bence (2006), p. 539) record the practical rule that the radius is in general the distance to the nearest singularity of the equation. Coddington (Coddington (1989), Ch. 4, §1, p. 144) states the shape of the answer in advance: near a regular singular point the solutions are

(x-x_0)^{r}\sigma(x) + (x-x_0)^{s}\rho(x)\log(x-x_0)

with $\sigma, \rho$ analytic at $x_0$ — which is the content of the trichotomy above.

The hypothesis of a regular singular point is not a technicality. Coddington (Coddington (1989), Ch. 4, §1, Eqs. (1.4)-(1.6), pp. 144-145) exhibits $x^{2}y'' - y' - \tfrac34y = 0$ , where the coefficient $-1$ of $y'$ is not of the form $x\,b_1(x)$ with $b_1$ analytic at $0$ , so the origin is a singular point that is not regular. A formal series $\sum_{k\ge0}c_kx^{k}$ satisfies the recursion $(k+1)c_{k+1} = \bigl(k^{2}-k-\tfrac34\bigr)c_k$ , and since the ratio of consecutive terms is $\bigl|(k^{2}-k-\tfrac34)/(k+1)\bigr||x| \to \infty$ for every $x \ne 0$ , the series converges nowhere but at the origin: a perfectly determinate formal solution that represents nothing. At an irregular singular point one needs asymptotic rather than convergent expansions (Classifying singular points, including infinity).

Method(The Frobenius method)AdaptedAdapted from Riley, Hobson & Bence (2006), §16.3, pp. 538-543; Arfken, Weber & Harris (2013), §7.5, pp. 347-348

Classify the point. Normalise to $y''+py'+qy=0$ and check that $s = xp$ and $t = x^{2}q$ are analytic at $0$ (translate first if the point of interest is not the origin). If they are not, stop: this method need not work.
Multiply through by $x^{2}$ (or by whatever polynomial clears the denominators) and substitute $y = \sum_{n\ge0}a_nx^{n+\sigma}$ , keeping $\sigma$ symbolic. Do not divide by the leading coefficient if that creates infinite series where polynomials would do.
Indicial equation. Set the coefficient of the lowest power of $x$ to zero. With $a_0 \ne 0$ this is $I(\sigma) = \sigma(\sigma-1)+s_0\sigma+t_0 = 0$ ; solve for $\sigma_1, \sigma_2$ and compute $N = \sigma_1-\sigma_2$ (taking $\operatorname{Re}\sigma_1 \ge \operatorname{Re}\sigma_2$ ).
Recurrence. Set the coefficient of the general power to zero, keeping $\sigma$ symbolic. This single recurrence serves both exponents.
First solution. Put $\sigma = \sigma_1$ and generate $a_1, a_2, \dots$ ; this never fails.
Second solution.
- If $N \notin \{0, 1, 2, \dots\}$ : put $\sigma = \sigma_2$ and repeat.
- If $N$ is a positive integer: put $\sigma = \sigma_2$ and run the recurrence; if it survives step $n = N$ , you have a second Frobenius series, otherwise go to the next item.
- If $N = 0$ , or if the recurrence failed: get the second solution either by reduction of order ((Riley, Hobson & Bence (2006), §16.4.1, p. 544) calls this the Wronskian method), or by substituting the ansatz $y_2 = y_1\log x + x^{\sigma_2}\sum_{n\ge1}b_nx^{n}$ and matching coefficients, or by the derivative method of (Riley, Hobson & Bence (2006), §16.4.2, pp. 545-547): keep $\sigma$ as a parameter in the recurrence, so that $L[y(x,\sigma)] = a_0(\sigma-\sigma_1)(\sigma-\sigma_2)x^{\sigma}$ , and take $\partial/\partial\sigma$ at $\sigma = \sigma_1$ .
Identify and check. Try to recognise the series; verify the first few coefficients by substitution; and remember the radius is limited by the next singularity of the equation.

Example(Exponents differing by one half)Original

Solve $2x\,y'' + y' + y = 0$ for $x > 0$ .

Classify. Dividing by $2x$ : $p = 1/(2x)$ , $q = 1/(2x)$ . Then $s = xp = \tfrac12$ and $t = x^{2}q = \tfrac{x}{2}$ , both analytic everywhere, so $x = 0$ is a regular singular point and the series will converge for all $x$ . Here
$s_0 = \tfrac12,\quad s_n = 0\ (n\ge1); \qquad t_0 = 0,\quad t_1 = \tfrac12,\quad t_n = 0\ (n\ge2) .$
Indicial equation. $I(\sigma) = \sigma(\sigma-1) + \tfrac12\sigma = \sigma^{2} - \tfrac12\sigma = \sigma\bigl(\sigma - \tfrac12\bigr)$ , so $\sigma_1 = \tfrac12$ , $\sigma_2 = 0$ . The difference $\tfrac12$ is not an integer, so two Frobenius solutions exist.
Recurrence. Only $s_1 = 0$ and $t_1 = \tfrac12$ are nonzero among the higher coefficients, so the sum on the right has one surviving term, $k = n-1$ :
$I(n+\sigma)\,a_n = -\Bigl[(n-1+\sigma)s_1 + t_1\Bigr]a_{n-1} = -\tfrac12a_{n-1} .$
Larger exponent $\sigma_1 = \tfrac12$ . $I(n+\tfrac12) = (n+\tfrac12)n = n(2n+1)/2$ , so $a_n = -a_{n-1}/\bigl(n(2n+1)\bigr)$ . With $a_0 = 1$ :
$a_1 = -\tfrac13, \quad a_2 = \tfrac{1}{30}, \quad a_3 = -\tfrac{1}{630}, \qquad y_1 = \sqrt{x}\left(1 - \frac{x}{3} + \frac{x^{2}}{30} - \frac{x^{3}}{630} + \cdots\right) .$
Smaller exponent $\sigma_2 = 0$ . $I(n) = n(2n-1)/2$ , so $b_n = -b_{n-1}/\bigl(n(2n-1)\bigr)$ . With $b_0 = 1$ :
$b_1 = -1,\quad b_2 = \tfrac16,\quad b_3 = -\tfrac{1}{90}, \qquad y_2 = 1 - x + \frac{x^{2}}{6} - \frac{x^{3}}{90} + \cdots .$
Closed forms. Put $u = \sqrt{2x}$ . Then $\cos u = 1 - \tfrac{u^{2}}{2} + \tfrac{u^{4}}{24} - \tfrac{u^{6}}{720} = 1 - x + \tfrac{x^{2}}{6} - \tfrac{x^{3}}{90}$ , matching $y_2$ ; and $\tfrac{1}{\sqrt2}\sin u = \sqrt{x}\bigl(1 - \tfrac{x}{3} + \tfrac{x^{2}}{30} - \cdots\bigr)$ , matching $y_1$ . So
$y_2 = \cos\sqrt{2x} , \qquad y_1 = \tfrac{1}{\sqrt2}\sin\sqrt{2x} ,$
and both satisfy $2xy''+y'+y=0$ by direct substitution. The general solution on $x > 0$ is $y = c_1\cos\sqrt{2x} + c_2\sin\sqrt{2x}$ .
Numbers. At $x = \tfrac12$ : $y_2 = \cos 1 = 0.540\,302$ and $y_1 = \sin 1/\sqrt2 = 0.595\,010$ . The four-term partial sums give $1 - 0.5 + 0.041\,667 - 0.001\,389 = 0.540\,278$ and $0.707\,107\,(1 - 0.166\,667 + 0.008\,333 - 0.000\,198) = 0.595\,008$ .
A warning the exponents give for free. Both solutions are bounded at $x = 0$ here, because both exponents are $\ge 0$ . Negative exponents are the usual case and signal a solution that blows up at the singular point; the indicial equation tells you that before any series is computed.

Example(A double exponent and the logarithm)Original

Solve $x\,y'' + y' - y = 0$ for $x > 0$ .

Classify. $p = 1/x$ , $q = -1/x$ , so $s = 1$ and $t = -x$ : a regular singular point at $0$ , with $s_0 = 1$ , $t_0 = 0$ , $t_1 = -1$ and all other coefficients zero.
Indicial equation. $I(\sigma) = \sigma(\sigma-1) + \sigma = \sigma^{2}$ : a double root $\sigma_1 = \sigma_2 = 0$ . By part 3 of the statement, there is exactly one Frobenius solution and the second involves a logarithm.
First solution. The recurrence is $I(n)a_n = -t_1a_{n-1}$ , that is $n^{2}a_n = a_{n-1}$ , so $a_n = a_{n-1}/n^{2}$ and, with $a_0 = 1$ ,
$a_n = \frac{1}{(n!)^{2}} , \qquad y_1(x) = \sum_{n=0}^{\infty}\frac{x^{n}}{(n!)^{2}} = 1 + x + \frac{x^{2}}{4} + \frac{x^{3}}{36} + \cdots ,$
convergent for all $x$ (the ratio of consecutive terms is $x/n^{2} \to 0$ ). This is $I_0(2\sqrt x)$ , the modified Bessel function of order zero.
Second solution: set up the ansatz. Try $y_2 = y_1\log x + v$ with $v = \sum_{n\ge1}b_nx^{n}$ . Writing $L(y) = xy''+y'-y$ and using $L(y_1) = 0$ ,
$L(y_1\log x) = \log x\,L(y_1) + 2y_1' + \frac{y_1}{x} - \frac{y_1}{x} = 2y_1' ,$
because $(y_1\log x)' = y_1'\log x + y_1/x$ and $(y_1\log x)'' = y_1''\log x + 2y_1'/x - y_1/x^{2}$ . So we need $L(v) = -2y_1'$ .
Match coefficients. With $b_0 = 0$ ,
$L(v) = \sum_{n\ge1}n^{2}b_nx^{n-1} - \sum_{n\ge1}b_nx^{n} = \sum_{m\ge0}\Bigl[(m+1)^{2}b_{m+1} - b_m\Bigr]x^{m} , \qquad -2y_1' = -\sum_{m\ge0}2(m+1)a_{m+1}x^{m} ,$
so $(m+1)^{2}b_{m+1} = b_m - 2(m+1)a_{m+1}$ . Hence
$b_1 = -2 , \qquad b_2 = \frac{b_1 - 4a_2}{4} = \frac{-2-1}{4} = -\frac34 , \qquad b_3 = \frac{b_2 - 6a_3}{9} = \frac{-\tfrac34 - \tfrac16}{9} = -\frac{11}{108} ,$
and then $b_4 = -\tfrac{25}{3456}$ . So
$y_2(x) = \left(\sum_{n\ge0}\frac{x^{n}}{(n!)^{2}}\right)\log x \;-\; 2x - \frac34x^{2} - \frac{11}{108}x^{3} - \frac{25}{3456}x^{4} - \cdots .$
Checks. Substituting the truncation of $y_2$ through $x^{7}$ into $xy''+y'-y$ leaves a residue of order $x^{7}$ only, as it must. At $x = 1$ the first solution is $y_1(1) = \sum_n 1/(n!)^{2} = 2.279\,585\,3$ , which is $I_0(2)$ to eight digits. As $x \to 0^{+}$ , $y_1 \to 1$ while $y_2 \to -\infty$ like $\log x$ : the second solution is the one a boundedness condition at the singular point discards.
Why the logarithm had to appear. Step 3 of the proof, applied here: $s_0 = 1$ and $\sigma = 0$ , so $2\sigma + s_0 = 1$ and the reduction-of-order integrand is $y_1^{-2}e^{-\int dx/x} = 1/(x\,y_1^{2})$ , whose residue at $0$ is $1/y_1(0)^{2} = 1$ . Integrating gives $\log x$ , with coefficient $1$ — matching the ansatz of step 4, where $y_1\log x$ enters with coefficient $1$ .

ApplicationOriginal

Bessel's equation and half-integer order. For $w'' + w'/z + (1 - \alpha^{2}/z^{2})w = 0$ one has $s_0 = 1$ , $t_0 = -\alpha^{2}$ , so $I(\nu) = \nu^{2}-\alpha^{2}$ and the exponents are $\pm\alpha$ . Hassani (Hassani (2013), Example 15.2.4(a), p. 466) notes the consequence: there are two Frobenius solutions unless $2\alpha$ is an integer, that is unless $\alpha$ is an integer or a half-integer. For integer $\alpha = n$ the obstruction is genuine and $Y_n$ carries a $\log$ ; for half-integer $\alpha$ the obstruction vanishes and both series survive — which is why the spherical Bessel functions of half-integer order are elementary, $j_0(z) = \sin z/z$ and so on, while $Y_0$ is not.

The Coulomb problem and the two radial behaviours. The radial equation for a $\beta/r$ potential has, in Hassani's (Hassani (2013), Example 15.2.4(b), p. 466) normalisation, exponents $-\tfrac12 \pm \tfrac12\sqrt{1+4\alpha}$ with $\alpha = \ell(\ell+1)$ , hence $\ell$ and $-\ell-1$ , differing by the integer $2\ell+1$ — so only one solution is obtained by series, and the other is the irregular one that diverges at the origin. The indicial equation is thus the formal reason that a bound-state wavefunction behaves like $r^{\ell}$ near the nucleus: the other exponent is excluded by normalisability. The exponents of this radial equation are those of the Euler equation of Euler-Cauchy equations, as they must be.

The hypergeometric equation, and why so many functions are one function. The equation with regular singular points at $0$ , $1$ and $\infty$ and no others is, up to a change of variable, the hypergeometric equation, whose exponents at $0$ are $0$ and $1-\gamma$ ((Hassani (2013), Example 15.2.4(c), p. 467)). Legendre's, Chebyshev's and Jacobi's equations, and in the confluent limit Bessel's, Laguerre's and Hermite's, are all special cases, which is why the handbook tables of special functions are so interrelated: the classification of equations by their singular points and exponents, continued in Classifying singular points, including infinity, explains the coincidences.

Boundary conditions at a singular point. In practice the exponents are often all one needs. A solution of the radial equation behaving like $r^{-\ell-1}$ is rejected at the origin; a solution of Legendre's equation behaving like $\log(1-x)$ is rejected at the pole of a sphere; a solution of the Euler-Bernoulli or diffusion equation growing at a corner is rejected on physical grounds. In each case the indicial equation at a regular singular point delivers the admissible exponents in two lines, long before any series is summed, and that is the usual way the method is used.

Definition

Classifying singular points, including infinity

Adapted

Add to PDF

A point is ordinary, regular singular or irregular singular according to how fast $p$ and $q$ blow up there; the point at infinity is classified by substituting $x = 1/t$ and looking at $t = 0$ ; an equation all of whose singular points in the extended plane are regular is Fuchsian, and with three such points it is the Riemann equation, whose six exponents sum to one.

Prerequisites: The Frobenius method, Euler-Cauchy equations

Definition(Classification of a finite point)AdaptedAdapted from Arfken, Weber & Harris (2013), §7.4, p. 344; Riley, Hobson & Bence (2006), §16.1.1, pp. 533-534; Hassani (2013), §15.2, Def. 15.2.2, p. 464

Write the second-order equation in normalised form

y'' + P(x)\,y' + Q(x)\,y = 0 ,

and let $x_0$ be a point of the (complex) plane. Then $x_0$ is

an ordinary point if $P$ and $Q$ are both analytic at $x_0$ ;
a regular singular point if not, but $(x-x_0)P(x)$ and $(x-x_0)^{2}Q(x)$ are both analytic at $x_0$ — equivalently, if $P$ has at worst a simple pole and $Q$ at worst a pole of order $2$ ;
an irregular singular point, also called an essential singularity of the equation, otherwise.

Arfken, Weber and Harris (Arfken, Weber & Harris (2013), p. 344) phrase the last case by the rates of divergence: $x_0$ is irregular if $P$ diverges faster than $(x-x_0)^{-1}$ , so that $(x-x_0)P(x) \to \infty$ , or if $Q$ diverges faster than $(x-x_0)^{-2}$ .

The classification controls what the solutions look like. At an ordinary point both solutions are analytic (Power-series solutions at ordinary points). At a regular singular point they have the form

(x-x_0)^{r}\sigma(x) + (x-x_0)^{s}\rho(x)\log(x-x_0)

with $\sigma, \rho$ analytic at $x_0$ , and the exponents are the roots of an indicial equation (The Frobenius method). At an irregular singular point no such description holds in general, and the series the method produces may diverge everywhere.

Definition(The point at infinity)AdaptedAdapted from Coddington (1989), Ch. 4, §9, Eq. (9.2) and the criterion, pp. 180-181; Arfken, Weber & Harris (2013), §7.4, Eqs. (7.20)-(7.22), p. 344; Riley, Hobson & Bence (2006), §16.1.1, p. 534

Substitute $x = 1/t$ and set $v(t) = y(1/t)$ . Since $dt/dx = -t^{2}$ ,

\frac{dy}{dx} = -t^{2}\frac{dv}{dt} , \qquad \frac{d^{2}y}{dx^{2}} = t^{4}\frac{d^{2}v}{dt^{2}} + 2t^{3}\frac{dv}{dt} ,

so the equation becomes the induced equation

t^{4}\,v'' + \Bigl[2t^{3} - t^{2}P\bigl(\tfrac1t\bigr)\Bigr]v' + Q\bigl(\tfrac1t\bigr)\,v = 0 ,

or, in normalised form,

v'' + \left[\frac{2}{t} - \frac{P(1/t)}{t^{2}}\right]v' + \frac{Q(1/t)}{t^{4}}\,v = 0 .

The point at infinity is called an ordinary point, a regular singular point or an irregular singular point of the original equation according to what $t = 0$ is for the induced equation. Explicitly, infinity is a regular singular point exactly when

P(x) = \sum_{k\ge1}\frac{a_k}{x^{k}} , \qquad Q(x) = \sum_{k\ge2}\frac{b_k}{x^{k}}

with both series convergent for $|x| > r_0$ ; equivalently, by Coddington (Coddington (1989), Ch. 4, §9, p. 181), exactly when the equation can be written

x^{2}y'' + a(x)\,x\,y' + b(x)\,y = 0

with $a$ and $b$ given by convergent power series in $1/x$ for $|x| > r_0$ . In that case the indicial polynomial at infinity is

I_\infty(\alpha) \;=\; \alpha(\alpha-1) + (2 - a_1)\,\alpha + b_2 ,

and a solution exists of the form

y(x) = x^{-\alpha}\left(1 + \sum_{k\ge1}\frac{C_k}{x^{k}}\right) ,

$\alpha$ a root of $I_\infty$ . Note the sign: an exponent $\alpha$ at infinity describes a solution decaying (or growing) like $x^{-\alpha}$ .

Definition(Fuchsian equations)AdaptedAdapted from Hassani (2013), §15.3, Def. 15.3.1, p. 470

A homogeneous linear differential equation with single-valued analytic coefficients is a Fuchsian differential equation if all of its singular points in the extended complex plane — the plane together with the point at infinity — are regular singular points.

Because a function whose only singularities in the extended plane are poles is a rational function, the coefficients of a Fuchsian equation are rational. The number and position of its regular singular points, together with the exponents at each, then determine the equation almost completely; that is what makes the class worth naming.

Definition(Fuchsian equations with few singular points)AdaptedAdapted from Hassani (2013), §15.3, Thms 15.3.2 and 15.3.3 with Eq. (15.19), pp. 471-472, and Eq. (15.21), p. 473; Coddington (1989), Ch. 4, §9, pp. 181-182

Let the equation $y'' + P y' + Qy = 0$ be Fuchsian.

Two regular singular points. Mapping them to $0$ and $\infty$ by $\xi = (z-z_1)/(z-z_2)$ forces $\Phi(\xi) = a_1/\xi$ and $\Theta(\xi) = b_2/\xi^{2}$ , so the equation becomes
$\xi^{2}u'' + a_1\,\xi u' + b_2\,u = 0 ,$
an Euler equation, equivalent to one with constant coefficients and therefore producing nothing new. (Coddington (Coddington (1989), p. 181) records the same fact from the other side: the Euler equation has $0$ and $\infty$ as its only singular points, both regular.)
Three regular singular points. Mapping them to $0$ , $1$ , $\infty$ by a Mobius transformation, the most general such equation is the Riemann differential equation
$w'' + \left[\frac{A_1}{z} + \frac{B_1}{z-1}\right]w' + \left[\frac{A_2}{z^{2}} + \frac{B_2}{(z-1)^{2}} - \frac{A_3}{z(z-1)}\right]w = 0 ,$
with five constants. Writing the exponent pairs at $0$ , $1$ and $\infty$ as $(\lambda_1,\lambda_2)$ , $(\mu_1,\mu_2)$ and $(\nu_1,\nu_2)$ , the equation is determined by them,
$w'' + \left[\frac{1-\lambda_1-\lambda_2}{z} + \frac{1-\mu_1-\mu_2}{z-1}\right]w' + \left[\frac{\lambda_1\lambda_2}{z^{2}} + \frac{\mu_1\mu_2}{(z-1)^{2}} + \frac{\nu_1\nu_2 - \lambda_1\lambda_2 - \mu_1\mu_2}{z(z-1)}\right]w = 0 ,$
and the six exponents obey the Riemann identity
$\lambda_1 + \lambda_2 + \mu_1 + \mu_2 + \nu_1 + \nu_2 \;=\; 1 .$
Multiplying $w$ by $z^{\lambda}(z-1)^{\mu}$ shifts the exponent pairs by $(\lambda,\lambda)$ , $(\mu,\mu)$ and $(-\lambda-\mu,-\lambda-\mu)$ , which normalises three of the six parameters away and brings the equation to the hypergeometric form
$z(1-z)\,w'' + \bigl[\gamma - (1+\alpha+\beta)z\bigr]w' - \alpha\beta\,w = 0 ,$
whose exponent pairs are $(0, 1-\gamma)$ at $0$ , $(0, \gamma-\alpha-\beta)$ at $1$ and $(\alpha,\beta)$ at $\infty$ .

Remark(How to use the classification)Original

The count of parameters. A Fuchsian equation with three regular singular points has five constants, and three of them can be absorbed; so there is a three-parameter family of genuinely different equations, which is the hypergeometric family $(\alpha,\beta,\gamma)$ . This is why so many of the classical functions are hypergeometric: there is simply nothing else with three regular singular points.
The Riemann identity is a cheap check. Compute the exponents at each singular point independently and add them; if the total is not $1$ , something is wrong — either the equation is not Fuchsian with exactly those three points, or an indicial equation has been miscomputed. The example below uses it this way.
Exponents at infinity carry a sign. A root $\alpha$ of $I_\infty$ corresponds to a solution behaving like $x^{-\alpha}$ at large $|x|$ . A negative exponent at infinity therefore means a growing solution. Confusion between the two conventions is the commonest error in using this material.
Irregular points are not a defect of the method. They are where the interesting asymptotics live. Bessel's equation has a regular singular point at $0$ and an irregular one at infinity, and correspondingly its solutions have convergent Frobenius series near $0$ but only asymptotic expansions, with the oscillatory factor $e^{\pm ix}$ , at large $x$ ((Arfken, Weber & Harris (2013), §7.5, p. 355)). An essential singularity at infinity is exactly what allows exponential behaviour there, which no Fuchsian equation can have.
Confluence. Letting two regular singular points merge produces an irregular one: the confluent hypergeometric equation $zw'' + (c-z)w' - aw = 0$ arises from the hypergeometric equation as $1$ and $\infty$ coalesce, and it has a regular singular point at $0$ and an irregular one at $\infty$ — the entry in (Riley, Hobson & Bence (2006), Table 16.1, p. 535). Bessel's and Laguerre's equations are in the same class, and Hermite's and the oscillator equation have no finite singular point at all and an irregular one at infinity.

Example(Classifying Chebyshev's equation)Original

Take

(1-x^{2})\,y'' - x\,y' + \nu^{2}y = 0 ,

with $\nu$ a constant. Normalised, $P(x) = -x/(1-x^{2})$ and $Q(x) = \nu^{2}/(1-x^{2})$ .

Finite singular points. $P$ and $Q$ are analytic except where $1-x^{2} = 0$ , so the only finite singular points are $x = \pm1$ , and every other point is ordinary.
The point $x = 1$ .
$s(x) = (x-1)P(x) = \frac{-x(x-1)}{(1-x)(1+x)} = \frac{x}{1+x} , \qquad t(x) = (x-1)^{2}Q(x) = \frac{-\nu^{2}(x-1)}{1+x} ,$
both analytic at $1$ , with $s_0 = \tfrac12$ and $t_0 = 0$ . So $x = 1$ is a regular singular point with indicial polynomial $\sigma(\sigma-1) + \tfrac12\sigma = \sigma(\sigma-\tfrac12)$ and exponents
$(\lambda_1, \lambda_2) = \bigl(0, \tfrac12\bigr) .$
The point $x = -1$ . The same computation with $(x+1)$ gives $s_0 = \tfrac12$ , $t_0 = 0$ again, so $x = -1$ is a regular singular point with exponents $(\mu_1,\mu_2) = (0, \tfrac12)$ . (The symmetry $x \mapsto -x$ leaves the equation invariant, which is the quick reason.)
The point at infinity. Expand in powers of $1/x$ :
$P(x) = \frac{x}{x^{2}-1} = \frac{1}{x} + \frac{1}{x^{3}} + \cdots , \qquad Q(x) = \frac{-\nu^{2}}{x^{2}-1} = -\frac{\nu^{2}}{x^{2}} - \frac{\nu^{2}}{x^{4}} - \cdots ,$
both of the required shape, so infinity is a regular singular point with $a_1 = 1$ and $b_2 = -\nu^{2}$ . Hence
$I_\infty(\alpha) = \alpha(\alpha-1) + (2-1)\alpha - \nu^{2} = \alpha^{2} - \nu^{2} , \qquad (\nu_1,\nu_2) = (\nu, -\nu) ,$
and the two solutions behave like $x^{-\nu}$ and $x^{\nu}$ at large $|x|$ .
Fuchsian, with three regular singular points. So Chebyshev's equation is Fuchsian with singular points $-1, 1, \infty$ , all regular — the entry for it in (Riley, Hobson & Bence (2006), Table 16.1, p. 535) — and by the statement above it is a Riemann equation, hence a disguised hypergeometric equation.
The Riemann identity as a check.
$\underbrace{0 + \tfrac12}_{x=1} + \underbrace{0 + \tfrac12}_{x=-1} + \underbrace{\nu + (-\nu)}_{x=\infty} \;=\; 1 ,$
as required, for every $\nu$ .
What this predicts, before any series. The exponent $0$ at each of $x = \pm1$ means one solution is bounded there and (since the exponents differ by $\tfrac12$ , not an integer) the other behaves like $(1\mp x)^{1/2}$ : no logarithms. The exponent $-\nu$ at infinity means a solution growing like $x^{\nu}$ , so for $\nu = n$ a non-negative integer one expects a solution that is a polynomial of degree $n$ — the Chebyshev polynomial $T_n$ . For $n = 2$ , $T_2 = 2x^{2}-1$ and indeed $(1-x^{2})\cdot4 - x\cdot 4x + 4(2x^{2}-1) = 4 - 4x^{2} - 4x^{2} + 8x^{2} - 4 = 0$ , and $T_2 \sim 2x^{2}$ at infinity, matching the exponent $\alpha = -2$ .

Example(Two equations that are not Fuchsian)Original

Hermite's equation $y'' - 2xy' + 2\nu y = 0$ has $P = -2x$ and $Q = 2\nu$ , both entire: there is no finite singular point at all. At infinity, $P$ must be a series in $1/x$ starting at $1/x$ ; instead it grows, and the induced equation has
$\frac{2}{t} - \frac{P(1/t)}{t^{2}} = \frac{2}{t} + \frac{2}{t^{3}} ,$
whose product with $t$ is $2 + 2/t^{2}$ , not analytic at $0$ . So infinity is an irregular singular point, Hermite's equation is not Fuchsian, and nothing prevents its solutions from growing like $e^{x^{2}}$ — which the non-terminating series solutions do, and which is why normalisability quantises $\nu$ .
Bessel's equation $x^{2}y'' + xy' + (x^{2}-\nu^{2})y = 0$ has $P = 1/x$ and $Q = 1 - \nu^{2}/x^{2}$ . At $x = 0$ : $xP = 1$ and $x^{2}Q = x^{2}-\nu^{2}$ are analytic, so the origin is a regular singular point with $s_0 = 1$ , $t_0 = -\nu^{2}$ , indicial polynomial $\sigma^{2}-\nu^{2}$ and exponents $\pm\nu$ . At infinity: $Q$ must be a series in $1/x$ starting at $1/x^{2}$ , but $Q \to 1 \ne 0$ , so infinity is irregular. Bessel's equation is therefore not Fuchsian, and its large- $x$ behaviour $x^{-1/2}e^{\pm ix}$ — oscillation, which no power can produce — is exactly what an irregular singular point permits. Compare Arfken's (Arfken, Weber & Harris (2013), Example 7.4.1, pp. 344-345) classification of the same equation and the Bessel row of (Riley, Hobson & Bence (2006), Table 16.1, p. 535).

ApplicationOriginal

A map of the special functions. The table of standard equations in (Riley, Hobson & Bence (2006), Table 16.1, p. 535) is a classification by singular points: hypergeometric, Legendre, associated Legendre and Chebyshev have regular singular points at $-1, 1, \infty$ (or $0, 1, \infty$ ) and nothing else, and so are all Riemann equations; confluent hypergeometric, Bessel, Laguerre and associated Laguerre have one regular singular point and an irregular one at infinity; Hermite and the oscillator have only the irregular point at infinity. The first group are therefore all the same equation after a change of variable and a factor $z^{\lambda}(z-1)^{\mu}$ , which explains the web of identities relating Legendre, Gegenbauer, Chebyshev and Jacobi functions: they are different normalisations of ${}_2F_1$ .

Choosing boundary conditions. In a physical problem the admissible solution is usually picked out by its behaviour at a singular point, and the exponents settle the matter without computing anything. The radial Schrodinger equation has exponents $\ell$ and $-\ell-1$ at the origin, so regularity selects $r^{\ell}$ ; Legendre's equation has exponents $0$ and $0$ at $x = \pm1$ , so one solution is bounded at the poles of the sphere and the other has the logarithm that forces $Q_\ell$ to be discarded; and the exponent at infinity decides whether a solution is normalisable.

Asymptotics and WKB. The exponent at infinity tells you the power-law behaviour of a solution at large distance when infinity is a regular singular point, which is precisely the information needed to match an interior solution onto a multipole expansion. When infinity is irregular the solutions instead behave like $e^{S(x)}$ with $S$ a power of $x$ , and the leading asymptotics are found by WKB or by the method of dominant balance rather than by series. The classification tells you in advance which of the two situations you are in.

Monodromy and the global theory. For a Fuchsian equation the analytic continuation of a solution around each singular point is a linear substitution, and the group these generate — the monodromy group — is a complete global invariant. This is the beginning of the Riemann-Hilbert problem, and it is why Fuchsian equations, and the hypergeometric equation in particular, recur in conformal field theory, in the theory of Painleve equations and in the study of periods of algebraic varieties.