Basic Finite Differences

Spatial and temporal grid discretization of differential operators for numerical physics simulations.

1. Domain Discretization and Differential Operators

In continuum physics, conservation laws (mass, momentum, thermal energy, electromagnetism) are formulated via Partial Differential Equations (PDEs). Unlike ODEs, PDEs govern scalar or vector fields $u(\vec{x}, t)$ distributed over both space and time.

The Finite Difference Method (FDM) replaces continuous derivatives with discrete difference quotients evaluated over a structured grid. For a one-dimensional spatial domain of length $L$ over duration $T$:

$x_i = i \Delta x, \quad i = 0, 1, \dots, N_x, \quad \Delta x = \frac{L}{N_x}$ $t_n = n \Delta t, \quad n = 0, 1, \dots, N_t, \quad \Delta t = \frac{T}{N_t}$

The discrete grid approximation to the true continuous solution at node $i$ and time step $n$ is written: $u_i^n \approx u(x_i, t_n)$

2. Finite Difference Formulas via Taylor Expansions

Expanding $u(x \pm \Delta x)$ in Taylor series:

$u(x + \Delta x) = u(x) + \Delta x \frac{\partial u}{\partial x} + \frac{\Delta x^2}{2} \frac{\partial^2 u}{\partial x^2} + \frac{\Delta x^3}{6} \frac{\partial^3 u}{\partial x^3} + \mathcal{O}(\Delta x^4)$ $u(x - \Delta x) = u(x) - \Delta x \frac{\partial u}{\partial x} + \frac{\Delta x^2}{2} \frac{\partial^2 u}{\partial x^2} - \frac{\Delta x^3}{6} \frac{\partial^3 u}{\partial x^3} + \mathcal{O}(\Delta x^4)$

Yields standard finite-difference stencils:

1. Forward Difference: $\left.\frac{\partial u}{\partial x}\right|_i = \frac{u_{i+1} - u_i}{\Delta x} + \mathcal{O}(\Delta x)$
2. Backward Difference: $\left.\frac{\partial u}{\partial x}\right|_i = \frac{u_i - u_{i-1}}{\Delta x} + \mathcal{O}(\Delta x)$
3. Central Difference: $\left.\frac{\partial u}{\partial x}\right|_i = \frac{u_{i+1} - u_{i-1}}{2 \Delta x} + \mathcal{O}(\Delta x^2)$
4. Second-Order Central Spatial Difference: $\left.\frac{\partial^2 u}{\partial x^2}\right|_i = \frac{u_{i+1} - 2u_i + u_{i-1}}{\Delta x^2} + \mathcal{O}(\Delta x^2)$

3. The 1D Diffusion / Heat Equation and the Explicit FTCS Scheme

The canonical homogeneous heat conduction equation is: $\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}$ where $\alpha$ denotes thermal diffusivity ($[\text{m}^2/\text{s}]$).

Discretizing using forward differences in time and central differences in space yields the FTCS (Forward-Time Central-Space) scheme: $\frac{u_i^{n+1} - u_i^n}{\Delta t} = \alpha \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{\Delta x^2}$

Solving explicitly for the future state $u_i^{n+1}$: $u_i^{n+1} = u_i^n + r \left(u_{i+1}^n - 2u_i^n + u_{i-1}^n\right)$ with dimensionless diffusion parameter: $r = \frac{\alpha \Delta t}{\Delta x^2}$

4. Von Neumann Stability Analysis

Von Neumann stability examines the temporal amplification of individual Fourier spatial harmonics: $u_i^n = G^n e^{i k x_i} = G^n e^{i k (i \Delta x)}$ where $G$ is the amplification factor and $k$ is the spatial wave number. Scheme stability requires that rounding and truncation perturbations do not grow exponentially: $|G| \le 1, \quad \forall k$

Substituting into the FTCS equation: $G^{n+1} e^{i k i \Delta x} = G^n e^{i k i \Delta x} \left[1 + r \left(e^{i k \Delta x} - 2 + e^{-i k \Delta x}\right)\right]$ $G = 1 + 2r(\cos(k \Delta x) - 1) = 1 - 4r \sin^2\left(\frac{k \Delta x}{2}\right)$

