Partial Differential Equations

Solve du/dx = 0 for a function of one variable and you get u = C, one arbitrary constant. Now ask for all functions u(x, y) with u_x = 0, where the subscript means the partial derivative. The answer is every function of y alone: u = y^3 works, so does u = sin(y), so does u = e^(y^2) - 7. The solution set is not a one-parameter family. It is indexed by an entire arbitrary function.

That single observation drives everything in this course. An ordinary differential equation of order n needs n numbers to pin down a solution, and you supply them at a point. A partial differential equation needs a function's worth of information, and you supply it along a curve or a surface. Where that data may be placed, and whether placing it there gives one solution or none or infinitely many, is the central question of the subject.

Vocabulary, fixed once

A partial differential equation relates an unknown function of two or more variables to its partial derivatives. Its order is the highest derivative appearing. The equation is linear if the unknown and its derivatives appear only to the first power and are not multiplied by one another, with coefficients depending on the independent variables alone. So

a(x,y) u_xx + b(x,y) u_xy + c(x,y) u_yy + d(x,y) u_x + e(x,y) u_y + f(x,y) u = g(x,y)

is the general second-order linear equation in two variables. It is homogeneous when g is identically zero. If the top-order coefficients depend on u or its lower derivatives, the equation is quasilinear; if they involve derivatives of the same order nonlinearly, it is fully nonlinear.

The reason linearity is worth naming first is that it buys a structural theorem, and the proof is two lines.

Proposition (superposition). If u_1 and u_2 both solve a linear homogeneous PDE, so does c_1 u_1 + c_2 u_2 for any constants.

Proof. Write the equation as L[u] = 0, where L is the operator assembled from the derivative terms. Differentiation is linear, so L[c_1 u_1 + c_2 u_2] = c_1 L[u_1] + c_2 L[u_2] = 0 + 0 = 0. This completes the proof.

Everything in Modules 2 through 5 depends on that proposition. Separation of variables produces infinitely many simple solutions and then adds them up; Fourier series, eigenfunction expansions and Green's functions are all superposition carried to a limit. For a nonlinear equation none of it is available, which is why Burgers' equation in the next module has to be handled by an entirely different technique.

General solutions carry arbitrary functions

Work two examples completely, because the pattern they show is what side conditions have to overcome.

Example 1. Solve u_xy = 0. Read it as (u_x)_y = 0: the function u_x has zero y-derivative, so it depends on x alone, say u_x = p(x). Integrating in x with y held fixed gives u = F(x) + G(y), where F is an antiderivative of p and G is the constant of integration, which may depend on y. Two arbitrary functions, and you can check directly that any such u satisfies the equation.

Example 2. Solve the transport equation u_t + c u_x = 0 for a constant c. Claim: every solution has the form u(x, t) = f(x - ct). Verification is immediate by the chain rule: u_t = -c f'(x - ct) and u_x = f'(x - ct), so u_t + c u_x = 0. That this is the general solution follows from a change of variables. Put s = x - ct and r = t, so that u(x,t) = v(s, r). Then u_t = -c v_s + v_r and u_x = v_s, so the equation becomes v_r = 0, which says v depends on s alone.

The picture is worth fixing now, because it recurs. The graph of f slides to the right at speed c without changing shape. Information travels along the lines x - ct = constant, and the value of u is constant along each such line. Those lines are the characteristics of the equation, and Module 2 is built on them.

Side conditions are the whole game

Since a general solution contains an arbitrary function, a problem is not posed until you say what that function is. Two kinds of extra information appear.

Initial conditions prescribe the solution, and possibly some of its time derivatives, at one instant: u(x, 0) = phi(x) for the heat equation, and both u(x, 0) = phi(x) and u_t(x, 0) = psi(x) for the wave equation. Boundary conditions prescribe behaviour at the edges of the spatial region for all time: u(0, t) = u(L, t) = 0 for a string clamped at both ends.

How much data is right depends on the equation, and getting it wrong is not a small error. The wave equation needs two initial conditions because it is second order in time; the heat equation needs one. Laplace's equation has no time variable at all and needs data all the way around the boundary of a region. Ask for the wrong amount and you get either no solution or a whole family, and Lesson 16 makes this precise under the heading of well-posedness.

Key idea: a PDE by itself almost never determines a solution. The equation plus the right side conditions on the right set does.

The three equations that set the agenda

These three differ by the placement of a few symbols and behave in almost opposite ways. A corner in the initial shape of a plucked string travels along the string forever; the same corner in an initial temperature profile is gone the instant after t = 0. A wave can be run backwards to recover its past; a temperature distribution cannot, which is the subject of Lesson 8. Explaining why such small differences in the symbols produce such different physics is the job of the classification in the next lesson.

Where these equations come from

One derivation covers two of them. Let u(x, t) be a density, meaning that the amount of stuff in the interval from a to b is the integral of u over that interval. Let q(x, t) be the flux, the rate at which stuff crosses the point x to the right. If nothing is created or destroyed, then the amount in the interval changes only through the two ends:

d/dt of the integral from a to b of u dx = q(a, t) - q(b, t).

Write the right side as minus the integral of q_x from a to b, by the fundamental theorem of calculus, and move the time derivative inside the integral on the left. Then the integral of u_t + q_x over [a, b] vanishes for every choice of a and b. A continuous function whose integral over every interval is zero is identically zero, so

u_t + q_x = 0.

That is the conservation law, and now the physics enters through a constitutive relation for the flux. If the stuff is simply carried along at speed c, then q = c u and the law becomes the transport equation u_t + c u_x = 0. If instead the stuff diffuses down its own gradient, Fourier's law says q = -k u_x with k positive, and the law becomes

u_t - k u_xx = 0,

the heat equation. The minus sign in Fourier's law is doing something specific: heat flows from hot to cold, so the flux points opposite to the temperature gradient. Reverse that sign and you get the backward heat equation, which Lesson 8 shows is catastrophically ill posed.

Common misconceptions

• "A general solution is the goal, as it is for ODEs." For most PDEs a general solution is either unavailable or useless, because it involves arbitrary functions that the side conditions must determine anyway. The working goal is a solution of a specific initial or boundary value problem.
• "More initial data is safer." Prescribing u and u_t at t = 0 for the heat equation generally admits no solution at all, because the equation already determines u_t from u_xx. Too much data is as fatal as too little.
• "Superposition always works." Only for linear equations. Adding two solutions of u_t + u u_x = 0 gives a function that solves nothing, since the nonlinear term produces a cross product.
• "Order tells you how many conditions to impose." It tells you the count in the time direction for evolution equations, but Laplace's equation is second order and needs one condition around a closed boundary rather than two at an initial time. The type, not the order, decides.

Try it: pin down the arbitrary functions

Exercise. Find the solution of u_xy = 0 on the square 0 <= x, y <= 1 satisfying u(x, 0) = x^2 and u(0, y) = sin(y).

Solution. The general solution is u = F(x) + G(y). Setting y = 0 gives F(x) + G(0) = x^2, so F(x) = x^2 - G(0). Setting x = 0 gives F(0) + G(y) = sin(y), so G(y) = sin(y) - F(0). Now F(0) = -G(0), so G(y) = sin(y) + G(0) and

u(x, y) = x^2 - G(0) + sin(y) + G(0) = x^2 + sin(y).

The unknown constant cancels, as it must, since F and G are only determined up to adding and subtracting the same constant. Check the two conditions: u(x, 0) = x^2 + 0 and u(0, y) = 0 + sin(y). Both hold.

Change one input. Try instead u(x, 0) = x^2 and u(0, y) = sin(y) + 1. Now the two conditions disagree at the corner: the first says u(0,0) = 0 and the second says u(0,0) = 1. No function has two values at a point, so the problem has no solution. Compatibility of the data at corners is a real constraint, not a technicality.

What you now know

A PDE relates a multivariable unknown to its partial derivatives, and its general solution carries arbitrary functions rather than arbitrary constants, so side conditions do the work of selection. Linear homogeneous equations obey superposition, which is the licence for every series method in this course. A conservation law in one dimension is u_t + q_x = 0, and choosing the flux gives either transport or diffusion. The wave, heat and Laplace equations are the three benchmarks, and they behave in almost opposite ways.

Looking ahead

Why should u_tt - c^2 u_xx = 0 and u_t - k u_xx = 0 behave so differently when they differ by one derivative? The next lesson answers that with an algebraic test on the second-order coefficients: compute b^2 - 4ac and read off whether the equation is hyperbolic, parabolic or elliptic. The test looks like a bookkeeping exercise and is nothing of the sort. It predicts whether signals travel at finite speed, whether rough data is instantly smoothed, and where you are allowed to put your data.

Sources

1. Wikipedia contributors. (n.d.). Partial differential equation. Wikipedia. en.wikipedia.org
2. Weisstein, E. W. (n.d.). Partial differential equation. Wolfram MathWorld. mathworld.wolfram.com
3. Dawkins, P. (n.d.). The heat equation. Paul's Online Notes, Differential Equations. tutorial.math.lamar.edu
4. Lebl, J. (2026). First order PDEs (Section 1.9). Notes on Diffy Qs: Differential Equations for Engineers. jirka.org
5. Strauss, W. A. (2008). Partial differential equations: An introduction (2nd ed.), Chapter 1. John Wiley and Sons.

Here are two equations that differ by a single sign:

u_tt - c^2 u_xx = 0 and u_tt + c^2 u_xx = 0.

The first is the wave equation. Give it an initial displacement and an initial velocity and it produces a unique solution that depends continuously on the data, for all time, forwards and backwards. The second is Laplace's equation with t playing the role of a second space coordinate. Give it the same two pieces of initial data and, as Lesson 16 will show in detail, arbitrarily small changes in the data produce arbitrarily large changes in the solution after any positive time. One sign, and a physically sensible problem becomes a numerically hopeless one.

There is a test that separates these cases, it takes ten seconds to apply, and it is the most useful piece of bookkeeping in the subject.

The discriminant

Consider the general second-order linear equation in two variables,

a u_xx + b u_xy + c u_yy + (lower order terms) = 0,

where a, b, c may depend on x and y. Form the discriminant D = b^2 - 4ac at a point. The equation is called

• hyperbolic at that point if D > 0,
• parabolic if D = 0,
• elliptic if D < 0.

The names come from conic sections: replacing u_xx by X^2, u_xy by XY and u_yy by Y^2 turns the top-order part into a quadratic form, and b^2 - 4ac is exactly the quantity that decides whether aX^2 + bXY + cY^2 = 1 is a hyperbola, a parabola or an ellipse. The analogy is a mnemonic and nothing more; the content is in what the sign predicts.

Read the heat equation entry carefully. Written as u_t - k u_xx = 0 with independent variables t and x, there is no u_tt, so a = 0; there is no mixed term, so b = 0; the coefficient of u_xx is -k. Hence D = 0 and the equation is parabolic. The absence of a second time derivative is the whole reason.

Why the test is really about characteristics

The discriminant is not arbitrary. It counts the families of curves along which the equation carries information.

A curve y = y(x) is characteristic for the equation above if

a (dy/dx)^2 - b (dy/dx) + c = 0.

This quadratic in dy/dx has discriminant b^2 - 4ac, and its roots are dy/dx = (b plus or minus sqrt(b^2 - 4ac))/(2a). So a hyperbolic equation has two distinct real families of characteristics, a parabolic one has exactly one, and an elliptic one has none at all in the real plane.

Characteristics are the curves along which the solution can be non-smooth. A wave equation has two of them, moving left and right at speed c, and a kink in the initial data slides along one of them forever. Laplace's equation has none, and correspondingly its solutions are infinitely differentiable no matter how rough the boundary data. Why this matters: the sign of one number tells you whether the equation supports travelling discontinuities.

Reduction to canonical form, worked

Choosing coordinates along the characteristics turns a hyperbolic equation into the simplest possible shape. Take

u_xx + 4 u_xy + 3 u_yy = 0.

Here a = 1, b = 4, c = 3, so D = 16 - 12 = 4 and the equation is hyperbolic. The characteristic equation (dy/dx)^2 - 4(dy/dx) + 3 = 0 factors as (dy/dx - 3)(dy/dx - 1) = 0, giving dy/dx = 3 and dy/dx = 1. Integrating, the two families are y - 3x = constant and y - x = constant. So set

s = y - 3x, r = y - x, and write u(x,y) = v(s,r).

Now grind out the chain rule. Since s_x = -3, s_y = 1, r_x = -1, r_y = 1,

u_x = -3 v_s - v_r and u_y = v_s + v_r.

Differentiating again,

u_xx = 9 v_ss + 6 v_sr + v_rr, u_xy = -3 v_ss - 4 v_sr - v_rr, u_yy = v_ss + 2 v_sr + v_rr.

Substituting into the equation, the v_ss coefficient is 9 - 12 + 3 = 0, the v_rr coefficient is 1 - 4 + 3 = 0, and the v_sr coefficient is 6 - 16 + 6 = -4. The equation collapses to

v_sr = 0,

whose general solution, from the previous lesson, is v = F(s) + G(r). Translating back,

u(x, y) = F(y - 3x) + G(y - x).

Check it. For u = F(y - 3x) we have u_xx = 9F'', u_xy = -3F'', u_yy = F'', and 9 - 12 + 3 = 0. For u = G(y - x) we get 1 - 4 + 3 = 0. Both are solutions, and by superposition so is the sum. The canonical form did not simplify the equation by luck: the characteristic coordinates are precisely the ones that annihilate the pure second derivatives.

The same procedure applied to u_tt - c^2 u_xx = 0 gives characteristics x - ct = constant and x + ct = constant, and hence u = F(x - ct) + G(x + ct). That is d'Alembert's solution, and Lesson 4 derives it again with the initial data attached.

What each type predicts

The row on infinite speed deserves a comment, since it sounds physically absurd. Put a unit of heat at the origin of an infinite bar at t = 0. The solution, computed in Lesson 14, is a Gaussian that is strictly positive at every point for every t > 0. So a thermometer a light year away registers a nonzero, if unimaginably small, temperature rise instantly. That is a defect of the model, not of the mathematics, and it is the price of the diffusion approximation.

An equation that changes its mind

Types need not be global. The Euler-Tricomi equation

u_xx + x u_yy = 0

has a = 1, b = 0, c = x, so D = -4x. It is elliptic where x > 0, parabolic on the line x = 0, and hyperbolic where x < 0. It is not a curiosity: it models transonic flow, where the same physical region contains subsonic parts, governed by an elliptic equation, and supersonic parts, governed by a hyperbolic one. The mathematics is telling you that a shock can form on one side of the sonic line and not the other.

The upshot: classify before you choose a method. Separation of variables and Fourier series suit parabolic and elliptic problems on bounded domains; characteristics suit hyperbolic ones; and a scheme that is stable for one type may be useless for another, as Lesson 17 shows.

Common misconceptions

• "The classification is about the shape of the solutions." It is computed from the top-order coefficients only, and it predicts qualitative behaviour. Lower-order terms change solutions but never change the type.
• "The heat equation is parabolic because of the word diffusion." It is parabolic because a = b = 0 in the t, x variables, so D = 0. Any equation missing the second derivative in one variable and having a second derivative in the other lands in the same box.
• "An elliptic equation has complex characteristics, so you can use them." Formally the roots are complex, and complex characteristic coordinates do reduce elliptic equations to Laplace's equation. But there are no real curves along which information travels, which is exactly why elliptic solutions are so smooth.
• "Type is a global property of an equation." The coefficients may vary, so the type may vary, and the Euler-Tricomi equation changes type across a line.

Try it: classify and reduce

Exercise. Classify u_xx - 5 u_xy + 6 u_yy = 0 and find its general solution.

Solution. Here a = 1, b = -5, c = 6, so D = 25 - 24 = 1 > 0 and the equation is hyperbolic. The characteristic equation (dy/dx)^2 + 5(dy/dx) + 6 = 0 factors as (dy/dx + 2)(dy/dx + 3) = 0, giving dy/dx = -2 and dy/dx = -3, so the characteristics are y + 2x = constant and y + 3x = constant. Taking s = y + 2x and r = y + 3x reduces the equation to v_sr = 0 by the same computation as above, and the general solution is

u(x, y) = F(y + 2x) + G(y + 3x).

Verify with u = F(y + 2x): u_xx = 4F'', u_xy = 2F'', u_yy = F'', and 4 - 10 + 6 = 0. Correct.

Change one input. Replace the middle coefficient by -4: the equation u_xx - 4u_xy + 6u_yy = 0 has D = 16 - 24 = -8 < 0 and is elliptic. There are no real characteristics, no travelling-wave general solution, and no sensible initial value problem. Two coefficients apart, and the entire method changes.

The short version

Compute b^2 - 4ac from the second-order coefficients. Positive means hyperbolic, two real characteristic families, finite propagation speed, kinks that persist, and initial data on a line. Zero means parabolic, one family, instantaneous smoothing, and a one-way time direction. Negative means elliptic, no real characteristics, analytic interiors, and data around a closed boundary. Characteristic coordinates reduce a hyperbolic equation to v_sr = 0, whose solution is a sum of two travelling profiles. Types can vary from point to point.

Looking ahead

Characteristics were introduced here as the curves the classification counts. In the next lesson they become a solution method in their own right, for equations of first order. Along a characteristic curve a PDE degenerates into an ODE, which you can integrate; the difficulty, and the interest, is that for nonlinear equations the characteristics can run into each other, at which point the solution stops being a function and a shock is born.

Sources

1. Weisstein, E. W. (n.d.). Elliptic partial differential equation. Wolfram MathWorld. mathworld.wolfram.com
2. Weisstein, E. W. (n.d.). Hyperbolic partial differential equation. Wolfram MathWorld. mathworld.wolfram.com
3. Wikipedia contributors. (n.d.). Euler-Tricomi equation. Wikipedia. en.wikipedia.org
4. Dawkins, P. (n.d.). Terminology for partial differential equations. Paul's Online Notes, Differential Equations. tutorial.math.lamar.edu
5. Evans, L. C. (2010). Partial differential equations (2nd ed.), Chapter 1. American Mathematical Society.

Stand on a bridge over a single-lane road and watch the traffic. Where cars are packed close together they crawl; where the road is empty they move at the speed limit. So the flux of cars past a point depends on the density of cars there, and the conservation law from Lesson 1 becomes nonlinear. The consequence is visible from the bridge: the back edge of a traffic jam is sharp. Cars approaching at 30 metres per second brake in a few car lengths, and that near-discontinuity propagates backwards along the road even though every individual car keeps moving forward.

That sharp edge is a shock, and no linear equation can produce one. This lesson develops the method that solves first-order equations, applies it to linear cases where it works cleanly, and then follows it into the nonlinear case where it predicts its own failure.

Along the right curve, a PDE is an ODE

Consider a(x, t) u_x + u_t = c(x, t, u). Suppose you travel through the x, t plane along a curve x = X(t) and record the value of u as you go, so that U(t) = u(X(t), t). By the chain rule,

dU/dt = u_x (dX/dt) + u_t.

Compare that with the equation. If you choose the curve so that dX/dt = a(X, t), the right-hand side becomes exactly the left side of the PDE, so

dU/dt = c(X, t, U).

The partial differential equation has become a pair of ordinary differential equations: one for where the curve goes, one for how the solution changes along it. Those curves are the characteristics, and this is the method of characteristics. Key idea: a first-order PDE is a statement about how u changes along a specific family of curves, and nothing else.

Three linear examples, worked completely

Example 1: constant speed. Solve u_t + 2u_x = 0 with u(x, 0) = e^(-x^2).

The characteristic equation is dX/dt = 2, so X(t) = x_0 + 2t, and along it dU/dt = 0, so U keeps its initial value u(x_0, 0) = e^(-x_0^2). To find u at a point (x, t), solve for the characteristic passing through it: x_0 = x - 2t. Hence

