Diffusion equation

A diffusion equation describes smoothing caused by a flux down spatial gradients.

A spatially uniform time-dependent reaction coefficient can be removed from a linear diffusion equation by a scalar integrating factor. If $A^{\prime }(t)=a(t)$ and $v=e^{A(t)}u$, substitution gives $u_t=u_{xx}$. Initial data and boundary values must also be divided by $e^{A(t)}$. A nonzero boundary profile can then be removed by boundary lifting before applying separation of variables. This reduction distinguishes growth of the whole amplitude from decay of the spatial diffusion modes.

For constant $\nu$, variation of this positive functional gives $({\partial }_t+\nu {\nabla }^2)({\partial }_t-\nu {\nabla }^2)h=0$, or $h_{tt}-{\nu }^2{\nabla }^{4}h=0$. The diffusion equation $h_t=\nu {\nabla }^{2}h$ gives zero-residual minima, but is only a subset of the stationary equations of the squared functional. Each nonzero spatial Fourier mode has both growing and decaying temporal solutions in the stationary equation. Boundary data decide which stationary paths are admissible. In two spatial dimensions, the bulk action is invariant under the anisotropic scale invariance $x\mapsto bx$, $t\mapsto b^{2}t$, $h\mapsto h$.

For a real scalar height field, integration by parts givesThe bulk expression is even under time reversal. The original density is not pointwise even, but changes by a total divergence; its stationary equation is time-reversal invariant. The first-order diffusion equation selected by zero residual is not invariant under ordinary time reversal at fixed positive $\nu$. Boundary conditions or the probabilistic interpretation can distinguish the two time directions even when the higher-order stationary equation does not.

These Robin boundary conditions give the spatial Sturm-Liouville operator $-d^2/d{x}^2$ one negative eigenvalue, $-1$, with eigenfunction $e^x$. The other eigenvalues are $n^2{\pi }^2$, with eigenfunctions $V_n=n\pi \cos (n\pi x)+\sin (n\pi x)$, $n\ge 1$; there is no zero eigenvalue. Hence the diffusion equation has a growing contribution $c_0{e}^{t+x}$ in addition to decaying modes.

There is no contradiction with ordinary diffusive smoothing. The outward normal derivative is $-u_x$ at the left endpoint and $u_x$ at the right endpoint, so imposing the same coordinate derivative at both ends does not impose dissipative outward-normal data at both ends. The energy identity is $\frac{1}{2}{\partial }_t\int u^2=[u{u}_x{]}_0^1-\int u_x^2$, whose boundary term is $u(1{)}^2-u(0{)}^2$ here and can supply energy. Orthogonal projection of initial data $f$ gives $c_0=2{\int }_0^{1}f{e}^{x}dx/(e^2-1)$; it vanishes precisely when the growing mode is absent.

The reflected boundary conditions $u_x(0,t)=-u(0,t)$ and $u_x(1,t)=-u(1,t)$ have the same spectrum, because reflection $x\mapsto 1-x$ reverses the coordinate derivative. Here the growing eigenfunction is $e^{-x}$ and the decaying functions are $W_n=n\pi \cos (n\pi x)-\sin (n\pi x)$, givingFor square-integrable initial data $f$, orthogonal projection givesThe norms are $(1-e^{-2})/2$ and $(n^2{\pi }^2+1)/2$. The spatial Sturm-Liouville problem has no zero eigenvalue and exactly one negative eigenvalue for $-D^2$, namely $-1$. The energy boundary term now reads $u(0{)}^2-u(1{)}^2$, so the left endpoint can supply energy. The growing mode disappears exactly when ${\int }_0^1{e}^{-x}f(x)dx=0$.

For positive diffusivity $D(x)$, a diffusion equation with potential-driven drift has diffusive flux $J=-D(c_x+\beta V^{\prime }c)$. On a finite interval with zero flux at both endpoints, the total amount is conserved. The substitution $c=e^{-\lambda t}{e}^{-\beta V}y$ converts the decay eigenvalue problem to the Sturm-Liouville problemThe Rayleigh quotient is $\lambda =\int p|y^{\prime }{|}^2/\int w|y{|}^2\ge 0$. Its zero eigenfunction is constant, giving the stationary solution $c_{\mathrm{eq}}=C{e}^{-\beta V}$. This is a spatial Boltzmann distribution; $D$ changes relaxation rates but cancels out of the zero-flux equilibrium profile. For a sufficiently regular positive $D$ and bounded $V$, the constant eigenfunction is the unique zero mode up to scale.

