Part 0 · The Toolkit — Chapter 0.8

Differential Equations and the Oscillator

One equation, solved completely. Then found hiding inside everything else in the book.

Tools you'll need  — Chapter 0.1: the derivative, the chain rule, $\ee$ defined by $\dv{}{t}\ee^{t}=\ee^{t}$, and the licence to separate differentials. Chapter 0.3: Taylor expansion, the small-oscillations result, and Euler's formula $\ee^{\ii\theta}=\cos\theta+\ii\sin\theta$, which we use constantly. Chapter 0.4: vector spaces, basis, dimension, linear maps. Chapter 0.5: eigenvalues and eigenvectors, the spectral theorem for symmetric matrices, and the warning about defective matrices. Sections §3 and §7 each cash one of those. Chapter 0.6: the Hessian and its eigenvalues as principal curvatures.

1 · What a differential equation is

An ordinary differential equation (ODE) is a relation among an unknown function $y(t)$, its derivatives, and $t$ itself. "Ordinary" means one independent variable. When there are several, as for a field depending on $x,y,z,t$, you get a partial differential equation instead. That is what §7 will build for you at the end.

Three words classify almost everything we will do.

Order is the highest derivative appearing. Newton's law is second order, because force determines acceleration. That one fact is why you always get two constants of integration in mechanics. It is also why "where is it and how fast is it going" is the right amount of information to specify a trajectory.

Linear means the unknown and its derivatives appear only to the first power, never multiplied together, never inside a $\sin$ or a square root. The general linear ODE of order $n$ is

The coefficients $a_i(t)$ may do whatever they like. What matters is that $y$ enters linearly. The point of writing it as $L[y]$ is that $L$ is a linear map in the sense of Chapter 0.4, meaning $L[\alpha y_1 + \beta y_2] = \alpha L[y_1] + \beta L[y_2]$. That holds because differentiation is linear, and so is multiplication by a fixed function. Everything in §3 follows from that single observation.

Homogeneous means $f\equiv0$, so there is no source term and nothing driving the system from outside. An inhomogeneous equation has a source, and §2 and §6 are about what a source does.

1.1 · Initial conditions, and the mathematical content of determinism

An ODE by itself has infinitely many solutions. To pick one you supply initial conditions: the value of $y$ and of its first $n-1$ derivatives at one instant $t_0$. That count is not a convention. It is exactly what the following theorem says is needed, and exactly what that same theorem says is enough.

Now let's read the theorem as a physical statement rather than an analytic one. It says: the state now determines the state later, uniquely. That is determinism. It is not a philosophical posture bolted onto physics. It is a theorem about the equations, given that the equations are of this form.

The "state" of a classical system is exactly the list of numbers the theorem requires, which for mechanics is position and velocity. That list is what Chapter 1.3 will call a point of phase space.

There are two places where this gets interesting. Both are flagged now and paid later.

• In §8 we meet systems that are perfectly deterministic and still unpredictable. The solution exists and is unique, but neighbouring initial conditions separate exponentially. A finite-precision measurement of the present therefore buys you only a finite time of usable future. That is chaos, and it is a statement about sensitivity rather than about determinism.
• In quantum mechanics the determinism survives intact, in a place people rarely look. The Schrödinger equation $\ii\hbar\,\partial_t\ket\psi = \hat H\ket\psi$ is a first-order linear ODE for the state, so the state evolves as deterministically as a planet. The randomness enters only at measurement, and Chapter 4.20 is about exactly where that seam is.

2 · First-order linear equations

Let's start where the mathematics is easiest and the physics is already familiar. Everything in this section is about equations with a single derivative in them, and the whole point is to see how a source term gets absorbed.

2.1 · Separation of variables, done honestly

Suppose the equation has the form $\dot y = g(y)\,h(t)$, so that the right-hand side factorises into a piece depending only on $y$ and a piece depending only on $t$. Then it can be integrated directly.

The usual manipulation is to write $\dd y/g(y) = h(t)\,\dd t$ and integrate both sides, which is the "deliberate abuse" that Chapter 0.1 warned you was legitimate. Here is why it is legitimate, in one line. Let $G$ be an antiderivative of $1/g$. Then by the chain rule,

so $G(y(t)) = \int h(t)\,\dd t + C$. No differentials were harmed. Separation of variables is the chain rule read from right to left.

Here is the example you already know. The equation $\dot y = -ky$ has $g(y)=y$ and $h(t)=-k$, so $G(y)=\ln\abs y$ and $\ln\abs y = -kt + C$, giving $y = y_0\ee^{-kt}$. Chapter 0.1 got this by recognising the defining property of $\ee$. We now have it as a special case of a method.

2.2 · The integrating factor, derived

Separation fails the moment there is a source. The equation $\dot y = -ky + f(t)$ does not factorise, so the method of §2.1 has nothing to bite on. The trick that works instead is one of the most quietly clever manoeuvres in elementary mathematics, and it is worth deriving rather than remembering.

Take the general first-order linear equation

We would like the left-hand side to be a single derivative, because then we could integrate it. As it stands, it is not one. So suppose we multiply the whole equation by some function $\mu(t)$, and then ask what $\mu$ would have to be for the result to come out as exactly $\dv{}{t}(\mu y)$. To answer that, expand the target with the product rule and set it beside what multiplying actually gives us:

The two agree for every $y$ precisely when $\dot\mu = p\,\mu$, and that is the separable equation we solved a moment ago. So $\mu$ is not chosen. It is forced on us:

Nothing was guessed. We asked for the function that makes the left side a total derivative and the demand determined it. Multiplying (0.8.3) by $\mu$ and integrating from $t_0$ to $t$ now gives the general solution outright:

2.3 · Exponential decay with a source — and a Green's function in embryo

Now let's take the case that matters most in practice, which is a decay rate that does not change with time. Specialise to constant $p=k$, so that $\mu=\ee^{kt}$, and set $t_0=0$. The general solution collapses to something worth staring at:

Stare at the integral. It says: the source deposited an amount $f(s)\,\dd s$ into the system at time $s$, that deposit has been decaying ever since for a duration $t-s$, and the present value is the sum over all past deposits of what is left of each. The system is a memory with a fading kernel, and the kernel is the response to a single sharp kick.

This is called the Duhamel form, and the function

is a Green's function. It is the response at time $t$ to a unit impulse at time $s$, and the second clause is what enforces that causes precede their effects.

Chapter 0.9 will give the impulse a name ($\delta$), prove that (0.8.7) is a convolution, and show that convolutions become ordinary products under Fourier transform. At that point, solving a linear differential equation becomes division. Everything in this book that computes a response to a disturbance, up to and including the propagators of Chapter 5.4, is (0.8.7) with more indices.

Grind box — checking the solution, and one nonlinear equation solved exactly

Verification. Never trust a derived formula you have not differentiated back. Differentiate (0.8.7), using the fundamental theorem of calculus on the upper limit and the chain rule on the $t$ inside the integrand:

$$\begin{aligned} \dot y &= -k\,y(0)\ee^{-kt} + \underbrace{\ee^{-k(t-t)}f(t)}_{\text{upper limit}} \;-\; k\int_{0}^{t}\ee^{-k(t-s)}f(s)\,\dd s\\[3pt] &= f(t) - k\left[y(0)\ee^{-kt} + \int_{0}^{t}\ee^{-k(t-s)}f(s)\dd s\right] \;=\; f(t)-k\,y. \end{aligned}$$