Setting $|G| \le 1$: $-1 \le 1 - 4r \sin^2\left(\frac{k \Delta x}{2}\right) \le 1$ At the worst-case wavenumber mode where $\sin^2(k \Delta x / 2) = 1$: $-1 \le 1 - 4r \implies 4r \le 2 \implies r = \frac{\alpha \Delta t}{\Delta x^2} \le \frac{1}{2}$

Severe Physical Implication:

Refining spatial grid resolution by half ($\Delta x \to \Delta x/2$) forces a four-fold reduction in maximum allowable timestep $\Delta t$: $\Delta t \le \frac{\Delta x^2}{2\alpha}$

5. 1D Wave Equation and the CFL Condition

For mechanical or electromagnetic waves traveling in 1D: $\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2}$ Approximating both second derivatives using central stencils: $\frac{u_i^{n+1} - 2u_i^n + u_i^{n-1}}{\Delta t^2} = c^2 \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{\Delta x^2}$

Isolating $u_i^{n+1}$: $u_i^{n+1} = 2u_i^n - u_i^{n-1} + C^2 \left(u_{i+1}^n - 2u_i^n + u_{i-1}^n\right)$ where $C$ is the dimensionless Courant Number: $C = \frac{c \Delta t}{\Delta x}$

Courant-Friedrichs-Lewy (CFL) Condition:

Stability holds if and only if: $C = \frac{c \Delta t}{\Delta x} \le 1$ Geometric / Physical Insight: The numerical numerical domain of dependence must fully enclose the physical characteristic domain of dependence (the wave light-cone). When $C > 1$, physical signals propagate across cells faster than the discrete grid transfers data, triggering catastrophic blow-up.

6. Implicit Schemes and Crank-Nicolson

To eliminate the quadratic timestep barrier $\Delta t \sim \Delta x^2$ in diffusive physics, implicit algorithms are deployed. The Crank-Nicolson scheme averages spatial second-derivatives between time levels $n$ and $n+1$:

$\frac{u_i^{n+1} - u_i^n}{\Delta t} = \frac{\alpha}{2} \left[\frac{u_{i+1}^{n+1} - 2u_i^{n+1} + u_{i-1}^{n+1}}{\Delta x^2} + \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{\Delta x^2}\right]$

• Accuracy: Second-order in both time and space $\mathcal{O}(\Delta t^2 + \Delta x^2)$.
• Stability: Unconditionally stable for all $\Delta t > 0$ ($|G| \le 1$ universally).
• Computational Cost: Requires solving a tridiagonal linear system at each step via the Thomas algorithm.

Protocol for Discretizing and Solving PDEs with Finite Differences

Systematic Numerical Methodology

Solving an Initial Boundary Value Problem (IBVP) modeled by PDEs using finite differences requires a structured five-stage protocol:

Step 1: Computational Domain Discretization

1. Define the spatial domain $[0, L]$ and target simulation window $[0, T]$.
2. Choose spatial discretization count $N_x$, yielding cell width $\Delta x = L / N_x$.
3. Establish coordinate array $x_i = i \Delta x$ across indices $i = 0, 1, \dots, N_x$.

Step 2: Timestep Selection Governed by Stability Constraints

• Parabolic PDEs (Heat / Diffusion): $\Delta t \le \frac{\Delta x^2}{2 \alpha} \implies r = \frac{\alpha \Delta t}{\Delta x^2} \le 0.5$ Standard operational practice adopts a safety factor, e.g., $r = 0.4$.
• Hyperbolic PDEs (Wave Equations): $\Delta t \le \frac{\Delta x}{c} \implies C = \frac{c \Delta t}{\Delta x} \le 1.0$ Operational practice adopts $C = 0.8 - 0.9$.

Step 3: Imposing Initial Conditions (t = 0)