A concentration varying over length $L$ has Laplacian of order its amplitude divided by $L^2$. Balancing the time derivative with $\kappa \Delta \chi$ in the diffusion equation gives the displayed relaxation time. A multiple-scale expansion with spatial scale $\epsilon x$ therefore uses slow time ${\epsilon }^{2}t$ when the effective diffusivity is of order one.

With positive diffusion coefficient $a(x)$, this diffusion equation expresses conservation with flux $-a(x)u_x$. For a differentiable coefficient, its spatial operator is $a{u}_{xx}+a^{\prime }{u}_x$, rather than either $a{u}_x$ or $a{u}_{xx}$ alone. A symmetric half-grid diffusion consistency discretization uses the difference of the two adjacent interface fluxes, preserving that divergence structure. Suitable Dirichlet boundary conditions make the spatial operator dissipative in the L2 norm.

Infinite propagation speed means that a disturbance initially confined to a bounded region can affect arbitrarily distant points after any positive time. The strictly positive heat kernel gives this behavior for a nontrivial nonnegative diffusion solution. It contrasts with the finite propagation speed of the wave equation; a very small distant response is still mathematically nonzero.

For $u_t={\partial }_x[a(x/ϵ)u_x]$ with positive periodic $a$, the one-dimensional homogenized diffusivity is the harmonic mean $a_{\mathrm{eff}}=⟨a^{-1}{⟩}^{-1}$. Microscopic flux continuity produces this harmonic rather than arithmetic averaging.

The porous medium equation $u_t=\Delta (u^m)$ with $m>1$ is a nonlinear partial differential equation whose diffusivity vanishes at $u=0$. Unlike the heat equation, it has compactly supported solutions with finite propagation speed.

For the nonlinear diffusion equation $u_t=(u{u}_x{)}_x+u/(2t)$ on $x>0,t>0$, the fixed boundary value $u(0,t)=1$ and zero slope at $x=0$ admit the displayed compact-support similarity solution, with front $x_F=2\sqrt{t}$. Substitution $u=f(x/\sqrt{t})$ reduces the equation to $(f{f}^{\prime }+\eta f/2{)}^{\prime }=0$. The boundary data make the integration constant zero. While $f>0$, division gives $f^{\prime }=-\eta /2$, hence $f=1-{\eta }^2/4$. At its first zero the diffusive flux $u{u}_x$ vanishes, so extension by zero creates no singular flux term. The solution is smooth inside its positive region and is a weak solution at the moving front.

For $C_t=k\nabla \cdot (C\nabla C)$ in three dimensions and total mass $4\pi M$, the nonnegative radial source solution is $C=M^{2/5}(kt{)}^{-3/5}f(r/(Mkt{)}^{1/5})$, with $f(\xi )=(75^{2/5}-{\xi }^2{)}_{+}/10$. Integrating the similarity equation once gives $f{f}^{\prime }=-\xi f/5$, and normalization ${\int }_0^{\infty }{\xi }^{2}fd\xi =1$ fixes the compact-support radius. The flux vanishes at the moving front, so the profile solves the equation weakly and approaches the prescribed point mass as time tends to zero.

For the planar diffusion equation $n_t=D_0\nabla \cdot [(n/n_0{)}^2\nabla n]$, a mass-preserving radial similarity solution has $n=n_0{\lambda }^{-2}F(r/(r_0\lambda ))$. Substitution reduces the equation to $(\xi F^2{F}^{\prime }{)}^{\prime }=-({\xi }^{2}F{)}^{\prime }$ when ${\lambda }^5\dot{\lambda }=D_0/r_0^2$. Zero central flux and $F(0)=1$ give the compact profile shown above, with ${\lambda }^6=6D_{0}t/r_0^2$. Its height decays as $t^{-1/3}$ and its radius grows as $t^{1/6}$; its total mass is $2\pi n_0{r}_0^2/3$.