So $\dot y + ky = f$, and $y(0)=y(0)$ because the integral is empty at $t=0$. ✓

A constant infusion. Put $f(t)=R$, a steady rate. Then $\int_0^t\ee^{-k(t-s)}R\,\dd s = \tfrac{R}{k}\big(1-\ee^{-kt}\big)$, so

$$y(t) = y(0)\ee^{-kt} + \frac{R}{k}\left(1-\ee^{-kt}\right) \;\longrightarrow\; \frac{R}{k}.$$

The steady state is rate in over rate out, approached with time constant $1/k$. So the time to reach steady state is set by the elimination constant and not by the infusion rate, which sets only the level. This is why a loading dose exists as a concept.

A nonlinear equation that still yields. Separation does not care about linearity. The logistic equation $\dot y = ry(1-y/K)$ separates:

$$\int\frac{\dd y}{y(1-y/K)} = \int r\,\dd t \quad\Longrightarrow\quad \ln\frac{y}{K-y} = rt + C,$$

using the partial fractions $\frac{1}{y(1-y/K)} = \frac1y + \frac{1/K}{1-y/K}$. Solving for $y$ gives $y(t) = K/\big(1+A\ee^{-rt}\big)$ with $A=(K-y_0)/y_0$. Note what we did not get: any right to add two solutions. Superposition is a linear privilege, and §8 is about how much that costs.

3 · A linear ODE is an eigenvalue problem

This section is the one to read twice. It contains no computation you could not do without it, and it changes what all the computations mean.

3.1 · The solution set is a vector space, and its dimension is the order

Take the homogeneous equation $L[y]=0$ with $L$ as in (0.8.1). Because $L$ is linear, if $L[y_1]=0$ and $L[y_2]=0$ then

Any combination of solutions is a solution, and $y\equiv0$ is a solution. Chapter 0.4 already identified what that makes it. The solution set of a homogeneous linear ODE is a vector space, namely the kernel of the linear map $L$. That is the structural reason superposition works. It is worth noticing that it has nothing to do with the physics of waves and everything to do with the algebra of $L$.

Now the dimension, which is where the quoted theorem of §1 gets spent. Consider the map that takes a solution to its initial data:

$\Phi$ is linear, because evaluating a derivative at a point is a linear operation. It is also both injective and surjective, and both halves of Picard–Lindelöf get used saying so.

• Injective. If two solutions share all $n$ initial values, then by the uniqueness half of the theorem they are the same solution, so $\Phi(y)=0$ forces $y=0$.
• Surjective. Given any list of $n$ numbers, the existence half produces a solution with those initial values.

So $\Phi$ is an isomorphism of vector spaces, and by Chapter 0.4 isomorphic spaces have equal dimension:

Read what has happened to a phrase you have used since school. "The general solution contains two arbitrary constants" is a vague and slightly apologetic statement. "The solution space is two-dimensional" is a precise one, and it tells you exactly what to do: find a basis. Two independent solutions is not a lucky haul, it is a complete answer, and the "arbitrary constants" are coordinates in that basis. Chapter 0.4 proved this for $\ddot y+\omega^{2}y=0$ by hand. We now have it for every linear ODE at once.

3.2 · The exponential is the eigenfunction of the derivative

Here is the punchline. Let $D=\dd/\dd t$, regarded as a linear map on the space of smooth functions. That space is an entirely legitimate vector space, as Chapter 0.4 insisted. Now put Chapter 0.5's question to it: which vectors does this map merely rescale? Written out, the question is

and we solved that equation back in §2. The answer is $y=C\ee^{\lambda t}$, and nothing else is an answer. That single fact deserves a box.

3.3 · The characteristic equation, from the eigenvalue view

Now specialise to constant coefficients, which is the case where diagonalising $D$ actually solves the problem. Write the equation using $D$:

where $p$ is an ordinary polynomial with the constants as coefficients. Feed an eigenfunction in. Since $D\ee^{\lambda t} = \lambda\ee^{\lambda t}$, applying $D$ again gives $D^{2}\ee^{\lambda t}=\lambda^{2}\ee^{\lambda t}$, and inductively $D^{m}\ee^{\lambda t}=\lambda^{m}\ee^{\lambda t}$. So every term of $p(D)$ acts on $\ee^{\lambda t}$ by multiplication, and

An exponential is a solution of $p(D)y=0$ exactly when $p(\lambda)=0$. That is the characteristic equation, and it did not have to be introduced as a recipe. It is (0.8.14), which is Chapter 0.5's statement that a polynomial in an operator acts on an eigenvector by the same polynomial in the eigenvalue. Compare (0.8.15) below with the matrix characteristic polynomial of Chapter 0.5, and you are looking at the same object.

Over $\C$ this degree-$n$ polynomial has $n$ roots counted with multiplicity, by the fact recorded in Chapter 0.4. Suppose for the moment that all $n$ of them, $\lambda_1,\dots,\lambda_n$, are distinct. Then the $n$ functions $\ee^{\lambda_i t}$ are independent, and by (0.8.11) they are a basis. That finishes the job. The general solution is $y=\sum_i c_i\ee^{\lambda_i t}$.

Real equations, complex roots. Physics equations have real coefficients, so complex roots come in conjugate pairs $\lambda = a\pm\ii b$. The corresponding solutions $\ee^{at}\ee^{\pm\ii bt}$ are complex, but their combinations

are real, by Euler's formula from Chapter 0.3. Since they span the same two-dimensional space, you may use either basis. A complex eigenvalue with non-zero imaginary part means oscillation, and its real part is growth or decay. That single sentence classifies the behaviour of every linear system in this book, and §5 is nothing but it applied once.

3.4 · Repeated roots: the defective case, and where $t\ee^{\lambda t}$ comes from

Suppose $p$ has a double root: $p(\lambda)=(\lambda-\lambda_0)^{2}q(\lambda)$ with $q(\lambda_0)\neq0$. Then $\ee^{\lambda_0 t}$ is a solution, but it is only one solution, and (0.8.11) insists there must be two associated with this root. The eigenvector supply has run short.

You have seen this before, and it was flagged as a warning. Chapter 0.5's warn box exhibited

a $2\times2$ matrix with a repeated eigenvalue and only a one-dimensional eigenspace. That is a defective matrix, with no basis of eigenvectors and no diagonal form. The correspondence about to appear is exact rather than an analogy.

Write $N = D-\lambda_0$, and restrict attention to the two-dimensional space $V=\ker\!\big((D-\lambda_0)^{2}\big)$ of solutions belonging to this root. On $V$ we have $N^{2}=0$ by construction, and $N\neq0$ (or the eigenspace would be two-dimensional and $\ee^{\lambda_0t}$ would not be alone). An operator with $N^2=0\neq N$ is nilpotent, and it cannot be diagonalised: a diagonal $N$ with $N^2=0$ would have all diagonal entries zero, hence be $0$.

What does $V$ contain? Take any $w\in V$ with $Nw\neq0$. Then $N(Nw)=N^{2}w=0$, so $Nw$ is an eigenvector, which means it is a multiple of $\ee^{\lambda_0t}$. What we are hunting for, then, is a function that $N$ maps to $\ee^{\lambda_0 t}$ rather than to zero. Try $w=t\ee^{\lambda_0t}$ and compute, using the product rule:

