The Heat Equation and the Diffusion Model

How does heat spread through a metal rod? Can we write an equation predicting temperature everywhere, for all future time, from the temperature right now? Joseph Fourier posed and answered this question in the early nineteenth century. His work produced one of the most important partial differential equations in mathematics and physics. The heat equation is the prototype of the parabolic class in the Classification of Second-Order Linear PDEs. It models not only literal heat flow but also chemical diffusion and Brownian motion. After transformation, it even applies to certain problems in finance. The big picture is this: physical conservation laws, combined with an empirical constitutive law relating flux to gradient, force a specific differential equation on us. We do not guess the heat equation — we derive it. This article carries out that derivation with full rigor. It states the precise theorem obtained, and proves it from first principles. It also previews the machinery (separation of variables, maximum principles) built on top of it.

Key takeaways
  • The heat equation arises from conservation of energy plus Fourier's law of heat conduction.
  • Its precise form is ut=α2Δuu_t = \alpha^2 \Delta u, a linear, second-order, parabolic PDE.
  • The derivation uses an arbitrary control volume and the divergence theorem — no symmetry is assumed.
  • Fourier's law q=−k∇u\mathbf{q} = -k\nabla u is an empirical constitutive assumption, not a logical necessity.
  • The same equation models diffusion of mass, with concentration replacing temperature.
  • Key idea: local conservation of a quantity, transported down its own gradient, always yields a diffusion-type PDE.

Where this lives today

The heat equation is far more than a model of literal heat flow. In probability theory, it is the forward Kolmogorov (Fokker–Planck) equation governing the transition density of Brownian motion, linking PDE theory directly to stochastic processes. In mathematical finance, a change of variables transforms the Black–Scholes equation for option pricing into the heat equation, making heat-kernel solutions directly usable for pricing formulas. In image processing and computer vision, applying the heat equation to an image (Gaussian blurring) is the basis of scale-space analysis and some denoising algorithms; anisotropic variants (Perona–Malik) modify Fourier's law itself to preserve edges. In differential geometry, the heat equation on manifolds underlies heat kernel proofs of the Atiyah–Singer index theorem, and Ricci flow (used in the proof of the Poincaré conjecture) is a nonlinear geometric analogue of the diffusion equation applied to a metric tensor. In numerical analysis and engineering, finite-difference and finite-element methods for the heat equation are the standard testbed for stability theory (e.g. the CFL-type conditions for parabolic problems), used in everything from semiconductor thermal design to climate modeling of heat diffusion in oceans and ice sheets.

Ask the tutor: How does separation of variables use the eigenfunctions of the Laplacian to actually solve this equation on a finite rod?

DefinitionThe heat/diffusion equation and its physical data

Let Ω⊆Rn\Omega \subseteq \mathbb{R}^n (n=1,2,3n = 1, 2, 3 in applications) be an open, connected region occupied by a conducting medium, and let T>0T > 0. Let

u:Ω×[0,T]→Ru : \Omega \times [0,T] \to \mathbb{R}

denote the temperature at position x∈Ωx \in \Omega and time t∈[0,T]t \in [0,T]. The medium is described by three material parameters, each a positive function of xx (constant in the homogeneous, isotropic case treated here):

  • ρ(x)>0\rho(x) > 0 — mass density,
  • c(x)>0c(x) > 0 — specific heat capacity,
  • k(x)>0k(x) > 0 — thermal conductivity.

The heat equation is the partial differential equation

ρ(x) c(x) ∂u∂t(x,t)  =  ∇⋅(k(x) ∇u(x,t))+f(x,t),\rho(x)\, c(x)\, \frac{\partial u}{\partial t}(x,t) \;=\; \nabla \cdot \big(k(x)\, \nabla u(x,t)\big) + f(x,t),

where f(x,t)f(x,t) is a volumetric heat source density (energy per unit volume per unit time). When ρ,c,k\rho, c, k are constants and f≡0f \equiv 0, this reduces to the standard heat equation