In a porous medium equation with exponent $m>1$, the effective diffusion coefficient vanishes at zero density. A nonnegative compactly supported source profile can therefore retain a finite front for every finite positive time. For example, the radial cubic-diffusion source profile has support radius proportional to $t^{1/6}$. This notion concerns a moving diffusion front; it does not require the constant characteristic cone used for a wave equation.

Consider a nonlinear diffusion equation with uniform supply, $\varphi h_t=D_m(h^m{h}_x{)}_x+R$, on $x>0$, with an absorbing boundary $h(0,t)=0$ and initially $h=0$. Far from the boundary, $h=Rt/\varphi$. Balancing the time derivative, supply, and diffusion givesThe similarity solution satisfiesThe outward boundary volume flux per unit width is $Q=R\ell c_m$, where $c_m={\lim }_{\eta \downarrow 0}{F}^m{F}^{\prime }$. Integrating the equation givesThus the growing region of depleted storage fixes the discharge prefactor, and $Q\propto t^{(m+1)/2}$.

For $\varphi h_t=D_m(h^m{h}_x{)}_x$ on $0<x<L$, with $h(0,t)=0$ and zero right-hand volume flux, a separation of variables giveswhere the positive profile satisfiesIts boundary volume flux per unit width isThese profiles are exact separated solutions and describe the leading long-time discharge for a broad class of positive initial data. The virtual origin $t_0$ depends on the initial profile; it does not make the separated solution an exact representation of every initial condition.

The Barenblatt solution is a mass-preserving similarity solution of the porous medium equation. For $u_t=\beta (u{u}_x{)}_x$ in one dimension it has the form

For the planar porous medium equation $n_t=[D_0/(3n_0^2)]\Delta (n^3)$, the mass-$Q$ point-source similarity solution is $n=n_0{\lambda }^{-2}\sqrt{[1-r^2/(r_0^2{\lambda }^2){]}_{+}}$, with $\lambda =(6D_{0}t/r_0^2{)}^{1/6}$ and $r_0^2=3Q/(2\pi n_0)$. Its front grows as $t^{1/6}$ and its central density falls as $t^{-1/3}$. The mass conservation integral uses the planar area element $2\pi rdr$. Although the density slope diverges at the front, its diffusion flux vanishes there, giving a compactly supported weak solution with initial measure $Q{\delta }^{(2)}$.

For $n,A>0$, the one-dimensional porous medium equation $h_t=A(h^n{h}_x{)}_x$ has a volume-preserving similarity solution $h=[c_n(R^2-x^2)/(At){]}_{+}^{1/n}$, where $c_n=n/[2(n+2)]$ and $R\propto t^{1/(n+2)}$. On a reflecting half-line with conserved area $M$, $R=[M(At/c_n{)}^{1/n}/I_n{]}^{n/(n+2)}$, where $I_n=\frac{1}{2}B(\frac{1}{2},1+1/n)$. For a whole symmetric release use half its total area as $M$. The Beta function fixes normalization, while finite-front zero volume flux fixes the edge condition.

The source-type similarity solution of $h_t=D(h^2{h}_x{)}_x$ with whole-line conserved area $\mathcal{V}$ is $h=H(t)[1-x^2/L(t{)}^2{]}_{+}^{1/2}$, where $L=2(\mathcal{V}/\pi {)}^{1/2}(Dt{)}^{1/4}$ and $H=(\mathcal{V}/\pi {)}^{1/2}(Dt{)}^{-1/4}$. Its volume flux vanishes at the finite front, while the velocity tends to $\dot{L}=L/(4t)$. A reflecting half-line uses twice its prescribed area in these formulas.

For the radially symmetric equationin two dimensions, the mass-$N$ similarity solution has the compactly supported formwhere mass normalization gives $x_0=(16/\pi {)}^{1/4}$. Its support radius grows as $t^{1/4}$ and its occupied area grows as $t^{1/2}$.

Fick's laws describe diffusion through a flux law and conservation of mass. The first law says that diffusive flux points down the concentration gradient:Combining it with local conservation gives
$C_t=\nabla \cdot (D(C)\nabla C)$.

The heat equationmodels diffusion with constant diffusivity $D$.

For thermal diffusivity $\kappa >0$ and thermal conductivity $k>0$, a semi-infinite conductor satisfies $u_t=\kappa u_{xx}$ for $x>0$. Take zero initial temperature and constant incoming heat flux $-k{u}_x(0,t)=Q$. Select the physical solution by boundedness at spatial infinity together with sufficient growth control for the time Laplace transform.