That is exactly what was needed. Applying $N$ once more kills it, so $(D-\lambda_0)^{2}\big(t\ee^{\lambda_0t}\big)=0$, and $t\ee^{\lambda_0t}$ is the second solution. It was derived rather than guessed.

Now let's see what $D$ looks like in the basis $\{u,w\}=\{\ee^{\lambda_0t},\,t\ee^{\lambda_0t}\}$. It acts by $Du=\lambda_0u$ and $Dw = \lambda_0 w + u$, whose matrix is

So (0.8.17) is not a curiosity a linear-algebra text invented to spoil the spectral theorem. It is what the derivative operator looks like at a repeated root, and the extra factor of $t$ in $t\ee^{\lambda_0 t}$ is the off-diagonal $1$. The vector $w$ has a name in linear algebra. It is a generalised eigenvector, and now you know what one is: a function that the operator almost, but not quite, merely rescales.

Grind box — higher multiplicity, reduction of order, and the confluent limit

Multiplicity $m$. The same argument iterates. Since $N\big(t^{j}\ee^{\lambda_0t}\big) = j\,t^{j-1}\ee^{\lambda_0t}$ by the identical product-rule computation, applying $N$ repeatedly walks down the ladder $t^{m-1}\to t^{m-2}\to\cdots\to 1\to 0$. Hence $N^{m}$ annihilates $t^{j}\ee^{\lambda_0t}$ for every $j\lt m$, and

$$\ee^{\lambda_0t},\ t\ee^{\lambda_0t},\ t^{2}\ee^{\lambda_0t},\ \dots,\ t^{m-1}\ee^{\lambda_0t}$$

are $m$ independent solutions. They are independent because a vanishing combination would make a polynomial of degree $\lt m$ identically zero. That is exactly the $m$ the dimension count demands, and the matrix of $D$ on this space is an $m\times m$ Jordan block.

Reduction of order, if you prefer no operators. Knowing one solution $y_1=\ee^{\lambda_0t}$ of $\ddot y - 2\lambda_0\dot y+\lambda_0^{2}y=0$, substitute $y = v(t)\,\ee^{\lambda_0t}$. Then $\dot y = (\dot v+\lambda_0v)\ee^{\lambda_0t}$ and $\ddot y = (\ddot v + 2\lambda_0\dot v + \lambda_0^{2}v)\ee^{\lambda_0t}$, and substituting, every term without a derivative of $v$ cancels:

$$\big(\ddot v + 2\lambda_0\dot v+\lambda_0^{2}v\big) - 2\lambda_0\big(\dot v + \lambda_0 v\big) + \lambda_0^{2}v = \ddot v = 0.$$

So $v = A+Bt$, recovering $y=(A+Bt)\ee^{\lambda_0t}$. The cancellation is not luck. It is (0.8.14) saying that the coefficient of $v$ is $p(\lambda_0)$ and the coefficient of $\dot v$ is $p'(\lambda_0)$, and both of those vanish at a double root.

The confluent limit, which is the most physical derivation. Take two distinct roots $\lambda_1\neq\lambda_2$. The pair $\{\ee^{\lambda_1t},\ee^{\lambda_2t}\}$ is a basis, but so is the pair $\{\ee^{\lambda_1t},\,(\ee^{\lambda_2t}-\ee^{\lambda_1t})/(\lambda_2-\lambda_1)\}$, since the second is a combination of the first two with non-zero coefficient. Now let $\lambda_2\to\lambda_1$. The second basis vector is a difference quotient in $\lambda$, and by the definition of the derivative (Chapter 0.1),

$$\frac{\ee^{\lambda_2t}-\ee^{\lambda_1t}}{\lambda_2-\lambda_1} \;\longrightarrow\; \left.\pdv{}{\lambda}\ee^{\lambda t}\right|_{\lambda_1} \;=\; t\,\ee^{\lambda_1 t}.$$

The second solution does not appear from nowhere when the roots collide. It is the derivative with respect to the eigenvalue, and it was there all along as a finite difference of the two exponentials. This is the same manoeuvre that will produce the secular term $t\sin\omega_0t$ at exact resonance in §6, and it happens for the same reason. Two things became equal, and the difference between them had to be rescaled to survive.

4 · The harmonic oscillator

Now we reach the equation the whole chapter is built around. Take a mass $m$ on a spring of stiffness $k$, with no friction and no driving. Hooke's law says the restoring force is $-kx$, so Newton gives $m\ddot x = -kx$, which we tidy into

By Chapter 0.3 this is not a statement about springs. Expand any potential $V$ about a stable minimum at $x_0$. The constant is irrelevant, the linear term vanishes because $V'(x_0)=0$, and what survives is $V\approx V(x_0)+\tfrac12V''(x_0)(x-x_0)^{2}$. That is a spring with $k=V''(x_0)$.

So (0.8.20) governs a pendulum, a diatomic molecule, an atom in a crystal, a mode of the electromagnetic field, and a small perturbation of a black hole. It is worth solving once, thoroughly, and that is what this section does.

4.1 · First way: real exponentials that are not real

Let's apply §3 straight off. The characteristic equation is $\lambda^{2}+\omega^{2}=0$, so $\lambda=\pm\ii\omega$. That is a conjugate pair with zero real part, which by the rule at the end of §3.3 means oscillation with no growth or decay attached. By (0.8.16) a real basis is $\{\cos\omega t,\ \sin\omega t\}$, so

The second form is the first one rewritten. Expanding $A\cos(\omega t+\varphi) = A\cos\varphi\cos\omega t - A\sin\varphi\sin\omega t$ identifies $C_1 = A\cos\varphi$ and $C_2=-A\sin\varphi$, and so

There are two constants either way, as (0.8.11) demands. It is one two-dimensional space seen in two coordinate systems. Fitting initial conditions $x(0)=x_0$ and $\dot x(0)=v_0$ gives $C_1=x_0$ and $C_2=v_0/\omega$.

Now let's look at what is not in (0.8.21). The frequency does not depend on the amplitude. A pendulum clock keeps time whether it swings wide or narrow, and §8 is about the fact that this is only approximately true.

4.2 · Second way: the complex amplitude

Physicists rarely write (0.8.21). They write

which is the same thing, because $\mathrm{Re}[A\ee^{\ii\varphi}\ee^{\ii\omega t}] = A\cos(\omega t+\varphi)$ by Euler. The gain is bookkeeping, and it is large. The single complex number $Z$ carries both the amplitude (its modulus) and the phase (its argument). Three operations get cheaper at once.

• Adding two oscillations of the same frequency becomes adding two complex numbers.
• Shifting a phase becomes multiplying by $\ee^{\ii\theta}$, which Chapter 0.3 showed is a rotation.
• Differentiating becomes multiplying by $\ii\omega$.

The last of those is the one that matters. The derivative has been turned into an algebraic operation, which is exactly (0.8.12) being used as a labour-saving device.

The move is legal for one reason, and the reason is worth stating out loud. The equation is linear with real coefficients, so if $z(t)$ solves it then so does $\bar z$, and hence so does $\mathrm{Re}\,z=(z+\bar z)/2$. You may solve the complex problem and take the real part at the end.

The moment the equation is nonlinear this fails, because $\mathrm{Re}(z^{2})\neq(\mathrm{Re}\,z)^{2}$. That is one more entry on §8's list.

