The 1D heat (diffusion) equation governs temperature \(u(x,t)\) along a rod of length \(L\) with thermal diffusivity \(\alpha > 0\):
\[\frac{\partial u}{\partial t} = \alpha\,\frac{\partial^2 u}{\partial x^2}, \qquad x \in [0,L],\; t > 0\]
Physically: the local rate of temperature change equals the Laplacian (curvature) of \(u\) scaled by \(\alpha\). Peaks lose heat; valleys gain it — the equation is a smoothing operator.
2. Boundary & Initial Conditions
\[u(0,t) = 0,\quad u(L,t) = 0\quad\text{(Dirichlet — rod ends held at zero)}\]
\[u(x,0) = f(x)\quad\text{(arbitrary initial profile)}\]
3. Separation of Variables
Writing \(u(x,t) = X(x)\,T(t)\) and substituting yields the separated equations with constant \(-\lambda\):
This app computes \(B_n\) via a trapezoidal quadrature rule \(\bigl(B_n \approx \frac{2\,\Delta x}{L}\sum_i f(x_i)\sin(n\pi x_i/L)\bigr)\), which is exact for the eigenfunction basis on the grid.
5a. FTCS Numerical Scheme (explicit)
Discretise with \(\Delta x = L/(N-1)\) and timestep \(\Delta t\). The Forward-Time Centered-Space update is:
Stability requires the CFL condition \(r \leq 0.5\). This app enforces \(r = 0.49\) automatically, giving first-order accuracy in time, second-order in space.
5b. Crank-Nicolson Scheme (implicit)
Averaging the spatial operator at time levels \(n\) and \(n+1\) gives the unconditionally stable implicit scheme:
Each step solves a tridiagonal linear system via the Thomas algorithm — \(O(N)\) cost. C-N is second-order in both time and space and remains stable for any \(r\), allowing timesteps 4× larger than FTCS here. The stability region covers the entire left half-plane, unlike FTCS which has a disk of radius 1.
6. Key Mathematical Properties
Infinite smoothing: \(u(\cdot,t) \in C^\infty\) for any \(t > 0\), even if \(f\) is discontinuous.
Maximum principle: \(\min f \leq u(x,t) \leq \max f\) for all \((x,t)\). No new extrema are created.
Energy dissipation: \(\tfrac{d}{dt}\|u\|^2 = -2\alpha\|\nabla u\|^2 \leq 0\) — energy is monotonically lost.
Mode decay ordering: The \(n=1\) half-life is \(T_{1/2} = L^2\ln 2\,/\,(\alpha\pi^2)\); every higher mode decays at least \(4\times\) faster.
Self-adjointness: The operator \(-\alpha\,\partial^2/\partial x^2\) is self-adjoint on \(L^2[0,L]\) with Dirichlet BCs, guaranteeing a complete real orthonormal eigenfunction basis.
Gibbs phenomenon: Truncating the Fourier series at \(N\) modes near a jump discontinuity causes ~8.9% overshoot/undershoot that does not vanish as \(N \to \infty\).