1. Sample initial continuous profiles at node locations: $u_i^0 = f(x_i), \quad \forall i = 0, 1, \dots, N_x$
2. For second-order time systems (wave equation) requiring $\dot{u}(x, 0) = g(x)$, employ a Taylor fictitious backward step to launch step $n = 1$: $u_i^1 = u_i^0 + \Delta t \, g(x_i) + \frac{C^2}{2} (u_{i+1}^0 - 2u_i^0 + u_{i-1}^0)$

Step 4: Imposing Boundary Conditions (BCs)

At each time step $n$, evaluate endpoints $i = 0$ and $i = N_x$:

• Dirichlet Boundary Conditions (Clamped values): $u_0^n = U_{\text{left}}, \quad u_{N_x}^n = U_{\text{right}}$
• Neumann Boundary Conditions (Prescribed flux / Adiabatic walls $\partial u / \partial x = 0$): Deploy a ghost node $u_{-1}^n$: $\frac{u_1^n - u_{-1}^n}{2 \Delta x} = 0 \implies u_{-1}^n = u_1^n$ Injecting into the internal stencil at $i = 0$: $u_0^{n+1} = u_0^n + 2r (u_1^n - u_0^n)$

Step 5: Time Marching and Field Update

1. Across iterations $n = 0, 1, \dots, N_t - 1$:
   • Apply the numerical stencil to all internal nodes $i = 1, 2, \dots, N_x - 1$.
   • Overwrite boundary points according to active BCs.
   • Advance array buffers ($u^n \leftarrow u^{n+1}$).
2. Perform periodic telemetry checks (total thermal energy, acoustic energy, mass conservation).

Worked Examples: Heat Conduction and 1D Wave Propagation

Example 1: Transient Thermal Diffusion Along a Metallic Rod

Consider a conductive metal rod of length $L = 1.0\text{ m}$ with thermal diffusivity $\alpha = 1.0\text{ m}^2/\text{s}$. Boundary ends are held at constant zero temperature: $u(0, t) = 0, \quad u(1, t) = 0$ Initial sinusoidal temperature distribution: $u(x, 0) = 100 \sin(\pi x)$ The exact analytical Fourier series solution is: $u_{\text{exact}}(x, t) = 100 e^{-\pi^2 \alpha t} \sin(\pi x)$

Numerical Grid Discretization:

• Subdivide rod into $N_x = 4$ uniform cells $\implies \Delta x = 0.25\text{ m}$.
• Spatial nodes:
  • $x_0 = 0.0$, $x_1 = 0.25$, $x_2 = 0.50$, $x_3 = 0.75$, $x_4 = 1.0$.
• Initial nodal state vector at $t = 0$:
  • $u_0^0 = 0$
  • $u_1^0 = 100 \sin(\pi/4) = 70.7107$
  • $u_2^0 = 100 \sin(\pi/2) = 100.0000$
  • $u_3^0 = 100 \sin(3\pi/4) = 70.7107$
  • $u_4^0 = 0$

Select stable diffusion parameter $r = \frac{\alpha \Delta t}{\Delta x^2} = 0.40 \le 0.5$. Timestep: $\Delta t = r \frac{\Delta x^2}{\alpha} = 0.40 \times (0.25)^2 / 1.0 = 0.025\text{ s}$.

First Step Computation ($n = 0 \to n = 1$, $t = 0.025\text{ s}$):

Applying FTCS update $u_i^1 = u_i^0 + 0.4 (u_{i+1}^0 - 2u_i^0 + u_{i-1}^0)$:

• Node 1 ($x_1 = 0.25$): $u_1^1 = 70.7107 + 0.4 [100.0000 - 2(70.7107) + 0.0] = 54.1421$

• Node 2 ($x_2 = 0.50$): $u_2^1 = 100.0000 + 0.4 [70.7107 - 2(100.0000) + 70.7107] = 76.5685$

• Node 3 ($x_3 = 0.75$): By physical symmetry: $u_3^1 = u_1^1 = 54.1421$.

Comparison with Exact Benchmark at $t = 0.025\text{ s}$:

$u_{\text{exact}}(0.5, 0.025) = 100 e^{-\pi^2 (1.0)(0.025)} \sin(\pi / 2) = 100 e^{-0.246740} = 78.1343$ Even on an ultra-coarse grid with only 4 cells, the discrete estimate ($76.5685$) achieves a relative error of only $2.0\%$.

Example 2: Pulse Propagation Along a Taut String (1D Wave Equation)

Consider a vibrating string with length $L = 1\text{ m}$, propagation velocity $c = 1\text{ m/s}$, and pinned ends $u(0, t) = u(1, t) = 0$. Initial triangular plucked state starting from rest ($g(x) = 0$):

• $N_x = 4 \implies \Delta x = 0.25\text{ m}$.
• Courant number set to the magic stability limit $C = \frac{c \Delta t}{\Delta x} = 1.0 \implies \Delta t = 0.25\text{ s}$.

Base state at $t = 0$: $u^0 = [0.0, \, 0.5, \, 1.0, \, 0.5, \, 0.0]$

Initial launch step ($g = 0$): $u_i^1 = u_i^0 + \frac{C^2}{2}(u_{i+1}^0 - 2u_i^0 + u_{i-1}^0) = \frac{1}{2}(u_{i+1}^0 + u_{i-1}^0)$

• $u_1^1 = \frac{1}{2}(1.0 + 0.0) = 0.5$
• $u_2^1 = \frac{1}{2}(0.5 + 0.5) = 0.5$
• $u_3^1 = \frac{1}{2}(0.0 + 1.0) = 0.5$ Array at $t = 0.25\text{ s}$: $u^1 = [0.0, \, 0.5, \, 0.5, \, 0.5, \, 0.0]$.

Second step ($t = 0.50\text{ s}$): $u_i^2 = 2u_i^1 - u_i^0 + (u_{i+1}^1 - 2u_i^1 + u_{i-1}^1) = u_{i+1}^1 + u_{i-1}^1 - u_i^0$

• $u_1^2 = 0.5 + 0.0 - 0.5 = 0.0$
• $u_2^2 = 0.5 + 0.5 - 1.0 = 0.0$
• $u_3^2 = 0.0 + 0.5 - 0.5 = 0.0$ Array at $t = 0.50\text{ s}$: $u^2 = [0.0, \, 0.0, \, 0.0, \, 0.0, \, 0.0]$. Because $C = 1.0$ matches the characteristic grid rays, D'Alembert wave splitting occurs without numerical dispersion or artificial numerical damping.

Common Pitfalls in Finite Differences and PDE Solvers

1. Violating the CFL or von Neumann Stability Threshold

• Pitfall: Halving the spatial grid spacing $\Delta x$ to achieve finer spatial resolution while keeping timestep $\Delta t$ constant.
• Consequence: In diffusive problems, the ratio $r = \alpha \Delta t / \Delta x^2$ quadruples. Once $r > 0.5$, amplification factor $|G| > 1$, triggering rapid numerical blow-up with NaN or inf outputs.
• Remedy: Automatically adapt $\Delta t$ dynamically whenever the spatial grid is altered to enforce $r \le 0.5$ or $C \le 1.0$.

2. Using First-Order Approximations for Neumann Boundary Conditions

• Pitfall: Discretizing insulating wall boundary condition $\partial u / \partial x = 0$ as a one-sided forward difference $\frac{u_1 - u_0}{\Delta x} = 0 \implies u_0 = u_1$.
• Consequence: This injects an $\mathcal{O}(\Delta x)$ first-order truncation error at the boundary that diffuses inward, degrading the second-order accuracy of the interior grid.
• Remedy: Implement ghost nodes with second-order central stencils $\mathcal{O}(\Delta x^2)$: $u_{-1} = u_1$.

3. Conflating Physical Dissipation with Artificial Numerical Diffusion

• Pitfall: Adopting first-order upwind stencils for advection or acoustic wave transport.
• Consequence: Leading truncation error behaves mathematically as an artificial diffusion coefficient proportional to $\frac{\Delta x}{2}(1 - C)$, excessively blurring peaks and dispersing wavepackets.
• Remedy: Deploy second-order centered stencils (such as Lax-Wendroff) or high-resolution TVD (Total Variation Diminishing) limiters.