4.3 · Third way: energy, and the first integral

The third route does not solve the equation at all. It reduces the order, and it generalises further than either of the others. We are after a combination whose time derivative vanishes, so multiply (0.8.20) by $m\dot x$ and see what appears:

because each term is the derivative of a square by the chain rule. So the bracket is constant:

A quantity that is constant along solutions is a first integral, and finding one converts a second-order equation into a first-order one. Here it does more than that. Solve (0.8.25) for $\dot x$ and you get a separable equation, and integrating that reproduces (0.8.21).

Chapter 1.4 will explain where first integrals come from in general. Every one of them is a symmetry, and the symmetry behind this one is invariance under shifting the clock.

4.4 · Phase space: the trajectory is a closed ellipse

Plot the state $(x,p)$ with $p=m\dot x$ rather than plotting $x$ against $t$. Rewriting (0.8.25) in these variables,

which is an ellipse with semi-axes $x_{\max}=A=\sqrt{2E/(m\omega^{2})}$ and $p_{\max}=m\omega A$. Every solution runs around it, clockwise. Rescale to the coordinates $(\omega x,\ \dot x)$ and the ellipse becomes a circle traversed at uniform angular rate, so the motion is literally a rotation. That picture is (0.8.23) drawn on the page.

That the curve closes is the whole content of energy conservation. An orbit that failed to close would return to the same position with a different speed, and $E$ would have changed. So "the trajectory is a closed loop in phase space" and "energy is conserved" are the same statement, one drawn and one written.

Let's compute the area the ellipse encloses, because the answer will matter a great deal later. The area of an ellipse is $\pi$ times the product of the semi-axes:

using $E=\tfrac12m\omega^{2}A^{2}$ and $T=2\pi/\omega$. Phase-space area is energy times period. Hold onto that, because Chapter 4.8 will find that a quantum oscillator has $E_n = (n+\tfrac12)\hbar\omega$, whence

The allowed orbits are those enclosing an integer-and-a-half multiple of Planck's constant. The classical ellipse is not replaced in quantum mechanics. It is rationed.

Chapter 4.8 derives this properly with ladder operators. Chapter 1.3 will show that phase-space area is preserved by all of classical mechanics rather than just by this one example, and that is why it was the right thing to quantise.

5 · Damping

Nothing oscillates forever. Let's add a resistive force proportional to velocity. It might be viscous drag, an electrical resistance, or the radiation of energy away to infinity. Newton's law then reads $m\ddot x = -kx - b\dot x$. Divide by $m$ and name the two constants:

The factor of $2$ in the definition of $\gamma$ is chosen purely to make what follows tidy, and you should expect it to earn its keep. Section 3 does all the work from here. The characteristic equation is $\lambda^{2}+2\gamma\lambda+\omega_0^{2}=0$, so its roots are

Everything about damped motion is in the discriminant $\gamma^{2}-\omega_0^{2}$, and there are exactly three cases because a real number is negative, zero, or positive.

5.1 · Underdamped: $\gamma \lt \omega_0$

The discriminant is negative, so the roots are a complex conjugate pair $\lambda_\pm = -\gamma\pm\ii\omega_{\mathrm d}$ with

Read it off (0.8.16). The real part $-\gamma$ is the decay and the imaginary part $\omega_{\mathrm d}$ is the oscillation. The system rings, and the ringing dies away with time constant $1/\gamma$.

Two details here are worth noticing.

• Damping lowers the frequency, since $\omega_{\mathrm d}\lt\omega_0$. The system is being dragged, so it takes longer to come back.
• The effect is second-order small, because $\omega_{\mathrm d}\approx\omega_0\big(1-\tfrac12\gamma^{2}/\omega_0^{2}\big)$. A system that loses a per cent of its amplitude per radian shifts its frequency by only a part in $2\times10^{4}$.

So weak damping kills the amplitude long before it moves the pitch.

5.2 · Overdamped: $\gamma \gt \omega_0$

The discriminant is positive, both roots are real and negative, and the motion is a sum of two decaying exponentials with no oscillation at all:

Both are negative because $\sqrt{\gamma^{2}-\omega_0^{2}}\lt\gamma$. At long times the slower of the two dominates, and that is $\lambda_+$, the root closer to zero. Now comes the counter-intuitive part. We want to know how $\lambda_+$ behaves under heavy damping, $\gamma\gg\omega_0$, so expand the square root using Chapter 0.3:

More damping makes the return to equilibrium slower, not faster. A door closer filled with treacle takes forever. The rate is inversely proportional to the drag, exactly as in the overdamped limit of any relaxation problem.

5.3 · Critically damped: $\gamma = \omega_0$

The discriminant vanishes and $\lambda=-\gamma$ is a repeated root. This is precisely the defective case of §3.4, and the second solution is forced to be (0.8.18):

So the extra factor of $t$ in a critically damped shock absorber and the off-diagonal $1$ in the Jordan block (0.8.17) are the same fact.

Critical damping is the boundary between ringing and crawling, and it is the fastest possible approach to equilibrium without overshoot. That is why car suspensions, galvanometers and control loops are tuned to it.

You can read the "fastest" claim off (0.8.33). The decay exponent here is $\gamma=\omega_0$. It is larger in magnitude than the $\omega_0^{2}/2\gamma$ available on the overdamped side, and larger than the envelope $\ee^{-\gamma t}$ available on the underdamped side, where $\gamma\lt\omega_0$.

5.4 · The quality factor

The dimensionless combination that controls everything is

and its meaning is energetic. Energy in an oscillator is quadratic in amplitude, and the amplitude envelope decays as $\ee^{-\gamma t}$, so

In time $t$ the oscillator has advanced very nearly $\omega_0 t$ radians of phase, and the exponent is that number divided by $Q$. Hence, for any oscillator worth calling one:

$Q$ is the number of radians of oscillation in which the energy falls by a factor of $\ee$.

Two approximations are buried in that sentence, and both are harmless for any oscillator you would call resonant. The phase actually advances at $\omega_{\mathrm d}$ rather than $\omega_0$, so the true radian count is $Q\sqrt{1-1/4Q^{2}}$, and $E(t)=E_0\ee^{-2\gamma t}$ is the cycle-averaged energy rather than the instantaneous one, which ripples about that envelope. At $Q=100$ the first correction is one part in $40{,}000$. At $Q=1$ it is $13\%$, and by then you no longer have an oscillator.

Divide by $2\pi$ and you have the number of cycles instead, so $Q$ is that cycle count times $2\pi$. A system with $Q=100$ rings for about $16$ cycles before losing $63\%$ of its energy.

You will often see the statement that "the fractional energy lost per cycle is $2\pi/Q$". That is the small-loss expansion of (0.8.36). Over one period $T=2\pi/\omega_0$ the energy ratio is $\ee^{-2\pi/Q}\approx 1-2\pi/Q$, which is accurate to a per cent only for $Q\gtrsim 300$. Use the exponential instead. It is no harder.

Section 6 will produce a second and entirely independent meaning for the same number, namely that $Q$ is the sharpness of the resonance peak. The fact that one number does both jobs is the width–lifetime relation that runs through the rest of this book.

Grind box — all three cases with the same initial conditions, and the seam between them

Release from rest: $x(0)=x_0$, $\dot x(0)=0$. Impose these on each case.