Writing the time transform as $U$ gives $U_{xx}=(p/\kappa )U$. The decaying branch is $U=A(p)e^{-x\sqrt{p/\kappa }}$, where $\mathrm{Re}\sqrt{p}>0$. The transformed flux condition gives $A=Q\sqrt{\kappa }/(k{p}^{3/2})$. Since $\mathcal{L}[2\sqrt{t}]=\sqrt{\pi }{p}^{-3/2}$, inversion gives the displayed surface temperature. More generally the full profile isDifferentiation verifies the heat equation and the applied flux; for every $x>0$ it tends to zero at initial time and decays at spatial infinity. Its surface temperature grows as $\sqrt{t}$, because the depth significantly heated by diffusion has the same scale.

On the real line, let $T(x)=\max (1-|x|/2,0)$. The Fourier transform of a triangular function gives $\hat{T}(\omega )=2{\sin }^2\omega /{\omega }^2$, with value two at zero. The Gaussian heat kernel has transform $e^{-t{\omega }^2}$. Hence $u=q_t\ast T$ solves the heat equation with initial profile $T$, and the inverse Fourier transform gives the displayed integral. Differentiating for positive time is justified by Gaussian decay. The uniform continuity of $T$ and the vanishing Gaussian mass outside any fixed neighborhood of zero prove uniform recovery of $T$ as $t\downarrow 0$. Convolution preserves total integral two, keeps the solution nonnegative, and smooths the corners for every positive time. Unlike the finite-interval problem in heat flow from a triangular initial temperature, this evolution has no endpoints through which heat can escape.

Let a homogeneous plate $0<x<L$ begin at uniform temperature $T_0$, cool with zero-temperature Dirichlet boundary conditions until time $t_b$, then switch to insulated Neumann boundary conditions. For thermal diffusivity $\kappa >0$, the Fourier sine series during cooling is $T=(4T_0/\pi ){\sum }_{n\ge 1,n\text{ odd}}\sin (n\pi x/L)e^{-\kappa (n\pi /L{)}^{2}t}/n$. Integrating the insulated heat equation shows that its integral average is conserved. All nonconstant Fourier cosine series modes then decay, proving that the final uniform temperature is $T_{\infty }=(8T_0/{\pi }^2){\sum }_{n\ge 1,n\text{ odd}}{e}^{-\kappa (n\pi /L{)}^2{t}_b}/n^2$. Its limits are $T_0$ at zero bath time and zero at infinite bath time. Instantaneous changes of boundary condition can create initial-boundary corner incompatibilities without invalidating the positive-time solution.

Integrate $u_t=\alpha \Delta u+Q$ over a periodic cell. Opposite boundary contributions to the Laplacian cancel, so the spatial mean obeys ${\overline{u}}^{\prime }=\overline{Q}$. With no source the mean is conserved, although the local temperature changes. The mean is the zero Fourier mode in a normalized Fourier series convention.

For a spatially constant PDE source term $Q(t)$, subtract its time integral from $u$. The result obeys the homogeneous heat equation with the same initial condition. This proves the displayed formula and shows that uniform heating changes only the zero Fourier mode. In an unnormalised $N_x\times N_y$ discrete Fourier transform, a constant increment $q$ adds $N_x{N}_{y}q$ to the zero-frequency coefficient.

For constant volumetric heat generation rate $q_v$ and thermal conductivity $\kappa$, spherical steady thermal conduction obeys $\kappa r^{-2}(r^2{T}^{\prime }{)}^{\prime }+q_v=0$. Regularity at the centre and a specified surface temperature give the displayed parabolic profile. For a black surface radiating into negligible background radiation, global energy conservation gives ${\sigma }_{\mathrm{SB}}{T}_S^4=q_{v}R/3$. The central excess scales as $R^2$, while the surface temperature scales only as $R^{1/4}$.

For two length-$L$ slabs with zero exterior temperatures and perfect thermal contact, put ${\alpha }_i=\sqrt{{\rho }_i{c}_i/k_i}$, $u=\omega {\alpha }_{1}L$, and $v=\omega {\alpha }_{2}L$. A separated solution with time factor $e^{-{\omega }^{2}t}$ and $\omega \ne 0$ exists precisely when the displayed interface determinant vanishes. Dividing by $\cos u\cos v$ gives a tangent characteristic equation only away from its poles. Common zeros of both cosines are legitimate normal modes and must be retained. The zero-frequency case gives only the trivial solution of a differential equation with these positive material parameters and exterior boundary conditions.