ut=α2Δu,α2:=kρc>0,u_t = \alpha^2 \Delta u, \qquad \alpha^2 := \frac{k}{\rho c} > 0,

where α2\alpha^2 is the thermal diffusivity and Δ=∑i=1n∂xi2\Delta = \sum_{i=1}^n \partial_{x_i}^2 is the Laplacian. The same equation, with uu reinterpreted as a chemical concentration and α2\alpha^2 as a diffusion constant DD, is called the diffusion equation.

Intuition

Picture a thin metal rod insulated along its sides so heat can only flow along its length. If one end is hotter than the other, common experience says heat flows from hot to cold, and it flows faster when the temperature difference (per unit length) is steeper. That is the essence of Fourier's law: the flow of heat, called the heat flux, is proportional to the negative gradient of temperature. No flux exists where the temperature is level; a kinked, sharply-varying temperature profile drives a strong flux. The second ingredient is bookkeeping: heat energy is neither created nor destroyed except by explicit sources, so any change in the thermal energy inside a region must equal the net flux of heat crossing its boundary, plus whatever a source contributes. This is exactly the same style of reasoning used to derive conservation laws in fluid dynamics and electromagnetism (continuity equations). When you combine "energy is conserved locally" with "flux follows the negative gradient," calculus (specifically the divergence theorem) automatically produces a second-order PDE relating one time derivative to two spatial derivatives — the parabolic signature. The picture below shows a one-dimensional rod, its temperature profile at a fixed time, and the heat flux at two points computed from the local slope. Steeper downward slope (bigger ∣ux∣|u_x|) means larger rightward flux.

Figure
Fourier's law on a rod: flux is proportional to the negative slope of u(x)u(x)
x u(x) hot end cool end large flux q small flux q q(x,t) = −k u_x(x,t)

Central result

TheoremDerivation of the heat equation from conservation of energy and Fourier's law

Let Ω⊆Rn\Omega \subseteq \mathbb{R}^n be open and connected, let u∈C2(Ω×(0,T))u \in C^2(\Omega \times (0,T)) be the temperature field, and assume: 1. (Conservation of energy) For every bounded subregion V⊂ΩV \subset \Omega with piecewise-smooth boundary ∂V\partial V and outward unit normal n\mathbf{n}, the rate of change of thermal energy stored in VV equals the net rate at which heat enters across ∂V\partial V, plus energy generated inside VV by a source of density f(x,t)f(x,t):

ddt∫Vρ(x)c(x) u(x,t) dV  =  −∫∂Vq(x,t)⋅n dS  +  ∫Vf(x,t) dV,\frac{d}{dt}\int_V \rho(x) c(x)\, u(x,t)\, dV \;=\; -\int_{\partial V} \mathbf{q}(x,t)\cdot \mathbf{n}\, dS \;+\; \int_V f(x,t)\, dV,

