Formula & theorem reference

A linear differential operator L is built entirely from partial derivatives combined by addition and multiplication by coefficients that depend only on the independent variables — never on u itself — so L inherits linearity directly from the linearity of differentiation: L[au+bv] = aL[u] + bL[v] for any functions u,v and constants a,b. Apply this repeatedly to u = c_1u_1+\cdots+c_nu_n: L[u] = L[c_1u_1+\cdots+c_nu_n] = c_1L[u_1]+\cdots+c_nL[u_n] = c_1(0)+\cdots+c_n(0) = 0, using that each u_i already solves L[u_i]=0. The sum collapses to 0 only because every individual right side is already 0 — exactly why superposition is stated for a homogeneous equation: for L[u]=g with g\not\equiv0, the same computation gives L[c_1u_1+c_2u_2] = (c_1+c_2)g, which equals g only in the special case c_1+c_2=1, not for arbitrary constants.

The argument extends to an infinite sum \sum_{n=1}^\infty c_nu_n wherever the series converges and can be differentiated term by term — the tool later chapters build an entire solution method on.

Work in one spatial dimension first: a thin rod with constant density \rho, specific heat c, and thermal conductivity K, all per unit cross-sectional area. Isolate the segment between x and x+\Delta x. Its heat energy is \int_x^{x+\Delta x} c\rho\,u(s,t)\,ds, and with no source, that energy changes only because heat flux \phi(x,t) carries energy across the segment's two ends: \frac{d}{dt}\int_x^{x+\Delta x} c\rho\,u(s,t)\,ds = \phi(x,t) - \phi(x+\Delta x,t). The left side is c\rho\,u_t(x,t)\,\Delta x to leading order as \Delta x\to0; dividing both sides by \Delta x and taking that limit turns the right side into -\phi_x(x,t), giving the conservation law c\rho\,u_t = -\phi_x. Fourier's law supplies \phi itself: \phi = -Ku_x (heat flows from hot toward cold, at a rate proportional to how steep the temperature gradient is). Substituting, c\rho\,u_t = -\big(-Ku_x\big)_x = Ku_{xx}, and setting k = K/(c\rho) gives u_t = ku_{xx}.

The same argument runs over an arbitrary region R in two or three dimensions, with the divergence theorem standing in for "flux across the two endpoints": the total heat energy in R changes at rate -\oint_{\partial R}\boldsymbol\phi\cdot\mathbf n\,dA = -\int_R \nabla\cdot\boldsymbol\phi\,dV, and Fourier's law generalizes to \boldsymbol\phi = -K\nabla u. The pointwise version of the resulting balance, c\rho\,u_t = -\nabla\cdot\boldsymbol\phi, becomes c\rho\,u_t = K\nabla^2u once K is constant, and dividing by c\rho gives u_t = k\nabla^2u exactly as in one dimension.

Model a stretched string of constant linear density \rho (mass per unit length) under constant tension T, vibrating transversely with small displacement u(x,t). "Small" means the string's slope u_x stays small everywhere, so the angle \theta(x,t) the string makes with the horizontal satisfies \sin\theta\approx\tan\theta = u_x — the approximation the whole derivation leans on.

Isolate the segment between x and x+\Delta x. Because the displacement is purely transverse, the horizontal component of tension is the same constant T at both ends (no horizontal motion to explain otherwise); the net vertical force on the segment is the difference between the vertical components of tension pulling at each end, T\sin\theta(x+\Delta x,t) - T\sin\theta(x,t) \approx T\big(u_x(x+\Delta x,t) - u_x(x,t)\big). Newton's second law equates this to mass times vertical acceleration, \rho\,\Delta x\,u_{tt}(x,t), to leading order: \rho\,\Delta x\,u_{tt} \approx T\big(u_x(x+\Delta x,t)-u_x(x,t)\big). Dividing by \Delta x and letting \Delta x\to0 turns the right side into Tu_{xx}, giving \rho\,u_{tt} = Tu_{xx}. Setting c^2 = T/\rho — tension makes the string snap back harder, density makes it resist moving — gives u_{tt} = c^2u_{xx}.

The same balance of forces, carried out on a small patch of a uniformly tensioned membrane instead of a segment of string, produces the two-dimensional case u_{tt} = c^2\nabla^2u by the identical small-displacement argument, with the divergence theorem replacing the two-endpoint difference the way it does in the heat equation's own higher-dimensional derivation.

Start from the conservation law behind any diffusing quantity with flux \boldsymbol\phi and source f: u_t + \nabla\cdot\boldsymbol\phi = f. A steady state is one that has stopped changing in time, so u_t\equiv0; with no internal source (f\equiv0 — no heat being added or removed, no charge being created), the conservation law itself collapses to \nabla\cdot\boldsymbol\phi = 0.

For a diffusive flux, the constitutive law is Fourier's law, \boldsymbol\phi = -K\nabla u with K the (constant) conductivity. Substituting, \nabla\cdot(-K\nabla u) = -K\nabla^2u = 0, and since K\neq0, this is exactly \nabla^2u=0.

Nothing in this argument is specific to heat — the same three lines apply to any quantity whose flux is proportional to the negative gradient of some scalar, unchanging in time, with nothing being created or destroyed: steady-state diffusion of a dissolved concentration, or an electrostatic potential in a charge-free region, satisfy the identical equation for the identical reason.

Let \xi=\xi(x,y), \eta=\eta(x,y) be a smooth change of variables with nonzero Jacobian J = \xi_x\eta_y - \xi_y\eta_x. Expanding u_{xx}, u_{xy}, u_{yy} by the chain rule in terms of \xi,\eta-derivatives and collecting the new second-order coefficients gives \tilde A = A\xi_x^2 + B\xi_x\xi_y + C\xi_y^2, \quad \tilde C = A\eta_x^2 + B\eta_x\eta_y + C\eta_y^2, \quad \tilde B = 2A\xi_x\eta_x + B(\xi_x\eta_y+\xi_y\eta_x) + 2C\xi_y\eta_y. Squaring \tilde B and subtracting 4\tilde A\tilde C, every cross term that mixes \xi- and \eta-derivatives cancels, leaving \tilde B^2 - 4\tilde A\tilde C = (B^2-4AC)(\xi_x\eta_y-\xi_y\eta_x)^2 = D\cdot J^2. Since J\neq0, J^2>0, so \tilde B^2-4\tilde A\tilde C and D always share the same sign — the classification survives the change of variables even though A, B, C themselves do not.

The differential form u_t+f(u)_x=0 already assumes u is smooth enough to differentiate, so it says nothing by itself about a jump — the jump condition comes from going back to the integral form the differential equation was derived from. For any fixed interval [x_1,x_2] (with x_1<s(t)<x_2), conservation of the total amount of u reads \frac{d}{dt}\int_{x_1}^{x_2}u\,dx = f(u(x_1,t)) - f(u(x_2,t)), the flux entering at x_1 minus the flux leaving at x_2. Split the integral at the shock and differentiate each piece with the Leibniz rule, which picks up a boundary term wherever the limit of integration itself moves: \frac{d}{dt}\int_{x_1}^{x_2}u\,dx = \frac{d}{dt}\int_{x_1}^{s(t)}u\,dx + \frac{d}{dt}\int_{s(t)}^{x_2}u\,dx = \int_{x_1}^{s(t)}u_t\,dx + u_L\,s'(t) + \int_{s(t)}^{x_2}u_t\,dx - u_R\,s'(t). Now shrink the interval onto the shock, x_1\to s(t)^- and x_2\to s(t)^+. Both remaining integrals vanish in this limit (u_t stays bounded on either side, away from the shock itself, so each integral is over a shrinking interval of a bounded integrand), leaving \frac{d}{dt}\int_{x_1}^{x_2}u\,dx \;\longrightarrow\; \sigma\,u_L - \sigma\,u_R = -\sigma\,\llbracket u\rrbracket. The right side of the conservation law shrinks the same way: f(u(x_1,t))-f(u(x_2,t)) \to f(u_L)-f(u_R) = -\llbracket f(u)\rrbracket. Equating the two limits and canceling the shared minus sign gives \sigma\,\llbracket u\rrbracket=\llbracket f(u)\rrbracket — the shock speed is exactly the ratio of the flux jump to the density jump, because that is the only speed at which the total amount of u on each side of the shock is still conserved through the discontinuity.