A heat equation problem in adjacent materials with separate thermal conductivities and volumetric heat capacities. Perfect thermal contact imposes continuous temperature and continuous normal heat flux at the interface. Multiplying each material's heat equation by its temperature cancels the interface boundary terms in a weighted energy estimate.

Suppose $k>0$ and two sufficiently regular steady state temperature profiles solve the same heat equation with the same PDE source term $H$ and the same Dirichlet boundary condition $h$. Their difference $w$ satisfies $\nabla \cdot (k\nabla w)=0$ and $w=0$ on the boundary. The divergence theorem applied to $wk\nabla w$ givesThus $\nabla w=0$, and the zero boundary values force $w=0$ on each connected component. The argument proves uniqueness without asserting existence and without replacing the spatially varying thermal conductivity by a constant. For weak solutions, the same calculation is the energy argument obtained by choosing $w$ as the test function, under the usual boundedness and positive lower bound assumptions on $k$.

Integrating energy conservation over a fixed material volume gives storage equal to incoming conductive heat flux plus internally generated heat. Here $c_{\mathrm{vol}}$ is volumetric heat capacity, with units $\mathrm{J}{\mathrm{m}}^{-3}{\mathrm{K}}^{-1}$. For a mass-specific heat capacity, the density factor must be included. Constant isotropic conductivity reduces the divergence to $k{\nabla }^{2}T$; variable conductivity must remain inside the divergence.

For a steady state heat equation $\nabla \cdot (k\nabla T)+H=0$ in a bounded region with zero outward heat flux, integration and the divergence theorem give the necessary compatibility conditionIn particular, a strictly positive PDE source term throughout a region of positive volume makes a steady state impossible. The time-dependent conservation of energy law instead gives $dE/dt={\int }_{R}HdV>0$: stored thermal energy must grow. This is a compatibility condition on the total source, not a restriction requiring the source to vanish pointwise.

A sinusoidal surface temperature applied to a homogeneous semi-infinite conductor produces a periodic depth profile with penetration length $\delta =\sqrt{2\kappa /(\omega C)}$. It follows by substituting a complex exponential into $C{\partial }_{t}T=\kappa {\partial }_z^{2}T$ and taking the root that decays at depth. A depth $z$ attenuates the amplitude by $e^{-z/\delta }$ and delays the phase by $z/\delta$. This is diffusion with phase delay, rather than a nondispersive travelling wave.

The thermal penetration depth of a harmonic ground temperature disturbance scales as the inverse square root of its angular frequency. For yearly and daily forcing, a 365-day year gives ${\delta }_Y/{\delta }_D=\sqrt{365}\simeq 19.1$. At ten daily penetration depths the daily component is attenuated by $e^{-10}$, whereas the annual one retains about $59.2\%$ of its surface amplitude. At ten annual penetration depths, even the annual component is attenuated by $e^{-10}$. These are consequences of the linear thermal diffusion wave model, with homogeneous constant thermal diffusivity.

For positive thermal diffusivity $\nu$, thermal conductivity $k$ and length $L$, subtract the stationary solution $w(x)=Hx/k$. The remainder satisfies the heat equation with a homogeneous Dirichlet boundary condition at zero and a homogeneous Neumann boundary condition at $L$. The mixed Dirichlet-Neumann modes on an interval are $\sin (k_{n}x)$, where $k_n=(n+1/2)\pi /L$. Thus for an initial temperature profile $u_0$,For zero initial temperature, $b_n=-2H(-1{)}^n/(kL{k}_n^2)$. The series gives the initial profile in the usual Fourier series sense and a smooth solution for $t>0$. If $H\ne 0$, the zero initial profile does not satisfy the prescribed endpoint derivative at $t=0$: the temperature can be continuous there without its spatial derivative extending continuously to that corner. The sign $k{u}_x(L,t)=H>0$ describes incoming heat flux, since the outward conductive heat flux is $-k{u}_x$.