where q\mathbf{q} is the heat flux vector (energy per unit area per unit time). 2. (Fourier's law) The heat flux is proportional to the negative temperature gradient:

q(x,t)=−k(x) ∇u(x,t),\mathbf{q}(x,t) = -k(x)\,\nabla u(x,t),

for a scalar conductivity function k(x)>0k(x) > 0. 3. (Regularity) ρ,c,k∈C1(Ω)\rho, c, k \in C^1(\Omega) are strictly positive, ff is continuous on Ω×(0,T)\Omega \times (0,T), and q∈C1\mathbf{q} \in C^1 so that the divergence theorem applies on every admissible VV. Then uu satisfies, at every point of Ω×(0,T)\Omega \times (0,T), the partial differential equation

ρ(x)c(x) ut(x,t)=∇⋅(k(x)∇u(x,t))+f(x,t).\rho(x) c(x)\, u_t(x,t) = \nabla\cdot\big(k(x)\nabla u(x,t)\big) + f(x,t).

In particular, if ρ,c,k\rho, c, k are constants and f≡0f \equiv 0, then

ut=α2Δu,α2=kρc.u_t = \alpha^2 \Delta u, \qquad \alpha^2 = \frac{k}{\rho c}.
Proof

We use the localization technique: derive an integral identity valid on every subregion, convert the boundary integral to a volume integral via the divergence theorem, and then argue that an integral vanishing on every subregion forces the integrand itself to vanish.

Step 1 (Convert the boundary flux integral to a volume integral). Fix t∈(0,T)t \in (0,T) and let V⊂ΩV \subset \Omega be any bounded region with piecewise-smooth boundary. By hypothesis 3, q(⋅,t)∈C1(V‾)\mathbf{q}(\cdot, t) \in C^1(\overline{V}), so the divergence theorem applies:

∫∂Vq(x,t)⋅n(x) dS=∫V∇⋅q(x,t) dV.\int_{\partial V} \mathbf{q}(x,t)\cdot \mathbf{n}(x)\, dS = \int_V \nabla \cdot \mathbf{q}(x,t)\, dV.

Substituting Fourier's law q=−k∇u\mathbf{q} = -k\nabla u from hypothesis 2,

∫∂Vq⋅n dS=−∫V∇⋅(k(x)∇u(x,t)) dV.\int_{\partial V} \mathbf{q}\cdot \mathbf{n}\, dS = -\int_V \nabla\cdot\big(k(x)\nabla u(x,t)\big)\, dV.

Step 2 (Differentiate under the integral on the left side of conservation). Since u∈C2(Ω×(0,T))u \in C^2(\Omega \times (0,T)) and ρ,c\rho, c do not depend on tt, and VV is a fixed bounded region independent of tt, the Leibniz rule for differentiating under the integral sign (valid because ρc ut\rho c\, u_t is continuous on the compact set V‾×[t1,t2]\overline V \times [t_1,t_2] for any [t1,t2]⊂(0,T)[t_1,t_2]\subset(0,T), hence bounded there) gives

ddt∫Vρ(x)c(x) u(x,t) dV=∫Vρ(x)c(x) ut(x,t) dV.\frac{d}{dt}\int_V \rho(x) c(x)\, u(x,t)\, dV = \int_V \rho(x) c(x)\, u_t(x,t)\, dV.

Step 3 (Assemble the integral balance). Substitute Steps 1 and 2 into the conservation law of hypothesis 1:

∫Vρc ut dV=∫V∇⋅(k∇u) dV+∫Vf dV.\int_V \rho c\, u_t \, dV = \int_V \nabla\cdot(k\nabla u)\, dV + \int_V f\, dV.

Since all three integrals are over the same region VV, we may combine them:

∫V[ρ(x)c(x) ut(x,t)−∇⋅(k(x)∇u(x,t))−f(x,t)]dV=0.(∗)\int_V \Big[\rho(x) c(x)\, u_t(x,t) - \nabla\cdot\big(k(x)\nabla u(x,t)\big) - f(x,t)\Big] dV = 0. \tag{$\ast$}

Step 4 (Localization: the integrand must vanish pointwise). Define

g(x,t):=ρ(x)c(x) ut(x,t)−∇⋅(k(x)∇u(x,t))−f(x,t).g(x,t) := \rho(x) c(x)\, u_t(x,t) - \nabla\cdot(k(x)\nabla u(x,t)) - f(x,t).

By hypothesis 3 and u∈C2u \in C^2, g(⋅,t)g(\cdot, t) is continuous on Ω\Omega. Equation (∗)(\ast) states ∫Vg(x,t) dV=0\int_V g(x,t)\, dV = 0 for every admissible bounded subregion V⊂ΩV \subset \Omega, in particular for every open ball Br(x0)⊂ΩB_r(x_0) \subset \Omega. Suppose, for contradiction, that g(x0,t)≠0g(x_0, t) \neq 0 for some x0∈Ωx_0 \in \Omega; without loss of generality g(x0,t)=2ε>0g(x_0,t) = 2\varepsilon > 0. By continuity of g(⋅,t)g(\cdot,t) at x0x_0, there exists r>0r > 0 with Br(x0)⊂ΩB_r(x_0) \subset \Omega such that g(x,t)>εg(x,t) > \varepsilon for all x∈Br(x0)x \in B_r(x_0). Then

∫Br(x0)g(x,t) dV≥ε⋅vol(Br(x0))>0,\int_{B_r(x_0)} g(x,t)\, dV \geq \varepsilon \cdot \mathrm{vol}(B_r(x_0)) > 0,

contradicting (∗)(\ast) applied to V=Br(x0)V = B_r(x_0). The case g(x0,t)<0g(x_0,t) < 0 is symmetric. Hence g(x0,t)=0g(x_0,t) = 0 for every x0∈Ωx_0 \in \Omega, and since t∈(0,T)t \in (0,T) was arbitrary,

ρ(x)c(x) ut(x,t)=∇⋅(k(x)∇u(x,t))+f(x,t)for all (x,t)∈Ω×(0,T).\rho(x) c(x)\, u_t(x,t) = \nabla\cdot\big(k(x)\nabla u(x,t)\big) + f(x,t) \quad \text{for all } (x,t) \in \Omega\times(0,T).

Step 5 (Constant-coefficient case). If ρ,c,k\rho, c, k are constants and f≡0f \equiv 0, then ∇⋅(k∇u)=k∇⋅∇u=kΔu\nabla \cdot (k\nabla u) = k\nabla\cdot\nabla u = k\Delta u, so the equation becomes ρc ut=kΔu\rho c\, u_t = k\Delta u, i.e.

ut=kρcΔu=α2Δu,u_t = \frac{k}{\rho c}\Delta u = \alpha^2 \Delta u,

with α2=k/(ρc)>0\alpha^2 = k/(\rho c) > 0 since k,ρ,c>0k, \rho, c > 0 by hypothesis. ■\blacksquare

∎
CorollaryDiffusion equation for mass concentration

Let u(x,t)u(x,t) denote the concentration (mass per unit volume) of a diffusing chemical species in a stationary medium, with no chemical reactions occurring (f≡0f \equiv 0). Suppose mass is conserved and the mass flux obeys Fick's law J=−D∇u\mathbf{J} = -D\nabla u for a constant diffusion coefficient D>0D > 0. Then the identical argument (with ρc≡1\rho c \equiv 1, since uu is already a density) yields

ut=D Δu.u_t = D\, \Delta u.
Proof

Conservation of mass for a chemical species states that the rate of change of total mass in any region VV equals the negative net outward flux (no sources, no sinks):

ddt∫Vu dV=−∫∂VJ⋅n dS.\frac{d}{dt}\int_V u\, dV = -\int_{\partial V}\mathbf{J}\cdot\mathbf n\, dS.

This is precisely hypothesis 1 of the theorem with ρc≡1\rho c \equiv 1 and f≡0f \equiv 0. Fick's law J=−D∇u\mathbf J = -D\nabla u is precisely hypothesis 2 with k≡Dk \equiv D. All steps of the proof above go through verbatim, yielding ut=DΔuu_t = D\Delta u. ■\blacksquare

∎

Worked examples

Example 1 (Non-uniform rod, direct verification of the PDE from a candidate solution). Consider a rod occupying Ω=(0,π)\Omega = (0, \pi) modeling the constant-coefficient case with α2=3\alpha^2 = 3 (chosen asymmetrically, not the "nice" value 11). Suppose we propose the candidate temperature profile

u(x,t)=5 e−3⋅4tsin⁡(2x)  +  2 e−3⋅9tsin⁡(3x),u(x,t) = 5\, e^{-3\cdot 4 t}\sin(2x) \;+\; 2\, e^{-3\cdot 9 t}\sin(3x),

combining two spatial modes with unequal amplitudes 55 and 22 and unequal wavenumbers 2,32, 3. We verify directly that this uu satisfies ut=3uxxu_t = 3 u_{xx}, confirming the theorem's equation is the correct one to solve (this also previews Separation of Variables for the Heat Equation, where such sums of modes are exactly how solutions are built). Compute the time derivative term by term. For the first term, ∂t(5e−12tsin⁡2x)=5(−12)e−12tsin⁡2x=−60 e−12tsin⁡2x\partial_t \big(5 e^{-12t}\sin 2x\big) = 5(-12) e^{-12t}\sin 2x = -60\, e^{-12t}\sin 2x. For the second, ∂t(2e−27tsin⁡3x)=2(−27)e−27tsin⁡3x=−54 e−27tsin⁡3x\partial_t\big(2 e^{-27t}\sin 3x\big) = 2(-27)e^{-27t}\sin 3x = -54\, e^{-27t}\sin 3x. So

ut=−60 e−12tsin⁡2x−54 e−27tsin⁡3x.u_t = -60\, e^{-12t}\sin 2x - 54\, e^{-27t}\sin 3x.

Now compute uxxu_{xx}. For the first term, ∂x2(5e−12tsin⁡2x)=5e−12t⋅(−4)sin⁡2x=−20 e−12tsin⁡2x\partial_x^2\big(5e^{-12t}\sin 2x\big) = 5e^{-12t}\cdot(-4)\sin 2x = -20\, e^{-12t}\sin 2x. For the second, ∂x2(2e−27tsin⁡3x)=2e−27t⋅(−9)sin⁡3x=−18 e−27tsin⁡3x\partial_x^2\big(2e^{-27t}\sin 3x\big) = 2e^{-27t}\cdot(-9)\sin 3x = -18\, e^{-27t}\sin 3x. So

uxx=−20 e−12tsin⁡2x−18 e−27tsin⁡3x.u_{xx} = -20\, e^{-12t}\sin 2x - 18\, e^{-27t}\sin 3x.

Multiply by α2=3\alpha^2 = 3:

3uxx=−60 e−12tsin⁡2x−54 e−27tsin⁡3x.3 u_{xx} = -60\, e^{-12t}\sin 2x - 54\, e^{-27t}\sin 3x.

This matches utu_t exactly, term by term. Note the exponents −3⋅4=−12-3\cdot 4 = -12 and −3⋅9=−27-3\cdot 9 = -27 are precisely −α2n2-\alpha^2 n^2 for wavenumbers n=2,3n=2,3; this is the signature pattern that separation of variables will later derive systematically, and here it is confirmed to satisfy the derived equation, not merely assumed.

Example 2 (A rod with a heat source and variable conductivity — using the general PDE from the theorem). Consider a rod along Ω=(0,L)\Omega = (0, L) with L=2L = 2, non-constant conductivity k(x)=1+xk(x) = 1 + x (the rod is a composite material, denser/more conductive toward x=Lx = L), constant ρ=2\rho = 2, c=1c = 1, and an internal heat source

f(x,t)=6x e−t.f(x,t) = 6x\, e^{-t}.

We are given the candidate temperature field

u(x,t)=x2 e−t,u(x,t) = x^2\, e^{-t},

and we check it against the general balance law of the theorem,

ρc ut=∂x(k(x)ux)+f(x,t).\rho c\, u_t = \partial_x\big(k(x) u_x\big) + f(x,t).

Left side. ut=−x2e−tu_t = -x^2 e^{-t}, so ρc ut=2⋅1⋅(−x2e−t)=−2x2e−t\rho c\, u_t = 2\cdot 1\cdot(-x^2 e^{-t}) = -2x^2 e^{-t}. Right side, flux term. First ux=2xe−tu_x = 2x e^{-t}. Then

k(x)ux=(1+x)(2xe−t)=(2x+2x2)e−t.k(x) u_x = (1+x)(2x e^{-t}) = (2x + 2x^2) e^{-t}.

Differentiate in xx:

∂x(k(x)ux)=(2+4x) e−t.\partial_x\big(k(x)u_x\big) = (2 + 4x)\, e^{-t}.

Right side, total.

∂x(kux)+f=(2+4x)e−t+6xe−t=(2+10x) e−t.\partial_x(k u_x) + f = (2+4x)e^{-t} + 6x e^{-t} = (2 + 10x)\, e^{-t}.

Compare: left side is −2x2e−t-2x^2 e^{-t}, right side is (2+10x)e−t(2+10x)e^{-t}. These do not match (e.g. at x=1,t=0x=1, t=0: left =−2=-2, right =12=12). So u=x2e−tu = x^2 e^{-t} is not a solution of this particular variable-coefficient heat equation with this source — which is the correct and expected outcome, since we invented uu and ff independently to illustrate the mechanics of plugging into the general divergence-form equation. This demonstrates the necessary discipline: one cannot simply guess uu and ff together and expect consistency; the theorem tells us exactly which combinations of (ρ,c,k,f,u)(\rho,c,k,f,u) are physically consistent, namely those satisfying the PDE identically, and checking this is a concrete, mechanical verification of the divergence-form balance ρc ut=(kux)x+f\rho c\,u_t = (k u_x)_x + f rather than the constant-coefficient shortcut ut=α2uxxu_t = \alpha^2 u_{xx}. To see the mechanism produce an actual solution, fix f≡0f \equiv 0 instead and ask what source-free steady state (ut=0u_t = 0) looks like with k(x)=1+xk(x) = 1+x: we need (kux)x=0(k u_x)_x = 0, i.e. (1+x)ux=C1(1+x)u_x = C_1 constant, so ux=C1/(1+x)u_x = C_1/(1+x), giving u(x)=C1ln⁡(1+x)+C2u(x) = C_1 \ln(1+x) + C_2. For instance, with boundary temperatures u(0)=0u(0) = 0 and u(2)=10u(2) = 10, we solve C2=0C_2 = 0 and C1ln⁡3=10C_1 \ln 3 = 10, so C1=10/ln⁡3≈9.102C_1 = 10/\ln 3 \approx 9.102, giving the genuine steady-state solution

u(x)=10ln⁡3ln⁡(1+x),u(x) = \frac{10}{\ln 3}\ln(1+x),

which one can check directly satisfies (1+x)ux=10/ln⁡3=(1+x)u_x = 10/\ln 3 = constant, hence (kux)x=0=ρc ut(k u_x)_x = 0 = \rho c\, u_t (as uu is time-independent). This is the variable-conductivity analogue of the familiar linear steady profile obtained when kk is constant.

CounterexampleWhy Fourier's law (hypothesis 2) cannot be dropped

Conservation of energy (hypothesis 1) alone, without a constitutive law relating q\mathbf q to uu, does not produce the heat equation — it produces only the continuity equation

ρc ut+∇⋅q=f,\rho c\, u_t + \nabla\cdot \mathbf q = f,

which involves the unknown flux q\mathbf q as an independent field, not expressible in terms of uu alone. To see concretely that the conclusion of the theorem fails without Fourier's law, consider Ω=(0,1)\Omega = (0,1), ρ=c=1\rho = c = 1, f≡0f \equiv 0, and define

u(x,t)=sin⁡(πx)cos⁡(t),q(x,t)=q(x,t)=πcos⁡(πx)sin⁡(t).u(x,t) = \sin(\pi x)\cos(t), \qquad \mathbf q(x,t) = q(x,t) = \pi \cos(\pi x)\sin(t).

Check conservation of energy directly in the pointwise (localized) form ut+qx=0u_t + q_x = 0:

ut=−sin⁡(πx)sin⁡(t),qx=−π2sin⁡(πx)sin⁡(t).u_t = -\sin(\pi x)\sin(t), \qquad q_x = -\pi^2\sin(\pi x)\sin(t).

These are not proportional with the constant −1-1 needed for ut+qx=0u_t + q_x = 0 unless π2=1\pi^2 = 1, which is false. So this particular (u,q)(u,q) pair does not even satisfy conservation as stated — showing that conservation alone is a genuine constraint linking uu and qq, but it is a single scalar equation relating two unknown functions uu and qq, hence underdetermined: infinitely many pairs (u,q)(u,q) satisfy ut+qx=0u_t + q_x = 0 pointwise (e.g. take any uu and set q(x,t)=−∫0xut(ξ,t) dξq(x,t) = -\int_0^x u_t(\xi,t)\,d\xi), and uu is not pinned down by conservation alone. Only after imposing a constitutive law such as Fourier's q=−kuxq = -k u_x does the system close into a well-posed PDE purely in uu. This is exactly why hypothesis 2 is listed as a separate, independent, empirical assumption in the theorem — it is not a consequence of conservation, but extra physical input (analogous to how Classification of Second-Order Linear PDEs separates the equation's algebraic type from the physical law that produced it).