Underdamped. Write (0.8.31) as $x=\ee^{-\gamma t}(A\cos\omega_{\mathrm d}t + B\sin\omega_{\mathrm d}t)$. At $t=0$ that gives $A=x_0$. Differentiating, $\dot x(0) = -\gamma A + \omega_{\mathrm d}B = 0$ gives $B = \gamma x_0/\omega_{\mathrm d}$, so

$$x(t) = x_0\,\ee^{-\gamma t}\left[\cos\omega_{\mathrm d}t + \frac{\gamma}{\omega_{\mathrm d}}\sin\omega_{\mathrm d}t\right].$$

Critically damped. From (0.8.34) we get $A=x_0$, and $\dot x(0) = B-\gamma A=0$ gives $B=\gamma x_0$, so $x(t) = x_0\ee^{-\gamma t}(1+\gamma t)$.

The seam is continuous. Let $\omega_{\mathrm d}\to0$ in the underdamped answer. Then $\cos\omega_{\mathrm d}t\to1$ and $\sin(\omega_{\mathrm d}t)/\omega_{\mathrm d}\to t$, so the bracket tends to $1+\gamma t$, which is exactly the critical solution. Nothing discontinuous happens at $\gamma=\omega_0$. What happens is that the oscillation has its first zero crossing pushed out to $t=\pi/\omega_{\mathrm d}\to\infty$. This is the confluent limit of the grind box in §3 seen in physical clothing, and it is why the repeated-root case is not a separate phenomenon.

Overdamped. Same algebra with $\omega_{\mathrm d}\to\ii\kappa$, $\kappa=\sqrt{\gamma^{2}-\omega_0^{2}}$, turns the trigonometric functions into hyperbolic ones:

$$x(t) = x_0\,\ee^{-\gamma t}\left[\cosh\kappa t + \frac{\gamma}{\kappa}\sinh\kappa t\right],$$

which is positive for all $t\gt0$: an overdamped system released from rest never crosses zero. The critically damped case is the last one for which that is true, and the underdamped case crosses first at $\omega_{\mathrm d}t = \pi-\arctan(\gamma/\omega_{\mathrm d})$.

6 · Driving, resonance, and the shape of every spectral line

Now push on it. Take a sinusoidal drive of angular frequency $\omega$, which need not be $\omega_0$, and write the force per unit mass as $F_0\cos\omega t$:

This is inhomogeneous, and Chapter 0.4 already told us the structure of the answer: general solution = one particular solution + the whole vector space of homogeneous solutions. The homogeneous part is §5, and it decays as $\ee^{-\gamma t}$. So after a few times $1/\gamma$ the system forgets how it was started and settles onto the particular solution. That survivor is the steady state. Hold on to one condition attached to it, because it will matter in a moment: it exists as a well-defined attractor only because $\gamma\gt0$.

6.1 · Solving it with a complex amplitude

Let's solve the complex equation $\ddot z + 2\gamma\dot z+\omega_0^{2}z = F_0\ee^{\ii\omega t}$ and take the real part at the end. That is legal because (0.8.37) is linear with real coefficients.

Try $z = Z\ee^{\ii\omega t}$ with $Z\in\C$ constant. Each derivative brings down a factor $\ii\omega$, which is (0.8.12) earning its keep, since $\ee^{\ii\omega t}$ is an eigenfunction of $D$. So the differential equation collapses to an algebraic one:

That is the entire solution. To read it in physical terms, write $Z = A\ee^{-\ii\delta}$, so that $x(t)=A\cos(\omega t-\delta)$ and $\delta$ is the amount by which the response lags the drive. Now take the modulus and the argument of (0.8.38):

The branch is fixed by requiring $\delta$ to be continuous and to start at $0$: the imaginary part of the denominator of $Z$ is $2\gamma\omega\gt0$, so $\delta$ is always in the upper half range.

6.2 · Where the amplitude peaks — and it is not at $\omega_0$

We want the frequency at which $A$ is largest, so we minimise the quantity under the root. Write $s=\omega^{2}$, which makes the job minimising $g(s) = (\omega_0^{2}-s)^{2}+4\gamma^{2}s$. Setting $g'(s) = -2(\omega_0^{2}-s)+4\gamma^{2}=0$ gives $s = \omega_0^{2}-2\gamma^{2}$, and so

with the amplitude following when we substitute $\omega_0^{2}-\omega_{\text{peak}}^{2}=2\gamma^{2}$ back into (0.8.39). Three things fall out of that.

• The peak sits at neither the natural frequency $\omega_0$ nor the damped frequency $\omega_{\mathrm d}=\sqrt{\omega_0^{2}-\gamma^{2}}$. It sits at a third value, lower than both.
• It exists at all only when $\omega_0^{2}\gt2\gamma^{2}$, which is to say $Q\gt1/\sqrt2$. Below that the response falls monotonically and there is no resonance.
• As $\gamma\to0$ the peak height diverges like $F_0/(2\gamma\omega_0)$. That divergence is what the warning box below is about.

6.3 · Power absorbed, and the exact width

Amplitude is not the physically cleanest quantity. Power is, because power is what a spectrometer measures. The instantaneous power delivered by the drive is force times velocity. We have $x = A\cos(\omega t-\delta)$, so $\dot x = -A\omega\sin(\omega t-\delta)$. Average their product over a cycle, using $\avg{\cos^{2}}=\tfrac12$ and $\avg{\sin\cos}=0$:

Look at what survived that average. Only the component of the response in phase with the velocity absorbs energy, and that component is the $\sin\delta$. This is why phase matters physically rather than only as bookkeeping. From (0.8.38) we have $\sin\delta = 2\gamma\omega A/F_0$, so

Now let's find where this curve peaks. The way to see it is to get $\omega$ into as few places as possible, so divide numerator and denominator by $\omega^{2}$:

The denominator is smallest exactly when the first bracket vanishes, and that happens at $\omega=\omega_0$ precisely, with $\bar P(\omega_0) = F_0^{2}/4\gamma$. The power resonance sits exactly at the natural frequency, even though the amplitude resonance does not.

The half-power points come out just as clean. Setting $\bar P = \tfrac12\bar P(\omega_0)$ requires the bracket to equal $\pm2\gamma$, so

taking the positive roots. We want the width of the peak, which is the distance between those two points, so subtract them. The square roots cancel, and

That is not approximately $2\gamma$. It is exactly $2\gamma$, for every $\gamma$, and the resonant frequency is recoverable as the geometric mean of the half-power points. Combining this with (0.8.35),

The same $Q$ that counted radians of ringdown in (0.8.36) now measures the sharpness of a peak in frequency. One number, two apparently unrelated jobs. That coincidence is not a coincidence, and it is the point of this section.

6.4 · The Lorentzian

Near resonance the exact curve (0.8.42) takes a universal form. To see it, factor $\omega_0^{2}-\omega^{2}=(\omega_0-\omega)(\omega_0+\omega)$, then set $\omega_0+\omega\approx2\omega_0$ and $\omega\approx\omega_0$ everywhere else. Those replacements are valid when $\abs{\omega-\omega_0}\ll\omega_0$, which is to say for a sharp resonance. Then

That shape is the Lorentzian: a peak of half-width $\gamma$, falling off as the inverse square of the detuning. It is worth being clear about what has and has not been assumed. Nothing about springs entered. Only two ingredients did, and they were one mode with a natural frequency $\omega_0$, and a decay rate $\gamma$ for its amplitude. Any system with those two properties, driven near resonance, has this lineshape.