Let S_N(x)=\dfrac{a_0}2+\displaystyle\sum_{n=1}^{N}\left(a_n\cos\dfrac{n\pi x}{L}+b_n\sin\dfrac{n\pi x}{L}\right) be the Nth partial sum, and expand the mean-square error directly: \frac1L\int_{-L}^{L}(f-S_N)^2\,dx = \frac1L\int_{-L}^{L}f^2\,dx - \frac2L\int_{-L}^{L}fS_N\,dx + \frac1L\int_{-L}^{L}S_N^2\,dx. \quad(\ast)

The cross term. Multiplying S_N's own definition by f and integrating term by term, every term collapses to the coefficient that defined it: \frac1L\int_{-L}^{L}fS_N\,dx = \frac{a_0}2\cdot\underbrace{\frac1L\int_{-L}^Lf\,dx}_{=\,a_0} + \sum_{n=1}^{N}\left(a_n\cdot\underbrace{\frac1L\int_{-L}^Lf\cos\tfrac{n\pi x}L\,dx}_{=\,a_n} + b_n\cdot\underbrace{\frac1L\int_{-L}^Lf\sin\tfrac{n\pi x}L\,dx}_{=\,b_n}\right) = \frac{a_0^2}2+\sum_{n=1}^{N}(a_n^2+b_n^2).

The S_N^2 term. Expanding S_N^2 produces a constant-squared term, N squared-cosine and squared-sine terms, and every cross product of two distinct terms in the sum defining S_N. Every cross product integrates to 0 over [-L,L] by the orthogonality of the trigonometric functions (a cosine against a different-frequency cosine, a sine against a different-frequency sine, a cosine against any sine, and the constant against any cosine or sine, are each orthogonal on [-L,L]), leaving only the squared terms: \frac1L\int_{-L}^{L}S_N^2\,dx = \left(\frac{a_0}2\right)^2\cdot2 + \sum_{n=1}^{N}\left(a_n^2\cdot\frac1L\!\int_{-L}^L\!\cos^2\tfrac{n\pi x}L\,dx + b_n^2\cdot\frac1L\!\int_{-L}^L\!\sin^2\tfrac{n\pi x}L\,dx\right) = \frac{a_0^2}2+\sum_{n=1}^{N}(a_n^2+b_n^2), using \frac1L\int_{-L}^L\cos^2(n\pi x/L)\,dx=\frac1L\int_{-L}^L\sin^2(n\pi x/L)\,dx=1 for n\geq1.

Substituting back into (\ast). Both bracketed sums above are identical, so \frac1L\int_{-L}^{L}(f-S_N)^2\,dx = \frac1L\int_{-L}^{L}f^2\,dx - \frac{a_0^2}2-\sum_{n=1}^{N}(a_n^2+b_n^2). \quad(\ast\ast) The left side is a mean-square average of a real-valued function, so it is never negative, for every N. Rearranging (\ast\ast) gives Bessel's inequality, \frac{a_0^2}2+\sum_{n=1}^{N}(a_n^2+b_n^2)\;\leq\;\frac1L\int_{-L}^{L}f^2\,dx, for every partial-sum length N — the partial sums of the squared coefficients never overshoot the function's own mean-square size.

From inequality to equality. The trigonometric functions form a complete system on [-L,L]: for a piecewise continuous f, the partial sums S_N converge to f in the mean-square sense, \frac1L\int_{-L}^L(f-S_N)^2\,dx\to0 as N\to\infty (this completeness fact is taken as given at this course's level, the same way it's assumed rather than re-derived for the abstract inner-product-space version of this same statement). Letting N\to\infty in (\ast\ast) then turns the left side to 0, and Bessel's inequality becomes the equality \frac1L\int_{-L}^{L}f(x)^2\,dx = \frac{a_0^2}2+\sum_{n=1}^{\infty}(a_n^2+b_n^2).

Sanity check. For f(x)=x on [-\pi,\pi] (L=\pi), a_n=0 for all n (odd function) and b_n=\dfrac{2(-1)^{n+1}}n, so Parseval's equality reads \dfrac1\pi\displaystyle\int_{-\pi}^{\pi}x^2\,dx=\sum_{n=1}^\infty\dfrac4{n^2}, i.e. \dfrac{2\pi^2}3=4\displaystyle\sum_{n=1}^\infty\dfrac1{n^2} — which rearranges to the classical \sum_{n=1}^\infty 1/n^2=\pi^2/6, recovering a famous numerical series sum from nothing but a single Fourier coefficient computation.