Scope, edge cases & common traps

SituationWhat happensWhy
k,ρ,ck, \rho, c constant, f=0f=0Reduces to ut=α2Δuu_t = \alpha^2 \Delta uDivergence of a constant times gradient is constant times Laplacian
kk depends on xxMust keep divergence form ∇⋅(k∇u)\nabla\cdot(k\nabla u)∇⋅(k∇u)≠kΔu\nabla\cdot(k\nabla u) \neq k\Delta u in general; product rule adds ∇k⋅∇u\nabla k \cdot \nabla u
Anisotropic mediumkk becomes a matrix K(x)K(x), flux q=−K∇u\mathbf q = -K\nabla uHeat need not flow parallel to −∇u-\nabla u; same proof technique applies with KK symmetric positive-definite
No boundary/initial data givenEquation alone has infinitely many solutionsThe PDE is derived from local physics; global uniqueness needs The Maximum Principle for Parabolic and Elliptic Equations
α2<0\alpha^2 < 0 mistakenly usedEquation becomes backward-parabolic, ill-posedSign of α2=k/(ρc)\alpha^2 = k/(\rho c) is forced positive by physical positivity of k,ρ,ck,\rho,c

The mathematician behind it

Origins · 1 minJean-Baptiste Joseph Fourier and how this was discovered
JJ
Jean-Baptiste Joseph Fourier
1768–1830 · Napoleonic-era France; École Polytechnique and Grenoble