6.5 · The phase, and why its sign structure matters

Let's look again at $\tan\delta = 2\gamma\omega/(\omega_0^{2}-\omega^{2})$, and take it in three regimes.

• $\omega\ll\omega_0$: the denominator is positive and large, so $\delta\to0$. The mass follows the force. Push slowly and it goes where you push. The response is stiffness-controlled.
• $\omega=\omega_0$: the denominator vanishes and $\delta=\pi/2$ exactly, for any $\gamma$ whatsoever. The displacement lags the force by a quarter cycle, which puts the velocity exactly in phase with the force. That is why (0.8.41) is maximal here, and why the power resonance is exactly at $\omega_0$ while the amplitude resonance is not.
• $\omega\gg\omega_0$: the denominator is negative and large, so $\delta\to\pi$. The mass moves opposite to the applied force. Shake a heavy object fast and it stays put while your hand goes back and forth around it. The response is mass-controlled.

That a resonance flips the sign of the response as you pass through it is one of the most consequential facts in physics. It is why the refractive index of glass falls with wavelength in the visible and rises again beyond an absorption line, the effect called anomalous dispersion.

There is a mathematical statement underneath that. The real and imaginary parts of a response function, here $\mathrm{Re}\,Z$ and $\mathrm{Im}\,Z$, are not independent and in fact determine each other. That is the Kramers–Kronig relation, and it is causality written in the frequency domain. The Green's function of (0.8.8) vanishing for $t\lt s$ is the same statement in the time domain.

quality Q

Q = 6.00 γ = 0.08333

measured FWHM = 0.1666 (2γ = 0.1667)

energy lifetime τ = 1/2γ = 6.00 → FWHM × τ = 1.000

7 · Coupled oscillators, normal modes, and the bridge to fields

Everything so far has had one degree of freedom. Almost nothing in nature does. This section adds a second one, discovers that the added complexity is an illusion, and then adds infinitely many.

7.1 · Two masses, three springs

Two equal masses $m$ sit in a line. Each is tied to a wall by a spring of stiffness $k$, and they are tied to each other by a spring of stiffness $\kappa$. Let $x_1,x_2$ be their displacements from equilibrium. The coupling spring is stretched by $x_1-x_2$, so it pulls mass 1 back with force $-\kappa(x_1-x_2)$ and mass 2 with the opposite sign:

Neither equation can be solved alone: each contains the other's unknown. That is what "coupled" means, and it is why the problem looks harder than §4. Write it as one vector equation with $\vv x = (x_1,x_2)^{\mathsf T}$:

$M$ is symmetric, and that is not an accident of this example. The forces came from a potential energy

so $mM_{ij} = \partial_i\partial_j V$ is the Hessian of the potential, and Chapter 0.6 proved, via Clairaut, that a Hessian is always symmetric.

Any system of masses moving in a potential near a stable equilibrium has this same form, because Chapter 0.6's second-order Taylor expansion says $V\approx V_0+\tfrac12\vv x^{\mathsf T}H\vv x$ at a critical point. So the analysis below is not about this apparatus. It is about every stable multi-degree-of-freedom system there is.

7.2 · The spectral theorem does the whole job

By Chapter 0.5, a real symmetric matrix has real eigenvalues and an orthonormal basis of eigenvectors. Call them $\vv u_1,\vv u_2$ with $M\vv u_k = \omega_k^{2}\vv u_k$.

Writing the eigenvalue as $\omega_k^{2}$ presumes it is positive. It is positive exactly when the equilibrium is a minimum rather than a saddle, which is Chapter 0.6's classification now doing physical work. A negative eigenvalue would give $\omega$ imaginary and $\ee^{\abs\omega t}$ growth, which is an instability, and an instability is what a saddle is.

Define the normal coordinates $q_k = \vv u_k\cdot\vv x$, the components of the displacement along the eigenvectors. Differentiate twice and use (0.8.49), remembering that $\vv u_k$ is a constant vector and that $M$ is symmetric so it can be moved across the dot product:

Let's look at what that line is saying. Each $q_k$ obeys $\ddot q_k = -\omega_k^{2}q_k$, which is the equation of §4, standing on its own, with no reference to the other coordinate. The coupled system has been turned into two independent harmonic oscillators by a rotation of axes.

The step that made it work is the one in the middle of (0.8.51), and it needed $M^{\mathsf T}=M$. Without symmetry there is no orthonormal eigenbasis, and the whole manoeuvre fails.

7.3 · The two modes, explicitly

So let's diagonalise $M$ from (0.8.49). It has the form $aI + b\sigma$ with $\sigma=\left(\begin{smallmatrix}0&1\\1&0\end{smallmatrix}\right)$, whose eigenvectors are $(1,1)$ and $(1,-1)$ by inspection. Swap the entries of either one and it comes back unchanged up to a sign. So

Let's verify the eigenvalues directly. We have $M(1,1)^{\mathsf T} = \tfrac1m(k+\kappa-\kappa)(1,1)^{\mathsf T}$ and $M(1,-1)^{\mathsf T} = \tfrac1m(k+\kappa+\kappa)(1,-1)^{\mathsf T}$. Note also that $\vv u_1\cdot\vv u_2 = \tfrac12(1-1)=0$. The two are perpendicular, as Chapter 0.5 promised, without our having arranged it.

The physical reading is the pleasure of the calculation.

• In mode 1 the two masses move together, in phase, with equal amplitude. The coupling spring is then never stretched, so it exerts no force at all. The frequency is $\sqrt{k/m}$, which is independent of $\kappa$, exactly as if the coupling spring were not there.
• In mode 2 they move exactly out of phase. The midpoint stays fixed, the coupling spring is stretched by $2x_1$, and the restoring force on each mass is $-(k+2\kappa)x_1$. That is stiffer, and hence faster.

Those are the only two ways this system can move sinusoidally, and everything else is a superposition of them.

Putting the two modes together, the general motion of the system is

with the coefficients read off by orthonormality. That is Chapter 0.5's coordinates have become inner products, used here for the first time on a physical problem. There are four constants for a four-dimensional solution space, since we have two second-order equations, and that is consistent with (0.8.11).

7.4 · Beats

Now let's pick the initial condition that makes the structure visible. Pull mass 1 aside by $d$, hold mass 2 at rest at its equilibrium, and release both from rest, so that $\vv x(0)=(d,0)$ and $\dot{\vv x}(0)=\vv 0$. Then $a_1 = d/\sqrt2$, $a_2=d/\sqrt2$ and $b_k=0$, and (0.8.53) gives

We want the slow timescale and the fast one separated out, so apply the sum-to-product identities, which are Euler's formula rearranged, with $\bar\omega=\tfrac12(\omega_1+\omega_2)$ and $\Delta=\omega_2-\omega_1$:

Each mass oscillates fast at the mean frequency $\bar\omega$, inside a slowly changing envelope of frequency $\Delta/2$. The two envelopes are in quadrature, so when one is at its largest the other is zero. The energy sloshes completely from one mass to the other and back, with the envelopes exchanging roles after a time $\pi/\Delta$ and returning after