u(x, t) = e^(-(x - 2t)^2).

Check: u_t = 4(x - 2t) e^(-(x-2t)^2) and u_x = -2(x - 2t) e^(-(x-2t)^2), so u_t + 2u_x = 0. The Gaussian bump slides right at speed 2 without changing shape.

Example 2: variable speed. Solve u_t + x u_x = 0 with u(x, 0) = f(x).

Now dX/dt = X, so X(t) = x_0 e^t, and dU/dt = 0 again. Inverting, x_0 = x e^(-t), so

u(x, t) = f(x e^(-t)).

Check: u_t = -x e^(-t) f'(x e^(-t)) and x u_x = x e^(-t) f'(x e^(-t)), and they cancel. Geometrically the characteristics fan out exponentially from the origin, so an initial profile is stretched rather than merely translated. Note that no characteristic ever crosses another, since they are the curves x e^(-t) = constant, and the solution stays single valued forever.

Example 3: decay along the way. Solve u_t + c u_x = -a u with u(x, 0) = f(x), where a is a positive constant.

The characteristics are still X(t) = x_0 + ct, but now dU/dt = -aU, so U(t) = U(0) e^(-at) = f(x_0) e^(-at). Substituting x_0 = x - ct,

u(x, t) = f(x - ct) e^(-at).

A pollutant carried downstream at speed c while decaying at rate a. The travelling and the decaying are entirely independent, which is exactly what the two ODEs say.

The quasilinear case: speed depends on the solution

Now take the inviscid Burgers equation

u_t + u u_x = 0, u(x, 0) = f(x).

The characteristic system is dX/dt = U and dU/dt = 0. Since U is constant along the curve, and its value there is f(x_0), the curve has constant speed f(x_0) and is therefore a straight line:

X(t) = x_0 + f(x_0) t, with u = f(x_0) along it.

That gives the solution implicitly: u(x, t) = f(x - u t). It is a genuine solution wherever you can solve that relation for u, and the trouble is that you cannot always do so.

The lines have different slopes. A characteristic carrying a large value moves faster than one carrying a small value. If f is increasing, the fast lines start ahead and the family spreads apart harmlessly. If f is decreasing anywhere, the fast lines start behind slower ones and must eventually catch up. Where two characteristics meet, they deliver two different values of u to the same point, and the solution stops being a function.

Breaking time, computed

Take the explicit profile

f(x) = 1 for x <= 0, f(x) = 1 - x for 0 <= x <= 1, f(x) = 0 for x >= 1.

For a starting point x_0 in [0, 1] the characteristic is X(t) = x_0 + (1 - x_0)t. Set t = 1 and every one of them gives X(1) = x_0 + 1 - x_0 = 1. The entire family, carrying every value of u from 0 to 1, arrives at the single point x = 1 at time t = 1. Before that time the solution is a well-defined ramp that steepens; at t = 1 the ramp becomes vertical; after it, the implicit relation has three solutions for u at some points and the classical solution has ceased to exist.

The general formula follows the same reasoning. Two nearby characteristics from x_0 and x_0 + h collide when x_0 + f(x_0) t = x_0 + h + f(x_0 + h) t, which as h goes to zero gives 1 + f'(x_0) t = 0. So the first collision happens at the breaking time

t_b = -1 / min f'(x_0),

taken over points where f' is negative, and never if f is nondecreasing everywhere. In the worked profile f' = -1 on the ramp, so t_b = 1, matching the direct computation.

After breaking: the shock condition

The equation does not stop being physically meaningful at t_b; the traffic on the road keeps flowing. What must be relaxed is the demand that u be differentiable. Return to the conservation form

u_t + (u^2/2)_x = 0,

which is equivalent to Burgers' equation wherever u is smooth, and require only that the integral balance from Lesson 1 hold. A jump discontinuity moving at speed s, with value u_L on the left and u_R on the right, is admissible provided the flux balance across it is satisfied. That is the Rankine-Hugoniot condition, which for a general flux q(u) reads

s = (q(u_L) - q(u_R)) / (u_L - u_R),

and for Burgers, with q = u^2/2, simplifies to s = (u_L + u_R)/2. The shock travels at the average of the two states it separates. In the worked example the shock forms at x = 1 with u_L = 1 and u_R = 0, so it moves right at speed 1/2.

The point: the method of characteristics does not merely fail for nonlinear equations. It tells you precisely when it will fail, where, and what has to replace it.

Common misconceptions

• "Characteristics are always straight lines." Only when the speed is constant along them, as for constant-coefficient linear equations and for Burgers. In u_t + x u_x = 0 they are exponential curves.
• "A shock means the model has broken down." The classical solution has broken down. The conservation law is still valid in integral form, and the shock is its correct solution; a real traffic jam has a genuinely sharp back edge.
• "Any decreasing initial profile shocks immediately." Breaking happens at t_b = -1/min f', which is finite but positive. A gently decreasing profile with min f' = -0.01 lasts until t = 100.
• "The implicit formula u = f(x - ut) is useless." It is a complete solution up to the breaking time, and it can be evaluated numerically at any point by solving one scalar equation.

Try it: a variable-coefficient problem

Exercise. Solve u_t + t u_x = 0 with u(x, 0) = sin(x).

Solution. The characteristic equation is dX/dt = t, whose solution is X(t) = x_0 + t^2/2. Along it dU/dt = 0, so u keeps the value sin(x_0). Inverting, x_0 = x - t^2/2, so

u(x, t) = sin(x - t^2/2).

Check: u_t = -t cos(x - t^2/2) and t u_x = t cos(x - t^2/2), and they cancel. The profile accelerates to the right but never distorts, because the speed depends on t alone and so is the same for every characteristic at a given instant.

Change one input. Replace the coefficient t by u, giving u_t + u u_x = 0 with the same initial data sin(x). Now the speed varies from point to point, f' = cos(x) attains the minimum -1, and the breaking time is t_b = 1. The same-looking problem now produces a shock at a computable moment.

Where this leaves us

Choosing a curve with dX/dt equal to the coefficient of u_x turns a first-order PDE into an ODE for the solution along that curve. For linear equations the characteristics are determined in advance and never cross, so the method gives a formula outright, whether the coefficient is constant, variable, or accompanied by a decay term. For quasilinear equations the characteristic speed is the solution itself, the lines can converge, and the classical solution breaks down at t_b = -1/min f'. Past that time the correct object is a discontinuity travelling at the Rankine-Hugoniot speed.

Looking ahead

The next lesson takes the second-order hyperbolic equation, u_tt = c^2 u_xx, and factors it into two first-order operators, one carrying signals to the right and one to the left. Running the method of characteristics through that factorisation produces d'Alembert's formula, a closed-form solution of the vibrating string on an infinite line, and with it the notions of domain of dependence and range of influence that make finite propagation speed precise.

Sources

1. Wikipedia contributors. (n.d.). Method of characteristics. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Burgers' equation. Wikipedia. en.wikipedia.org
3. Wikipedia contributors. (n.d.). Rankine-Hugoniot conditions. Wikipedia. en.wikipedia.org
4. Lebl, J. (2026). First order PDEs (Section 1.9). Notes on Diffy Qs: Differential Equations for Engineers. jirka.org
5. Evans, L. C. (2010). Partial differential equations (2nd ed.), Chapter 3. American Mathematical Society.

Take a long rope, lie it out straight, pinch it at the middle and lift it into a narrow triangle. Let go. You do not see one triangle. You see two, each half the height of the original, running away from each other in opposite directions at the same speed.

That is the entire content of d'Alembert's formula in one observation, and he published it in 1747, in an article on vibrating strings that contains the first appearance of the wave equation in print. This lesson derives the equation from Newton's second law, solves it in closed form on the infinite line, and reads off from the formula exactly which part of the initial data can affect a given point at a given time.

Where the equation comes from

Model a string as a curve u(x, t) under tension T with linear mass density rho, and assume displacements and slopes are small. Consider the piece of string between x and x + dx. The tension pulls at each end along the string, and for small slopes the horizontal components cancel while the vertical components are approximately T times the slope. So the net upward force is

T u_x(x + dx, t) - T u_x(x, t), which is approximately T u_xx(x, t) dx.

The mass of the piece is rho dx and its acceleration is u_tt. Newton's second law gives rho dx u_tt = T u_xx dx, and cancelling dx leaves

u_tt = c^2 u_xx with c = sqrt(T/rho).

The constant c has units of speed, and it is a speed: a guitar string tightened raises T and so raises c and so raises the pitch, while a thicker string raises rho and lowers it. Both effects are audible on any stringed instrument.

Factoring the operator

Write the equation as u_tt - c^2 u_xx = 0 and notice that the operator factors exactly as a difference of squares:

u_tt - c^2 u_xx = (d/dt - c d/dx)(d/dt + c d/dx) u,

which you can verify by expanding, using equality of mixed partials to cancel the cross terms. So set

v = u_t + c u_x.

The wave equation says v_t - c v_x = 0, a transport equation moving left at speed c. By Lesson 3 its solution is v(x, t) = h(x + ct) for some function h. Now solve u_t + c u_x = h(x + ct), again by characteristics: along X(t) = x_0 + ct we get dU/dt = h(x_0 + 2ct), and integrating in t produces a function of x + ct plus a constant of integration depending on x_0 = x - ct. The upshot is

u(x, t) = F(x - ct) + G(x + ct),

a right-moving profile plus a left-moving one. This is the same conclusion the canonical-form computation of Lesson 2 reached, arrived at here by a route that shows the mechanism.

Attaching the initial data

Now impose u(x, 0) = phi(x) and u_t(x, 0) = psi(x). Setting t = 0 in the general solution and in its time derivative,

F(x) + G(x) = phi(x) and -c F'(x) + c G'(x) = psi(x).

Divide the second by c and integrate from 0 to x:

G(x) - F(x) = (1/c) times the integral from 0 to x of psi(s) ds, plus a constant K.

Adding and subtracting the two relations,

G(x) = phi(x)/2 + (1/(2c)) times the integral from 0 to x of psi, plus K/2,

F(x) = phi(x)/2 - (1/(2c)) times the integral from 0 to x of psi, minus K/2.

Substitute x - ct into F and x + ct into G and add. The constants K/2 cancel, and the two integrals combine into a single integral over the interval between the arguments. The result is d'Alembert's formula:

u(x, t) = [phi(x - ct) + phi(x + ct)]/2 + (1/(2c)) times the integral from x - ct to x + ct of psi(s) ds.

Check the two conditions. At t = 0 the bracket is 2 phi(x)/2 = phi(x) and the integral is over an empty interval, so u(x, 0) = phi(x). Differentiating in t and using the fundamental theorem of calculus on the integral gives, at t = 0, the value [-c phi' + c phi']/2 + (1/(2c))(c psi(x) + c psi(x)) = psi(x). Both conditions hold, and since the derivation was forced at every step, the solution is unique.

Two solutions, computed

The plucked rope. Take psi = 0, so the string starts at rest with shape phi. Then

u(x, t) = [phi(x - ct) + phi(x + ct)]/2,

which is exactly the observation from the opening: two copies of the initial shape, each of half the height, moving in opposite directions. If phi is a triangle of height 1 supported on [-1, 1], then at t = 2/c the two half-height triangles occupy [-3, -1] and [1, 3] and the middle of the rope is flat and momentarily motionless in shape, though not in velocity.

The struck string. Take phi = 0 and psi(x) = 1 for |x| <= 1, zero otherwise, as if a hammer gave the segment a uniform kick. Then

u(x, t) = (1/(2c)) times the length of the overlap of [x - ct, x + ct] with [-1, 1].

At the origin, once ct >= 1, the overlap is the whole interval of length 2, so u(0, t) = 1/c and stays there forever: a struck string does not return to its original position. That permanent displacement is a real feature of the one-dimensional wave equation and is worth contrasting with the plucked case, where the string does return.

Domain of dependence and range of influence

Read the formula again. The value u(x_0, t_0) depends on phi at exactly two points, x_0 - c t_0 and x_0 + c t_0, and on psi throughout the interval between them. Nothing outside that interval matters at all. The interval

[x_0 - c t_0, x_0 + c t_0]

is the domain of dependence of the point. Reversing the picture, initial data at a point x_0 can affect the solution only within the wedge |x - x_0| <= ct, the range of influence. Signals travel at exactly the speed c and no faster.

Contrast this with the heat equation of Lesson 8, whose solution at any point depends on the initial data everywhere. That is the concrete meaning of hyperbolic versus parabolic, and here you can see it in a formula rather than in a slogan.

Energy is conserved

For data vanishing outside a bounded set, define

E(t) = (1/2) times the integral over all x of (u_t^2 + c^2 u_x^2) dx.

Proposition. E is constant in time.

Proof. Differentiate under the integral sign: E'(t) = the integral of (u_t u_tt + c^2 u_x u_xt) dx. Replace u_tt by c^2 u_xx using the equation, giving c^2 times the integral of (u_t u_xx + u_x u_xt). That integrand is exactly the x-derivative of the product u_t u_x. So the integral equals c^2 times the value of u_t u_x at the two ends, which is zero because the data has bounded support and signals travel at finite speed. Hence E'(t) = 0. This completes the proof.

The two terms are kinetic and potential energy, and the proposition says a vibrating string neither gains nor loses energy. It also gives uniqueness free of charge: if two solutions share the same data, their difference starts with zero energy, so it has zero energy forever, so both u_t and u_x vanish identically and the difference is a constant, which must be zero. Why this matters: an energy identity proves uniqueness without producing any formula, which is why it still works for equations that cannot be solved explicitly.

Common misconceptions

• "The two travelling waves are the two initial conditions." They are not. Both F and G depend on phi and on psi. Only in the special case psi = 0 do the halves reduce to copies of the initial shape.
• "A wave equation smooths its data." It does the opposite of smoothing: it preserves roughness. A corner in phi produces two corners that travel forever along the characteristics x plus or minus ct.
• "The struck string returns to rest position." On the infinite line it does not. After the disturbance passes, the centre of a struck string sits at a permanent displacement 1/c times the total impulse, because the velocity integral never cancels.
• "Finite propagation speed is an approximation." It is exact for the wave equation, and the formula shows it: change the data outside the domain of dependence, however violently, and u(x_0, t_0) does not move at all.

Try it: apply the formula

Exercise. Solve u_tt = 4 u_xx with u(x, 0) = 0 and u_t(x, 0) = x.

Solution. Here c = 2, phi = 0 and psi(s) = s. D'Alembert's formula gives

u(x, t) = (1/4) times the integral from x - 2t to x + 2t of s ds = (1/4)[(x + 2t)^2 - (x - 2t)^2]/2.

Expanding, (x + 2t)^2 - (x - 2t)^2 = 8xt, so u(x, t) = 8xt/8 = xt. Check: u_tt = 0 and u_xx = 0, so the equation holds; u(x, 0) = 0 and u_t(x, 0) = x. Both conditions are satisfied.

Change one input. Keep the same equation but take phi(x) = x and psi = 0. Now u = [(x - 2t) + (x + 2t)]/2 = x, a string that never moves at all. Two problems with the same functions in different slots, and one produces motion while the other does not.

Putting it together

Newton's second law on a string element gives u_tt = c^2 u_xx with c^2 = T/rho. The operator factors into a left-moving and a right-moving transport operator, so the general solution is F(x - ct) + G(x + ct), and matching the two initial conditions produces d'Alembert's formula: half the sum of the displaced initial profiles, plus the average of the initial velocity over the interval between them. That formula displays finite propagation speed as a domain of dependence, and the energy integral is constant, which proves uniqueness without any formula at all.

Looking ahead

A real guitar string is not infinite; it is clamped at both ends, and a wave reaching a clamp reflects. D'Alembert's formula can be adapted by extending the data oddly about each endpoint, but the more powerful approach, and the one that generalises to the heat equation and beyond, is to look for solutions that are products of a function of x and a function of t. That is separation of variables, and it turns a PDE into two ODEs plus an eigenvalue problem.

Sources

1. Weisstein, E. W. (n.d.). d'Alembert's solution. Wolfram MathWorld. mathworld.wolfram.com
2. Wikipedia contributors. (n.d.). D'Alembert's formula. Wikipedia. en.wikipedia.org
3. Lebl, J. (2026). D'Alembert solution of the wave equation (Section 4.8). Notes on Diffy Qs: Differential Equations for Engineers. jirka.org
4. O'Connor, J. J., & Robertson, E. F. (1998). Jean Le Rond d'Alembert. MacTutor History of Mathematics Archive. mathshistory.st-andrews.ac.uk
5. Strauss, W. A. (2008). Partial differential equations: An introduction (2nd ed.), Chapter 2. John Wiley and Sons.

Rest a finger lightly at the exact midpoint of a guitar string and pluck. You do not hear a muffled version of the open note. You hear a clean note an octave higher. Move the finger to a third of the way along and you get the note an octave and a fifth above. The string will vibrate at the fundamental frequency, or at twice it, or three times it, and at no frequency in between.

That discreteness is not a fact about guitars. It comes out of the wave equation the moment you clamp both ends, and it comes out as an eigenvalue problem: a differential equation that has nonzero solutions only for a discrete list of values of a parameter. This lesson derives that list and then builds the general solution out of it.

The problem

Solve

u_tt = c^2 u_xx for 0 < x < L and t > 0,

with boundary conditions u(0, t) = u(L, t) = 0 for all t, and initial conditions u(x, 0) = phi(x), u_t(x, 0) = psi(x). The boundary conditions say the string is held fixed at both ends.

The separation ansatz

Look for solutions of the special form u(x, t) = X(x) T(t). There is no reason a general solution should look like this, and it does not; the point is that such solutions are easy to find and superposition will assemble them into general ones later. Substituting,

X(x) T''(t) = c^2 X''(x) T(t).

Divide by c^2 X T, which is legitimate wherever neither factor vanishes:

T''/(c^2 T) = X''/X.

The left side depends only on t and the right only on x, and they are equal. Fix x and vary t: the right side does not move, so the left side is constant. Fix t and vary x: the same argument on the other side. So both equal a single constant, which we name -lambda with the minus sign chosen for later convenience. That gives two ordinary differential equations:

X'' + lambda X = 0 and T'' + c^2 lambda T = 0.

The boundary conditions transfer to X. Since u(0, t) = X(0)T(t) = 0 for all t, and T is not identically zero for a nontrivial solution, we need X(0) = 0, and likewise X(L) = 0.

The eigenvalue problem, all three cases

We must find every lambda for which X'' + lambda X = 0 with X(0) = X(L) = 0 has a solution other than the zero function. Split by sign.

Case lambda < 0. Write lambda = -mu^2 with mu > 0. The general solution is X = A cosh(mu x) + B sinh(mu x). From X(0) = 0 we get A = 0. From X(L) = B sinh(mu L) = 0, and since sinh is nonzero for positive arguments, B = 0. Only the zero solution.

Case lambda = 0. The equation is X'' = 0, so X = A + Bx. Then X(0) = A = 0 and X(L) = BL = 0 force B = 0. Again only the zero solution.

Case lambda > 0. Write lambda = mu^2. Now X = A cos(mu x) + B sin(mu x). From X(0) = A = 0. From X(L) = B sin(mu L) = 0, either B = 0, which is the zero solution, or sin(mu L) = 0, which happens exactly when mu L = n pi for a positive integer n.

So the eigenvalues and eigenfunctions are

lambda_n = (n pi / L)^2 and X_n(x) = sin(n pi x / L), for n = 1, 2, 3, ....

Key idea: the discreteness of musical harmonics is the statement that this boundary value problem has solutions only for a discrete list of lambda. Nothing about music went into the derivation.

The time equation now reads T'' + (n pi c / L)^2 T = 0, with general solution

T_n(t) = a_n cos(n pi c t / L) + b_n sin(n pi c t / L).