On $0<x<L$ with zero-temperature endpoints, let the initial temperature rise linearly from zero to height $H$ at $L/2$ and then fall linearly to zero. Integration by parts on the two halves gives Fourier sine series coefficients $8H\sin (n\pi /2)/({\pi }^2{n}^2)$. Under the heat equation $T_t=D{T}_{xx}$ with $D>0$, each coefficient decays by $e^{-D(n\pi /L{)}^{2}t}$, giving the displayed solution. Only odd spatial modes occur; absolute coefficient summability proves uniform convergence for $t\ge 0$ via the Weierstrass M-test. This spatial profile also occurs in a midpoint-plucked fixed-end string, where the time factors oscillate instead of decaying.

On a bar of length $L$, let the initial temperature increase linearly from zero at either endpoint to height $H$ at the midpoint. With thermal diffusivity $D>0$ and insulated Neumann boundary conditions, a Fourier cosine series gives the displayed solution. Its constant term is the initial integral average $H/2$. Integration by parts on the two straight segments gives zero coefficients for odd modes and modes divisible by four; the remaining coefficients are $-4H/[{\pi }^2(2j+1{)}^2]$. Nonconstant modes decay, leaving the uniform temperature $H/2$, consistent with conservation of the total heat. This contrasts with the zero-temperature boundary solution in heat flow from a triangular initial temperature, where heat can leave the bar and the limit is zero. The initial endpoint slopes do not satisfy the insulated condition, so derivatives need not converge to their initial endpoint values as $t\downarrow 0$; the series satisfies the Neumann condition for every $t>0$ and recovers the continuous initial profile.

For a square of side $\ell$ initially at constant temperature $T_0$, held at zero boundary temperature for $t>0$, the heat equation $T_t=\kappa \Delta T$ with $\kappa >0$ has solutionThe double Fourier sine series of the initial constant gives the coefficients. The spatial modes are Dirichlet Laplacian eigenfunctions, so separation of variables gives their decay rates. At late times the $(1,1)$ mode dominates, giving the displayed profile; the next nonzero decay exponent is $10\kappa {\pi }^2/{\ell }^2$. At early times the series approaches the constant at each interior point and in L2 space, while the boundary remains at zero. An instantaneous change to the boundary condition therefore need not produce a solution continuous at the initial boundary in space-time.

For a circle of circumference $2\pi$, the heat equation $u_t=D{u}_{xx}$ has periodic matching of both $u$ and $u_x$ at $x=\pm \pi$. For $D>0$, the displayed Fourier series solves the initial-value problem, with $\bar{f}=(2\pi {)}^{-1}{\int }_{-\pi }^{\pi }f(x)dx$ and the usual real Fourier coefficients. Every nonconstant normal mode decays, leaving the conserved mean. In particular $f(x)=(\pi -|x|{)}^2$ gives $\bar{f}={\pi }^2/3$, $a_n=4/n^2$ and $b_n=0$. The initial cusp is smoothed for every positive time.

For $L>0$, extend $f_L(x)=\sin x-x\sin L/L$ from $[-L,L]$ periodically. The linear term makes both endpoint values zero and the endpoint first derivatives match. The function is odd, so its Fourier series has only sine terms, with the displayed Fourier coefficients. Indeed, for $k_n=n\pi /L\ne 1$,subtracting and multiplying by $2/L$ proves the formula. The coefficients decay as $O(n^{-3})$. At $L=m\pi$, the single coefficient $b_m$ has removable limit $1$, with all others zero. Multiplying $b_n$ by $e^{-\kappa k_n^{2}t}$ solves the corresponding heat equation on the circle.

On $0<x<L$, let $u_t=\nu u_{xx}+G\sin (n\pi x/L)$, with $\nu >0$, zero Dirichlet boundary conditions, and zero initial data. Write $k_n=n\pi /L$. Its stationary solution is $u_s=G\sin (k_{n}x)/(\nu k_n^2)$. Subtracting it removes the forcing from the partial differential equation and leaves a heat equation with initial data $-u_s$.

Because the forcing is a single spatial eigenfunction, separation of variables gives the entire solution without mixing modes:Its integral is $G[1-(-1{)}^n](1-e^{-\nu k_n^{2}t})/(\nu k_n^3)$; even modes have zero signed net integral. Removing forcing from the evolution equation does not remove its amplitude from the initial data for the transient.