For weak coupling, expand $\omega_2 = \omega_1\sqrt{1+2\kappa/k}\approx\omega_1(1+\kappa/k)$ by Chapter 0.3. Then $\Delta\approx\omega_1\kappa/k$ and $T_{\text{beat}}\approx 2\pi k/(\omega_1\kappa)$, so the weaker the coupling, the longer the beat. It diverges as $\kappa\to0$, which is right, since uncoupled masses never exchange anything.

This is the mechanism behind two pendulums on a shared rail, and behind a photon oscillating between two cavities. With the same algebra and different words, it is also neutrino oscillation, where the two "masses" are two mass eigenstates with slightly different frequencies and the beat is measured in kilometres of flight path.

preset mode 1 (in phase) mode 2 (out of phase) one mass only

coupling κ

κ = 0.140

start x₁

start x₂

pause

ω₁ = 1.0000 ω₂ = 1.1314 Δω = 0.1314 T_beat = 47.8 s

a₁ = 0.500 a₂ = 0.500 mode energies 43.9% / 56.1%

7.5 · Normal modes are principal axes

Pause on what the eigenvectors are. We observed after (0.8.50) that $mM$ is the Hessian of $V$. So the normal-mode directions are the eigenvectors of the Hessian, which Chapter 0.6 identified as the principal axes of curvature of the potential surface, with eigenvalues equal to the curvatures along them. Hence

squared frequency $=$ curvature of the potential $/$ mass,

direction by direction. That gives three cases to keep straight.

• A stiff direction of the potential valley is a fast mode.
• A shallow direction is a slow mode.
• A flat direction is a zero-frequency mode, which does not oscillate at all but drifts. You will meet one in Problem 2.

So "normal modes" and "principal axes" are not two ideas that happen to use the same theorem. They are the same idea, and Chapter 6.6 will use the flat-direction version of it to make a massless particle out of a symmetry.

Grind box — $N$ masses in a chain, solved exactly

Nothing above used $n=2$. Take $N$ equal masses in a line, each connected to its neighbours by springs of stiffness $k$, with the chain pinned to walls at both ends. Writing $u_j$ for the displacement of mass $j$ and setting $u_0=u_{N+1}=0$:

$$m\ddot u_j = k\left(u_{j+1}-2u_j+u_{j-1}\right), \qquad j=1,\dots,N.$$

The matrix is symmetric, with $2$ on the diagonal, $-1$ next to it and zeros elsewhere. So the spectral theorem guarantees $N$ real frequencies and $N$ orthogonal modes before we compute anything. To find them, try $u_j = \sin(j\theta)$ for a constant $\theta$. Then

$$u_{j+1}+u_{j-1} = \sin(j\theta+\theta)+\sin(j\theta-\theta) = 2\cos\theta\,\sin j\theta,$$

by the sum-to-product identity, so the right-hand side is $k(2\cos\theta-2)u_j$. That is proportional to $u_j$ itself, which is exactly the eigenvector condition. Hence $-m\omega^{2} = 2k(\cos\theta-1)$, and using $1-\cos\theta = 2\sin^{2}(\theta/2)$,

$$\omega(\theta) = 2\sqrt{k/m}\ \left|\sin\frac{\theta}{2}\right|.$$

The boundary condition $u_0=0$ is automatic. The other one, $u_{N+1}=0$, requires $\sin\big((N+1)\theta\big)=0$, so $\theta_n = n\pi/(N+1)$ with $n=1,\dots,N$, since larger $n$ repeats the same modes. Therefore

$$\omega_n = 2\sqrt{\frac{k}{m}}\,\sin\!\left(\frac{n\pi}{2(N+1)}\right), \qquad u_j^{(n)} = \sin\!\left(\frac{n\pi j}{N+1}\right).$$

Check against §7.3. Put $N=2$: $\omega_1 = 2\sqrt{k/m}\sin(\pi/6)=\sqrt{k/m}$ and $\omega_2 = 2\sqrt{k/m}\sin(\pi/3)=\sqrt{3k/m}$, matching (0.8.52) at $\kappa=k$. ✓ And $u^{(1)}=(\sin\tfrac\pi3,\sin\tfrac{2\pi}3)\propto(1,1)$, $u^{(2)}=(\sin\tfrac{2\pi}3,\sin\tfrac{4\pi}3)\propto(1,-1)$. ✓

What to notice. $N$ masses give exactly $N$ modes, one per degree of freedom, which is (0.8.11) again. The mode shapes are sampled sine waves, with $n$ counting half-wavelengths across the chain. And the frequencies are not evenly spaced. They crowd together near the top of the band at $2\sqrt{k/m}$, because $\sin$ flattens. That crowding is the first hint that a discrete chain and a continuous string are not the same thing at short wavelengths, which is the subject of §7.6.

7.6 · The limit that makes a field

Now let the chain become continuous. Give the masses a spacing $a$ and let position along the chain be $x=ja$, so that $u_j(t)$ becomes a function $u(x,t)$ of two variables. Expand the neighbours by Taylor's theorem (Chapter 0.3):

Now add the two expansions together. Every odd term cancels, which is the payoff of using both neighbours at once, and what is left is

That second difference is exactly the combination appearing in the chain equation, so substitute it in and divide through by $m$:

Now take $a\to0$, holding the physical quantities fixed. Those are the mass per unit length $\rho = m/a$ and the tension $T = ka$, the tension being the force needed for unit fractional stretch of a spring of length $a$. Then $ka^{2}/m = (ka)(a/m) = T/\rho$, which is a finite limit, and the correction term dies. What survives is the wave equation:

Let's take stock of what that limit did. The discrete chain has become a continuous medium. The list of coordinates $u_1,\dots,u_N$ has become a field $u(x,t)$, one number for every point of space. And the matrix $M$ has become the operator $-v^{2}\partial_x^{2}$. Nothing else changed.

The modes survive the limit too. A normal mode is a solution in which every point oscillates at the same frequency, so let's look for one of those by putting $u(x,t)=X(x)\cos\omega t$ into (0.8.60). Substituting and cancelling the time factor gives $-\omega^{2}X = v^{2}X''$, which is $X'' = -(\omega/v)^{2}X$. That is the oscillator equation again, now living in space rather than in time. Its solutions are $\sin$ and $\cos$ of $\omega x/v$, and pinning the ends at $x=0$ and $x=L$ kills the cosine and quantises the rest:

These are the standing waves, and they are the $N\to\infty$ limit of the chain's modes. Setting $x=ja$ and $L=(N+1)a$ turns $u_j^{(n)}=\sin\big(n\pi j/(N+1)\big)$ into $\sin(n\pi x/L)$ exactly, while the grind box's frequency $2\sqrt{k/m}\sin\big(n\pi/2(N+1)\big)$ tends to $\sqrt{k/m}\,n\pi a/L = n\pi v/L$ for $n\ll N$.

The finite chain has $N$ modes. The field has one for every $n$, or equivalently one for every wavenumber $k_n=n\pi/L$. And the argument of §7.2 never used the dimension, so it carries over untouched: each mode is an independent harmonic oscillator, with its own frequency $\omega_n$, its own amplitude and its own conserved energy.

8 · Nonlinearity, briefly and honestly

Everything above rested on linearity. It is worth one section to see exactly what is lost when it goes, because most of the physics in the second half of this book is nonlinear.

8.1 · The pendulum, exactly