Fourier served as a mathematics instructor and later as a scientific administrator in Napoleon's Egyptian campaign before returning to France to serve as prefect of the Isère department. It was there, largely isolated from the Paris mathematical establishment, that he developed his theory of heat conduction, submitting his prize memoir to the French Academy of Sciences in 1807. His central innovation was twofold: the empirical flux law now bearing his name, and the audacious claim that arbitrary periodic functions could be represented as infinite sums of sines and cosines — the germinal idea of Fourier series, developed specifically to solve the heat equation on bounded domains.

History
- 1804–1807 — Fourier develops his theory of heat propagation and submits a memoir to the French Academy; the committee (including Lagrange and Laplace) is skeptical of the trigonometric series claims.
- 1822 — Fourier publishes Théorie analytique de la chaleur, presenting the heat equation, Fourier's law, and Fourier series in unified, polished form.
- 1828 — Poisson and others begin extending and rigorizing the analysis, addressing convergence questions Fourier had left informal.
- 1870s–1890s — Rigorous foundations of Fourier series (Dirichlet, Riemann, later Cantor's set theory motivated partly by convergence questions) solidify the analytic tools needed to justify Fourier's formal solutions.
- 20th century — The heat equation becomes the canonical model for parabolic PDEs generally, with maximum principles, semigroup theory, and probabilistic (Brownian motion) interpretations all built on Fourier's original derivation.
Why is the heat equation called "parabolic"?

The classification of second-order linear PDEs Auxx+2Buxt+Cutt+⋯=0A u_{xx} + 2B u_{xt} + C u_{tt} + \dots = 0 is by the sign of the discriminant B2−ACB^2 - AC, in analogy with conic sections. For ut=α2uxxu_t = \alpha^2 u_{xx}, treating tt as the second variable, the second-order terms are only uxxu_{xx} (coefficient α2\alpha^2) with no uttu_{tt} or mixed term, giving discriminant zero — the borderline case, exactly as a parabola is the borderline conic. This classification, covered in Classification of Second-Order Linear PDEs, predicts qualitative behavior: parabolic equations have one-directional (irreversible) time evolution and smoothing properties, unlike the two-directional wave equation (hyperbolic) or the equilibrium Laplace equation (elliptic).

Is Fourier's law exact, or an approximation?

Fourier's law is an empirical constitutive law, not a fundamental physical principle like conservation of energy. It is an excellent approximation for conduction in most solids and stationary fluids under normal conditions, but it is known to fail at extremely short time or length scales (e.g. it predicts infinite propagation speed for thermal disturbances, which is unphysical). More refined models like the Cattaneo equation introduce a relaxation time to fix this. For everyday engineering purposes at macroscopic scales, however, Fourier's law is highly accurate and is the standard starting point.

Does the heat equation apply to non-uniform materials?

Yes — this is exactly why the theorem is stated for spatially varying ρ(x),c(x),k(x)\rho(x), c(x), k(x) rather than only constants. The derivation goes through unchanged because conservation of energy and Fourier's law are stated locally, for an arbitrary control volume, and never assumed the material to be uniform. The only difference is that the equation retains the divergence form ∇⋅(k(x)∇u)\nabla\cdot(k(x)\nabla u) instead of simplifying to kΔuk\Delta u; the extra term ∇k⋅∇u\nabla k \cdot \nabla u appears whenever conductivity varies in space, as shown explicitly in the second worked example.

What is the difference between the heat equation and the diffusion equation?

Mathematically, none — both are ut=DΔuu_t = D\Delta u for a positive constant DD (thermal diffusivity or diffusion coefficient, respectively). The difference is purely in physical interpretation: uu is temperature governed by Fourier's law in one case, and concentration governed by Fick's law in the other. The corollary in this article shows the derivations are structurally identical, since both rest on the same mathematical skeleton — a conserved density transported down its own gradient — despite describing entirely different physical phenomena.

Why do we need boundary and initial conditions if the PDE is already derived?

The theorem derives a local differential relationship that any physically consistent temperature field must satisfy at each interior point and instant. But a PDE alone, like an ODE alone, admits infinitely many solutions (e.g. shifting or scaling initial data yields different valid solutions of the same equation). To pin down a unique, physically meaningful temperature history, one must additionally specify an initial temperature distribution u(x,0)u(x,0) and boundary behavior (fixed temperature, insulated ends, or radiative conditions). Establishing that such data indeed produces a unique, well-behaved solution is the role of results like The Maximum Principle for Parabolic and Elliptic Equations.

How does separation of variables fit in with this derivation?

This article establishes which equation governs heat flow; it does not yet solve it. Separation of Variables for the Heat Equation is the standard technique for actually solving ut=α2uxxu_t = \alpha^2 u_{xx} on a bounded domain with given boundary and initial conditions, by seeking solutions of the product form u(x,t)=X(x)T(t)u(x,t) = X(x)T(t) and expanding the initial data in a Fourier series of eigenfunctions of the spatial operator. The first worked example above previews this: the solution was built from exponentially decaying sine modes, which is exactly the output format that separation of variables produces.

What role does the divergence theorem play in the proof?

The divergence theorem converts a boundary integral (net flux of heat through the surface of a region) into a volume integral (of the divergence of the flux, over the interior of the region). This step is essential because conservation of energy is naturally stated as a boundary flux balance, while a pointwise differential equation requires everything expressed as integrals over the same region. Without the divergence theorem, one could not localize the conservation law into a PDE at all; the technique is a direct application of vector calculus that anyone using this derivation should verify applies (i.e., that the flux field is C1C^1 and the boundary is sufficiently regular), as flagged explicitly in the regularity hypothesis.

Before this

Where this leads