4. Believing Stability Guarantees Convergence Without Consistency

• Pitfall: Assuming that because an algorithm executes without numeric overflow, it converges to the true physical solution.
• Consequence: By Lax's Equivalence Theorem, stability guarantees convergence if and only if the discrete stencil is consistent (local truncation error vanishes as $\Delta x, \Delta t \to 0$). Inconsistent stencils compute mathematically invalid states.
• Remedy: Prove Taylor-expansion consistency analytically before programming numerical equations.

5. Mismanaging Discontinuities and Shock Fronts

• Pitfall: Applying linear centered difference stencils across shock waves or sharp contact discontinuities.
• Consequence: Severe Gibbs-type numerical ringing and spurious oscillations pollute the solutions upstream and downstream of the shock.
• Remedy: Formulate equations in conservative divergence form and apply flux limiters or Godunov-type Riemann solvers.

Physical Applications of Basic Finite Differences

1. Transient Thermal Conduction in Materials and Cryogenics

In mechanical engineering and solid-state physics, spatial temperature distributions across structural components, nuclear fuel rods, and semiconductor dies are computed via 2D/3D finite differences: $\rho c_p \frac{\partial T}{\partial t} = \vec{\nabla} \cdot (k \vec{\nabla} T) + Q(\vec{x}, t)$ Discretization captures transient thermal stress in spacecraft heat shields, microprocessor heatsinks, and liquid-helium cryostats cooling superconducting MRI magnets.

2. Room Acoustics and Computational Seismology

In geophysics, tectonic earthquakes propagate elastic $P$-waves and $S$-waves through heterogeneous geological strata, simulated using time-domain finite differences (FDTD): $\rho \frac{\partial^2 \vec{u}}{\partial t^2} = (\lambda + \mu) \vec{\nabla}(\vec{\nabla} \cdot \vec{u}) + \mu \nabla^2 \vec{u}$ These models forecast ground shake intensity for civil engineering and optimize acoustic reflection in architectural concert hall designs.

3. Computational Electromagnetics (The Yee FDTD Algorithm)

Kane Yee's FDTD formulation discretizes Maxwell's curl equations by staggering electric $\vec{E}$ and magnetic $\vec{H}$ field components across interlocking spatial grids: $\vec{\nabla} \times \vec{E} = -\frac{\partial \vec{B}}{\partial t}, \quad \vec{\nabla} \times \vec{H} = \vec{J} + \frac{\partial \vec{D}}{\partial t}$ This algorithm serves as the primary tool for modeling 5G/6G antennas, radar cross-sections, dielectric waveguides, and photonic crystal bandgaps.

4. Time-Dependent Schrödinger Equation in Quantum Systems

In condensed matter physics and quantum device modeling, electron wavepacket tunneling across potential barriers is computed via: $i \hbar \frac{\partial \psi}{\partial t} = -\frac{\hbar^2}{2m} \frac{\partial^2 \psi}{\partial x^2} + V(x) \psi$ The implicit Crank-Nicolson method is universally adopted here because its unitary Cayley operator strictly conserves total probability $\int |\psi|^2 dx = 1$.

5. Hydrodynamic Contaminant Advection-Diffusion Transport

In environmental fluid dynamics and oceanography, chemical dispersion through groundwater aquifers, rivers, or the atmosphere couples velocity advection with turbulent diffusion: $\frac{\partial C}{\partial t} + \vec{v} \cdot \vec{\nabla} C = D \nabla^2 C + S$ Finite difference models predict contaminant plume migration following industrial spills, oil leaks, or urban particulate pollution dispersion.

Practice Problems: Finite Differences and Grid Stability

Problem 1: Catastrophic Instability from Diffusion Boundary Violation

A 1D thermal simulation with diffusivity $\alpha = 2.0 \times 10^{-4}\text{ m}^2/\text{s}$ uses cell size $\Delta x = 0.01\text{ m}$.

