ODE: The Heat Equation

The heat equation is the fundamental governing law for thermal diffusion, describing how temperature evolves in a solid object over time. Mastering it is not just an academic exercise; it is essential for predicting thermal stress in engine components, designing efficient heat sinks for electronics, and understanding processes from metal forging to climate modeling. As a cornerstone of partial differential equations (PDEs), its solution techniques form the bedrock for analyzing a vast array of diffusion-type phenomena in engineering and physics.

Derivation from Energy Conservation and Fourier's Law

We begin by deriving the heat equation from first principles, combining a physical conservation law with an empirical constitutive relation. Consider a one-dimensional rod of constant cross-sectional area $A$, material density $\rho$, and specific heat capacity $c$. The fundamental principle is conservation of energy: the rate of change of thermal energy within a small segment $[x,x+\Delta x]$ equals the net rate of heat flow into the segment plus any internal heat generation.

The thermal energy $E$ in the segment is $E=\rho cA\Delta xu(x,t)$, where $u(x,t)$ is the temperature. Thus, the rate of change is $\rho cA\Delta x\frac{\partial u}{\partial t}$.

Heat flow is governed by Fourier's law of heat conduction, which states that the heat flux $q$ (heat flow per unit area per unit time) is proportional to the negative temperature gradient: $q=-k\frac{\partial u}{\partial x}$. The constant $k$ is the thermal conductivity, a material property. The net heat flow into the segment is $A[q(x,t)-q(x+\Delta x,t)]$.

Assuming no internal heat generation for the standard derivation, energy conservation gives: $$\rho cA\Delta x\frac{\partial u}{\partial t}\approx A\left[-k\frac{\partial u}{\partial x}{|}_x+k\frac{\partial u}{\partial x}{|}_{x+\Delta x}\right].$$ Dividing by $\rho cA\Delta x$ and taking the limit as $\Delta x\to 0$ yields the classic one-dimensional heat equation: $$\frac{\partial u}{\partial t}=\alpha \frac{{\partial }^{2}u}{\partial x^2}.$$ Here, $\alpha =k/(\rho c)$ is the thermal diffusivity. This parameter dictates how quickly heat spreads through a material. A high $\alpha$ (like in metals) means rapid temperature equalization; a low $\alpha$ (like in insulators) means slow diffusion.

Initial and Boundary Conditions

A PDE alone has infinitely many solutions. To model a specific physical scenario, we must impose an initial condition and boundary conditions. The initial condition specifies the temperature distribution throughout the rod at time zero: $u(x,0)=f(x)$.

Boundary conditions specify what happens at the ends of the rod ($x=0$ and $x=L$, for a rod of length $L$). Two primary types are essential:

1. Dirichlet Condition: This specifies the temperature at the boundary. For example, $u(0,t)=T_0$ models an end held at a constant fixed temperature by contact with a reservoir.
2. Neumann Condition: This specifies the temperature gradient (and thus the heat flux via Fourier's law) at the boundary. For example, $\frac{\partial u}{\partial x}(L,t)=0$ models a perfectly insulated end where no heat flows out.

Mixed or Robin conditions (a linear combination of $u$ and its gradient) are also possible, modeling convective heat transfer. The choice of boundary conditions profoundly influences the solution's form and behavior.

Solution via Separation of Variables

The primary analytical tool for solving the linear heat equation on a finite domain with homogeneous boundary conditions is the method of separation of variables. We assume a solution of the form $u(x,t)=X(x)T(t)$. Substituting into the heat equation and separating variables leads to: $$\frac{T^{\prime }(t)}{\alpha T(t)}=\frac{X^{\prime \prime }(x)}{X(x)}=-\lambda .$$ Since one side depends only on $t$ and the other only on $x$, both must equal a constant, which we call the separation constant $-\lambda$.

This yields two ordinary differential equations:

1. Temporal: $T^{\prime }+\alpha \lambda T=0$, with solution $T(t)=e^{-\alpha \lambda t}$.
2. Spatial: $X^{\prime \prime }+\lambda X=0$.

The spatial ODE is solved subject to the homogeneous boundary conditions. For instance, with Dirichlet conditions $u(0,t)=u(L,t)=0$, we require $X(0)=X(L)=0$. This turns into an eigenvalue problem. Non-trivial solutions ($X(x)\not\equiv 0$) exist only for specific eigenvalues ${\lambda }_n=(\frac{n\pi }{L}{)}^2$, with corresponding eigenfunctions (spatial modes) $X_n(x)=\sin (\frac{n\pi x}{L})$, for $n=1,2,3,...$.

Fourier Series Solution and Interpretation

The general solution is a superposition (infinite sum) of all possible separated solutions: $$u(x,t)=\sum _{n=1}^{\infty }{b}_n\sin \left(\frac{n\pi x}{L}\right)e^{-\alpha (n\pi /L{)}^{2}t}.$$ The coefficients $b_n$ are determined by the initial condition. At $t=0$, the solution must equal $f(x)$: $$f(x)=\sum _{n=1}^{\infty }{b}_n\sin \left(\frac{n\pi x}{L}\right).$$ This is a Fourier sine series expansion of $f(x)$. The coefficients are found using orthogonality: $$b_n=\frac{2}{L}{\int }_0^{L}f(x)\sin \left(\frac{n\pi x}{L}\right)dx.$$

This solution reveals deep physical insight. The temperature distribution is decomposed into a sum of spatial sine waves (modes). Each mode decays exponentially in time at a rate proportional to $\alpha {\lambda }_n=\alpha (n\pi /L{)}^2$. Higher-frequency modes (larger $n$) decay much faster because $(n\pi /L{)}^2$ grows rapidly. This means sharp temperature spikes (composed of many high-n modes) smooth out quickly, while the overall shape (dominated by low-n modes) evolves more slowly.

Steady-State and Transient Temperature Distributions

For problems with non-homogeneous boundary conditions (e.g., one end held at $T_1$, the other at $T_2$), it is powerful to split the solution into steady-state and transient parts: $u(x,t)=u_{ss}(x)+u_{tr}(x,t)$.

The steady-state solution $u_{ss}(x)$ is the temperature distribution reached as $t\to \infty$, when all time derivatives vanish. It satisfies $\frac{d^2{u}_{ss}}{d{x}^2}=0$, subject to the original non-homogeneous boundary conditions. The solution is simply a straight line: $u_{ss}(x)=T_1+\frac{(T_2-T_1)}{L}x$.

The transient solution $u_{tr}(x,t)$ then satisfies the homogeneous heat equation with homogeneous boundary conditions (obtained by subtracting $u_{ss}$ from the original problem). Its initial condition is $u_{tr}(x,0)=f(x)-u_{ss}(x)$. We solve for $u_{tr}$ using separation of variables as before, and it will decay to zero as $t\to \infty$, leaving only the steady-state profile. This decomposition simplifies the math and clarifies the physics: the system evolves from an initial state toward a final equilibrium linear profile.

Common Pitfalls

1. Misapplying Boundary Conditions to the Spatial ODE: A frequent error is applying the original non-homogeneous boundary conditions directly to the ODE $X^{\prime \prime }+\lambda X=0$ during separation of variables. This yields no useful solution. You must first use the decomposition into steady-state and transient parts, or use a different technique, to handle non-homogeneous conditions. Separation of variables applies directly only when the boundary conditions for the entire PDE are homogeneous.

1. Incorrect Interpretation of the Separation Constant: Setting the separated equation equal to $+\lambda$ instead of $-\lambda$ is a common sign error. This leads to spatial solutions like $\sinh$ and $\cosh$ and a temporal solution $e^{+\alpha \lambda t}$, which grows exponentially—a physical impossibility for heat diffusion without an internal source. The negative sign is crucial for obtaining decaying exponentials in time.

1. Neglecting the Role of Thermal Diffusivity ($\alpha$): Treating $\alpha$ as merely a constant factor misses its critical physical meaning. When comparing solutions or scaling time, remember that the decay rate for each mode is $\alpha {\lambda }_n$. A problem solved for a copper rod ($\alpha \approx 1.1\times 10^{-4}{m}^2/s$) will evolve over seconds, while the same geometry in brick ($\alpha \approx 5\times 10^{-7}{m}^2/s$) will take hours.

1. Confusing Heat Flux with Temperature Gradient: A Neumann boundary condition specifies the gradient, $\frac{\partial u}{\partial x}$. Since Fourier's law states $q=-k\frac{\partial u}{\partial x}$, a condition of $\frac{\partial u}{\partial x}=0$ means zero heat flux (insulation). A condition like $\frac{\partial u}{\partial x}=C$ (a constant) means a constant inward heat flux of $-kC$. Always connect the mathematical condition back to the physical heat flow.

Summary

• The one-dimensional heat equation, $\frac{\partial u}{\partial t}=\alpha \frac{{\partial }^{2}u}{\partial x^2}$, is derived from energy conservation and Fourier's Law ($q=-k\frac{\partial u}{\partial x}$). The key material property is thermal diffusivity $\alpha =k/(\rho c)$.
• Physically meaningful solutions require an initial condition $u(x,0)=f(x)$ and boundary conditions (Dirichlet for fixed temperature, Neumann for specified heat flux/insulation).
• The primary solution method for homogeneous problems is separation of variables, which leads to an eigenvalue problem for spatial modes and exponential decay in time. The full solution is a Fourier series whose coefficients are determined by the initial temperature distribution.
• The solution can be understood as a sum of spatial modes, where higher-frequency components decay much faster, explaining the rapid smoothing of sharp temperature variations.
• For non-homogeneous boundary conditions, splitting the solution into a steady-state part (a time-independent linear function) and a decaying transient part simplifies analysis and clarifies the approach to equilibrium.