Write L[y]=(py')'+qy for the Sturm–Liouville operator, so the equation reads L[y]=-\lambda\sigma y. The whole proof rests on one algebraic identity, true for any two twice-differentiable u,v on [a,b]: vL[u]-uL[v] = v(pu')'-u(pv')' = \frac{d}{dx}\Big[p\,(vu'-uv')\Big], which expands by the product rule and the two terms pv'u'-pu'v' cancelling. Integrating both sides from a to b gives Lagrange's identity: \int_a^b\big(vL[u]-uL[v]\big)\,dx = \Big[p\,(vu'-uv')\Big]_a^b. \quad(\ast)

The boundary term vanishes for any u,v satisfying the same separated conditions. At x=a, both (u(a),u'(a)) and (v(a),v'(a)) solve the same linear equation \alpha_1w+\alpha_2w'=0 in the plane — a single line through the origin, since (\alpha_1,\alpha_2)\neq(0,0). Two points on the same line through the origin are proportional (or one is the origin), so the 2\times2 determinant built from them vanishes: v(a)u'(a)-u(a)v'(a)=0. The identical argument at x=b (using \beta_1,\beta_2) kills the term there too, so \big[p(vu'-uv')\big]_a^b=0 whenever u and v both satisfy the problem's boundary conditions.

Reality. Suppose \lambda is an eigenvalue, possibly complex, with eigenfunction y (possibly complex-valued): L[y]=-\lambda\sigma y. Since p,q,\sigma are real-valued, conjugating the whole equation gives L[\bar y]=-\bar\lambda\sigma\bar y. The boundary conditions have real coefficients \alpha_1,\alpha_2,\beta_1,\beta_2, so \bar y satisfies exactly the same conditions y does — the vanishing argument above applies to the pair u=y, v=\bar y. Apply (\ast): \int_a^b\big(\bar yL[y]-yL[\bar y]\big)\,dx = 0 \quad\Longrightarrow\quad \int_a^b\big(-\lambda\,\sigma|y|^2+\bar\lambda\,\sigma|y|^2\big)\,dx = 0 \quad\Longrightarrow\quad (\bar\lambda-\lambda)\int_a^b\sigma|y|^2\,dx = 0. An eigenfunction is not identically zero and \sigma(x)>0 on [a,b], so \int_a^b\sigma|y|^2\,dx>0. Dividing forces \bar\lambda=\lambda: every eigenvalue is real. (Once \lambda is known real, the real and imaginary parts of y each solve the same real equation with the same real \lambda, so a real eigenfunction can always be chosen — used freely from here on.)

Orthogonality. Now take y_m,y_n to be real eigenfunctions for distinct real eigenvalues \lambda_m\neq\lambda_n: L[y_m]=-\lambda_m\sigma y_m and L[y_n]=-\lambda_n\sigma y_n. Apply (\ast) to u=y_m, v=y_n — both satisfy the problem's boundary conditions, so the right-hand side is again 0: \int_a^b\big(y_nL[y_m]-y_mL[y_n]\big)\,dx = 0 \quad\Longrightarrow\quad (\lambda_n-\lambda_m)\int_a^b\sigma\,y_my_n\,dx = 0. Since \lambda_m\neq\lambda_n, dividing gives \int_a^b\sigma(x)y_m(x)y_n(x)\,dx=0 — orthogonality, weighted by \sigma.

Sanity check. Take p\equiv1, q\equiv0, \sigma\equiv1, [a,b]=[0,L], with y(0)=y(L)=0 — this is exactly the eigenvalue problem y''+\lambda y=0 behind the Fourier sine series, whose eigenfunctions y_n=\sin(n\pi x/L) are already known to be real and mutually orthogonal on [0,L] with weight 1: this theorem is what makes that a special case of one general fact, not a coincidence specific to sines.

Let w=u_1-u_2, the difference of two solutions sharing the same initial and boundary data — so w itself solves w_t=kw_{xx} with w(x,0)=0 and, at each end, either w=0 (Dirichlet) or w_x=0 (Neumann). Define the energy E(t) = \int_0^L w(x,t)^2\,dx. Differentiating under the integral sign and using the PDE, E'(t) = 2\int_0^Lww_t\,dx = 2k\int_0^Lww_{xx}\,dx. Integrate by parts: \int_0^Lww_{xx}\,dx = \Big[ww_x\Big]_0^L - \int_0^Lw_x^2\,dx. At each endpoint, the boundary term vanishes: if the condition there is Dirichlet, w=0 at that end kills the product outright; if Neumann, w_x=0 at that end does the same — any combination of the two types at the two ends works identically. So E'(t) = -2k\int_0^Lw_x^2\,dx \le 0, and since E(0)=0 (from w(x,0)=0) and E(t)\ge0 always (it's an integral of a square), E(t)\equiv0 for every t\in[0,T] — forcing w\equiv0, i.e. u_1=u_2. This is uniqueness.

Continuous dependence follows the same computation with w=u_1-u_2 now the difference of two solutions sharing the same boundary data but different initial data f_1,f_2: the identical argument shows E(t)\le E(0)=\int_0^L(f_1-f_2)^2\,dx for every t (since E is non-increasing, not merely zero), so a small initial difference in the mean-square sense stays small for all later time.

Differentiate E under the integral sign (justified since u is smooth): E'(t) = \int_0^L\Big(u_tu_{tt} + c^2u_xu_{xt}\Big)\,dx. Integrate the second term by parts in x, moving one derivative off u_{xt} and onto u_x: \int_0^L c^2u_xu_{xt}\,dx = \Big[c^2u_xu_t\Big]_0^L - \int_0^L c^2u_{xx}u_t\,dx.

The boundary term vanishes. At each endpoint, the stated condition kills one of the two factors in c^2u_xu_t: at a Dirichlet end, u(0,t)=0 (or u(L,t)=0) holds for every t, so differentiating in t gives u_t=0 there too; at a Neumann end, u_x=0 there directly by hypothesis. Either way the product u_xu_t is zero at that endpoint, for every t -- and this holds independently at x=0 and at x=L, so any combination of the two condition types across the two ends still kills \big[c^2u_xu_t\big]_0^L.

Substituting back: E'(t) = \int_0^L u_tu_{tt}\,dx - \int_0^L c^2u_{xx}u_t\,dx = \int_0^L u_t\big(u_{tt}-c^2u_{xx}\big)\,dx. \quad(\ast) Since u solves the wave equation, u_{tt}-c^2u_{xx}\equiv0 everywhere on [0,L], so the integrand in (\ast) is identically zero and E'(t)=0 for every t>0. A function with zero derivative on an interval is constant there, so E(t)=E(0).

Corollary (uniqueness). Suppose u_1,u_2 both solve the same wave equation on [0,L] with the same boundary conditions (each end Dirichlet or Neumann) and the same initial data u(x,0)=f(x), u_t(x,0)=g(x). The difference w=u_1-u_2 solves the same wave equation with the same type of homogeneous boundary conditions at each end (linearity: a Dirichlet end has w=u_1-u_2=0-0=0; a Neumann end has w_x=0-0=0) and zero initial data, w(x,0)=0, w_t(x,0)=0. The energy of w is therefore constant by the theorem just proved, and its value at t=0 is E_w(0) = \frac12\int_0^L\big(w_t(x,0)^2+c^2w_x(x,0)^2\big)\,dx = \frac12\int_0^L\big(0^2+c^2\cdot0^2\big)\,dx = 0 (using w_t(x,0)=0 from the zero initial velocity, and w_x(x,0)=0 since w(x,0)\equiv0 on the whole interval makes its x-derivative zero too). So E_w(t)\equiv0 for all t. Since E_w(t) is a sum of two squared terms integrated over [0,L], E_w(t)=0 forces w_t(x,t)\equiv0 and w_x(x,t)\equiv0 on [0,L] for every t, so w is constant in both x and t -- and w(x,0)=0 pins that constant at 0. Hence w\equiv0 and u_1\equiv u_2: the wave equation's initial-boundary value problem has at most one solution.

Sanity check. For a string fixed at both ends (u(0,t)=u(L,t)=0) vibrating in its lowest mode, u(x,t)=\sin(\pi x/L)\cos(\pi ct/L): at t=0 the string is at its extreme displacement, momentarily at rest (u_t\equiv0), so all of E(0) is potential energy; a quarter-period later, at t=L/(2c), the string passes through u\equiv0 moving at its fastest, so all of E is kinetic. The two values of E computed directly from the formula agree -- energy has moved between the two forms without changing its total, exactly what the theorem guarantees.

Start from the disk's own series solution (separation of variables in polar coordinates, matched to f's full Fourier series): u(r,\theta) = \frac{a_0}2+\sum_{n=1}^{\infty}r^n\big(a_n\cos n\theta+b_n\sin n\theta\big), \qquad a_n=\frac1\pi\int_{-\pi}^{\pi}f(\varphi)\cos n\varphi\,d\varphi,\ \ b_n=\frac1\pi\int_{-\pi}^{\pi}f(\varphi)\sin n\varphi\,d\varphi (and a_0=\frac1\pi\int_{-\pi}^\pi f\,d\varphi). Substitute the coefficient integrals directly into the series and pull the sum inside the integral -- justified for r<1, since \sum r^n converges absolutely there: u(r,\theta) = \frac1{2\pi}\int_{-\pi}^{\pi}f(\varphi)\,d\varphi + \frac1\pi\sum_{n=1}^{\infty}r^n\int_{-\pi}^{\pi}f(\varphi)\big[\cos n\varphi\cos n\theta+\sin n\varphi\sin n\theta\big]\,d\varphi. The bracket is \cos n(\theta-\varphi) (the cosine angle-subtraction identity), so swapping the (now finite, for each fixed r<1) sum back inside the integral gives u(r,\theta) = \frac1\pi\int_{-\pi}^{\pi}f(\varphi)\left[\frac12+\sum_{n=1}^{\infty}r^n\cos n(\theta-\varphi)\right]d\varphi.

Summing the bracket. Write \psi=\theta-\varphi and sum the geometric series in complex form: with z=re^{i\psi} (|z|=r<1), \sum_{n=1}^{\infty}r^n\cos n\psi = \operatorname{Re}\sum_{n=1}^\infty z^n = \operatorname{Re}\frac{z}{1-z} = \operatorname{Re}\frac{re^{i\psi}(1-re^{-i\psi})}{|1-re^{i\psi}|^2} = \frac{r\cos\psi-r^2}{1-2r\cos\psi+r^2} (using |1-re^{i\psi}|^2=(1-r\cos\psi)^2+(r\sin\psi)^2=1-2r\cos\psi+r^2). Adding the \frac12: \frac12+\frac{r\cos\psi-r^2}{1-2r\cos\psi+r^2} = \frac{\tfrac12(1-2r\cos\psi+r^2)+r\cos\psi-r^2}{1-2r\cos\psi+r^2} = \frac{\tfrac12-\tfrac12r^2}{1-2r\cos\psi+r^2} = \frac12\,P_r(\psi). Substituting back gives u(r,\theta)=\frac1{2\pi}\int_{-\pi}^\pi f(\varphi)P_r(\theta-\varphi)\,d\varphi, the claimed formula.

P_r\geq0 and its average is 1. The denominator 1-2r\cos\psi+r^2=(1-r\cos\psi)^2+(r\sin\psi)^2\geq(1-r)^2>0 for r<1, and the numerator 1-r^2>0, so P_r(\psi)>0 everywhere. Setting f\equiv1 (so every a_n,b_n=0 except a_0=2) makes the series solution the constant function u\equiv1 for every r, and the formula above then forces \frac1{2\pi}\int_{-\pi}^\pi P_r(\psi)\,d\psi=1.

Sanity check. At r=0, P_0(\psi)=\frac{1-0}{1-0+0}=1 for every \psi, so u(0,\theta)=\frac1{2\pi}\int_{-\pi}^\pi f(\varphi)\,d\varphi -- the disk's center value is the plain average of the boundary data, exactly the mean value property AMVC's own harmonic-function theory already states for a general domain, recovered here as the r=0 special case of the explicit disk formula.

The wave operator \partial_t^2-c^2\partial_x^2 factors as (\partial_t-c\partial_x)(\partial_t+c\partial_x), the product of two first-order transport operators (ch.02's own method of characteristics, run on each factor in turn). Equivalently, and more directly, change to characteristic coordinates \xi=x+ct, \eta=x-ct. The chain rule gives \partial_x=\partial_\xi+\partial_\eta and \partial_t=c(\partial_\xi-\partial_\eta), so u_{tt}-c^2u_{xx} = c^2\big[(\partial_\xi-\partial_\eta)^2-(\partial_\xi+\partial_\eta)^2\big]u = -4c^2u_{\xi\eta}. The equation becomes u_{\xi\eta}=0: integrating in \eta first shows u_\xi doesn't depend on \eta, then integrating in \xi gives u(\xi,\eta) = F(\xi)+G(\eta) = F(x+ct)+G(x-ct) for some single-variable functions F,G -- every solution of the wave equation on the whole line is a right-traveling profile plus a left-traveling profile, with no other freedom.

Matching the initial data. u(x,0)=F(x)+G(x)=f(x). Differentiating u in t and setting t=0: u_t(x,0)=cF'(x)-cG'(x)=g(x), so F'(x)-G'(x)=g(x)/c; integrating from 0 to x, F(x)-G(x)=\frac1c\int_0^xg(s)\,ds+C for a constant C. Solving the two equations for F and G: F(x)=\frac12\Big[f(x)+\frac1c\int_0^xg(s)\,ds+C\Big],\qquad G(x)=\frac12\Big[f(x)-\frac1c\int_0^xg(s)\,ds-C\Big]. The constant C cancels in the final substitution: with u(x,t)=F(x+ct)+G(x-ct), u(x,t) = \frac{f(x+ct)+f(x-ct)}2 + \frac1{2c}\left[\int_0^{x+ct}g(s)\,ds-\int_0^{x-ct}g(s)\,ds\right] = \frac{f(x-ct)+f(x+ct)}2 + \frac1{2c}\int_{x-ct}^{x+ct}g(s)\,ds, the claimed formula.

The heat equation has a scaling symmetry: if u(x,t) solves u_t=ku_{xx}, so does u(\lambda x,\lambda^2t) for any \lambda>0 (the \lambda^2 on t exactly compensates the two x-derivatives). A solution built from data concentrated at a single point has no length scale of its own to break this symmetry, so it should be expressible as a function of the single combined variable \eta=x/\sqrt t alone, up to an overall power of t fixing units. Try u(x,t) = \frac1{\sqrt t}\,\varphi(\eta),\qquad \eta=\frac x{\sqrt t}.

Compute the two derivatives the equation needs, using \partial\eta/\partial t=-\tfrac12\eta t^{-1} and \partial\eta/\partial x=t^{-1/2}: u_t = -\tfrac12t^{-3/2}\big[\varphi(\eta)+\eta\varphi'(\eta)\big],\qquad u_{xx} = t^{-3/2}\varphi''(\eta). Substituting into u_t=ku_{xx} and cancelling the common t^{-3/2}: k\varphi''+\tfrac12\eta\varphi'+\tfrac12\varphi=0,\qquad\text{i.e.}\qquad k\varphi'' + \tfrac12(\eta\varphi)' = 0 (the last step using (\eta\varphi)'=\varphi+\eta\varphi') -- the self-similarity ansatz has turned the PDE into an ODE in \eta alone.

Solving the ODE. Integrating once, k\varphi'+\tfrac12\eta\varphi=A for a constant A; a fundamental solution and all its derivatives decay to 0 as \eta\to\pm\infty, which forces A=0. What remains is separable: \varphi'/\varphi=-\eta/(2k), so \ln\varphi=-\eta^2/(4k)+\text{const}, i.e. \varphi(\eta)=Be^{-\eta^2/(4k)} for a constant B. Substituting back, u(x,t) = \frac Bt{\vphantom{1}}^{1/2}e^{-x^2/(4kt)}.

Fixing B. A fundamental solution's defining property is that it represents unit total heat for every t: \int_{-\infty}^\infty u(x,t)\,dx=1. Substituting s=x/\sqrt{4kt}, \int_{-\infty}^\infty Bt^{-1/2}e^{-x^2/(4kt)}\,dx = B\sqrt{4\pi k} (a standard Gaussian integral), so B=1/\sqrt{4\pi k}, giving S(x,t)=\frac1{\sqrt{4\pi kt}}e^{-x^2/(4kt)}, the claimed formula.

Heat equation. For u_t-ku_{xx}=F(x,t), u(x,0)=0, let v(x,t;s) solve the HOMOGENEOUS heat equation v_t=kv_{xx} for t>s with v(x,s;s)=F(x,s) (the source at time s, injected as an initial condition at that instant). Define u(x,t) = \int_0^t v(x,t;s)\,ds. Differentiating under the integral sign (Leibniz's rule, the boundary term at t=s contributing exactly v(x,t;t)=F(x,t)): u_t = v(x,t;t) + \int_0^tv_t(x,t;s)\,ds = F(x,t) + \int_0^tkv_{xx}(x,t;s)\,ds = F(x,t)+ku_{xx}, using v_t=kv_{xx} inside the integral -- exactly u_t-ku_{xx}=F. And u(x,0)=\int_0^0(\cdots)\,ds=0 automatically. By linearity, adding this to the solution of the homogeneous problem with the ACTUAL initial data (§4's whole-line heat solution) handles a heat equation with both nonzero initial data and a source term at once.

Wave equation. For u_{tt}-c^2u_{xx}=F(x,t), u(x,0)=u_t(x,0)=0, the same idea applies ONE DERIVATIVE HIGHER, since the wave equation is second order in t: let w(x,t;s) solve the homogeneous wave equation for t>s with w(x,s;s)=0 and w_t(x,s;s)=F(x,s) -- the source injected as an initial VELOCITY, not an initial position, at time s. Define u(x,t)=\int_0^tw(x,t;s)\,ds. Differentiating once, u_t = w(x,t;t) + \int_0^tw_t(x,t;s)\,ds = 0 + \int_0^tw_t(x,t;s)\,ds (the boundary term vanishes because w(x,s;s)=0 is exactly how each slice was built). Differentiating again, u_{tt} = w_t(x,t;t) + \int_0^tw_{tt}(x,t;s)\,ds = F(x,t) + \int_0^tc^2w_{xx}(x,t;s)\,ds = F(x,t)+c^2u_{xx} (this time the boundary term is w_t(x,t;t)=F(x,t), the injected velocity, and w_{tt}=c^2w_{xx} is used inside the integral) -- exactly u_{tt}-c^2u_{xx}=F, with u(x,0)=0 and u_t(x,0)=0 both automatic from the empty-interval integral.

Both cases are the SAME construction: differentiate the accumulated integral, let the boundary term at t=s (which derivative order it lands on depends on the equation's own order in t) absorb the source, and let the homogeneous equation inside the integral absorb everything else.

Heat equation. For u_t-ku_{xx}=F(x,t), u(x,0)=0, let v(x,t;s) solve the HOMOGENEOUS heat equation v_t=kv_{xx} for t>s with v(x,s;s)=F(x,s) (the source at time s, injected as an initial condition at that instant). Define u(x,t) = \int_0^t v(x,t;s)\,ds. Differentiating under the integral sign (Leibniz's rule, the boundary term at t=s contributing exactly v(x,t;t)=F(x,t)): u_t = v(x,t;t) + \int_0^tv_t(x,t;s)\,ds = F(x,t) + \int_0^tkv_{xx}(x,t;s)\,ds = F(x,t)+ku_{xx}, using v_t=kv_{xx} inside the integral -- exactly u_t-ku_{xx}=F. And u(x,0)=\int_0^0(\cdots)\,ds=0 automatically. By linearity, adding this to the solution of the homogeneous problem with the ACTUAL initial data (§4's whole-line heat solution) handles a heat equation with both nonzero initial data and a source term at once.

Wave equation. For u_{tt}-c^2u_{xx}=F(x,t), u(x,0)=u_t(x,0)=0, the same idea applies ONE DERIVATIVE HIGHER, since the wave equation is second order in t: let w(x,t;s) solve the homogeneous wave equation for t>s with w(x,s;s)=0 and w_t(x,s;s)=F(x,s) -- the source injected as an initial VELOCITY, not an initial position, at time s. Define u(x,t)=\int_0^tw(x,t;s)\,ds. Differentiating once, u_t = w(x,t;t) + \int_0^tw_t(x,t;s)\,ds = 0 + \int_0^tw_t(x,t;s)\,ds (the boundary term vanishes because w(x,s;s)=0 is exactly how each slice was built). Differentiating again, u_{tt} = w_t(x,t;t) + \int_0^tw_{tt}(x,t;s)\,ds = F(x,t) + \int_0^tc^2w_{xx}(x,t;s)\,ds = F(x,t)+c^2u_{xx} (this time the boundary term is w_t(x,t;t)=F(x,t), the injected velocity, and w_{tt}=c^2w_{xx} is used inside the integral) -- exactly u_{tt}-c^2u_{xx}=F, with u(x,0)=0 and u_t(x,0)=0 both automatic from the empty-interval integral.

Both cases are the SAME construction: differentiate the accumulated integral, let the boundary term at t=s (which derivative order it lands on depends on the equation's own order in t) absorb the source, and let the homogeneous equation inside the integral absorb everything else.

Integrate the defining transform of f' by parts, with u = e^{-st}, dv = f'(t)\,dt: \mathcal{L}\{f'(t)\} = \int_0^\infty f'(t)e^{-st}\,dt = \left[f(t)e^{-st}\right]_0^\infty + s\int_0^\infty f(t)e^{-st}\,dt. The boundary term at t=0 contributes -f(0); the term at t\to\infty vanishes for every s in the transform's region of convergence, since e^{-st} decays faster than f can grow there. What remains is sF(s) - f(0).

For f'', apply this same rule to g = f': \mathcal{L}\{g'(t)\} = s\mathcal{L}\{g(t)\} - g(0) = s\big(sF(s)-f(0)\big) - f'(0) = s^2F(s) - sf(0) - f'(0). Each further derivative repeats the same integration-by-parts step once more, producing one more factor of s and one more initial-condition term — the pattern the general nth- order statement packages. This is exactly what makes the transform solve an IVP directly (§12): every derivative in the equation becomes a power of s times F(s), with the initial conditions built in from the start rather than applied afterward.

Fix (x_0,t_0) and let B_\tau=\{x:|x-x_0|\le c(t_0-\tau)\} be the ball obtained by shrinking the backward light cone from x_0,t_0 down to a point as \tau runs from 0 to t_0. Define the energy carried inside this shrinking ball, E(\tau) = \frac12\int_{B_\tau}\Big(u_t^2+c^2|\nabla u|^2\Big)\,dV. Differentiating under the integral sign, the boundary of B_\tau itself moves inward at speed c, which contributes a boundary term in addition to the usual energy-identity term. Using the divergence theorem exactly as in the bounded-string energy argument, the interior contribution vanishes by the PDE itself; the moving-boundary contribution can be shown (by choosing the ball's shrink rate to match the wave speed exactly) to be \le0 throughout. Hence E(\tau) is non-increasing. Since E(0)=0 whenever f,g vanish identically on the starting ball B_0=\{x:|x-x_0|\le ct_0\}, and E cannot decrease below 0, E(\tau)\equiv0 on all of [0,t_0] — forcing u\equiv const inside every shrinking ball, in particular at the final point (x_0,t_0). So data vanishing on B_0 produces u(x_0,t_0)= the value that vanishing data would give (zero, after accounting for the homogeneous case), which is exactly the statement that only data inside B_0 can affect u(x_0,t_0).

The range-of-influence statement is the same fact read in the opposite direction: if (x,t) is not in the solid forward cone from (x_0,0), i.e. |x-x_0|>ct, then x_0 does not lie in (x,t)'s own domain of dependence \{y:|y-x|\le ct\}, so data at x_0 alone cannot affect u(x,t).

It suffices to solve the "special" case g=0 first: define v(x,t)=t\,M_f(x,ct) and check directly (by expanding the spherical average and differentiating in t using the divergence theorem on the sphere) that v_{tt}=c^2\Delta v with v(x,0)=0, v_t(x,0)=f(x). The general initial-velocity case u(x,0)=0, u_t(x,0)=g is then this same v with f replaced by g; the general case with both f and g present follows by superposition, and differentiating the velocity-only solution built from f in t recovers the position-only case (a reduction lemma the same in spirit to how the 1D theory reduces the general case to two simpler ones). Summing the two pieces gives the boxed formula.

The method of descent. If the initial data f,g happen not to depend on the third coordinate x_3 (a genuinely 2D initial state, viewed as a special 3D one), the sphere-average integral above can be converted into an integral over the disk of radius ct in the x_1x_2-plane (each vertical line through the disk meets the sphere twice, and the two contributions combine into one planar integral) — producing the 2D wave equation's own solution formula as a special case of this one. Applying the same reduction a second time (data independent of x_2 as well) collapses the disk integral to a single line integral, recovering d'Alembert's formula exactly.

Bessel's equation of order m, written in Sturm-Liouville form (ch.04 §8), is (rR')'-\frac{m^2}rR+\lambda rR=0 with p(r)=r, q(r)=-m^2/r, \sigma(r)=r on [0,1] -- singular at r=0, since p(0)=\sigma(0)=0. R_k(r)=J_m(j_{m,k}r) solves this equation with \lambda_k=j_{m,k}^2 (substituting s=j_{m,k}r turns the equation into Bessel's equation in s, solved by J_m(s)), stays bounded as r\to0 (the singular endpoint's own admissibility condition, ch.04 §8), and satisfies R_k(1)=J_m(j_{m,k})=0 by definition of j_{m,k} as a zero.

Ch.04 §2's own orthogonality argument applies unchanged: for two eigenfunctions R_k,R_l with \lambda_k\neq\lambda_l, (\lambda_k-\lambda_l)\int_0^1\sigma R_kR_l\,dr=\big[p(R_k'R_l-R_kR_l')\big]_0^1. At r=1, both R_k(1)=R_l(1)=0, killing that boundary contribution directly (no separated condition needed, just the shared zero). At r=0, p(0)=0 kills the boundary term regardless of R_k,R_l's own values there, exactly ch.04 §8's singular-endpoint mechanism. Both ends contribute zero, so \int_0^1rR_k(r)R_l(r)\,dr=0 whenever \lambda_k\neq\lambda_l, i.e. whenever k\neq l (the zeros are strictly increasing, so distinct k give distinct \lambda_k=j_{m,k}^2).

Route 1: Sturm-Liouville. Legendre's equation in Sturm-Liouville form (ch.04 §8) is [(1-x^2)y']'+\lambda y=0 with p(x)=1-x^2, q=0, \sigma=1 on [-1,1] -- singular at BOTH endpoints, since p(\pm1)=0. P_n solves this with \lambda_n=n(n+1) and stays bounded (in fact, finite-degree polynomials are bounded everywhere) at both x=\pm1. Ch.04 §2's boundary-term argument, (\lambda_n-\lambda_k)\int_{-1}^1P_nP_k\,dx=\big[(1-x^2)(P_n'P_k-P_nP_k')\big]_{-1}^1, has its right side killed at BOTH endpoints by p(\pm1)=0 alone (ch.04 §8's singular mechanism, needing no separated condition at either end). Since n\neq k gives \lambda_n=n(n+1)\neq k(k+1)=\lambda_k (true for nonnegative integers n\neq k), the integral is forced to 0.

Route 2: Rodrigues' formula. P_n(x)=\dfrac1{2^nn!}\dfrac{d^n}{dx^n}(x^2-1)^n (direct differentiation reproduces P_0=1,P_1=x,P_2=\frac{3x^2-1}2,\ldots). Suppose k<n. Integrate \int_{-1}^1P_n(x)P_k(x)\,dx by parts n times, moving all n derivatives off (x^2-1)^n and onto P_k(x) (a degree-k polynomial) at each step. Every boundary term vanishes: (x^2-1)^n and each of its first n-1 derivatives all still have a factor of (x^2-1) surviving (differentiating (x^2-1)^n fewer than n times never fully removes that factor), which is 0 at x=\pm1. After n integrations by parts, \int_{-1}^1P_n(x)P_k(x)\,dx=\frac{(-1)^n}{2^nn!}\int_{-1}^1(x^2-1)^n\,\frac{d^nP_k}{dx^n}\,dx. P_k has degree k<n, so its n-th derivative is identically 0 -- the integral vanishes for exactly that reason, independent of the Sturm-Liouville argument above.

\nabla^2\Phi=0 away from the origin is a direct computation (the curvilinear Laplacian formula, chapter 04, applied in spherical coordinates to \Phi=-1/(4\pi\rho), a function of \rho alone). Concentrated entirely AT the origin: integrate \nabla^2\Phi over a small ball of radius \varepsilon and apply the Divergence Theorem — \iiint_{B_\varepsilon}\nabla^2\Phi\,dV = \oiint_{\partial B_\varepsilon}\nabla\Phi\cdot\hat{\mathbf n}\,dS. Direct computation gives \nabla\Phi = \dfrac{\mathbf x}{4\pi\|\mathbf x\|^3}, so on the sphere of radius \varepsilon, \nabla\Phi\cdot\hat{\mathbf n} = \dfrac{1}{4\pi\varepsilon^2} (constant), and the flux integral is \dfrac{1}{4\pi\varepsilon^2}\times(\text{area }4\pi\varepsilon^2)=1, independent of \varepsilon — exactly the defining property of a unit point source, \delta(\mathbf x).

The integral formula follows from Green's second identity applied to u and G(\cdot,\mathbf x_0) on V minus a small ball around \mathbf x_0 (to avoid G's own singularity), then shrinking that ball to zero — the singular contribution converges to exactly u(\mathbf x_0) by the same point-source argument above, while G=0 on \partial V kills the other boundary term.

The heat equation has a scaling symmetry: if u(x,t) solves u_t=ku_{xx}, so does u(\lambda x,\lambda^2t) for any \lambda>0 (the \lambda^2 on t exactly compensates the two x-derivatives). A solution built from data concentrated at a single point has no length scale of its own to break this symmetry, so it should be expressible as a function of the single combined variable \eta=x/\sqrt t alone, up to an overall power of t fixing units. Try u(x,t) = \frac1{\sqrt t}\,\varphi(\eta),\qquad \eta=\frac x{\sqrt t}.

Compute the two derivatives the equation needs, using \partial\eta/\partial t=-\tfrac12\eta t^{-1} and \partial\eta/\partial x=t^{-1/2}: u_t = -\tfrac12t^{-3/2}\big[\varphi(\eta)+\eta\varphi'(\eta)\big],\qquad u_{xx} = t^{-3/2}\varphi''(\eta). Substituting into u_t=ku_{xx} and cancelling the common t^{-3/2}: k\varphi''+\tfrac12\eta\varphi'+\tfrac12\varphi=0,\qquad\text{i.e.}\qquad k\varphi'' + \tfrac12(\eta\varphi)' = 0 (the last step using (\eta\varphi)'=\varphi+\eta\varphi') -- the self-similarity ansatz has turned the PDE into an ODE in \eta alone.

Solving the ODE. Integrating once, k\varphi'+\tfrac12\eta\varphi=A for a constant A; a fundamental solution and all its derivatives decay to 0 as \eta\to\pm\infty, which forces A=0. What remains is separable: \varphi'/\varphi=-\eta/(2k), so \ln\varphi=-\eta^2/(4k)+\text{const}, i.e. \varphi(\eta)=Be^{-\eta^2/(4k)} for a constant B. Substituting back, u(x,t) = \frac Bt{\vphantom{1}}^{1/2}e^{-x^2/(4kt)}.

Fixing B. A fundamental solution's defining property is that it represents unit total heat for every t: \int_{-\infty}^\infty u(x,t)\,dx=1. Substituting s=x/\sqrt{4kt}, \int_{-\infty}^\infty Bt^{-1/2}e^{-x^2/(4kt)}\,dx = B\sqrt{4\pi k} (a standard Gaussian integral), so B=1/\sqrt{4\pi k}, giving S(x,t)=\frac1{\sqrt{4\pi kt}}e^{-x^2/(4kt)}, the claimed formula.

The wave operator \partial_t^2-c^2\partial_x^2 factors as (\partial_t-c\partial_x)(\partial_t+c\partial_x), the product of two first-order transport operators (ch.02's own method of characteristics, run on each factor in turn). Equivalently, and more directly, change to characteristic coordinates \xi=x+ct, \eta=x-ct. The chain rule gives \partial_x=\partial_\xi+\partial_\eta and \partial_t=c(\partial_\xi-\partial_\eta), so u_{tt}-c^2u_{xx} = c^2\big[(\partial_\xi-\partial_\eta)^2-(\partial_\xi+\partial_\eta)^2\big]u = -4c^2u_{\xi\eta}. The equation becomes u_{\xi\eta}=0: integrating in \eta first shows u_\xi doesn't depend on \eta, then integrating in \xi gives u(\xi,\eta) = F(\xi)+G(\eta) = F(x+ct)+G(x-ct) for some single-variable functions F,G -- every solution of the wave equation on the whole line is a right-traveling profile plus a left-traveling profile, with no other freedom.

Matching the initial data. u(x,0)=F(x)+G(x)=f(x). Differentiating u in t and setting t=0: u_t(x,0)=cF'(x)-cG'(x)=g(x), so F'(x)-G'(x)=g(x)/c; integrating from 0 to x, F(x)-G(x)=\frac1c\int_0^xg(s)\,ds+C for a constant C. Solving the two equations for F and G: F(x)=\frac12\Big[f(x)+\frac1c\int_0^xg(s)\,ds+C\Big],\qquad G(x)=\frac12\Big[f(x)-\frac1c\int_0^xg(s)\,ds-C\Big]. The constant C cancels in the final substitution: with u(x,t)=F(x+ct)+G(x-ct), u(x,t) = \frac{f(x+ct)+f(x-ct)}2 + \frac1{2c}\left[\int_0^{x+ct}g(s)\,ds-\int_0^{x-ct}g(s)\,ds\right] = \frac{f(x-ct)+f(x+ct)}2 + \frac1{2c}\int_{x-ct}^{x+ct}g(s)\,ds, the claimed formula.

Heat equation. For u_t-ku_{xx}=F(x,t), u(x,0)=0, let v(x,t;s) solve the HOMOGENEOUS heat equation v_t=kv_{xx} for t>s with v(x,s;s)=F(x,s) (the source at time s, injected as an initial condition at that instant). Define u(x,t) = \int_0^t v(x,t;s)\,ds. Differentiating under the integral sign (Leibniz's rule, the boundary term at t=s contributing exactly v(x,t;t)=F(x,t)): u_t = v(x,t;t) + \int_0^tv_t(x,t;s)\,ds = F(x,t) + \int_0^tkv_{xx}(x,t;s)\,ds = F(x,t)+ku_{xx}, using v_t=kv_{xx} inside the integral -- exactly u_t-ku_{xx}=F. And u(x,0)=\int_0^0(\cdots)\,ds=0 automatically. By linearity, adding this to the solution of the homogeneous problem with the ACTUAL initial data (§4's whole-line heat solution) handles a heat equation with both nonzero initial data and a source term at once.

Wave equation. For u_{tt}-c^2u_{xx}=F(x,t), u(x,0)=u_t(x,0)=0, the same idea applies ONE DERIVATIVE HIGHER, since the wave equation is second order in t: let w(x,t;s) solve the homogeneous wave equation for t>s with w(x,s;s)=0 and w_t(x,s;s)=F(x,s) -- the source injected as an initial VELOCITY, not an initial position, at time s. Define u(x,t)=\int_0^tw(x,t;s)\,ds. Differentiating once, u_t = w(x,t;t) + \int_0^tw_t(x,t;s)\,ds = 0 + \int_0^tw_t(x,t;s)\,ds (the boundary term vanishes because w(x,s;s)=0 is exactly how each slice was built). Differentiating again, u_{tt} = w_t(x,t;t) + \int_0^tw_{tt}(x,t;s)\,ds = F(x,t) + \int_0^tc^2w_{xx}(x,t;s)\,ds = F(x,t)+c^2u_{xx} (this time the boundary term is w_t(x,t;t)=F(x,t), the injected velocity, and w_{tt}=c^2w_{xx} is used inside the integral) -- exactly u_{tt}-c^2u_{xx}=F, with u(x,0)=0 and u_t(x,0)=0 both automatic from the empty-interval integral.

Both cases are the SAME construction: differentiate the accumulated integral, let the boundary term at t=s (which derivative order it lands on depends on the equation's own order in t) absorb the source, and let the homogeneous equation inside the integral absorb everything else.

Argue by contradiction, and suppose the set is linearly dependent. Among all dependence relations c_1 v_1 + c_2 v_2 + \cdots + c_k v_k = \mathbf{0} with not every c_i zero, pick one using the fewest nonzero vectors — say it uses exactly m of them, relabeled so they're v_1, \ldots, v_m with every coefficient c_1, \ldots, c_m nonzero: c_1 v_1 + c_2 v_2 + \cdots + c_m v_m = \mathbf{0}. \qquad (\star) Since eigenvectors are nonzero by definition, m \geq 2 — a single nonzero vector can't satisfy c_1 v_1 = \mathbf{0} with c_1 \neq 0.

Apply A to both sides of (\star) and use Av_i = \lambda_i v_i: c_1 \lambda_1 v_1 + c_2 \lambda_2 v_2 + \cdots + c_m \lambda_m v_m = \mathbf{0}. \qquad (\star\star) Now multiply (\star) by \lambda_m and subtract (\star\star) from it, so the v_m terms cancel: c_1(\lambda_m - \lambda_1) v_1 + c_2(\lambda_m - \lambda_2) v_2 + \cdots + c_{m-1}(\lambda_m - \lambda_{m-1}) v_{m-1} = \mathbf{0}. This is a dependence relation on only m - 1 of the vectors — fewer than the minimal count m assumed above. For it to be genuinely shorter (not just relabeled with a zero coefficient tacked on), at least one of its coefficients must be nonzero. Since the eigenvalues are distinct, \lambda_m - \lambda_i \neq 0 for every i < m, and c_i \neq 0 by choice of (\star), so each coefficient c_i(\lambda_m - \lambda_i) is nonzero for i = 1, \ldots, m-1. That contradicts m being the fewest vectors any dependence relation among v_1, \ldots, v_k can use.

The contradiction shows no dependence relation exists, so \{v_1, \ldots, v_k\} is linearly independent.

Sanity check: take k = 2. The argument reduces to: if c_1 v_1 + c_2 v_2 = \mathbf{0} with c_1 \neq 0, applying A gives c_1\lambda_1 v_1 + c_2\lambda_2 v_2 = \mathbf{0}; subtracting \lambda_2 times the first equation leaves c_1(\lambda_1 - \lambda_2)v_1 = \mathbf{0}, and since \lambda_1 \neq \lambda_2 and v_1 \neq \mathbf{0}, this forces c_1 = 0 — contradicting c_1 \neq 0. So two eigenvectors for distinct eigenvalues are always independent, exactly the base case the induction-style argument above builds on.

(\Leftarrow) n independent eigenvectors give A = PDP^{-1}. Suppose v_1, \dots, v_n are linearly independent eigenvectors of A, with Av_i = \lambda_i v_i. Let P be the matrix with these vectors as its columns, and let D = \operatorname{diag}(\lambda_1, \dots, \lambda_n). Multiplying A into P column by column, AP = A\begin{bmatrix} v_1 & \cdots & v_n \end{bmatrix} = \begin{bmatrix} Av_1 & \cdots & Av_n \end{bmatrix} = \begin{bmatrix} \lambda_1 v_1 & \cdots & \lambda_n v_n \end{bmatrix}. Multiplying P by D on the right scales each column of P by the matching diagonal entry of D, producing the very same matrix: PD = \begin{bmatrix} v_1 & \cdots & v_n \end{bmatrix} \operatorname{diag}(\lambda_1,\dots,\lambda_n) = \begin{bmatrix} \lambda_1 v_1 & \cdots & \lambda_n v_n \end{bmatrix}. So AP = PD. Because P's columns are n linearly independent vectors in \mathbb{R}^n, P is invertible, and right-multiplying both sides by P^{-1} gives A = PDP^{-1}.

(\Rightarrow) A = PDP^{-1} gives n independent eigenvectors. Suppose A = PDP^{-1} with P invertible and D = \operatorname{diag}(\lambda_1, \dots, \lambda_n). Then AP = PD, and reading that equation column by column exactly as above shows Av_i = \lambda_i v_i, where v_i is the i-th column of P. Each v_i is nonzero — a zero column would make P singular — so v_1,\dots,v_n are eigenvectors of A. Since P is invertible, its columns are linearly independent, which produces the needed n independent eigenvectors.

Sanity check: take A = \begin{bmatrix} 4 & 1 \\ 0 & 3 \end{bmatrix}, upper triangular, so its eigenvalues sit on the diagonal: \lambda_1 = 4, \lambda_2 = 3. Solving (A - 4I)v = 0 gives eigenvector v_1 = (1,0), and (A-3I)v = 0 gives v_2 = (1,-1). These two are not scalar multiples of each other, so they're independent, and the theorem guarantees A = PDP^{-1} with P = \begin{bmatrix} 1 & 1 \\ 0 & -1 \end{bmatrix} and D = \begin{bmatrix} 4 & 0 \\ 0 & 3 \end{bmatrix} — which a direct multiplication PDP^{-1} confirms reproduces A.

Conservation of mass and momentum, with p_x=c(\rho)^2\rho_x from the isentropic pressure law, give the quasilinear system \rho_t+u\rho_x+\rho u_x=0,\qquad u_t+uu_x+\frac{c(\rho)^2}\rho\rho_x=0. Seek r=u+\phi(\rho) constant along dx/dt=u+c, i.e. r_t+(u+c)r_x=0. Expanding and substituting the two PDEs above for u_t and \rho_t collapses the equation to u_x\big[c-\rho\phi'(\rho)\big]+\rho_x\Big[\phi'(\rho)c-\frac{c^2}\rho\Big]=0. Since u_x,\rho_x are otherwise free, both bracketed coefficients must vanish, and both vanish for the same choice \phi'(\rho)=c(\rho)/\rho. For p=K\rho^\gamma, c(\rho)/\rho=\sqrt{K\gamma}\,\rho^{(\gamma-3)/2}, and integrating gives \phi(\rho)=2c(\rho)/(\gamma-1), so r=u+2c/(\gamma-1) is constant along C^+. The identical computation with u-c in place of u+c gives s=u-2c/(\gamma-1) constant along C^-.

The differential form u_t+f(u)_x=0 already assumes u is smooth enough to differentiate, so it says nothing by itself about a jump — the jump condition comes from going back to the integral form the differential equation was derived from. For any fixed interval [x_1,x_2] (with x_1<s(t)<x_2), conservation of the total amount of u reads \frac{d}{dt}\int_{x_1}^{x_2}u\,dx = f(u(x_1,t)) - f(u(x_2,t)), the flux entering at x_1 minus the flux leaving at x_2. Split the integral at the shock and differentiate each piece with the Leibniz rule, which picks up a boundary term wherever the limit of integration itself moves: \frac{d}{dt}\int_{x_1}^{x_2}u\,dx = \frac{d}{dt}\int_{x_1}^{s(t)}u\,dx + \frac{d}{dt}\int_{s(t)}^{x_2}u\,dx = \int_{x_1}^{s(t)}u_t\,dx + u_L\,s'(t) + \int_{s(t)}^{x_2}u_t\,dx - u_R\,s'(t). Now shrink the interval onto the shock, x_1\to s(t)^- and x_2\to s(t)^+. Both remaining integrals vanish in this limit (u_t stays bounded on either side, away from the shock itself, so each integral is over a shrinking interval of a bounded integrand), leaving \frac{d}{dt}\int_{x_1}^{x_2}u\,dx \;\longrightarrow\; \sigma\,u_L - \sigma\,u_R = -\sigma\,\llbracket u\rrbracket. The right side of the conservation law shrinks the same way: f(u(x_1,t))-f(u(x_2,t)) \to f(u_L)-f(u_R) = -\llbracket f(u)\rrbracket. Equating the two limits and canceling the shared minus sign gives \sigma\,\llbracket u\rrbracket=\llbracket f(u)\rrbracket — the shock speed is exactly the ratio of the flux jump to the density jump, because that is the only speed at which the total amount of u on each side of the shock is still conserved through the discontinuity.

Add a small diffusive term, u_t+f(u)_x=\varepsilon u_{xx} (\varepsilon>0 small), and seek a traveling-wave profile u(x,t)=U(\xi), \xi=x-\sigma t, connecting U(-\infty)=u_L to U(\infty)=u_R with U'(\pm\infty)=0 -- a smooth stand-in for the sharp shock, which the true (inviscid) shock is understood as the limit of as \varepsilon\to0^+. Substituting into the viscous equation: -\sigma U'+f(U)'=\varepsilon U''. Integrating once in \xi, -\sigma U+f(U) = \varepsilon U' + C for a constant C. Evaluating as \xi\to\pm\infty (where U'\to0) gives C=f(u_L)-\sigma u_L=f(u_R)-\sigma u_R -- the two expressions must be equal, which is exactly the Rankine-Hugoniot condition, consistently recovered.

Define g(U)=f(U)-\sigma U, so \varepsilon U'=g(U)-g(u_L) and, by the above, g(u_L)=g(u_R)=C. Say u_L>u_R (the profile decreases from u_L to u_R as \xi increases). For U'<0 throughout, g(U)-g(u_L) must stay negative for every U strictly between u_R and u_L -- i.e. g dips BELOW its own shared boundary value C in the interior. Near \xi\to-\infty (U just below u_L), this forces g to be non-decreasing as U increases toward u_L: g'(u_L)\ge0, i.e. f'(u_L)\ge\sigma. Near \xi\to\infty (U just above u_R), the same requirement forces g'(u_R)\le0, i.e. f'(u_R)\le\sigma. Together, f'(u_L)\ge\sigma\ge f'(u_R) -- strict, generically, since equality throughout would make the viscous profile fail to connect the two states in finite \xi.

A jump where this FAILS admits no such viscous traveling-wave profile approximating it -- the vanishing-viscosity limit simply does not produce that jump, which is exactly the sense in which it is not the physically correct weak solution.

Let g(r) = \frac{1}{2\pi r}\oint_{\|\mathbf x-\mathbf a\|=r}u\,ds = \frac{1}{2\pi}\int_0^{2\pi}u(\mathbf a+r(\cos\theta,\sin\theta))\, d\theta. Differentiate under the integral sign (chapter 08) with respect to r: g'(r) = \frac{1}{2\pi}\int_0^{2\pi}\nabla u\cdot(\cos\theta,\sin\theta)\,d\theta = \frac{1}{2\pi r}\oint \nabla u\cdot\hat{\mathbf n}\,ds (recognizing (\cos\theta,\sin\theta) as the outward unit normal). By the Divergence Theorem (chapter 00), this boundary flux equals \frac{1}{2\pi r}\iint_{\text{disk}} \nabla^2u\,dA = 0, since u is harmonic. So g'(r)=0 for every r — g is CONSTANT — and \lim_{r\to0}g(r)=u(\mathbf a) (the average over a shrinking circle approaches the center value by continuity), giving g(r)=u(\mathbf a) for every r.

Suppose u attains an interior maximum at \mathbf a\in D, with u(\mathbf a) = M the largest value u takes anywhere on the closure. By the Mean Value Property, u(\mathbf a) equals the average of u over any small circle centered at \mathbf a. But every value on that circle is \le M (since M is the overall maximum) — an average of values \le M equaling exactly M forces EVERY value on the circle to equal M too (otherwise the average would be strictly less). This holds for circles of every radius small enough to stay in D, so u \equiv M on an entire neighborhood of \mathbf a. The set where u=M is therefore open; it's also closed (as the preimage of a point under a continuous function); on a connected region, a nonempty set that's both open and closed must be the WHOLE region — so u\equiv M everywhere, the constant case. If u isn't constant, no interior maximum can occur, forcing the maximum onto the boundary.