1. Find the critical timestep $\Delta t_{\text{crit}}$ above which explicit FTCS becomes unstable.
2. If an operator mistakenly runs the code with $\Delta t = 0.3\text{ s}$, demonstrate using von Neumann amplification what happens to the highest spatial frequency mode ($k \Delta x = \pi$).

Step-by-Step Solution:

1. Stability requires: $r = \frac{\alpha \Delta t}{\Delta x^2} \le 0.5$ $\Delta t_{\text{crit}} = \frac{\Delta x^2}{2 \alpha} = \frac{(0.01)^2}{2(2.0 \times 10^{-4})} = 0.25\text{ s}$
2. At $\Delta t = 0.3\text{ s}$: $r = \frac{(2.0 \times 10^{-4})(0.3)}{10^{-4}} = 0.60 > 0.50$ The amplification factor for the grid Nyquist mode ($k \Delta x = \pi$) is: $G = 1 - 4r \sin^2(\pi / 2) = 1 - 4(0.60)(1) = -1.40$ Since $|G| = 1.40 > 1$, this mode grows in amplitude by $40\%$ per timestep with alternating signs:
   • After 10 steps: factor of $(1.40)^{10} \approx 28.9$.
   • After 50 steps: factor of $(1.40)^{50} \approx 1.8 \times 10^7$. The field degenerates into an unphysical high-frequency zigzag blowup.

Problem 2: Insulated Wall Implementation via Ghost Nodes

On a thermal rod of length $L = 0.5\text{ m}$, the left boundary $x = 0$ is thermally insulated (zero heat flux): $\left.\frac{\partial u}{\partial x}\right|_{x=0} = 0$

1. Utilizing a second-order central difference stencil with fictitious ghost node $x_{-1} = -\Delta x$, derive the update formula for boundary temperature $u_0^{n+1}$.
2. Given $\Delta x = 0.1\text{ m}$, $\alpha = 1\text{ m}^2/\text{s}$, $\Delta t = 0.004\text{ s}$, with current values $u_0^n = 100^\circ\text{C}$ and $u_1^n = 90^\circ\text{C}$, compute $u_0^{n+1}$.

Step-by-Step Solution:

1. Centered Neumann boundary condition: $\frac{u_1^n - u_{-1}^n}{2 \Delta x} = 0 \implies u_{-1}^n = u_1^n$ Standard internal stencil applied at $i = 0$: $u_0^{n+1} = u_0^n + r (u_1^n - 2u_0^n + u_{-1}^n)$ Substituting ghost node identity: $u_0^{n+1} = u_0^n + 2r (u_1^n - u_0^n)$
2. Parameter $r$: $r = \frac{(1)(0.004)}{(0.1)^2} = 0.40$ Update: $u_0^{n+1} = 100 + 2(0.40)(90 - 100) = 100 - 8.0 = 92.0^\circ\text{C}$ Boundary temperature cools as heat conducts into the colder adjacent internal node.

Problem 3: Numerical Dispersion in Acoustic Wavefields

An acoustic solver with speed of sound $c = 340\text{ m/s}$ runs with $\Delta x = 0.1\text{ m}$ and $\Delta t = 0.0002\text{ s}$.

1. Verify if the CFL condition is satisfied.
2. Explain why operating at a Courant number substantially below 1 (e.g., $C = 0.3$) induces artificial numerical dispersion that distorts sharp wavefronts.

Step-by-Step Solution:

1. Courant number: $C = \frac{c \Delta t}{\Delta x} = \frac{(340)(0.0002)}{0.1} = 0.68 \le 1.0$ The solver is unconditionally bounded and stable.
2. When $C < 1$, the discrete numerical dispersion relation deviates from the linear physical law $\omega = ck$: $\sin\left(\frac{\omega_{\text{num}} \Delta t}{2}\right) = C \sin\left(\frac{k \Delta x}{2}\right)$ Consequently, numerical phase velocity $v_{\text{phase}} = \omega_{\text{num}} / k$ becomes wavenumber-dependent. Short wavelengths lag behind long wavelengths, generating spurious oscillations trailing shock or impulse fronts. Only at $C = 1.0$ is numerical dispersion identically zero in 1D.