A bob of mass $m$ on a rigid massless rod of length $L$, at angle $\theta$ from vertical. The restoring torque is $-mgL\sin\theta$ and the moment of inertia is $mL^{2}$, so

This is not linear, since $\sin(\theta_1+\theta_2)\neq\sin\theta_1+\sin\theta_2$. Chapter 0.3's move is to expand, $\sin\theta = \theta - \theta^{3}/6+\cdots$, and keep the first term, giving $\ddot\theta=-\omega_0^{2}\theta$ and $T_0=2\pi\sqrt{L/g}$. That this is an approximation is not in question. The questions worth asking are how good an approximation it is, and where it stops being good at all.

The energy first integral of §4.3 still works here, because multiplying by $\dot\theta$ does not care about linearity. Multiply (0.8.62) by $\dot\theta$ and integrate:

Releasing from rest at $\theta_0$ fixes $\mathcal E = -\omega_0^{2}\cos\theta_0$, and the grind box turns this into an exact period. The answer is $T = T_0\cdot\tfrac2\pi K\!\left(\sin\tfrac{\theta_0}{2}\right)$ with $K$ the complete elliptic integral, and expanding it,

The period depends on amplitude. That single fact is the death of superposition. Two solutions of different sizes have different frequencies, so their sum is not a solution of anything.

Numerically, the correction is $+0.19\%$ at $\theta_0=10^{\circ}$, $+1.74\%$ at $30^{\circ}$, and $+18.0\%$ at $90^{\circ}$. That is why a pendulum clock is built to swing through a small angle, and why Huygens had to think hard about escapements.

8.2 · The phase portrait and the separatrix

Let's plot the state $(\theta,\dot\theta)$, as in §4.4. From (0.8.63) the trajectories are the level curves of $\mathcal E$, and there are two qualitatively different families.

• For $\mathcal E\lt\omega_0^{2}$ the curve closes. The bob swings back and forth, and near the origin the closed curve is the ellipse of §4.4, slightly distorted.
• For $\mathcal E\gt\omega_0^{2}$ there is enough energy to pass over the top. Now $\dot\theta$ never vanishes and the trajectory runs off to $\theta\to\infty$, which is to say the pendulum whirls.

Between them sits the single curve $\mathcal E = \omega_0^{2}$, on which the bob arrives at the inverted position with exactly zero speed. Substituting into (0.8.63) and using $1+\cos\theta = 2\cos^{2}(\theta/2)$:

This is the separatrix, and it does something the linear oscillator cannot. It takes infinite time to complete. Near $\theta=\pi$ the right-hand side vanishes linearly in $(\pi-\theta)$, so the approach is exponential and the destination is never reached.

The inverted position is an unstable equilibrium. It is a maximum of the potential, a negative Hessian eigenvalue in Chapter 0.6's classification, and hence an $\ee^{+\omega_0 t}$ in §3's language rather than an oscillation. Linearising about the top of the potential gives you exponential growth, and that is the correct local answer. There is nothing wrong with linearisation except that it is local.

8.3 · Chaos, in one paragraph

Add a drive and damping to (0.8.62) and it becomes possible for trajectories to be bounded, aperiodic, and exponentially sensitive to initial conditions. Two starting points a distance $\epsilon$ apart then separate like $\epsilon\,\ee^{\Lambda t}$ for some positive $\Lambda$.

Picard–Lindelöf is untouched by any of that. The solution still exists and is still unique. What changes is that predicting the state at time $t$ to fixed accuracy now requires knowing the initial state to accuracy $\ee^{-\Lambda t}$, so every factor of $\ee$ in your instruments buys one more unit of $1/\Lambda$ in forecast horizon.

That is the entire content of "chaos". Determinism and predictability are different properties, and only the first is guaranteed. It cannot happen in a linear system, where solutions are sums of $\ee^{\lambda_i t}$ and neighbouring trajectories either converge, or diverge in a fixed direction, or wind around each other forever, but never explore.

Grind box — the exact pendulum period

From (0.8.63) with $\mathcal E=-\omega_0^{2}\cos\theta_0$,

$$\dot\theta = \omega_0\sqrt{2\left(\cos\theta-\cos\theta_0\right)}.$$

A quarter period is the time from $\theta=0$ to $\theta=\theta_0$, so separating variables,

$$\frac{T}{4} = \int_0^{\theta_0}\frac{\dd\theta}{\omega_0\sqrt{2(\cos\theta-\cos\theta_0)}}.$$

Use $\cos\theta = 1-2\sin^{2}(\theta/2)$ to write $\cos\theta-\cos\theta_0 = 2\big(s^{2}-\sin^{2}(\theta/2)\big)$ where $s\equiv\sin(\theta_0/2)$. Then substitute $\sin(\theta/2) = s\sin\varphi$, which maps $\theta:0\to\theta_0$ onto $\varphi:0\to\pi/2$ and gives $\tfrac12\cos(\theta/2)\,\dd\theta = s\cos\varphi\,\dd\varphi$ with $\cos(\theta/2)=\sqrt{1-s^{2}\sin^{2}\varphi}$. The square root in the denominator becomes $2s\cos\varphi$, and it cancels against the numerator exactly:

$$\frac{T}{4} = \frac{1}{\omega_0}\int_0^{\pi/2}\frac{\dd\varphi}{\sqrt{1-s^{2}\sin^{2}\varphi}} \;\equiv\; \frac{K(s)}{\omega_0}.$$

$K$ is the complete elliptic integral of the first kind, a function defined by this integral, which cannot be reduced to elementary ones. Since $T_0 = 2\pi/\omega_0$, $T/T_0 = 2K(s)/\pi$, and $K(0)=\pi/2$ recovers $T=T_0$. ✓

The expansion. Expand the integrand by the binomial series (Chapter 0.3), $(1-u)^{-1/2} = 1+\tfrac12u+\tfrac38u^{2}+\cdots$ with $u=s^{2}\sin^{2}\varphi$, and integrate term by term using $\int_0^{\pi/2}\sin^{2}\varphi\,\dd\varphi = \pi/4$ and $\int_0^{\pi/2}\sin^{4}\varphi\,\dd\varphi = 3\pi/16$:

$$K(s) = \frac\pi2\left(1+\frac{s^{2}}{4}+\frac{9s^{4}}{64}+\cdots\right).$$

Finally put $s = \sin(\theta_0/2) = \tfrac{\theta_0}{2}-\tfrac{\theta_0^{3}}{48}+\cdots$, so $s^{2} = \tfrac{\theta_0^{2}}{4}-\tfrac{\theta_0^{4}}{48}+\cdots$ and $s^{4}=\tfrac{\theta_0^{4}}{16}+\cdots$, and collect:

$$\frac{T}{T_0} = 1 + \frac{\theta_0^{2}}{16} + \left(\frac{9}{64\cdot16}-\frac{1}{4\cdot48}\right)\theta_0^{4} + \cdots = 1+\frac{\theta_0^{2}}{16}+\frac{11\theta_0^{4}}{3072}+\cdots$$

since $\tfrac{9}{1024}-\tfrac{1}{192} = \tfrac{27-16}{3072} = \tfrac{11}{3072}$, confirming (0.8.64). At $\theta_0=1$ rad the series gives $1.06608$ against the exact $1.06633$. The next term is therefore doing $2\times10^{-4}$ of work, which is what an $O(\theta_0^{6})$ term should do.

9 · Worked examples

10 · Your turn