Each product X_n T_n is a standing wave: a fixed spatial shape whose amplitude oscillates. Its frequency is n c / (2L) cycles per unit time, so the modes are the fundamental and its integer multiples. Mode n has n - 1 interior points where sin(n pi x/L) vanishes; those are the nodes, and the finger at the midpoint in the opening paragraph is suppressing every mode that does not have a node there, which is every odd one, leaving the octave.

Assembling the solution

By superposition, any sum of these is a solution, and we take the full series

u(x, t) = the sum over n from 1 to infinity of [a_n cos(n pi c t/L) + b_n sin(n pi c t/L)] sin(n pi x/L).

Setting t = 0 gives phi(x) = the sum of a_n sin(n pi x/L). Differentiating in t and then setting t = 0 gives psi(x) = the sum of b_n (n pi c/L) sin(n pi x/L). So both coefficient families come from expanding a given function in sines. The tool for that is orthogonality: for positive integers m and n,

the integral from 0 to L of sin(m pi x/L) sin(n pi x/L) dx = 0 when m is not n, and = L/2 when m = n.

Multiplying the series for phi by sin(m pi x/L) and integrating kills every term but one, leaving

a_n = (2/L) times the integral from 0 to L of phi(x) sin(n pi x/L) dx,

b_n = (2/(n pi c)) times the integral from 0 to L of psi(x) sin(n pi x/L) dx.

The next two lessons examine whether the resulting series really converges to phi, and what happens if phi has a jump. For now, take the expansion on faith and compute one.

A worked example

Release a string from rest in the parabolic shape phi(x) = x(L - x), with psi = 0. Then every b_n is zero and

a_n = (2/L) times the integral from 0 to L of (Lx - x^2) sin(n pi x/L) dx.

Write a = n pi / L. Two integrations by parts give the integral from 0 to L of x sin(ax) dx = -L(-1)^n / a and the integral from 0 to L of x^2 sin(ax) dx = -L^2(-1)^n/a + 2((-1)^n - 1)/a^3. Combining,

the integral from 0 to L of (Lx - x^2) sin(ax) dx = 2(1 - (-1)^n)/a^3 = 2(1 - (-1)^n) L^3/(n^3 pi^3).

Multiplying by 2/L,

a_n = 4L^2 (1 - (-1)^n) / (n^3 pi^3),

which is 8L^2/(n^3 pi^3) for odd n and zero for even n. So

u(x, t) = the sum over odd n of (8 L^2/(n^3 pi^3)) sin(n pi x/L) cos(n pi c t/L).

Two things are worth reading off. First, the even modes are absent, because the parabola is symmetric about the midpoint and the even sine modes are antisymmetric there; the integral of a symmetric function against an antisymmetric one vanishes. Second, the coefficients decay like 1/n^3, so the series converges quickly: taking three terms already reproduces the shape to within about one percent. Both observations generalise, and the rate of decay of Fourier coefficients turns out to be a precise measure of the smoothness of the function, as the next lesson shows.

Common misconceptions

• "Assuming u = X(x)T(t) loses solutions." It would, if we stopped there. The separated solutions are a basis, and the series recovers everything the boundary conditions permit. The assumption is a device for generating candidates, not a claim about the general solution.
• "The separation constant could be any function." It cannot. Two expressions in different variables that are always equal must both be constant, and the two-step argument of fixing one variable and varying the other proves it.
• "Negative eigenvalues are excluded by convention." They are excluded by computation. With lambda < 0 the solutions are hyperbolic sines and cosines, and no nonzero combination vanishes at two distinct points.
• "Every initial shape excites every harmonic." A shape symmetric about the midpoint excites no even harmonic at all, as the parabola example shows, and a shape that happens to be a single sine excites exactly one mode.

Try it: one mode only

Exercise. Solve the clamped string problem on [0, pi] with c = 1, phi(x) = 3 sin(2x) - sin(5x), and psi = 0.

Solution. With L = pi, the eigenfunctions are sin(nx). The initial shape is already a finite combination of them, so no integrals are needed: a_2 = 3, a_5 = -1, and every other coefficient is zero. Since psi = 0, all b_n vanish. Hence

u(x, t) = 3 sin(2x) cos(2t) - sin(5x) cos(5t).

Verify: each term satisfies u_tt = u_xx because the same integer appears in both factors, both terms vanish at x = 0 and x = pi, at t = 0 the sum is phi, and the time derivative at t = 0 is zero since each factor is a cosine.

Change one input. Take psi(x) = sin(2x) instead, with phi = 0. Now b_2 (2)(1) = 1, so b_2 = 1/2 and u = (1/2) sin(2x) sin(2t). The same mode, now started by a velocity rather than a displacement, and the time factor is a sine instead of a cosine, so the string passes through the flat position rather than starting there.

What to carry forward

Substituting a product into a linear PDE and dividing forces each side to be constant, converting the equation into two ODEs joined by a separation constant. The boundary conditions turn the spatial ODE into an eigenvalue problem, and for a clamped string only positive eigenvalues survive, giving lambda_n = (n pi/L)^2 and eigenfunctions sin(n pi x/L). Superposing the resulting standing waves and matching the initial data determines the coefficients through orthogonality integrals. A parabolic release excites only odd harmonics, with coefficients decaying like 1/n^3.

Looking ahead

Everything above assumed that an arbitrary initial shape can be written as a sum of sines, and that assumption did all the work. The next module examines it. What class of functions have such expansions, in what sense does the series converge, and what does it converge to at a point where the function jumps? The answers are more interesting than a simple yes, and the phenomenon at a jump has a name and a specific numerical size.

Sources

1. Dawkins, P. (n.d.). Separation of variables. Paul's Online Notes, Differential Equations. tutorial.math.lamar.edu
2. Dawkins, P. (n.d.). Vibrating string. Paul's Online Notes, Differential Equations. tutorial.math.lamar.edu
3. Lebl, J. (2026). One-dimensional wave equation (Section 4.7). Notes on Diffy Qs: Differential Equations for Engineers. jirka.org
4. Wikipedia contributors. (n.d.). Separation of variables. Wikipedia. en.wikipedia.org
5. Wikipedia contributors. (n.d.). Standing wave. Wikipedia. en.wikipedia.org

On 21 December 1807, Joseph Fourier read a memoir on the propagation of heat in solid bodies to the Paris Institute. The committee appointed to assess it consisted of Lagrange, Laplace, Monge and Lacroix, and two of them refused to accept its central device. Lagrange and Laplace objected in 1808 to Fourier's claim that a function could be expanded in a series of sines and cosines, and further explanation did not convince them. They were right to be uneasy: the claim is true, but only in a sense that took another sixty years to state correctly.

Here is the phenomenon that troubled them. Add up sin(x) + sin(3x)/3 + sin(5x)/5 + sin(7x)/7 and plot it. Each term is a smooth wave, and the sum is a wobbly approximation to a square wave: flat at height near 0.8, then a sharp drop, then flat at about -0.8. Keep adding odd terms and the flats flatten and the drop steepens. Infinitely many perfectly smooth curves, adding up to something with a corner in it.

The series and its coefficients

Let f be defined on [-L, L]. Its Fourier series is

a_0/2 + the sum over n of [a_n cos(n pi x/L) + b_n sin(n pi x/L)],

where the coefficients are given by the Euler-Fourier formulas

a_n = (1/L) times the integral from -L to L of f(x) cos(n pi x/L) dx, for n = 0, 1, 2, ...,

b_n = (1/L) times the integral from -L to L of f(x) sin(n pi x/L) dx, for n = 1, 2, ....

The a_0/2 rather than a_0 is a convention that lets one formula cover n = 0; the constant term is the average value of f over the interval. Where do these formulas come from? From one computation, done next.

Orthogonality, proved

Proposition. For non-negative integers m and n, and writing p = pi/L:

• The integral from -L to L of sin(mpx) cos(npx) dx is 0 for every m and n.
• The integral of cos(mpx) cos(npx) dx is 0 when m is not n, is L when m = n > 0, and is 2L when m = n = 0.
• The integral of sin(mpx) sin(npx) dx is 0 when m is not n, and L when m = n > 0.

Proof of the second family. Use the product-to-sum identity cos A cos B = [cos(A - B) + cos(A + B)]/2. Then

the integral from -L to L of cos(mpx)cos(npx) dx = (1/2) the integral of cos((m-n)px) dx + (1/2) the integral of cos((m+n)px) dx.

For an integer k not zero, the integral of cos(kpx) from -L to L is [sin(kpx)/(kp)] evaluated at the ends, which is [sin(k pi) - sin(-k pi)]/(kp) = 0. So when m is not n, both m - n and m + n are nonzero (unless both are zero, which cannot happen for distinct non-negative integers), and the whole thing vanishes. When m = n > 0, the first integrand is cos(0) = 1, contributing (1/2)(2L) = L, and the second still vanishes. When m = n = 0, both integrands are 1 and the total is 2L. The other two families are proved the same way with the identities for sin A sin B and sin A cos B. This completes the proof.

Now the coefficient formulas fall out. Suppose f equals its series and that term-by-term integration is legitimate. Multiply both sides by cos(m pi x/L) and integrate from -L to L. Every term on the right vanishes except the one with n = m, which contributes a_m L. Solving gives the stated formula. The same move with sin(m pi x/L) gives b_m. Why this matters: orthogonality is what makes the coefficients computable one at a time, rather than by solving an infinite linear system.

Three series, computed

The square wave. On [-pi, pi] let f(x) = -1 for x < 0 and f(x) = 1 for x > 0. This is odd, so every a_n vanishes: the integral of an odd function times a cosine is odd and integrates to zero over a symmetric interval. For the sines,

b_n = (2/pi) times the integral from 0 to pi of sin(nx) dx = (2/pi) [1 - cos(n pi)]/n = 2[1 - (-1)^n]/(n pi),

which is 4/(n pi) for odd n and zero for even. So

f(x) = (4/pi)[sin(x) + sin(3x)/3 + sin(5x)/5 + ...],

the series from the opening paragraph. Evaluating at x = pi/2, where f = 1 and sin(n pi/2) cycles through 1, 0, -1, 0, gives 1 = (4/pi)(1 - 1/3 + 1/5 - 1/7 + ...), which is Leibniz's series for pi/4.

The sawtooth. On [-pi, pi] let f(x) = x. Again odd, so the cosines vanish. Integrating by parts,

b_n = (2/pi) times the integral from 0 to pi of x sin(nx) dx = (2/pi)(-pi cos(n pi)/n) = 2(-1)^(n+1)/n.

So x = 2[sin(x) - sin(2x)/2 + sin(3x)/3 - ...] on the open interval. Note the coefficients decay only like 1/n, much more slowly than the 1/n^3 of the parabola in the previous lesson. The reason is that the periodic extension of x jumps by 2 pi at every odd multiple of pi, and jumps are expensive.

The parabola. On [-pi, pi] let f(x) = x^2. This is even, so every b_n vanishes. The constant term is

a_0/2 = (1/(2 pi)) times the integral from -pi to pi of x^2 dx = (1/(2 pi))(2 pi^3/3) = pi^2/3.

Two integrations by parts give the integral from 0 to pi of x^2 cos(nx) dx = 2 pi (-1)^n/n^2, so a_n = 4(-1)^n/n^2 and

x^2 = pi^2/3 + 4 times the sum over n of (-1)^n cos(nx)/n^2.

Set x = pi. The left side is pi^2, and cos(n pi) = (-1)^n cancels the other sign, leaving pi^2 = pi^2/3 + 4 times the sum of 1/n^2. Rearranging,

the sum over n from 1 to infinity of 1/n^2 = pi^2/6,

Euler's answer to the Basel problem, obtained here as a by-product of expanding a parabola.

Sine series, cosine series, and boundary conditions

Separation of variables produced sines alone on [0, L], not the full series. The two are connected by extension. Given f on [0, L], extend it to [-L, L] as an odd function; its full Fourier series then contains only sines, and restricting to [0, L] gives the sine series

f(x) = the sum of b_n sin(n pi x/L) with b_n = (2/L) times the integral from 0 to L of f(x) sin(n pi x/L) dx.

Extend evenly instead and you get the cosine series, with only cosines and the analogous coefficient formula. Which one you want is dictated by the boundary conditions: sines vanish at both ends and so match a clamped string or a bar held at zero temperature; cosines have zero derivative at both ends and so match an insulated bar. Lesson 11 makes the connection precise, and Lesson 12 shows both are special cases of one eigenvalue framework.

Energy in the coefficients

One more identity is worth having, both as a check and as a tool. Parseval's theorem states that

(1/L) times the integral from -L to L of f(x)^2 dx = a_0^2/2 + the sum over n of (a_n^2 + b_n^2).

It says the total energy of a signal equals the sum of the energies of its harmonics, which is the natural statement once you view the sines and cosines as orthogonal directions and the coefficients as components. Apply it to the sawtooth: the left side is (1/pi)(2 pi^3/3) = 2 pi^2/3, and the right side is the sum of 4/n^2. Equating, the sum of 1/n^2 is pi^2/6, confirming the value obtained from the parabola by an entirely different route.

Common misconceptions

• "A Fourier series represents f everywhere it is defined." It represents the periodic extension of f. The sawtooth series equals x on the open interval and equals zero at x = pi, where the extension jumps. The next lesson makes this precise.
• "Odd functions have odd-numbered coefficients." Odd functions have only sine terms, and even functions only cosine terms. Whether the surviving indices are odd or even is a separate matter of symmetry about the midpoint, as the square wave and the parabola show.
• "Slow decay of coefficients means the function is badly behaved everywhere." The decay rate reports on the smoothness of the periodic extension. The function x is perfectly smooth on the open interval; it is the jump created at the ends that forces the 1/n decay.
• "Term-by-term integration is obviously allowed." It is a theorem, not an observation, and it was the point Lagrange and Laplace balked at. Uniform convergence, or square integrability plus Parseval, is what licenses the manipulation.

Try it: a cosine series

Exercise. Find the cosine series of f(x) = x on [0, pi].

Solution. Here the even extension is used, so a_n = (2/pi) times the integral from 0 to pi of x cos(nx) dx. The constant term is a_0/2 = (1/pi) times the integral from 0 to pi of x dx = pi/2. For n at least 1, integrating by parts,

the integral from 0 to pi of x cos(nx) dx = [x sin(nx)/n] from 0 to pi - (1/n) times the integral of sin(nx) dx = 0 + [cos(nx)/n^2] from 0 to pi = ((-1)^n - 1)/n^2.

So a_n = 2((-1)^n - 1)/(pi n^2), which is -4/(pi n^2) for odd n and zero for even. Hence

x = pi/2 - (4/pi)[cos(x) + cos(3x)/9 + cos(5x)/25 + ...] on [0, pi].

Change one input. Compare with the sine series of the same function, computed above, whose coefficients decay like 1/n. The cosine series decays like 1/n^2 and therefore converges faster, because the even extension of x is continuous everywhere while the odd extension jumps at pi. Same function, two expansions, and the choice of extension decides the convergence rate. Setting x = 0 in the cosine series also gives the bonus identity 1 + 1/9 + 1/25 + ... = pi^2/8.

Summing up

The sines and cosines on a symmetric interval are mutually orthogonal, proved with product-to-sum identities, and that orthogonality converts the problem of finding coefficients into a single integral apiece. The square wave gives Leibniz's series, the sawtooth gives coefficients decaying like 1/n, and the parabola gives the Basel sum pi^2/6, confirmed independently by Parseval applied to the sawtooth. Sine series come from odd extensions and cosine series from even ones, and which you use is decided by the boundary conditions.

Looking ahead

Every computation above assumed the series converges to the function. At x = pi the sawtooth series is exactly zero while the function is pi, so that assumption is false somewhere. The next lesson states Dirichlet's theorem, which says precisely what the series converges to at every point including the jumps, and examines the overshoot that partial sums exhibit near a discontinuity, a nine percent excess that does not go away no matter how many terms you take.

Sources

1. Weisstein, E. W. (n.d.). Fourier series. Wolfram MathWorld. mathworld.wolfram.com
2. Weisstein, E. W. (n.d.). Fourier series: square wave. Wolfram MathWorld. mathworld.wolfram.com
3. Dawkins, P. (n.d.). Fourier series. Paul's Online Notes, Differential Equations. tutorial.math.lamar.edu
4. O'Connor, J. J., & Robertson, E. F. (1997). Jean Baptiste Joseph Fourier. MacTutor History of Mathematics Archive. mathshistory.st-andrews.ac.uk
5. Lebl, J. (2026). The trigonometric series (Section 4.2). Notes on Diffy Qs: Differential Equations for Engineers. jirka.org

Henry Wilbraham noticed it first, in an 1848 paper that nobody read for sixty-six years. In 1898 Albert Michelson built a machine that computed and resynthesised Fourier series, and a correspondence in Nature between Michelson and A. E. H. Love about the square wave prompted J. Willard Gibbs to write a note that year. His first note got the limit wrong. In April 1899 he published a correction describing the overshoot correctly, and in 1906 Maxime Bocher gave the full analysis and named the effect after Gibbs.

The effect is this. Take the partial sums of the square wave series from the previous lesson. Near the jump they do not merely approximate the step; they overshoot it, by about nine percent of the jump height, and the overshoot does not shrink as you add terms. It only moves closer to the jump. Reconciling that with the claim that the series converges is the business of this lesson.

Three different questions

"Does the series converge to f" is three questions wearing one coat. Let S_N denote the N-th partial sum.

• Pointwise convergence: for each fixed x, S_N(x) approaches a limit as N grows. Different points may converge at wildly different rates.
• Uniform convergence: the largest error, the maximum over all x of |S_N(x) - f(x)|, approaches zero. This is much stronger and is what licenses term-by-term integration and differentiation.
• Mean-square convergence: the integral of (S_N - f)^2 approaches zero. This ignores what happens on small sets entirely.

The Gibbs phenomenon is precisely a case where the first and third hold and the second fails. Keeping the three apart is the whole resolution.

Dirichlet's theorem

Call f piecewise smooth on [-L, L] if the interval splits into finitely many pieces on each of which f and f' are continuous, with finite one-sided limits at the break points.

Theorem (Dirichlet). If f is piecewise smooth, its Fourier series converges at every point x, and the value it converges to is

[f(x+) + f(x-)]/2,

the average of the one-sided limits. In particular it converges to f(x) wherever f is continuous.

The theorem is doing something the naive statement cannot. A Fourier series is built from the values of f under an integral sign, and changing f at one point changes no integral, so the series cannot possibly know the value assigned at a jump. What it does know is both one-sided limits, and it splits the difference.

Worked check. The sawtooth series 2[sin x - sin(2x)/2 + sin(3x)/3 - ...] at x = pi: every term is sin(n pi) = 0, so the series sums to 0. Dirichlet predicts [f(pi-) + f(pi+)]/2, where the periodic extension approaches pi from the left and -pi from the right, giving (pi + (-pi))/2 = 0. Agreement, and the function value f(pi) = pi is simply not what the series reports.

Worked check. The square wave series at x = 0: every term is sin(0) = 0, so the sum is 0, and Dirichlet predicts (-1 + 1)/2 = 0. Again the midpoint.

When convergence is uniform

Uniform convergence needs continuity of the periodic extension, and a little more.

Theorem. If f is continuous on [-L, L], piecewise smooth, and satisfies f(-L) = f(L), then the Fourier series converges to f uniformly and absolutely on the whole interval.

The condition f(-L) = f(L) is not cosmetic. It is exactly the requirement that the periodic extension have no jump at the seam. The function f(x) = x on [-pi, pi] is continuous and smooth and fails only this condition, and its series is famously non-uniform near the endpoints.

There is a quantitative version of the same idea, and it is the most useful rule of thumb in the subject. If the periodic extension of f is continuous with k continuous derivatives, the coefficients decay faster than 1/n^k. Read backwards: measured decay tells you smoothness.

Remember: a jump costs you a factor of n, a corner costs you another, and each further degree of smoothness buys one more.

Mean-square convergence always works

The third mode of convergence is the most forgiving and the most useful in theory. If f is merely square integrable, meaning the integral of f^2 is finite, then the partial sums converge to f in mean square, with no continuity assumption whatever. Two facts are attached. Bessel's inequality says that for any N, the sum of the first N squared coefficients is at most the mean square of f; and equality in the limit is Parseval's theorem from the previous lesson, which is equivalent to saying that the trigonometric system is complete, so nothing is left over.

This is the sense in which Lagrange's and Laplace's objection is finally answered. Fourier's claim is not that the series converges to the function at every point, which is false; it is that the series converges to the function in mean square, which is true for every function of finite energy.

The Gibbs phenomenon, quantified

Return to the square wave, jumping from -1 to 1, so the jump has size 2. Compute the partial sum S_N and find its first maximum to the right of the jump. As N grows, that maximum does not approach 1. It approaches

(2/pi) times the integral from 0 to pi of sin(t)/t dt = 1.17897974....

The excess above 1 is about 0.17898, which is 0.0894898... of the jump size 2, or roughly 8.95 percent. That number, the Wilbraham-Gibbs constant, is universal: every jump discontinuity of every piecewise smooth function produces an overshoot of the same fraction of its own jump.

Where does the location go? The first maximum of S_N sits at a distance of order pi/N from the jump. So as N increases, the spike gets no shorter but it gets narrower and slides toward the discontinuity.

Now the reconciliation. Fix any point x_0 away from the jump. Once N is large enough that pi/N is smaller than the distance from x_0 to the jump, the spike has already passed x_0, and S_N(x_0) settles down to f(x_0). So pointwise convergence holds at every point, exactly as Dirichlet promises. What fails is uniform convergence: for every N there is some point, admittedly a moving one, where the error is about 0.09 times the jump. And mean-square convergence holds because the area under the spike shrinks like 1/N, so its contribution to the integral of the squared error vanishes.

Key idea: the graph of the limit and the limit of the graphs are different objects. That was exactly the distinction Gibbs drew in his 1899 correction.

What this costs in a PDE

Suppose you solve the heat equation on [0, pi] with a square-wave initial temperature. The Fourier coefficients start with a 1/n decay and are then multiplied by exp(-k n^2 t), which for any t > 0 crushes the high harmonics. The Gibbs wiggles are gone instantly, and the solution is smooth for all positive time. This is the smoothing property of Lesson 8, visible in the coefficients.

Now solve the wave equation with the same initial shape. The coefficients are multiplied by cos(n c t), which has modulus at most 1 and never decays. The wiggles are still there at every later time, moving along with the wave. That difference between multiplying by a decaying exponential and by a bounded oscillation is the difference between parabolic and hyperbolic, expressed in the language of Fourier coefficients.

Common misconceptions

• "The Gibbs overshoot shrinks if you take enough terms." Its height is asymptotically constant at about 8.95 percent of the jump. Only its width shrinks. Averaging the partial sums, a device called Cesaro summation, does remove it, at the cost of slower convergence elsewhere.
• "Gibbs contradicts the convergence theorem." It contradicts uniform convergence, which was never claimed for a discontinuous function. Pointwise convergence at every fixed point is untouched.
• "The series gets the value at a jump wrong." It returns the average of the one-sided limits, which is the only answer determined by the integrals defining the coefficients. Changing f at a single point changes no coefficient at all.
• "Michelson's machine discovered the phenomenon." The popular anecdote overstates it. The machine's graphs were not sharp enough to display the effect clearly, Michelson never mentioned it, and the correct analysis came from Gibbs in 1899 and Bocher in 1906, half a century after Wilbraham.

Try it: predict the limit at every point

Exercise. Let f(x) = 0 for -pi < x < 0 and f(x) = x for 0 <= x < pi. State what the Fourier series converges to at x = -pi/2, at x = 0, and at x = pi.

Solution. At x = -pi/2 the function is continuous with value 0, so the series converges to 0. At x = 0 the one-sided limits are f(0-) = 0 and f(0+) = 0, so the function is continuous there and the series converges to 0, with no jump at all despite the change of formula. At x = pi the periodic extension approaches pi from the left and f(-pi+) = 0 from the right, so the series converges to (pi + 0)/2 = pi/2, which is neither one-sided value.

Change one input. Redefine f to be 1 rather than 0 on the left half. Now x = 0 is a genuine jump from 1 to 0, the series converges there to 1/2, and partial sums overshoot to about 1 + 0.0895 just left of the origin and undershoot symmetrically to the right. One changed constant, and a smooth crossing becomes a Gibbs jump.

What to remember

Convergence comes in three strengths, and a discontinuous function separates them. Dirichlet's theorem gives pointwise convergence for piecewise smooth functions, to the average of the one-sided limits. Uniform convergence requires the periodic extension to be continuous, and the decay rate of the coefficients measures how smooth that extension is: 1/n for a jump, 1/n^2 for a corner, faster for smoother. Mean-square convergence needs nothing but finite energy. Near a jump, partial sums overshoot by about 8.95 percent of the jump forever, in a spike that narrows without shortening.

Looking ahead

The observation at the end, that heat multiplies coefficients by a decaying exponential while waves multiply by a bounded oscillation, is about to become the main event. The next module takes the heat equation seriously: it derives the solution on a bounded interval, extracts the smoothing property from the exponential factor, and then shows that trying to run the same formula backwards in time amplifies every high harmonic without limit.

Sources

1. Wikipedia contributors. (n.d.). Gibbs phenomenon. Wikipedia. en.wikipedia.org
2. Weisstein, E. W. (n.d.). Gibbs phenomenon. Wolfram MathWorld. mathworld.wolfram.com
3. Wikipedia contributors. (n.d.). Convergence of Fourier series. Wikipedia. en.wikipedia.org
4. Dawkins, P. (n.d.). Convergence of Fourier series. Paul's Online Notes, Differential Equations. tutorial.math.lamar.edu
5. Wikipedia contributors. (n.d.). Parseval's theorem. Wikipedia. en.wikipedia.org

Film a drop of ink spreading in a glass of still water, then run the film backwards. Everyone watching knows immediately which way is forward. Nobody has ever seen a uniformly grey glass spontaneously concentrate its ink into a sphere, and this is not merely improbable: the equation that governs the spreading is mathematically incapable of running the other way, in a sense this lesson makes exact.

The same equation is the one derived in Lesson 1 from conservation plus Fourier's law, u_t = k u_xx. We will solve it on a bar, read off why it smooths whatever you give it, and then watch the backward problem fail.

The bar with cold ends

Solve u_t = k u_xx on 0 < x < L, with u(0, t) = u(L, t) = 0 and u(x, 0) = phi(x). Physically: a bar with initial temperature profile phi, whose two ends are suddenly clamped to an ice bath.

Separating variables as in Lesson 5, put u = X(x)T(t). Substituting and dividing by kXT,

T'/(kT) = X''/X = -lambda.

The spatial problem X'' + lambda X = 0 with X(0) = X(L) = 0 is exactly the eigenvalue problem already solved: lambda_n = (n pi/L)^2 and X_n = sin(n pi x/L). What is new is the time equation, which is now first order:

T' = -k lambda_n T, so T_n(t) = exp(-k (n pi/L)^2 t).

Instead of the oscillation cos(n pi c t/L) of the wave equation, we get a decaying exponential, and the decay rate grows like n^2. That single difference is responsible for everything else in this lesson. Superposing,

u(x, t) = the sum over n of b_n exp(-k n^2 pi^2 t/L^2) sin(n pi x/L),

with b_n = (2/L) times the integral from 0 to L of phi(x) sin(n pi x/L) dx, the sine coefficients of the initial profile.

A worked case, with numbers

Take L = pi, k = 1, and a bar initially at the uniform temperature phi(x) = 100, with both ends dropped to zero at t = 0. From Lesson 6, the sine coefficients of a constant are

b_n = (2/pi) times the integral from 0 to pi of 100 sin(nx) dx = 200[1 - (-1)^n]/(n pi),

which is 400/(n pi) for odd n and zero for even. So

u(x, t) = the sum over odd n of (400/(n pi)) exp(-n^2 t) sin(nx).

Now watch the exponentials. At t = 1 the first mode carries the factor e^(-1) = 0.368, the third mode e^(-9) = 1.2 times 10^(-4), and the fifth e^(-25) = 1.4 times 10^(-11). Relative to the first mode, the third contributes about three parts in ten thousand and the fifth about four parts in a hundred billion. After one time unit the temperature profile is a pure half-sine arch to within a fraction of a percent, whatever the initial profile happened to be, and thereafter it just decays.

That is the long-time behaviour of every solution of this problem: the slowest mode wins, and it decays at rate k pi^2/L^2. Notice the L^2: doubling the length of a bar quadruples the time it takes to cool. Diffusion times scale with the square of the distance, which is why a thick steak takes disproportionately longer to cook than a thin one.

Smoothing, proved

Proposition. For any bounded initial profile phi and any t > 0, the solution u(., t) is infinitely differentiable in x.

Proof. The coefficients b_n are bounded, since |b_n| is at most (2/L) times the integral of |phi|, a fixed number M. Differentiating the series m times in x multiplies the n-th term by (n pi/L)^m up to sign, so the n-th term of the differentiated series is bounded in absolute value by

M (n pi/L)^m exp(-k n^2 pi^2 t/L^2).

For fixed t > 0 this decays faster than any power of 1/n, because an exponential in n^2 beats any polynomial in n. So the differentiated series converges absolutely and uniformly for every m, which justifies differentiating term by term as often as we like. This completes the proof.

Read what that says. Give the bar a sawtooth initial temperature with corners, or a square-wave profile with jumps, and one instant later the temperature is a perfectly smooth function of position. No hyperbolic equation does this; the wave equation multiplies its coefficients by cos(n c t), which never decays and so preserves every corner forever. Why this matters: smoothing is not a side effect of diffusion, it is diffusion, and it is visible in a single exponential factor.

Uniqueness, by an energy argument

Proposition. The problem above has at most one solution.

Proof. Let u and v be solutions with the same data and put w = u - v, which solves the equation with zero initial and boundary data. Define E(t) = the integral from 0 to L of w(x, t)^2 dx, which is non-negative. Then

E'(t) = 2 times the integral of w w_t dx = 2k times the integral of w w_xx dx.

Integrate by parts: the boundary term 2k[w w_x] vanishes because w is zero at both ends, leaving E'(t) = -2k times the integral of (w_x)^2 dx, which is at most zero. So E is non-increasing, and E(0) = 0 because w starts at zero. A non-negative, non-increasing function starting at zero is identically zero, so w vanishes and u = v. This completes the proof.

Notice what the proof did not need: any formula for the solution. Energy arguments of exactly this kind are the standard tool for uniqueness in problems where no series is available.

The backward problem fails

Now reverse the question. You are handed the temperature profile of the bar at time T = 1 and asked for the profile at t = 0. Formally the answer is easy: expand the given profile in sines with coefficients c_n, and since forward evolution multiplies the n-th coefficient by e^(-n^2), the initial coefficients must have been b_n = c_n e^(n^2). Every step is algebraically legitimate. The problem is what it does to errors.

Take two candidate final profiles: g(x) = 0, and the perturbed g_n(x) = (1/n) sin(nx). As n grows, g_n approaches g uniformly, since its maximum height is 1/n. Their preimages at t = 0, however, are 0 and

(1/n) e^(n^2) sin(nx),

whose maximum height is e^(n^2)/n. At n = 5 that is already about 1.4 times 10^(10); at n = 10 it is beyond any physical scale. So two final states differing by less than a thousandth of a degree correspond to initial states differing by astronomical amounts. Arbitrarily small changes in the data produce arbitrarily large changes in the answer, which is exactly the failure of continuous dependence that Lesson 16 calls ill-posedness.

The physical reading is the ink drop. Diffusion multiplies fine spatial detail by a factor that shrinks like e^(-k n^2 t), so within moments the fine structure of the initial state is not merely small but has fallen below the noise floor of any measurement. The information is gone, and no amount of arithmetic recovers it. The upshot: the heat equation is not reversible not because reversing it is hard, but because reversing it is discontinuous.

Two footnotes keep this honest. First, the backward problem is not meaningless: it is the mathematical form of deblurring an image, and it is solved in practice by regularisation, which trades exactness for stability. Second, molecular dynamics is time reversible; irreversibility appears when you pass to the diffusion description, and the heat equation inherits it.

Insulated ends, and what is conserved

Replace the boundary conditions by u_x(0, t) = u_x(L, t) = 0, meaning no heat flux through the ends. Then

d/dt of the integral from 0 to L of u dx = k times the integral of u_xx dx = k[u_x] evaluated at the ends = 0.

Total heat is exactly conserved, and the solution approaches the constant equal to the average of the initial profile rather than to zero. The eigenfunctions become cosines, which is the cosine series of Lesson 6 arriving for a physical reason. Lesson 11 catalogues the boundary conditions systematically.

Common misconceptions

• "The heat equation is reversible because the algebra can be run backwards." The algebra can. The map from final data to initial data is unbounded, so the reversed problem has no stable solution, and that is a property of the equation and not a defect of the method.
• "Smoothing takes a while." It is instantaneous. For every t > 0, however small, the solution is infinitely differentiable, because the exponential factor is already crushing high harmonics.
• "The bar cools uniformly." It cools mode by mode, and the high modes go first. Whatever the initial shape, after a short time the profile is essentially the fundamental half-sine.
• "Doubling the length of a bar doubles the cooling time." It quadruples it, since the slowest decay rate is k pi^2/L^2. Diffusion time scales as the square of the distance.

Try it: a single-mode start

Exercise. Solve u_t = 3 u_xx on [0, pi] with u(0, t) = u(pi, t) = 0 and u(x, 0) = 5 sin(2x) - 2 sin(7x). Then find the time at which the second term contributes less than one percent of the first.

Solution. The initial profile is already a combination of eigenfunctions, so no integrals are needed: b_2 = 5, b_7 = -2, all others zero. With k = 3 and L = pi the decay factors are exp(-3 n^2 t), so

u(x, t) = 5 e^(-12t) sin(2x) - 2 e^(-147t) sin(7x).

The ratio of amplitudes is (2/5) e^(-135t). Setting that equal to 0.01 gives e^(-135t) = 0.025, so t = ln(40)/135, about 0.0273. After roughly three hundredths of a time unit the seventh mode is negligible and the bar is essentially oscillating in its second mode alone as it cools.

Change one input. Solve the wave equation u_tt = 3 u_xx with the same data and zero initial velocity. Now u = 5 cos(2 sqrt(3) t) sin(2x) - 2 cos(7 sqrt(3) t) sin(7x), and the ratio of amplitudes stays at 2/5 forever. The seventh mode never becomes negligible, because nothing damps it.

Where this leaves us

Separation of variables gives the same eigenfunctions as for the wave equation but a first-order time equation, so each coefficient is multiplied by exp(-k n^2 pi^2 t/L^2). That factor makes the solution infinitely smooth for every positive time, makes the slowest mode dominate quickly, makes cooling time scale like the square of the length, and makes the backward problem violate continuous dependence by amplifying high harmonics without bound. An energy integral proves uniqueness forward in time without using the series at all.

Looking ahead

Let a heat problem run until nothing changes and the time derivative drops out, leaving u_xx = 0 in one dimension, or u_xx + u_yy = 0 in two. That is Laplace's equation, the elliptic member of the trio. Its solutions turn out to satisfy a remarkable averaging property, and from that one property come the maximum principle, uniqueness for the Dirichlet problem, and the smoothness that makes elliptic equations the best behaved of the three.

Sources

1. Wikipedia contributors. (n.d.). Heat equation. Wikipedia. en.wikipedia.org
2. Dawkins, P. (n.d.). Solving the heat equation. Paul's Online Notes, Differential Equations. tutorial.math.lamar.edu
3. Lebl, J. (2026). The heat equation (Section 4.6). Notes on Diffy Qs: Differential Equations for Engineers. jirka.org
4. Weisstein, E. W. (n.d.). Heat conduction equation. Wolfram MathWorld. mathworld.wolfram.com
5. Evans, L. C. (2010). Partial differential equations (2nd ed.), Chapter 2, Section 3. American Mathematical Society.

Bend a loop of wire into any shape you like, dip it in soapy water, and pull it out. The film that spans the loop is stretched taut, and for small deflections its height satisfies Laplace's equation. Look at the film: there is no bump anywhere in the middle. Every point of the film sits below the highest point of the wire and above the lowest. You cannot make an interior peak, no matter how you bend the wire.

That observation is the maximum principle, and it is not a soap-film fact. It follows from the equation, it holds in every dimension, and it delivers uniqueness and stability for the Dirichlet problem in two lines each.

Laplace's equation and harmonic functions

A function u on a region of the plane is harmonic if

u_xx + u_yy = 0,

Laplace's equation, written compactly as Laplacian(u) = 0. With a right-hand side it becomes Poisson's equation Laplacian(u) = f, which describes an electrostatic potential with charge density proportional to f, or a soap film pushed by a distributed load.

Where does the equation come from? Set the time derivative in the heat equation to zero. A temperature distribution that has stopped changing satisfies u_xx + u_yy = 0, so harmonic functions are exactly the steady states of diffusion. That interpretation will keep paying off.

Some harmonic functions, worth memorising:

The pattern behind the first four is that the real and imaginary parts of any holomorphic function of a complex variable are harmonic, which is why complex analysis and two-dimensional potential theory keep meeting.

The mean value property, proved

Theorem. Let u be harmonic on a region containing the closed disc of radius R about a point P. Then u(P) equals the average of u over the circle of any radius r at most R centred at P.

Proof. Define

phi(r) = (1/(2 pi)) times the integral from 0 to 2 pi of u(P + r w(theta)) d theta,

where w(theta) is the unit vector at angle theta. This is the average of u over the circle of radius r. Differentiate under the integral sign:

phi'(r) = (1/(2 pi)) times the integral of (gradient of u) dotted with w d theta.

The dot product of the gradient with the outward unit normal is the normal derivative, and rewriting the angular integral as an integral over the circle of circumference 2 pi r gives

phi'(r) = (1/(2 pi r)) times the integral over the circle of the normal derivative of u.

By the divergence theorem, that integral over the circle equals the integral of Laplacian(u) over the enclosed disc, which is zero because u is harmonic. So phi'(r) = 0 for every r in (0, R], meaning phi is constant. Finally, as r approaches 0, continuity of u gives phi(r) approaching u(P). Hence phi(r) = u(P) for every r. This completes the proof.

Averaging the circle averages over r from 0 to R gives the same statement for the solid disc: a harmonic function equals its average over any disc centred at the point. Key idea: a harmonic function is precisely one that is always equal to its own local average, which is why it can have no local bumps.

The maximum principle

Theorem (strong maximum principle). Let u be harmonic on a connected open region D. If u attains a maximum at an interior point of D, then u is constant on D.

Proof. Suppose u(P) = M is the maximum value and P lies inside D. Choose a disc about P contained in D. By the mean value property, M = u(P) is the average of u over that disc. But u is at most M everywhere, and if u were strictly less than M anywhere on the disc, then by continuity it would be less than M on a small patch of positive area, and the average would come out strictly below M. That contradiction forces u = M throughout the disc. So the set where u = M is open. It is also closed in D, being the preimage of a point under a continuous function, and it is nonempty. In a connected region the only such set is the whole region, so u = M on D. This completes the proof.

Corollary (weak maximum principle). A function harmonic inside a bounded region and continuous up to the boundary attains its maximum and its minimum on the boundary.

The minimum statement comes free by applying the maximum principle to -u, which is harmonic whenever u is. This is the soap film: the highest and lowest points of the film are on the wire.

Two theorems for free

The Dirichlet problem asks for a function harmonic inside a region D and equal to a prescribed continuous function g on the boundary.

Theorem (uniqueness). The Dirichlet problem has at most one solution.

Proof. If u_1 and u_2 both solve it, then w = u_1 - u_2 is harmonic in D and zero on the boundary. By the weak maximum principle its maximum over the closed region is attained on the boundary and so equals 0, and its minimum likewise equals 0. Hence w is identically zero. This completes the proof.

Theorem (continuous dependence). If u_1 and u_2 solve the Dirichlet problem with boundary data g_1 and g_2, then the maximum of |u_1 - u_2| over the region is at most the maximum of |g_1 - g_2| over the boundary.

Proof. The difference is harmonic, so its maximum and minimum occur on the boundary, where it equals g_1 - g_2. This completes the proof.

Compare with the backward heat problem of the previous lesson, where an arbitrarily small change in the data produced an arbitrarily large change in the answer. Here a change of one millidegree in the boundary temperature changes the interior temperature by at most one millidegree, everywhere, with no constant and no loss. Why this matters: the Dirichlet problem is the model of a well-posed problem, and Lesson 16 uses that word precisely.

Two further facts, stated without proof but worth knowing. Harmonic functions are infinitely differentiable, indeed real analytic, in the interior, however rough the boundary data is; that is the elliptic version of the smoothing seen for the heat equation. And Liouville's theorem says a harmonic function bounded on the whole plane must be constant, which rules out any nontrivial global solution without boundaries.

Laplace's equation on a rectangle, solved

Solve u_xx + u_yy = 0 on the rectangle 0 < x < a, 0 < y < b, with u = 0 on the left, right and bottom sides, and u(x, b) = f(x) on the top.

Separating variables with u = X(x)Y(y) gives X''/X = -Y''/Y = -lambda. The conditions u(0, y) = u(a, y) = 0 force X(0) = X(a) = 0, which is the familiar eigenvalue problem: lambda_n = (n pi/a)^2 and X_n = sin(n pi x/a). The y equation is then Y'' = lambda_n Y, with real exponential solutions, so

Y_n(y) = C cosh(n pi y/a) + D sinh(n pi y/a).

The condition u(x, 0) = 0 forces Y_n(0) = 0, killing the cosh term. So

u(x, y) = the sum over n of A_n sin(n pi x/a) sinh(n pi y/a).

Matching the top edge, f(x) = the sum of A_n sinh(n pi b/a) sin(n pi x/a), so the products A_n sinh(n pi b/a) are the sine coefficients of f, giving

A_n = (2/(a sinh(n pi b/a))) times the integral from 0 to a of f(x) sin(n pi x/a) dx.

Worked instance. Take f(x) = sin(pi x/a). Then only the n = 1 coefficient survives, with A_1 sinh(pi b/a) = 1, so

u(x, y) = sin(pi x/a) sinh(pi y/a) / sinh(pi b/a).

Check: the second x-derivative brings down -(pi/a)^2 and the second y-derivative brings down +(pi/a)^2, and they cancel. On the bottom, sinh(0) = 0. On the sides, sin(0) = sin(pi) = 0. On the top the sinh factors cancel and the value is f. All four conditions hold, and the maximum principle is visible: the largest value in the rectangle is 1, attained on the top edge, and nothing exceeds it inside.

Common misconceptions

• "The maximum principle says the solution is bounded." It says something sharper: the extreme values occur on the boundary, so the interior can be read off as being trapped between boundary values. Boundedness is a consequence, not the content.
• "Harmonic functions have no critical points." They can, and x^2 - y^2 has one at the origin. What they cannot have is an interior local maximum or minimum; the origin here is a saddle.
• "Laplace's equation needs initial conditions." There is no time variable. It needs data on the whole boundary of the region, and prescribing data on only part of it generally destroys uniqueness.
• "Rough boundary data gives a rough solution." The solution is real analytic everywhere in the interior, no matter how jagged the boundary values are; the roughness is confined to the boundary itself.

Try it: read the maximum principle off a formula

Exercise. The function u(x, y) = x^2 - y^2 is harmonic on the square -1 <= x, y <= 1. Find its maximum and minimum on the square and check that both occur on the boundary.

Solution. On the square, x^2 is at most 1 and y^2 is at most 1, so u is at most 1 and at least -1. The value 1 is attained where x^2 = 1 and y = 0, that is at (1, 0) and (-1, 0), both on the boundary. The value -1 is attained at (0, 1) and (0, -1), also on the boundary. At the interior critical point (0, 0) the value is 0, which is neither the maximum nor the minimum: it is a saddle.

Change one input. Consider v(x, y) = x^2 + y^2, which is not harmonic since Laplacian(v) = 4. Its minimum is 0, attained at the interior point (0, 0). The maximum principle fails immediately once the Laplacian is nonzero, and the sign of that Laplacian tells you which half of the principle survives: functions with non-negative Laplacian, called subharmonic, still attain their maxima on the boundary.

Putting it together

Harmonic functions are steady states of diffusion, and the divergence theorem shows each one equals its average over any circle centred at a point. That mean value property forbids interior maxima unless the function is constant, so extremes live on the boundary. From the maximum principle, uniqueness and continuous dependence for the Dirichlet problem follow in two lines apiece, with a sharp constant of 1. On a rectangle, separation of variables gives sines in the direction with two zero edges and hyperbolic sines in the other.

Looking ahead

Rectangles are convenient because they match Cartesian coordinates. A disc does not, and forcing sines and cosines onto it produces nothing. The next lesson repeats the whole construction in polar coordinates, and the answer collapses into a single integral formula, the Poisson kernel, which writes the value at any interior point as a weighted average of the boundary data. That formula makes the mean value property a special case and proves existence for the Dirichlet problem on a disc outright.

Sources

1. Wikipedia contributors. (n.d.). Harmonic function. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Maximum principle. Wikipedia. en.wikipedia.org
3. Weisstein, E. W. (n.d.). Laplace's equation. Wolfram MathWorld. mathworld.wolfram.com
4. Dawkins, P. (n.d.). Laplace's equation. Paul's Online Notes, Differential Equations. tutorial.math.lamar.edu
5. Lebl, J. (2026). Steady state temperature and the Laplacian (Section 4.9). Notes on Diffy Qs: Differential Equations for Engineers. jirka.org

Take a thin copper disc of radius 1. Hold the right half of its rim at 100 degrees and the left half at 0, and wait until the temperature stops changing. What is the temperature at the exact centre?

Fifty degrees, and you already know why: the mean value property from the previous lesson says the centre value is the average of the boundary values, and the average of a function that is 100 on half the circle and 0 on the other half is 50. Now move halfway out toward the hot side, to the point at radius 1/2 on the horizontal axis. The mean value property says nothing about that point. This lesson builds the formula that does, and the answer turns out to be about 79.5 degrees.

Laplace's equation in polar coordinates

Writing x = r cos(theta) and y = r sin(theta) and grinding through the chain rule converts u_xx + u_yy = 0 into

u_rr + u_r/r + u_(theta theta)/r^2 = 0.

The extra u_r/r term is what a circle costs you, and it is the reason the radial equation below is not a constant-coefficient one. Multiplying through by r^2 gives the form we will separate:

r^2 u_rr + r u_r + u_(theta theta) = 0.

Separating in polar coordinates

Put u = R(r) Theta(theta). Substituting and dividing by R Theta,

(r^2 R'' + r R')/R = -Theta''/Theta = lambda.

The angular problem is Theta'' + lambda Theta = 0, but the boundary condition is unusual: there is no boundary in the angle, only the requirement that Theta be 2 pi-periodic, since theta and theta + 2 pi name the same point. Negative lambda gives exponentials, which are never periodic; lambda = 0 gives A + B theta, periodic only if B = 0; positive lambda = mu^2 gives cos(mu theta) and sin(mu theta), periodic exactly when mu is an integer. So

lambda_n = n^2 for n = 0, 1, 2, ..., with eigenfunctions 1, cos(n theta), sin(n theta).

The radial equation is r^2 R'' + r R' - n^2 R = 0, an Euler equation. Trying R = r^p gives p(p-1) + p - n^2 = p^2 - n^2 = 0, so p = n or p = -n. For n = 0 the two roots coincide and the second solution is ln(r). Now the physical condition enters: the temperature must be finite at the centre, and both r^(-n) and ln(r) blow up there. So only r^n survives.

Key idea: nothing in the differential equation excluded the singular solutions. Boundedness at the origin did, and it is a genuine side condition even though it is not written on a boundary.

The series solution

Superposing, for a disc of radius a with boundary data g(theta),

u(r, theta) = a_0/2 + the sum over n of (r/a)^n [a_n cos(n theta) + b_n sin(n theta)],

where a_n and b_n are the ordinary Fourier coefficients of g on [-pi, pi]. Setting r = a makes the factor 1 and recovers the Fourier series of g, so the boundary condition holds. Setting r = 0 kills every term but the first and gives u(0) = a_0/2, the average of g: the mean value property, recovered from the formula.

Worked example. Take a = 1 and g(theta) = 100 for |theta| < pi/2 and 0 otherwise, the half-hot disc from the opening. The data is even in theta, so every b_n vanishes. The constant is a_0/2 = 50. For n at least 1,

a_n = (100/pi) times the integral from -pi/2 to pi/2 of cos(n theta) d theta = (100/pi)(2 sin(n pi/2)/n) = 200 sin(n pi/2)/(n pi).

Now sin(n pi/2) is 0 for even n, +1 when n is 1, 5, 9, ..., and -1 when n is 3, 7, 11, .... So

u(r, theta) = 50 + (200/pi)[r cos(theta) - r^3 cos(3 theta)/3 + r^5 cos(5 theta)/5 - ...].

Evaluate at r = 1/2, theta = 0. The bracket is

0.5 - 0.0416667 + 0.00625 - 0.0011161 + 0.000217 - 0.0000444 + ... = 0.46384,

and 200/pi = 63.662, so u = 50 + 63.662 times 0.46384 = 79.53 degrees. The series converges quickly here because r = 1/2 makes each successive term smaller by a factor of about 8; near the rim, where r approaches 1, convergence slows markedly, which is a hint that a closed form would be more useful.

Summing the series: the Poisson kernel

Substitute the integral formulas for a_n and b_n into the series and exchange the sum with the integral. Using the identity cos(n theta) cos(n phi) + sin(n theta) sin(n phi) = cos(n(theta - phi)), the solution becomes

u(r, theta) = (1/(2 pi)) times the integral from -pi to pi of [1 + 2 times the sum over n of (r/a)^n cos(n(theta - phi))] g(phi) d phi.

The bracket is a summable series. With t = r/a and psi = theta - phi, the standard identity

1 + 2 times the sum over n from 1 of t^n cos(n psi) = (1 - t^2)/(1 - 2t cos(psi) + t^2)

holds for |t| < 1, and can be checked by writing cos(n psi) as the real part of a complex exponential and summing two geometric series. Substituting t = r/a and clearing denominators gives the Poisson integral formula:

u(r, theta) = (1/(2 pi)) times the integral from -pi to pi of [(a^2 - r^2)/(a^2 - 2ar cos(theta - phi) + r^2)] g(phi) d phi.

The bracketed function is the Poisson kernel. It expresses the value at any interior point as a weighted average of the boundary values, with the weights depending only on the geometry.

Reading the kernel

Three properties of the kernel carry the whole theory, and each is visible in the formula.

• It is positive. The numerator a^2 - r^2 is positive inside the disc, and the denominator is |a e^(i(theta-phi)) - r|^2 in disguise, hence positive. So the value at an interior point is a genuine weighted average, with no negative weights.
• It integrates to 1. Taking g identically 1 must give u identically 1, since constants are harmonic, so the kernel has total mass one. Combined with positivity, this gives the maximum principle again: an average of numbers between m and M lies between m and M.
• At the centre it is constant. Setting r = 0 makes the kernel equal a^2/a^2 = 1, so u(0) is the plain average of g. The mean value property is the special case of the Poisson formula at the centre.

As r approaches a, the numerator goes to zero while the denominator also goes to zero at phi = theta, so the kernel spikes near that angle and flattens elsewhere: the interior value is drawn almost entirely from the nearest piece of boundary. That concentration is what makes the boundary data get attained continuously at every point where g is continuous, which is the existence half of the Dirichlet problem on a disc. The upshot: the Poisson formula proves existence, uniqueness comes from the maximum principle, and together they settle the disc completely.

Common misconceptions

• "The singular solutions r^(-n) and ln(r) are wrong solutions." They are perfectly good harmonic functions away from the origin, and ln(r) is the fundamental solution used to build Green's functions in Lesson 13. They are excluded here only because this problem requires finiteness at the centre.
• "The Poisson kernel is a different method from the series." It is the same solution with the sum performed. Which form to use is a matter of convenience: the series near the centre, the integral near the rim.
• "An average of boundary values must be the mean value property." Only at the centre are the weights equal. Off centre the kernel weights nearby boundary points much more heavily, which is why the point at r = 1/2 on the hot side reads 79.5 rather than 50.
• "Periodicity in theta is a technicality." It is the eigenvalue condition. It is what forces lambda to be a perfect square and hence forces the exponents in the radial solutions to be integers.

Try it: a single harmonic on the rim

Exercise. Solve the Dirichlet problem on the unit disc with boundary data g(theta) = 3 + 2 cos(theta) - cos(3 theta), and evaluate the solution at r = 1/2, theta = 0.

Solution. The data is already a finite Fourier series, so no integrals are needed. Reading off a_0/2 = 3, a_1 = 2, a_3 = -1 and all else zero, the series solution is

u(r, theta) = 3 + 2 r cos(theta) - r^3 cos(3 theta).

At r = 1/2 and theta = 0 this gives 3 + 2(0.5) - (0.125) = 3.875. Check the boundary: at r = 1 the formula returns g exactly. Check harmonicity: in Cartesian terms r cos(theta) = x and r^3 cos(3 theta) = x^3 - 3xy^2, both of which have vanishing Laplacian.

Change one input. Replace cos(3 theta) by cos(3 theta)/r in the boundary data. On the rim, where r = 1, that is the same function, so the boundary data has not changed at all and the solution is unchanged. Boundary data is a function of angle alone, and any r appearing in it is evaluated at r = a before it means anything.

Recap

In polar coordinates Laplace's equation separates into an angular equation whose periodicity forces integer frequencies, and a radial Euler equation whose two solutions are r^n and r^(-n), of which boundedness at the centre keeps only the first. The resulting series multiplies the n-th Fourier coefficient of the boundary data by (r/a)^n. Summing the series gives the Poisson integral, whose kernel is positive, has total mass one, is constant at the centre, and concentrates near the boundary point being approached. Positivity plus unit mass gives the maximum principle; constancy at the centre gives the mean value property; concentration gives existence.

Looking ahead

Every problem so far has used one of two boundary conditions, prescribing the function or prescribing its normal derivative, and the choice has quietly decided whether the eigenfunctions were sines or cosines. The next lesson catalogues the conditions properly, says what each one means physically, explains the compatibility constraint that a pure Neumann problem must satisfy, and shows how a nonzero condition is reduced to a zero one by subtracting a steady state.

Sources

1. Lebl, J. (2026). Dirichlet problem in the circle and the Poisson kernel (Section 4.10). Notes on Diffy Qs: Differential Equations for Engineers. jirka.org
2. Wikipedia contributors. (n.d.). Poisson kernel. Wikipedia. en.wikipedia.org
3. Weisstein, E. W. (n.d.). Poisson kernel. Wolfram MathWorld. mathworld.wolfram.com
4. Wikipedia contributors. (n.d.). Dirichlet problem. Wikipedia. en.wikipedia.org
5. Strauss, W. A. (2008). Partial differential equations: An introduction (2nd ed.), Chapter 6. John Wiley and Sons.

Take a copper rod one metre long. Wrap the left end in a thick sleeve of glass wool so no heat escapes, and plunge the right end into an ice bath held at zero degrees. Two ends, two entirely different physical situations, and the mathematics has to say which is which. The ice bath fixes the value of the temperature at that end. The insulation fixes nothing about the value; what it fixes is the flux, which is zero.

Those are the two basic boundary conditions, and the difference between them changes the eigenfunctions, the long-time behaviour, and even whether a solution exists at all.

The catalogue

The Robin condition deserves a closer look, because it contains the other two. Newton's law of cooling says the heat flux out of a surface is proportional to the difference between the surface temperature and the ambient temperature, so at the right end of a rod

-k u_x(L, t) = h [u(L, t) - u_ambient],

with h the heat transfer coefficient. Let h go to zero, meaning perfect insulation, and the condition becomes u_x(L, t) = 0, which is Neumann. Let h grow without bound, meaning perfect contact with a reservoir, and dividing by h forces u(L, t) = u_ambient, which is Dirichlet. Key idea: Dirichlet and Neumann are the two extremes of one physical parameter, and real interfaces sit somewhere in between.

What each combination does to the eigenfunctions

Every separation of variables in this course produces the same equation X'' + lambda X = 0 on [0, L] and differs only in the conditions imposed at the two ends. Here are the four standard cases, each solved by the three-case argument of Lesson 5.

Check the third row, since it is the one from the opening scene. With X(0) = 0 the solution must be sin(mu x). Then X'(L) = mu cos(mu L) = 0 forces cos(mu L) = 0, so mu L is an odd multiple of pi/2, giving mu = (2n - 1) pi/(2L). The lowest mode is a quarter of a sine wave rather than a half, so the insulated-plus-cold rod cools at rate k pi^2/(4L^2), four times more slowly than a rod cold at both ends. Insulating one end really does double the effective length.

Note the second row too. The Neumann problem admits lambda = 0 with the constant eigenfunction, which the Dirichlet problem does not. That single extra mode is why an insulated rod settles to a nonzero average rather than to zero.

Compatibility: a Neumann problem can be impossible

Consider steady heat flow, so Laplacian(u) = 0 on a region D, with prescribed flux u_n = g all around the boundary. Integrate the equation over D and apply the divergence theorem:

0 = the integral over D of Laplacian(u) dA = the integral over the boundary of u_n ds = the integral over the boundary of g ds.

So a solution can exist only if g has total integral zero around the boundary. This is not a technicality but a statement of physics: at steady state the heat flowing in must exactly equal the heat flowing out, or the temperature could not have stopped changing. Prescribe a net inflow and there is no steady state.

The Neumann problem also fails uniqueness. If u is a solution then so is u + C for any constant, since the constant has zero normal derivative. So the pure Neumann problem is solvable exactly when the data balances, and then only up to an additive constant, which physically is the choice of reference temperature. Compare with the Dirichlet problem, which is always solvable and always unique. Why this matters: the type of boundary condition changes the existence and uniqueness theory, not just the formulas.

Turning inhomogeneous conditions into homogeneous ones

Separation of variables needs zero boundary conditions, because otherwise the superposition does not respect them. The standard fix is to subtract the steady state.

Worked example. Solve u_t = k u_xx on [0, 1] with u(0, t) = 0, u(1, t) = 100, and u(x, 0) = 0. Physically: a bar starting at zero degrees, with its right end suddenly clamped to 100.

First find the steady state, the solution with no time dependence: v'' = 0 with v(0) = 0 and v(1) = 100 gives the straight line v(x) = 100x. Now set w = u - v. Since v is independent of time and satisfies v_xx = 0, the function w satisfies the same heat equation, and it has homogeneous boundary conditions w(0, t) = w(1, t) = 0. Its initial data is

w(x, 0) = u(x, 0) - v(x) = -100x.

Expand that in sines. Using the integral from 0 to 1 of x sin(n pi x) dx = (-1)^(n+1)/(n pi),

b_n = 2 times the integral from 0 to 1 of (-100x) sin(n pi x) dx = 200(-1)^n/(n pi).

So the full solution is

u(x, t) = 100x + the sum over n of (200(-1)^n/(n pi)) exp(-k n^2 pi^2 t) sin(n pi x).

Two checks. As t grows every exponential dies and u approaches the linear profile 100x, which is the steady state we started by computing. At t = 0, evaluating at x = 1/2 gives 50 + (200/pi)(-1 + 1/3 - 1/5 + ...) = 50 + (200/pi)(-pi/4) = 50 - 50 = 0, matching the initial condition. The Leibniz series from Lesson 6 turns up as a consistency check.

How much data is the right amount

The type of the equation decides how conditions may be distributed, and Lesson 2 already sketched the rule.

• Parabolic. One initial condition, plus one boundary condition at each end for all time. Adding a second initial condition overdetermines the problem, since u_t is already fixed by u_xx.
• Hyperbolic. Two initial conditions, plus one boundary condition at each end. The extra initial condition is exactly what the second time derivative requires.
• Elliptic. One condition at every point of a closed boundary, and no initial data at all, because there is no distinguished variable to march in.

Prescribing both u and u_n on the boundary of an elliptic problem is asking for a Cauchy problem for Laplace's equation, and Lesson 16 shows that this is the classic ill-posed problem: solutions exist only for very special data, and they do not depend continuously on it.

Common misconceptions

• "Neumann conditions are just Dirichlet conditions on the derivative." Formally yes, structurally no. The Neumann eigenvalue problem has a zero eigenvalue and a constant eigenfunction, the pure Neumann problem needs a compatibility condition on the data, and its solution is unique only up to a constant.
• "An insulated end is a boundary condition on the temperature." It is a condition on the flux, and it says nothing at all about the temperature there, which is generally nonzero and time dependent.
• "Inhomogeneous boundary conditions block separation of variables." They block it directly, but subtracting a steady state converts the problem to a homogeneous one with modified initial data, at no cost.
• "More boundary conditions make a problem better determined." Only up to the right number. Prescribing both value and normal derivative on the boundary of an elliptic region is overdetermined and ill posed.

Try it: one end insulated

Exercise. Solve u_t = u_xx on [0, 1] with u_x(0, t) = 0, u(1, t) = 0, and u(x, 0) = cos(pi x/2).

Solution. The eigenvalue problem is X'' + lambda X = 0 with X'(0) = 0 and X(1) = 0. From X'(0) = 0 the solution must be cos(mu x), and X(1) = cos(mu) = 0 forces mu = (2n - 1) pi/2. So the eigenfunctions are cos((2n-1) pi x/2) with eigenvalues ((2n-1) pi/2)^2. The initial data is exactly the first eigenfunction, so only that mode is present:

u(x, t) = cos(pi x/2) exp(-pi^2 t/4).

Check: u_t = -(pi^2/4)u and u_xx = -(pi^2/4)u, so the equation holds; u_x(0, t) = 0 since the sine vanishes at zero; u(1, t) = cos(pi/2) exp(...) = 0. All conditions are met.

Change one input. Insulate the right end too, so u_x(1, t) = 0, keeping the same initial data. Now the eigenfunctions are cos(n pi x), including the constant, and the solution tends as t grows to the average value of the initial data, the integral from 0 to 1 of cos(pi x/2) dx = 2/pi, roughly 0.637, rather than to zero. Insulating the second end changes not just the rate but the destination.

What to carry forward

Dirichlet fixes the value, Neumann the flux, Robin a combination that reduces to each in a limit, and periodic conditions describe a loop. Each choice produces a different eigenvalue problem: sines with no zero mode, cosines with one, quarter-wave modes for a mixed pair. A pure Neumann problem requires its data to have zero net flux and determines the solution only up to a constant. An inhomogeneous condition is removed by subtracting the steady state, which shifts the work into the initial data. And the type of the equation dictates how many conditions are appropriate and where.

Looking ahead

Four sets of eigenfunctions have now appeared, all from the same differential equation with different conditions attached, and in each case they turned out to be mutually orthogonal and to support an expansion. That is not a coincidence. The next module identifies the general framework, Sturm-Liouville theory, which explains why orthogonality always holds, why the eigenvalues are real and increase without bound, and why the expansions converge. It also handles equations with variable coefficients, where the eigenfunctions are Bessel functions rather than sines.

Sources

1. Wikipedia contributors. (n.d.). Boundary value problem. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Robin boundary condition. Wikipedia. en.wikipedia.org
3. Weisstein, E. W. (n.d.). Neumann boundary conditions. Wolfram MathWorld. mathworld.wolfram.com
4. Dawkins, P. (n.d.). Heat equation with non-zero temperature boundaries. Paul's Online Notes, Differential Equations. tutorial.math.lamar.edu
5. Lebl, J. (2026). Boundary value problems (Section 4.1). Notes on Diffy Qs: Differential Equations for Engineers. jirka.org

Strike a guitar string and its overtones are exact integer multiples of the fundamental: 1, 2, 3, 4. Strike a circular drumhead and they are not. The first few ratios are approximately 1, 1.593, 2.136, 2.295, and no two of them are in a simple whole-number relation. That is why a plucked string sings a definite pitch and a drum, however carefully tuned, does not.

Those numbers are ratios of zeros of Bessel functions, and the reason they appear is that a drum's eigenvalue problem is not X'' + lambda X = 0. It has a variable coefficient, because a circle has more area near its rim than near its centre. This lesson identifies the general framework that covers both cases, and every other case in this course, and proves the two facts that all the expansion methods have been quietly assuming.

The general form

A regular Sturm-Liouville problem on an interval [a, b] is

(p(x) X')' + q(x) X + lambda w(x) X = 0,

with separated boundary conditions alpha X(a) + beta X'(a) = 0 and gamma X(b) + delta X'(b) = 0, where p, p', q and w are continuous, p > 0 and w > 0 on the closed interval, and neither pair of boundary constants is both zero. The function w is called the weight, and it is the one piece of the structure that has no analogue in the constant-coefficient case, where it is simply 1.

The last two rows are singular rather than regular, because p vanishes at an endpoint. The theory still applies with modifications, and boundedness at the singular endpoint replaces the boundary condition there, exactly as it did on the disc in Lesson 10.

The theorem the whole course has been using

Theorem. For a regular Sturm-Liouville problem:

• The eigenvalues are real, and they form an increasing sequence tending to infinity.
• Each eigenvalue has a one-dimensional space of eigenfunctions, and the eigenfunction for the n-th eigenvalue has exactly n - 1 zeros inside the interval.
• Eigenfunctions for different eigenvalues are orthogonal with respect to the weight: the integral of w X_m X_n over [a, b] is zero when m is not n.
• The eigenfunctions are complete: any function with finite weighted energy can be expanded as f = the sum of c_n X_n, converging in the weighted mean-square sense, with

c_n = [the integral of w f X_n] / [the integral of w X_n^2].

Two of these can be proved here in a few lines, and the two proofs share one identity.

The Lagrange identity, and what it gives

Write L[X] = (pX')' + qX, so the problem is L[X] = -lambda w X. For any two twice-differentiable functions u and v,

u L[v] - v L[u] = u(pv')' - v(pu')' = [p(u v' - v u')]',

as you can check by expanding the right side and cancelling the p u' v' terms; the q terms cancel between the two sides. Integrating from a to b,

the integral of (u L[v] - v L[u]) dx = [p(u v' - v u')] evaluated at b minus at a.

Lemma. If u and v both satisfy the separated boundary conditions, that bracket vanishes at both ends.

Proof. At x = a the two conditions alpha u(a) + beta u'(a) = 0 and alpha v(a) + beta v'(a) = 0 say the vector (alpha, beta), which is not the zero vector, is orthogonal to both (u(a), u'(a)) and (v(a), v'(a)). Two vectors in the plane orthogonal to the same nonzero vector are parallel, so the determinant u(a)v'(a) - v(a)u'(a) is zero. Multiplying by p(a) leaves zero. Same at b. This completes the proof.

Theorem (orthogonality). Eigenfunctions for distinct eigenvalues are orthogonal with weight w.

Proof. Let L[X_m] = -lambda_m w X_m and L[X_n] = -lambda_n w X_n. Substituting into the integrated identity with u = X_m and v = X_n, the left side becomes

the integral of (X_m (-lambda_n w X_n) - X_n (-lambda_m w X_m)) dx = (lambda_m - lambda_n) times the integral of w X_m X_n dx,

and the right side is zero by the lemma. Since lambda_m is not lambda_n, the integral vanishes. This completes the proof.

Theorem (reality). Every eigenvalue is real.

Proof. Suppose L[X] = -lambda w X with X not identically zero. Taking complex conjugates and using that p, q and w are real gives L[conjugate of X] = -(conjugate of lambda) w (conjugate of X). Running the previous computation with u = X and v = conjugate of X yields

(lambda - conjugate of lambda) times the integral of w |X|^2 dx = 0.

The integral is strictly positive, since w > 0 and X is not identically zero, so lambda equals its own conjugate and is therefore real. This completes the proof.

Key idea: orthogonality is not a lucky property of sines. It is forced by the self-adjoint structure of the operator together with the separated boundary conditions, and it holds for Bessel functions, Legendre polynomials and every other eigenfunction family in this subject.

A Robin problem with no closed-form eigenvalues

Solve X'' + lambda X = 0 on [0, 1] with X(0) = 0 and X'(1) + X(1) = 0, a mixture of a clamp at one end and Newton cooling at the other.

As before, only positive lambda = mu^2 can work, and X(0) = 0 forces X = sin(mu x). The second condition gives

mu cos(mu) + sin(mu) = 0, that is, tan(mu) = -mu.

This transcendental equation has no closed-form solution, but it has infinitely many roots, one just below each odd multiple of pi/2 beyond the first. Numerically the first three are

mu_1 = 2.0288, mu_2 = 4.9132, mu_3 = 7.9787,

so lambda_1 = 4.116, lambda_2 = 24.14, lambda_3 = 63.66. The eigenfunctions sin(mu_n x) are still orthogonal on [0, 1] with weight 1, by the theorem, even though the mu_n are not nice numbers, and expansions in them converge exactly as Fourier sine series do. That is what the general theory buys: you lose the closed forms and keep every structural property.

The drum

Separating the wave equation on a disc of radius a in polar coordinates, for solutions independent of angle, produces the radial problem

(r R')' + lambda r R = 0 on [0, a], with R bounded at 0 and R(a) = 0.

Here p(r) = r and the weight is w(r) = r, which is the area element of the plane: outer rings count for more because there is more of them. Substituting s = sqrt(lambda) r turns the equation into Bessel's equation of order zero, whose bounded solution is J_0(s). So R(r) = J_0(sqrt(lambda) r), and the condition R(a) = 0 requires

sqrt(lambda) a = j_(0,n), the n-th positive zero of J_0.

Those zeros begin 2.40483, 5.52008, 8.65373, and they are not evenly spaced. Including the modes that depend on angle brings in the zeros of J_1, J_2 and so on, and the lowest four frequencies of a drumhead are proportional to

j_(0,1) = 2.40483, j_(1,1) = 3.83171, j_(2,1) = 5.13562, j_(0,2) = 5.52008,

whose ratios to the first are 1, 1.593, 2.136, 2.295. There are the numbers from the opening paragraph. The orthogonality that lets you expand an initial drumhead shape in these functions is the weighted orthogonality of the theorem, with weight r.

Common misconceptions

• "Orthogonality is a property of trigonometric functions." It is a property of the operator and boundary conditions. Bessel functions are orthogonal with weight r, Legendre polynomials with weight 1 on [-1, 1], and the Robin eigenfunctions above with weight 1 despite having irrational frequencies.
• "The weight is a normalisation convention." It is part of the problem. Omitting the r in a disc expansion gives coefficients that are simply wrong, because the integrals are not the ones the theorem certifies.
• "Eigenvalues always have closed forms." The Robin example has none. This is the normal situation, and the theory is arranged so that nothing depends on having closed forms.
• "A drum can be tuned to a pitch like a string." Its overtones are ratios of Bessel zeros, which are not integer multiples, so a drumhead produces a sound with no clear harmonic series. Timpani are shaped and loaded specifically to push a few modes closer to harmonic ratios.

Try it: identify the Sturm-Liouville data

Exercise. Put x^2 X'' + 2x X' + lambda X = 0 into Sturm-Liouville form and identify p, q and w.

Solution. Notice that x^2 X'' + 2x X' = (x^2 X')', since differentiating the product gives 2x X' + x^2 X''. So the equation is already

(x^2 X')' + lambda X = 0,

which is Sturm-Liouville with p(x) = x^2, q(x) = 0 and w(x) = 1. On an interval away from the origin it is regular, and its eigenfunctions are orthogonal with weight 1.

Change one input. Consider x^2 X'' + x X' + lambda X = 0 instead, where the middle coefficient is x rather than 2x. Now the left side is not a derivative of a product. Dividing by x gives x X'' + X' + (lambda/x) X = 0, which is (x X')' + lambda (1/x) X = 0: Sturm-Liouville with p = x, q = 0 and weight w = 1/x. The same equation can hide a weight, and finding it is the whole content of putting a problem in standard form.

Summing up

Every boundary value problem in this course is a Sturm-Liouville problem, with a coefficient p, a potential q and a weight w. The Lagrange identity converts the self-adjointness of the operator into a boundary term that separated conditions annihilate, and from that one identity come both the reality of the eigenvalues and the weighted orthogonality of the eigenfunctions. Completeness then licenses the expansions used since Lesson 5. The framework survives without closed forms, as the Robin problem shows, and it handles the variable coefficient of a circular drum, whose non-integer overtone ratios are ratios of Bessel zeros.

Looking ahead

Eigenfunction expansions solve problems driven by their initial data. Many problems are driven instead by a source: a heater inside a bar, a load on a beam, a charge in a region. For those, the natural object is a function that answers the question of what happens if you put a unit source at one point and nothing anywhere else. Adding up those answers with weights solves the general problem, and the next lesson builds those Green's functions explicitly.

Sources

1. Wikipedia contributors. (n.d.). Sturm-Liouville theory. Wikipedia. en.wikipedia.org
2. Lebl, J. (2026). Sturm-Liouville problems (Section 5.1). Notes on Diffy Qs: Differential Equations for Engineers. jirka.org
3. National Institute of Standards and Technology. (n.d.). Chapter 10: Bessel functions. NIST Digital Library of Mathematical Functions. dlmf.nist.gov
4. Dawkins, P. (n.d.). Eigenvalues and eigenfunctions. Paul's Online Notes, Differential Equations. tutorial.math.lamar.edu
5. Wikipedia contributors. (n.d.). Bessel function. Wikipedia. en.wikipedia.org

George Green ran a windmill in Nottingham and taught himself mathematics on its top floor. On 14 December 1827 he placed an advertisement in the Nottingham Review announcing a book to be published by subscription, and in March 1828 An Essay on the Application of Mathematical Analysis to the Theories of Electricity and Magnetism appeared with 51 subscribers, most of whom cannot have read it. It contains the idea this lesson is about, and it is one of the most useful ideas in the subject.

The idea is this. Instead of solving a forced problem for every possible forcing, solve it once for a single point source, and then build any other answer by adding up point sources with weights. If the equation is linear, superposition guarantees the answer is right.

The simplest case, built from scratch

Consider the boundary value problem

-u'' = f(x) on [0, L], with u(0) = u(L) = 0.

Physically: a string under tension, held at both ends, sagging under a distributed load f. The Green's function G(x, s) is the sag at position x caused by a unit point load applied at position s. Once you have it,

u(x) = the integral from 0 to L of G(x, s) f(s) ds,

because the distributed load is a superposition of point loads of strength f(s) ds.

To find G, note what a point load does. Away from s, there is no force, so G satisfies G'' = 0 and is a straight line on each side. It vanishes at both ends. It is continuous at s, since a string does not tear. And its slope jumps at s, by an amount determined by the size of the load: integrating -G'' = delta(x - s) across the point gives

G'(s+) - G'(s-) = -1.

Now solve. Write G = A x for x below s, which vanishes at 0, and G = B(L - x) for x above s, which vanishes at L. Continuity at s gives A s = B(L - s). The jump condition gives -B - A = -1, that is, A + B = 1. Solving the pair,

A = (L - s)/L and B = s/L,

so

G(x, s) = x(L - s)/L for x <= s, and G(x, s) = s(L - x)/L for x >= s.

The two branches are the same expression with x and s interchanged, so G(x, s) = G(s, x). That symmetry is reciprocity: the sag at x from a load at s equals the sag at s from an equal load at x. It is a general feature of self-adjoint problems, and it is a genuine physical prediction that can be tested with a real string.

Using it, and checking

Worked example. Take L = 1 and f = 1, a uniformly loaded string. Then

u(x) = the integral from 0 to x of s(1 - x) ds + the integral from x to 1 of x(1 - s) ds,

splitting the integral at s = x so that each piece uses the correct branch. The first integral is (1 - x) x^2/2. The second is x times the integral from x to 1 of (1 - s) ds, which is x (1 - x)^2/2. Adding,

u(x) = (x(1 - x)/2)[x + (1 - x)] = x(1 - x)/2.

Check directly: u'' = -1, so -u'' = 1 = f; and u(0) = u(1) = 0. The parabola sags to a maximum of 1/8 at the midpoint, which is the classic answer for a uniformly loaded string.

Second example, no extra work. Take f(s) = sin(pi s) on [0, 1]. Rather than integrate the Green's function, note that -u'' = sin(pi x) with zero ends is solved by u = sin(pi x)/pi^2, and the Green's function integral must reproduce it. The point of G is that it handles every f at once, including ones for which no clever guess is available, such as a load concentrated on part of the interval.

What a delta function is, honestly

The symbol delta(x - s) is not a function. No function is zero everywhere except one point and still has integral 1. It is defined by what it does inside an integral:

the integral of delta(x - s) g(x) dx = g(s)

for every continuous g. Everything above can be phrased without it: G is characterised by being harmonic on each side of s, continuous at s, satisfying the boundary conditions, and having a slope jump of -1. That characterisation is what we actually solved, and the delta is bookkeeping for it. Worth holding on to: when a computation involving a delta feels shaky, restate it as a continuity condition plus a jump condition and the shakiness disappears.

The eigenfunction form

There is a second route to the same object, using the previous lesson. Expand both G and f in the eigenfunctions of the problem. For -u'' = f with zero ends on [0, L], the eigenfunctions are sin(n pi x/L) with eigenvalues (n pi/L)^2, and one finds

G(x, s) = (2/L) times the sum over n of sin(n pi x/L) sin(n pi s/L) / (n pi/L)^2.

The symmetry in x and s is manifest here, and so is the reason the series converges nicely: the eigenvalues appear in the denominator, so high modes are suppressed. This form generalises immediately to any Sturm-Liouville problem, replacing sines by the appropriate eigenfunctions and dividing by the corresponding eigenvalues. It also shows what goes wrong when zero is an eigenvalue, as it is for a pure Neumann problem: one denominator vanishes, and the Green's function does not exist, which is the compatibility failure of Lesson 11 seen from another angle.

Two dimensions, and the method of images

For Poisson's equation Laplacian(u) = f in the plane, the free-space Green's function is

(1/(2 pi)) ln|X - Y|,

the potential of a unit point source at Y. It is harmonic away from Y, as noted in Lesson 10 where ln(r) appeared and was discarded for blowing up at the origin. Here that blow-up is exactly the point.

To solve the Dirichlet problem on a region, one needs a Green's function vanishing on the boundary, and for simple geometries the method of images supplies it. Take the upper half plane {y > 0} and a source at Y = (a, b) with b > 0. Place an equal and opposite source at the mirror point Y* = (a, -b), which lies outside the region so it does not disturb the equation there. Then

G(X, Y) = (1/(2 pi))[ln|X - Y| - ln|X - Y*|]

vanishes on the axis y = 0, because a point on the axis is equidistant from Y and its mirror image, so the two logarithms cancel. The physical picture is a grounded conducting plate: the field above it is the same as that of the original charge plus an image charge of opposite sign below.

Differentiating this Green's function along the boundary normal and integrating against the boundary data produces the Poisson kernel for the half plane:

u(x, y) = (y/pi) times the integral over all t of g(t)/[(x - t)^2 + y^2] dt.

Check the structure against the disc kernel of Lesson 10: it is positive, it integrates to 1 in t for each fixed y > 0, and as y shrinks it concentrates near t = x. The same three properties, the same three consequences: weighted average, maximum principle, and continuous attainment of the data.

Common misconceptions

• "The delta function is a very tall narrow spike." That is a useful picture of approximating sequences, not a definition. The definition is the sampling property inside an integral, and every honest computation can be rewritten with continuity and jump conditions instead.
• "A Green's function always exists." It exists when the homogeneous problem has only the zero solution. If zero is an eigenvalue, as for a pure Neumann problem, the construction fails, and that failure is the compatibility condition reappearing.
• "Symmetry of G is a coincidence of the example." It is reciprocity, and it holds for every self-adjoint problem, as the eigenfunction form makes obvious.
• "The image charge is physically present." It is a construction that happens to satisfy the boundary condition. It lives outside the region, where the equation is not being solved, which is exactly why it is allowed to be there.

Try it: a Green's function with mixed conditions

Exercise. Build the Green's function for -u'' = f on [0, 1] with u(0) = 0 and u'(1) = 0.

Solution. On each side of s the function is linear. To the left it must vanish at 0, so G = A x. To the right it must have zero slope at 1, so it is constant: G = C. Continuity at s gives A s = C. The jump condition gives G'(s+) - G'(s-) = 0 - A = -1, so A = 1. Hence C = s and

G(x, s) = x for x <= s, and G(x, s) = s for x >= s,

compactly G(x, s) = min(x, s), which is symmetric as reciprocity requires. Test it with f = 1: u(x) = the integral from 0 to x of s ds + the integral from x to 1 of x ds = x^2/2 + x(1 - x) = x - x^2/2. Check: u'' = -1, u(0) = 0, and u'(1) = 1 - 1 = 0. All three conditions hold.

Change one input. Try Neumann conditions at both ends, u'(0) = u'(1) = 0. Now the left branch must be constant and so must the right, and continuity forces them equal, leaving no room for a slope jump. No Green's function exists, matching the fact that -u'' = f with two Neumann conditions is solvable only when f integrates to zero.

The takeaway

A Green's function is the response to a unit point source, and integrating it against a forcing solves the forced problem. In one dimension it is pinned down by four requirements: the homogeneous equation on each side, the boundary conditions, continuity at the source, and a prescribed jump in the derivative. It is symmetric in its two arguments, a genuine reciprocity statement. It can also be written as an eigenfunction series with the eigenvalues in the denominator, which shows immediately why a zero eigenvalue destroys it. In the plane the free-space solution is a logarithm, and for a half plane the method of images produces the boundary-vanishing version and, from it, the half-plane Poisson kernel.

Looking ahead

Everything so far has lived on a bounded interval or region, where eigenvalues are discrete. On an infinite line there is no quantisation: every frequency is allowed, and the sum over modes becomes an integral over frequencies. That is the Fourier transform, and it turns the heat equation on the whole line into an algebraic problem whose solution is a Gaussian. The Laplace transform does the analogous job in the time variable for semi-infinite problems.

Sources

1. Weisstein, E. W. (n.d.). Green's function. Wolfram MathWorld. mathworld.wolfram.com
2. Wikipedia contributors. (n.d.). Green's function. Wikipedia. en.wikipedia.org
3. O'Connor, J. J., & Robertson, E. F. (1998). George Green. MacTutor History of Mathematics Archive. mathshistory.st-andrews.ac.uk
4. Wikipedia contributors. (n.d.). Method of images. Wikipedia. en.wikipedia.org
5. Wikipedia contributors. (n.d.). Dirac delta function. Wikipedia. en.wikipedia.org

Release a single drop of dye at one point of a long, still channel of water and photograph the concentration profile at intervals. What you see is a bell curve that flattens and widens, and its width grows like the square root of the elapsed time, not proportionally to the time. Wait four times as long and the dye has spread twice as far, not four times.

That square root is not an empirical curiosity. It falls out of solving the heat equation on an infinite line, and the method that produces it is the one this lesson develops: on an unbounded domain, the discrete sum over modes becomes a continuous integral over frequencies.

From series to integrals

On a bounded interval, boundary conditions quantise the allowed frequencies to a discrete list. On the whole line there is no boundary and no quantisation, so every frequency is allowed and the sum becomes an integral. Define the Fourier transform of a function u(x) by

U(xi) = the integral over all x of u(x) e^(-i xi x) dx,

with inversion

u(x) = (1/(2 pi)) times the integral over all xi of U(xi) e^(i xi x) d xi.

Other conventions distribute the 2 pi differently; what matters is consistency. The single property that makes the transform useful for differential equations is what it does to derivatives. Integrating by parts, and assuming u decays at infinity,

the transform of u_x is i xi U(xi), and hence the transform of u_xx is -xi^2 U(xi).

Key idea: differentiating in x becomes multiplying by i xi. A differential operator in space becomes an algebraic one, and a PDE in two variables becomes an ODE in one.

The heat equation on the line

Solve u_t = k u_xx for all real x and t > 0, with u(x, 0) = phi(x). Transform in x, treating t as a parameter. The left side becomes U_t and the right becomes -k xi^2 U, so

U_t = -k xi^2 U, giving U(xi, t) = U(xi, 0) e^(-k xi^2 t) = Phi(xi) e^(-k xi^2 t),

where Phi is the transform of the initial data. So the entire effect of running the heat equation for time t is to multiply each frequency component by e^(-k xi^2 t). This is exactly the factor from Lesson 8, with the discrete n pi/L replaced by the continuous xi.

To recover u, invert. A product of transforms corresponds to a convolution, and the inverse transform of the Gaussian e^(-k xi^2 t) is another Gaussian:

S(x, t) = (1/sqrt(4 pi k t)) e^(-x^2/(4 k t)),

the heat kernel. Hence

u(x, t) = the integral over all y of S(x - y, t) phi(y) dy = (1/sqrt(4 pi k t)) times the integral of e^(-(x-y)^2/(4kt)) phi(y) dy.

The kernel is a normalised Gaussian: its integral over all x is 1 for every t, which says heat is conserved. Its standard deviation is sqrt(2 k t), so the profile widens like the square root of time, which is the observation from the opening paragraph and the reason diffusion over a distance d takes a time proportional to d^2.

Three readings of the kernel

• Infinite propagation speed. For any t > 0, the exponential is strictly positive at every x. So an initial profile confined to a small interval produces, immediately, a nonzero temperature everywhere, however far away. Lesson 2 asserted this; here is the formula that proves it.
• Smoothing. The kernel is infinitely differentiable in x, and differentiating under the integral sign transfers those derivatives to u without touching phi. So u is infinitely smooth for t > 0 even if phi is a step function.
• Point source. Taking phi to be a unit point release at the origin gives u(x, t) = S(x, t) exactly. The heat kernel is the Green's function of the heat equation on the line, in the sense of the previous lesson.

Worked example. Take phi(x) = 1 for |x| <= 1 and 0 otherwise: a slab of hot material one unit wide on each side of the origin. Then

u(x, t) = (1/sqrt(4 pi k t)) times the integral from -1 to 1 of e^(-(x-y)^2/(4kt)) dy.

Substituting z = (y - x)/sqrt(4kt) turns the integral into a standard error function computation, and

u(x, t) = (1/2)[erf((1 - x)/sqrt(4kt)) + erf((1 + x)/sqrt(4kt))],

where erf is the error function, the normalised integral of a Gaussian. Two checks. At the centre, u(0, t) = erf(1/sqrt(4kt)), which is close to 1 for small t and decays toward 0 as t grows, since the heat spreads out. And as t shrinks to zero with |x| < 1 fixed, both arguments grow without bound, both error functions approach 1, and u approaches 1, recovering the initial data.

The Laplace transform, in time

The Fourier transform handles unbounded space. For a problem on t > 0 with an initial condition and possibly time-dependent boundary data, the natural tool is the Laplace transform

U(s) = the integral from 0 to infinity of u(t) e^(-s t) dt,

whose defining property for differential equations is

the transform of u_t is s U(s) - u(0).

Notice that the initial condition is built into the transform rather than imposed afterwards, which is what makes the method suit initial value problems.

Worked example: a semi-infinite bar. Solve u_t = k u_xx on x > 0 and t > 0, with u(x, 0) = 0, u(0, t) = T_0 for t > 0, and u bounded as x grows. Physically: a cold half-space whose face is suddenly brought to temperature T_0 and held there.

Transform in t. Since u(x, 0) = 0, the equation becomes

s U = k U_xx, that is, U_xx - (s/k) U = 0.

The general solution is a combination of e^(x sqrt(s/k)) and e^(-x sqrt(s/k)), and boundedness discards the growing one. So U(x, s) = A(s) e^(-x sqrt(s/k)). The boundary condition transforms to U(0, s) = T_0/s, since the Laplace transform of a constant T_0 is T_0/s. So A = T_0/s and

U(x, s) = (T_0/s) e^(-x sqrt(s/k)).

Inverting this is a standard table entry, and the result is

u(x, t) = T_0 erfc(x/(2 sqrt(k t))),

where erfc = 1 - erf is the complementary error function. Check the two conditions: at x = 0, erfc(0) = 1 and u = T_0; for fixed x > 0 as t shrinks to zero the argument grows and erfc tends to 0, recovering the cold initial state.

The engineering content is in the argument. The temperature is appreciable only where x is not much larger than 2 sqrt(k t), which is called the penetration depth. Since erfc(1) = 0.157, the point at x = 2 sqrt(kt) has reached about 16 percent of the surface temperature. That single expression tells you how deep a cold snap reaches into soil, how thick a fire-resistant layer must be, and why the square-root law of the opening paragraph is unavoidable.

Which transform, and when

The last row is worth stating explicitly. Transforms and series are the same idea applied to different geometries: decompose into modes, solve for each mode separately because the operator acts on it by multiplication, and reassemble. Which modes exist is decided by the domain.

Common misconceptions

• "The Fourier transform requires periodicity." Fourier series require periodicity. The transform is designed for functions on the whole line that decay, and it replaces the discrete spectrum with a continuous one.
• "The heat kernel is an approximation of the true solution." It is the exact solution for a point release, and convolving it with any integrable initial profile gives the exact solution for that profile.
• "Penetration depth is where the heat stops." The solution is strictly positive at every depth for every positive time. The penetration depth is where the effect ceases to be appreciable, and erfc(1) = 0.157 quantifies what appreciable means.
• "The Laplace transform is only for ODEs." Applied in time to a PDE, it converts the time derivative into multiplication and leaves an ODE in the space variable, which is exactly what happened above.

Try it: transport with a transform

Exercise. Solve u_t + c u_x = 0 on the whole line with u(x, 0) = phi(x), using the Fourier transform, and confirm the answer against Lesson 3.

Solution. Transforming in x turns u_x into i xi U, so the equation becomes U_t + i c xi U = 0, an ODE in t for each xi. Its solution is

U(xi, t) = Phi(xi) e^(-i c xi t).

Multiplying a transform by e^(-i xi a) corresponds to translating the function by a, so with a = ct the inverse transform is u(x, t) = phi(x - ct). That is the travelling profile found by characteristics in Lesson 3, obtained here without a single characteristic curve.

Change one input. Add a diffusion term, giving u_t + c u_x = k u_xx. The transform now gives U_t = (-i c xi - k xi^2)U, so U = Phi(xi) e^(-i c xi t) e^(-k xi^2 t): a translation multiplied by the heat kernel factor. The solution is the initial profile shifted by ct and simultaneously smeared by a Gaussian of width sqrt(2kt), which is exactly what a pollutant in a river does.

Where this leaves us

On an unbounded domain the Fourier transform replaces the eigenfunction series, converting u_xx into multiplication by -xi^2 and the heat equation into a first-order ODE in time for each frequency. Inverting gives the heat kernel, a Gaussian whose width grows like the square root of time, which proves infinite propagation speed, instant smoothing and the conservation of total heat in one formula. The Laplace transform does the corresponding job in time, and applied to a suddenly heated half-space it yields the complementary error function solution and the penetration depth 2 sqrt(kt).

Looking ahead

The final module returns to waves and asks what changes in more than one space dimension. The answer is surprising: a sharp sound in three dimensions passes a listener and stops cleanly, but the same event in two dimensions leaves a decaying tail that never quite ends. That difference is Huygens' principle, and it is the reason a clap sounds like a clap in air and like a drawn-out rumble on the surface of a pond.

Sources

1. Weisstein, E. W. (n.d.). Fourier transform. Wolfram MathWorld. mathworld.wolfram.com
2. Weisstein, E. W. (n.d.). Laplace transform. Wolfram MathWorld. mathworld.wolfram.com
3. Wikipedia contributors. (n.d.). Heat kernel. Wikipedia. en.wikipedia.org
4. Lebl, J. (2026). Solving PDEs with the Laplace transform (Section 6.5). Notes on Diffy Qs: Differential Equations for Engineers. jirka.org
5. National Institute of Standards and Technology. (n.d.). Section 1.14: Integral transforms. NIST Digital Library of Mathematical Functions. dlmf.nist.gov

Clap your hands in an open field. A listener a hundred metres away hears the clap about 0.29 seconds later, since sound travels at roughly 343 metres per second in air at twenty degrees, and then hears nothing more. The signal arrives, and it stops.

Now stretch a large rubber membrane flat and tap it once at a point. A sensor elsewhere on the membrane does not record a single sharp pulse. It records an arrival followed by a tail that decays but never quite ends. The same equation, the same initial data, and the only difference is that the membrane is two dimensional and the air is three dimensional.

That difference has a name and a precise statement, and it is the reason speech works.

The equation in n dimensions

The wave equation in n space dimensions is

u_tt = c^2 (u_(x_1 x_1) + ... + u_(x_n x_n)),

written compactly as u_tt = c^2 Laplacian(u), with initial data u = phi and u_t = psi at t = 0. The derivation is the same as for the string: Newton's law applied to a small piece of a taut membrane or a small parcel of gas.

Spherical waves in three dimensions

Start with the case that admits a complete hand computation: a solution in three dimensions depending only on the distance r from the origin and on time. In spherical coordinates the radial part of the Laplacian is

u_rr + (2/r) u_r,

so the equation reads u_tt = c^2 (u_rr + (2/r) u_r). Now make the substitution v = r u, so u = v/r. Then

u_r = v_r/r - v/r^2 and u_rr = v_rr/r - 2 v_r/r^2 + 2v/r^3.

Adding (2/r)u_r = 2v_r/r^2 - 2v/r^3, the last two terms of each cancel in pairs and

u_rr + (2/r)u_r = v_rr/r.

Since u_tt = v_tt/r, multiplying through by r gives

v_tt = c^2 v_rr,

the one-dimensional wave equation. So v = F(r - ct) + G(r + ct) and

u(r, t) = [F(r - ct) + G(r + ct)]/r.

Two facts fall out. An outgoing spherical wave keeps its shape exactly, translating outward at speed c, so a sharp pulse stays sharp. And its amplitude falls off like 1/r, hence its intensity, proportional to the square, falls off like 1/r^2. The inverse square law of sound and light is this factor. Key idea: the three-dimensional wave equation is the one-dimensional one in disguise, once you multiply by the radius.

Kirchhoff and Poisson formulas

For general data the closed-form solutions are these. In three dimensions, Kirchhoff's formula gives

u(X, t) = (d/dt)[the average of phi over the sphere of radius ct about X] times t + t times [the average of psi over that sphere].

Everything on the right is computed from data on a sphere of radius exactly ct. In two dimensions, Poisson's formula gives

u(X, t) = (1/(2 pi c)) [ (d/dt) the integral over the disc |Y - X| < ct of phi(Y)/sqrt(c^2 t^2 - |Y - X|^2) dY + the integral over that disc of psi(Y)/sqrt(c^2 t^2 - |Y - X|^2) dY ].

Here the integration is over the whole disc of radius ct, not just its boundary circle. That single structural difference is everything.

Where does the two-dimensional formula come from? By descent: a two-dimensional problem is a three-dimensional problem whose data happens not to depend on the third coordinate. Integrating the three-dimensional sphere formula over that ignored direction converts the sphere integral into an integral over its shadow, which is a disc.

The strong Huygens principle

Compare the domains of dependence. In three dimensions, u(X, t) depends on the data only on the sphere |Y - X| = ct. In two dimensions it depends on the data throughout the disc |Y - X| <= ct. The strong Huygens principle is the statement that the first behaviour occurs, and it holds in odd space dimensions three and above but not in even dimensions.

Trace the consequence. Take initial data supported in a small ball around a point P. In three dimensions, a listener at distance d sees nothing until ct reaches d minus the source radius, hears the signal while the expanding sphere sweeps across the source region, and then hears nothing again, because from then on the sphere of radius ct misses the support entirely. In two dimensions the disc of radius ct, once it has reached the source, contains it forever, so the integral never returns to zero. The signal has a tail.

One dimension is worth checking against this, because it is odd and yet has a tail: d'Alembert's formula integrates psi over the whole interval [x - ct, x + ct], so an initial velocity kick leaves a permanent displacement, exactly as Lesson 4 found. The correct statement is that strong Huygens holds in odd dimensions three and above, and one dimension is an exception.

Why this matters: if air were two dimensional, every syllable you spoke would be smeared over the echoes of every syllable before it. Clean transmission of sharp signals is a property of three-dimensional space, and it is why acoustics and optics work at all.

A caution about the pond. Ripples on water do have a long tail after the leading edge, but the dominant reason is dispersion: water waves of different wavelengths travel at different speeds, and the governing equation is not the wave equation. The clean two-dimensional example is a stretched membrane, where the wave equation really does apply. Physical intuition and mathematical statement should be kept apart here.

A note on the optical principle

The name also attaches to an older and different idea. Christiaan Huygens proposed in the seventeenth century that every point of a wavefront acts as a source of secondary wavelets whose envelope forms the next wavefront, and Fresnel added the interference of those wavelets in the nineteenth. That construction, the Huygens-Fresnel principle, explains reflection, refraction and diffraction. It is related to but distinct from the sharp statement above about domains of dependence, and the two are often conflated.

Common misconceptions

• "Huygens' principle says waves spread out." Waves spread in every dimension. The principle is the much sharper statement that in three dimensions the solution depends only on data on a sphere, so signals arrive and then completely stop.
• "Odd dimensions are the good ones." Odd dimensions three and above are. One dimension has a tail from initial velocity, as d'Alembert's formula shows.
• "The 1/r decay means energy is lost." Energy is conserved: the amplitude falls like 1/r so the intensity falls like 1/r^2, exactly compensating the growth of the sphere's area, which is proportional to r^2.
• "Ripples on a pond are the two-dimensional wave equation." They are governed by a dispersive equation. The membrane is the honest example.

Try it: a spherical pulse

Exercise. A spherically symmetric outgoing wave in three dimensions has u(r, t) = F(r - ct)/r, where F(s) = 1 for 0 < s < 1 and 0 otherwise. Describe what a detector at r = 10 records, and give the amplitude.

Solution. The detector reads u(10, t) = F(10 - ct)/10, which is nonzero exactly when 0 < 10 - ct < 1, that is, when t lies between 9/c and 10/c. So it records a pulse of duration 1/c starting at time 9/c, of constant amplitude 1/10, and then silence. The pulse has the same shape and duration as at any other radius; only the amplitude has fallen, by the factor 1/r.

Change one input. Put the same initial disturbance on a two-dimensional membrane. Now Poisson's formula integrates over a disc, so once the expanding disc of radius ct swallows the source region it keeps it forever. The detector records an arrival at the same moment and then a tail decaying roughly like 1/t rather than dropping to zero. The onset is identical; the ending is not.

Pulling it together

In three dimensions the substitution v = ru turns the spherically symmetric wave equation into the one-dimensional one, so spherical pulses keep their shape and decay like 1/r in amplitude, giving the inverse square law for intensity. Kirchhoff's formula uses data only on the sphere of radius ct; Poisson's two-dimensional formula, obtained by descent, uses data throughout the disc. That contrast is the strong Huygens principle, valid in odd dimensions three and above, and it is why a clap ends cleanly in air and a tap on a membrane does not.

Looking ahead

The course has now solved the three model equations under the boundary and initial conditions that suit them, and has twice run into problems where the natural-looking question turns out to be unanswerable: the backward heat equation in Lesson 8, and the overdetermined elliptic problem mentioned in Lesson 11. The next lesson gives the framework that names what those failures have in common, states Hadamard's three conditions, and works his classic ill-posed example in full.

Sources

1. Wikipedia contributors. (n.d.). Wave equation. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Huygens-Fresnel principle. Wikipedia. en.wikipedia.org
3. Wikipedia contributors. (n.d.). Kirchhoff integral theorem. Wikipedia. en.wikipedia.org
4. Weisstein, E. W. (n.d.). Wave equation. Wolfram MathWorld. mathworld.wolfram.com
5. Evans, L. C. (2010). Partial differential equations (2nd ed.), Chapter 2, Section 4. American Mathematical Society.

In 1902 Jacques Hadamard proposed that a mathematical problem modelling something physical ought to satisfy three conditions: a solution should exist, it should be unique, and it should depend continuously on the data. The third one is the interesting one. It says that if you measure your inputs slightly wrong, the answer is only slightly wrong, and without it a problem may be perfectly well determined on paper and completely useless in practice.

Hadamard produced an example to show that the condition can fail spectacularly for an equation everyone regarded as the tamest in the subject. That example is worth working through in full, because it is short, it is explicit, and the numbers are startling.

The three conditions

A problem is well posed in the sense of Hadamard if

• Existence: a solution exists for every admissible choice of data;
• Uniqueness: there is at most one;
• Continuous dependence: small changes in the data produce small changes in the solution, in an appropriate measure of size.

The three are logically independent, and the course has proved instances of each separately. Uniqueness came from an energy integral for the wave and heat equations and from the maximum principle for Laplace's equation. Existence came from explicit formulas: d'Alembert on the line, the heat kernel, the Poisson integral. Continuous dependence came from the maximum principle for the Dirichlet problem, with the sharp constant 1, and from the decaying exponential factors for the forward heat equation.

Hadamard's example, worked

Consider Laplace's equation in the upper half plane,

u_xx + u_yy = 0 for y > 0,

with Cauchy data on the line y = 0: both the value and the normal derivative are prescribed, exactly as one would do for the wave equation. Take

u(x, 0) = 0 and u_y(x, 0) = (1/n) sin(nx),

for a positive integer n. This data is as small as you like: the first function is zero, and the second has maximum size 1/n, which goes to zero as n grows. Compare it with the data u(x, 0) = 0 and u_y(x, 0) = 0, whose solution is u identically zero.

The solution with the perturbed data is

u(x, y) = (1/n^2) sin(nx) sinh(ny).

Verify all three requirements. First, u_xx = -(n^2/n^2) sin(nx) sinh(ny) = -sin(nx) sinh(ny), while u_yy = (n^2/n^2) sin(nx) sinh(ny) = sin(nx) sinh(ny), so their sum is zero and Laplace's equation holds. Second, at y = 0 the hyperbolic sine vanishes, so u(x, 0) = 0. Third, u_y = (n/n^2) sin(nx) cosh(ny), which at y = 0 equals (1/n) sin(nx), matching the prescribed derivative.

Now measure the two solutions apart at a fixed height. At y = 1 the maximum of |u| is sinh(n)/n^2. Put in numbers.

The data shrinks toward zero while the solution one unit away grows without bound. There is no constant C for which the solution size is bounded by C times the data size, and therefore no continuous dependence, and therefore the Cauchy problem for Laplace's equation is ill posed. The point: the failure is not a defect of any solution method. It is a property of the equation together with the way the data was posed.

The classification of Lesson 2 predicted this. Laplace's equation is elliptic, it has no real characteristics, and it does not support propagation from a line of initial data. Prescribing Cauchy data on a curve is the natural thing to do for a hyperbolic equation and is a category error for an elliptic one.

Existence and uniqueness are not enough

It might be hoped that ill-posedness signals some deeper failure, such as non-existence. It does not, and the cleanest way to see this is the Cauchy-Kovalevskaya theorem, proved in the general form by Sofia Kovalevskaya in 1875. It says that a Cauchy problem with analytic coefficients and analytic data has a unique analytic solution in a neighbourhood of the initial surface.

Laplace's equation has analytic coefficients, and Hadamard's data (1/n) sin(nx) is analytic. So the theorem applies: a solution exists and is unique, and indeed we wrote it down. Every clause of the theorem holds, and the problem is still useless, because the map from data to solution is unbounded. Continuous dependence is an independent requirement, and no existence theorem supplies it.

The same lesson appeared in Lesson 8. The backward heat equation has a unique solution whenever it has one at all, since running it forward again recovers the data, and the coefficients are recovered by multiplying by e^(+k n^2 T). Uniqueness is not the issue; unbounded amplification is.

Ill-posed problems that people solve anyway

Ill-posed does not mean unimportant. Some of the most valuable problems in applied mathematics are ill posed, and they are all of the same shape: recover a cause from an observed effect.

• Deblurring an image. Blurring is a diffusion-like smoothing, so deblurring is a backward diffusion and amplifies high frequencies, which is why naive deblurring produces noise.
• Computed tomography. Reconstructing a density from line integrals is unstable to measurement noise at fine scales.
• Geophysical inversion. Inferring subsurface structure from surface measurements is a Cauchy problem for an elliptic equation, structurally identical to Hadamard's example.

The standard remedy is regularisation: replace the exact problem with a nearby one that is well posed, for example by adding a penalty term that discourages wild high-frequency content, and accept a slightly wrong answer that is stable instead of an exactly right answer that is unobtainable. Tikhonov regularisation is the classic form. Choosing how much to regularise is a genuine trade-off between fidelity and stability, and it is decided by how noisy the data is.

Common misconceptions

• "Ill posed means no solution exists." It usually means the third condition fails. Hadamard's example has an explicit solution for every n.
• "Ill-posedness is a numerical problem." It is an analytical property of the problem itself. A perfect computer with infinite precision would still amplify a perturbation in the data by sinh(n)/n^2.
• "Continuous dependence follows from uniqueness." It does not. Both ill-posed problems above have unique solutions, and both fail continuity.
• "Ill-posed problems should be abandoned." Medical imaging, seismology and image restoration are all ill posed. They are handled by regularisation, which trades exactness for stability rather than giving up.

Try it: which condition fails

Exercise. Consider the problem u_xx + u_yy = 0 on the unit disc with u = g on the boundary and, in addition, u_r = h on the boundary. For which pairs (g, h) does a solution exist, and what does that say about the problem?

Solution. Given g alone, the Dirichlet problem already has a unique solution, by Lesson 10, and that solution has a definite normal derivative on the boundary, say h_g. So a solution of the augmented problem exists only if h happens to equal h_g. For every other h, and that is almost all of them, no solution exists at all. The problem fails existence, because it is overdetermined: the elliptic equation admits one condition per boundary point, and two have been supplied.

Change one input. Drop the condition u = g and keep only u_r = h. Now the problem is pure Neumann, and Lesson 11 showed it fails differently: existence requires the compatibility condition that h integrate to zero around the boundary, and when a solution exists it is not unique, since any constant may be added. Three variants of one equation on one disc, failing three different Hadamard conditions.

The short version

Hadamard's three conditions are existence, uniqueness and continuous dependence, and the third is the one that decides whether a problem can be used. The Cauchy problem for Laplace's equation satisfies the first two and fails the third catastrophically: data of size 1/n produces a solution of size sinh(n)/n^2 one unit away. The Cauchy-Kovalevskaya theorem shows that existence and uniqueness are genuinely insufficient. Ill-posed problems are common and important, and the practical response is regularisation rather than abandonment.

Looking ahead

The final lesson takes the same concern into computation. A finite difference scheme replaces derivatives by differences and marches forward step by step, and the question of whether small errors grow or shrink from step to step turns out to be the numerical analogue of continuous dependence. For the explicit scheme for the heat equation there is a sharp threshold, expressed as a relation between the time step and the square of the space step, and crossing it destroys the computation completely.

Sources

1. Wikipedia contributors. (n.d.). Well-posed problem. Wikipedia. en.wikipedia.org
2. O'Connor, J. J., & Robertson, E. F. (2003). Jacques Salomon Hadamard. MacTutor History of Mathematics Archive. mathshistory.st-andrews.ac.uk
3. Wikipedia contributors. (n.d.). Cauchy-Kovalevskaya theorem. Wikipedia. en.wikipedia.org
4. Massachusetts Institute of Technology. (2011). 18.152 Introduction to Partial Differential Equations [Course materials]. MIT OpenCourseWare. ocw.mit.edu
5. Evans, L. C. (2010). Partial differential equations (2nd ed.), Chapter 2. American Mathematical Society.

Here is a numerical experiment you can run on paper. Discretise the heat equation u_t = u_xx on a grid, start with every value zero except a single 1 at the origin, and step forward using the obvious formula with mesh ratio 1. The centre values you get are

1, -1, 3, -7, ...

alternating in sign and roughly doubling in size each step. The true solution of the heat equation is positive everywhere, decreasing at the origin, and utterly smooth. The scheme is not approximating it badly; it is approximating something else entirely, and after twenty steps the numbers are meaningless.

Now change one number, the ratio of the time step to the square of the space step, from 1 to 0.4, and the same scheme produces a smooth spreading bump that matches the exact solution to within the truncation error. This lesson explains where the threshold sits and why.

Setting up the grid

Choose a space step h and a time step tau, and write U_j^m for the approximation to u(jh, m tau). Taylor expansion gives the standard approximations

u_t is approximately (U_j^(m+1) - U_j^m)/tau, with error of order tau,

u_xx is approximately (U_(j+1)^m - 2U_j^m + U_(j-1)^m)/h^2, with error of order h^2.

The second one is worth checking: expanding u(x plus or minus h) to fourth order and adding, the odd terms cancel and the result is 2u + h^2 u_xx + h^4 u_xxxx/12, so dividing by h^2 gives u_xx plus an error proportional to h^2. A scheme built from these is called consistent: as the steps shrink, it approximates the differential equation.

The explicit scheme

Substituting both approximations into u_t = k u_xx and solving for the new value gives the forward-time centred-space scheme,

U_j^(m+1) = U_j^m + r [U_(j+1)^m - 2U_j^m + U_(j-1)^m], where r = k tau/h^2.

The number r is the mesh ratio, and it is the only parameter that matters. Each new value is computed directly from three old ones, so the scheme is explicit: no equations need solving.

The instability, by hand. Set k = 1 and r = 1, with initial data 1 at j = 0 and zero elsewhere. Then

• Step 1: U_0 = 1 + (0 - 2 + 0) = -1, and U_(±1) = 0 + (0 - 0 + 1) = 1. Profile: 1, -1, 1.
• Step 2: U_0 = -1 + (1 + 2 + 1) = 3, U_(±1) = 1 + (0 - 2 - 1) = -2, U_(±2) = 1. Profile: 1, -2, 3, -2, 1.
• Step 3: U_0 = 3 + (-2 - 6 - 2) = -7, U_(±1) = -2 + (1 + 4 + 3) = 6, U_(±2) = 1 + (0 - 2 - 2) = -3. Profile: 1, -3, 6, -7, 6, -3, 1.

The peak values 1, 1, 3, 7 grow while the sign alternates from cell to cell. That sawtooth is the signature of this instability, and it is the reason a physically sensible computation can turn into noise within a dozen steps.

von Neumann stability analysis

The general test asks how the scheme treats a single Fourier mode. Substitute

U_j^m = G^m e^(i theta j),

where theta ranges over [-pi, pi] and G is the amplification factor for that mode. If |G| > 1 for any theta, that mode grows geometrically and the scheme is unstable.

Substituting into the explicit scheme and cancelling the common factor G^m e^(i theta j),

G = 1 + r[e^(i theta) - 2 + e^(-i theta)] = 1 + r[2cos(theta) - 2].

Using 2 - 2cos(theta) = 4 sin^2(theta/2), this is

G(theta) = 1 - 4r sin^2(theta/2).

Now impose |G| <= 1 for every theta. Since sin^2(theta/2) runs over [0, 1], the quantity G runs from 1 down to 1 - 4r. The upper end is fine. The lower end requires 1 - 4r >= -1, that is,

r <= 1/2.

That is the stability condition. Check it against the hand computation: with r = 1 the worst mode is theta = pi, giving G = 1 - 4 = -3. Each step multiplies that mode by -3, which is precisely the alternating, roughly tripling behaviour seen in the peak values. Key idea: instability is not accumulated rounding error. It is one Fourier mode being multiplied by a number larger than 1 in magnitude, over and over.

Written out, r = k tau/h^2 <= 1/2 means tau <= h^2/(2k). Halving the space step forces the time step down by a factor of four. Refining a heat computation is therefore expensive, and that expense is what motivates the alternatives.

Implicit schemes

Evaluate the space difference at the new time level instead:

U_j^(m+1) - U_j^m = r[U_(j+1)^(m+1) - 2U_j^(m+1) + U_(j-1)^(m+1)].

Now the unknowns appear on both sides, so each step requires solving a linear system. The system is tridiagonal, and tridiagonal systems are solved in time proportional to the number of unknowns, so the cost per step is modest. The payoff is in the amplification factor: repeating the analysis gives

G = 1/(1 + 4r sin^2(theta/2)),

which lies in (0, 1] for every r > 0. The scheme is unconditionally stable: any time step is permitted. It remains only first order accurate in time.

Averaging the explicit and implicit right-hand sides gives the Crank-Nicolson scheme, whose amplification factor is

G = (1 - 2r sin^2(theta/2))/(1 + 2r sin^2(theta/2)),

a ratio whose numerator is never larger in magnitude than its denominator, so again |G| <= 1 for all r. Crank-Nicolson is unconditionally stable and second order accurate in both time and space, which is why it is the default choice for parabolic problems.

The CFL condition for waves

For the wave equation u_tt = c^2 u_xx, the natural explicit scheme replaces both second derivatives by central differences, and von Neumann analysis gives the stability requirement

c tau/h <= 1.

Richard Courant, Kurt Friedrichs and Hans Lewy described this condition in a 1928 paper, and it carries their initials. Its interpretation is geometric and memorable. The exact solution at a grid point depends on the initial data throughout the interval [x - ct, x + ct], the domain of dependence of Lesson 4. The numerical solution depends only on the grid points reachable by the stencil, which after m steps is the interval [x - mh, x + mh]. If the true domain of dependence sticks out beyond the numerical one, the scheme is computing without access to data that the answer requires, and no amount of accuracy can rescue it. Requiring ct <= mh with t = m tau is exactly c tau <= h.

Notice the difference in cost. The heat condition ties tau to h^2; the wave condition ties it to h. Refining a hyperbolic computation is far cheaper than refining a parabolic one, which is a direct consequence of the classification in Lesson 2.

Consistency, stability, convergence

The three ideas are tied together by the Lax equivalence theorem: for a consistent finite difference scheme applied to a well-posed linear initial value problem, stability is necessary and sufficient for convergence. Read the parts. Consistency says the scheme approximates the right equation as the steps shrink, and it is checked by Taylor expansion. Stability says errors do not grow without bound, and it is checked by von Neumann analysis. Convergence says the computed answer approaches the true solution, and it is what you actually want.

The theorem is why numerical analysis is organised the way it is: nobody proves convergence directly, because consistency and stability are each straightforward to check and together they suffice. Bottom line: a consistent scheme that is unstable does not converge slowly, it does not converge at all, and the hand computation at the start of this lesson is what that looks like.

Common misconceptions

• "Instability comes from rounding error." It comes from the scheme multiplying a Fourier mode by a factor larger than 1 in magnitude at every step. Exact arithmetic would still produce the alternating growth, seeded by the initial data itself.
• "A smaller time step always helps." For the explicit heat scheme it does, since r shrinks. For a scheme that is unstable for structural reasons, it does not, and for the explicit heat scheme the required smallness is proportional to h^2, which becomes brutal quickly.
• "Unconditionally stable means accurate." Backward Euler is unconditionally stable and only first order in time. Stability is a licence to take large steps, not a promise that the answer is good.
• "The CFL condition is a rule of thumb." It is a necessary condition with a proof: if the numerical domain of dependence does not contain the analytical one, the scheme cannot converge, because refining the grid never gives it access to the missing data.

Try it: find the largest safe time step

Exercise. You are solving u_t = 0.5 u_xx on [0, 1] with the explicit scheme and 100 intervals, so h = 0.01. What is the largest stable time step, and how many steps are needed to reach t = 1?

Solution. Stability requires r = k tau/h^2 <= 1/2, so tau <= h^2/(2k) = 0.0001/1 = 0.0001. Reaching t = 1 therefore takes at least 10000 steps, each touching 101 grid points, for about a million arithmetic operations. If you refine to 200 intervals, h halves, the maximum tau drops by four, and the total work goes up by a factor of eight.

Change one input. Solve the wave equation u_tt = 0.25 u_xx on the same grid, so c = 0.5. Now CFL requires tau <= h/c = 0.01/0.5 = 0.02, and reaching t = 1 takes only 50 steps. Same grid, the same target time, two hundred times fewer steps, entirely because one condition scales with h and the other with h^2.

What you now know

Replacing derivatives by differences gives a consistent scheme, and consistency is only half of what is needed. Substituting a single Fourier mode produces an amplification factor, and requiring its magnitude to stay at most 1 for every mode gives the stability condition: r = k tau/h^2 <= 1/2 for the explicit heat scheme, and c tau/h <= 1 for the explicit wave scheme, the CFL condition of Courant, Friedrichs and Lewy. Implicit and Crank-Nicolson schemes trade a tridiagonal solve for unconditional stability, with Crank-Nicolson also second order accurate. The Lax equivalence theorem says consistency plus stability is exactly convergence.

Where the course leaves you

Seventeen lessons ago the subject was the observation that u_x = 0 has a whole function's worth of solutions, so side conditions decide everything. From there: characteristics that carry information until they collide, d'Alembert's two travelling profiles, Fourier's series and the exact size of its overshoot at a jump, the exponential factor that makes diffusion smooth and irreversible, the averaging property that traps a harmonic function between its boundary values, Sturm-Liouville theory explaining why every expansion in this subject works, Green's functions and transforms, Huygens' principle in three dimensions, Hadamard's three conditions, and finally the mesh ratio that decides whether a computation means anything. The three model equations were the spine throughout, and the classification of Lesson 2 predicted the behaviour of each one before any of it was solved.

Sources

1. Wikipedia contributors. (n.d.). FTCS scheme. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Von Neumann stability analysis. Wikipedia. en.wikipedia.org
3. Wikipedia contributors. (n.d.). Courant-Friedrichs-Lewy condition. Wikipedia. en.wikipedia.org
4. Wikipedia contributors. (n.d.). Crank-Nicolson method. Wikipedia. en.wikipedia.org
5. Massachusetts Institute of Technology. (2014). 18.303 Linear Partial Differential Equations: Analysis and Numerics [Course materials]. MIT OpenCourseWare. ocw.mit.edu