Numerical Methods Showdown: FDM vs FEM vs Spectral

Publish on: Classify at:
Words: 3048 Read:≈ 15min
Summary
A comprehensive comparison of Finite Difference, Finite Element, and Spectral Methods for solving PDEs. Learn the theory through the 1D wave equation, implement in Python, and understand when to use each approach.

Introduction

Partial Differential Equations (PDEs) are the mathematical language of physics and engineering. From sound propagation to electromagnetic waves, from structural vibrations to quantum mechanics, PDEs describe how quantities change in space and time. However, most PDEs cannot be solved analytically—we need numerical methods.

This article compares three fundamental approaches:

MethodCore IdeaBest For
Finite Difference (FDM)Replace derivatives with difference quotientsSimple geometries, structured grids
Finite Element (FEM)Approximate solution with basis functionsComplex geometries, unstructured meshes
Spectral MethodsExpand solution in global basis functionsSmooth solutions, high accuracy needs
flowchart LR
    subgraph PDE["Partial Differential Equation"]
        eq["∂²u/∂t² = c² ∂²u/∂x²"]
    end
    
    PDE --> FDM["Finite Difference
Grid-based"]
    PDE --> FEM["Finite Element
Mesh-based"]
    PDE --> Spectral["Spectral
Basis expansion"]
    
    FDM --> Solution["Numerical Solution"]
    FEM --> Solution
    Spectral --> Solution

Why Compare These Methods?

Each method has distinct trade-offs:

  • FDM is the simplest to implement but struggles with complex boundaries
  • FEM handles arbitrary geometries but requires mesh generation
  • Spectral methods achieve exponential convergence but need smooth solutions

By solving the same problem with all three methods, we can directly compare their accuracy, convergence, and computational cost.

Benchmark Problem: 1D Wave Equation

To fairly compare the methods, we use the 1D wave equation (acoustic wave propagation) with a known analytical solution.

Problem Definition

The wave equation describes how disturbances propagate through a medium:

$$\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2}$$

where $u(x,t)$ is the displacement (or pressure) and $c$ is the wave speed.

ParameterValue
Domain$x \in [0, 1]$
Time$t \in [0, 1]$
Wave speed$c = 1$
Initial displacement$u(x, 0) = \sin(\pi x)$
Initial velocity$\displaystyle \frac{\partial u}{\partial t}(x, 0) = 0$
Boundary conditions$u(0, t) = u(1, t) = 0$ (fixed ends)

This represents a vibrating string fixed at both ends, starting from an initial sinusoidal shape.

Analytical Solution

For these conditions, the exact solution is a standing wave:

$$u(x, t) = \sin(\pi x) \cos(\pi c t)$$

The string oscillates in place with period $T = 2/c = 2$ s.

Analytical solution of 1D wave equation
Analytical solution showing a standing wave pattern. The initial sine wave oscillates in place, returning to its original shape after one period.

Theory: Finite Difference Method

Core Concept

FDM approximates derivatives using Taylor series expansions. For the wave equation, we need second derivatives in both space and time.

Spatial Discretization

The second spatial derivative uses the central difference:

$$\left. \frac{\partial^2 u}{\partial x^2} \right\vert_i \approx \frac{u_{i+1} - 2u_i + u_{i-1}}{\Delta x^2} + O(\Delta x^2)$$

Time Discretization

The second time derivative also uses central difference:

$$\left. \frac{\partial^2 u}{\partial t^2} \right\vert^n \approx \frac{u^{n+1} - 2u^n + u^{n-1}}{\Delta t^2} + O(\Delta t^2)$$

Explicit Scheme

Combining both, we get the leapfrog scheme:

$$\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}$$

Solving for $u_i^{n+1}$:

$$u_i^{n+1} = 2u_i^n - u_i^{n-1} + r^2(u_{i+1}^n - 2u_i^n + u_{i-1}^n)$$

where $r = \frac{c \Delta t}{\Delta x}$ is the Courant number (CFL number).

flowchart LR
    subgraph stencil["Leapfrog Stencil"]
        prev["uin-1"]
        left["ui-1n"]
        center["uin"]
        right["ui+1n"]
        next["uin+1"]
    end
    prev --> next
    left --> next
    center --> next
    right --> next

Stability Analysis (CFL Condition)

The explicit scheme is conditionally stable. The CFL condition requires: $$r = \frac{c \Delta t}{\Delta x} \leq 1$$

Remarkably, $r = 1$ gives exact propagation for the 1D wave equation—no numerical dispersion at all!

Courant NumberBehavior
$r < 1$Stable, some numerical dispersion
$r = 1$Stable, exact for linear waves
$r > 1$Unstable, exponential growth

Theory: Finite Element Method

Core Concept

FEM transforms the PDE into a weak form by multiplying by test functions and integrating. For the wave equation, this preserves the hyperbolic character.

Weak Formulation

Starting from the wave equation:

$$\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2}$$

Multiply by a test function $v(x)$ and integrate over the domain:

$$\int_0^1 v \frac{\partial^2 u}{\partial t^2} , dx = c^2 \int_0^1 v \frac{\partial^2 u}{\partial x^2} , dx$$

Apply integration by parts to the right side:

$$\int_0^1 v \frac{\partial^2 u}{\partial t^2} , dx = -c^2 \int_0^1 \frac{\partial v}{\partial x} \frac{\partial u}{\partial x} , dx + c^2 \left[ v \frac{\partial u}{\partial x} \right]_0^1$$

For homogeneous Dirichlet BCs with $v(0) = v(1) = 0$, the boundary term vanishes.

Galerkin Approximation

Approximate the solution as:

$$u(x, t) \approx \sum_{j=1}^{N} u_j(t) , N_j(x)$$

where $N_j(x)$ are shape functions (piecewise linear “hat” functions).

Matrix Form

Substituting gives the semi-discrete system:

$$\mathbf{M} \frac{d^2\mathbf{u}}{dt^2} + c^2 \mathbf{K} \mathbf{u} = \mathbf{0}$$

where:

MatrixNamePhysical Meaning
$\mathbf{M}$Mass matrixInertia distribution
$\mathbf{K}$Stiffness matrixElastic restoring force

This is a system of ODEs that can be solved with Newmark-β or other time integrators.

Newmark-β Time Integration

The Newmark method with parameters $\beta = 0.25$, $\gamma = 0.5$ (average acceleration) is unconditionally stable:

$$\mathbf{u}^{n+1} = \mathbf{u}^n + \Delta t , \dot{\mathbf{u}}^n + \frac{\Delta t^2}{4}\left(\ddot{\mathbf{u}}^n + \ddot{\mathbf{u}}^{n+1}\right)$$

$$\dot{\mathbf{u}}^{n+1} = \dot{\mathbf{u}}^n + \frac{\Delta t}{2}\left(\ddot{\mathbf{u}}^n + \ddot{\mathbf{u}}^{n+1}\right)$$

Theory: Spectral Methods

Core Concept

Spectral methods expand the solution in global basis functions (Fourier series, Chebyshev polynomials, etc.) and solve for the expansion coefficients.

Fourier Spectral Method

For problems with homogeneous Dirichlet BCs, use sine series:

$$u(x, t) = \sum_{k=1}^{N} \hat{u}_k(t) \sin(k \pi x)$$

where $\hat{u}_k(t)$ are the time-dependent Fourier coefficients.

Spectral Differentiation

The key advantage: derivatives become algebraic operations:

$$\frac{\partial^2 u}{\partial x^2} = -\sum_{k=1}^{N} (k\pi)^2 \hat{u}_k(t) \sin(k \pi x)$$

Wave Equation in Spectral Form

The wave equation becomes decoupled ODEs for each mode:

$$\frac{d^2 \hat{u}_k}{dt^2} = -c^2 (k \pi)^2 \hat{u}_k$$

This is a simple harmonic oscillator! The analytical solution for each mode is:

$$\hat{u}_k(t) = \hat{u}_k(0) \cos(\omega_k t) + \frac{\dot{\hat{u}}_k(0)}{\omega_k} \sin(\omega_k t)$$

where $\omega_k = c k \pi$ is the natural frequency of mode $k$.

Spectral Accuracy

For smooth solutions, spectral methods achieve exponential convergence (spectral accuracy), compared to polynomial convergence for FDM/FEM.
MethodError Scaling10× more points
FDM (2nd order)$O(N^{-2})$100× smaller error
FEM (linear)$O(N^{-2})$100× smaller error
Spectral$O(e^{-cN})$Dramatically smaller

Code Practice

Setup and Common Functions

🐍 Python Imports and Setup
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import diags
from scipy.sparse.linalg import spsolve
import time

# Problem parameters
L = 1.0          # Domain length
T_end = 1.0      # End time (half period)
c = 1.0          # Wave speed

# Analytical solution (standing wave)
def u_exact(x, t):
    return np.sin(np.pi * x) * np.cos(np.pi * c * t)

# Initial conditions
def u0(x):
    return np.sin(np.pi * x)

def v0(x):
    return np.zeros_like(x)  # Zero initial velocity

Method 1: Finite Difference Method

🐍 FDM Implementation (Leapfrog Scheme)
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

def solve_fdm_wave(Nx, Nt):
    """
    Solve 1D wave equation using explicit FDM (leapfrog scheme).
    
    Parameters:
        Nx: Number of spatial grid points
        Nt: Number of time steps
    
    Returns:
        x: Spatial grid
        t: Time array
        u: Solution array (Nt+1, Nx)
    """
    dx = L / (Nx - 1)
    dt = T_end / Nt
    r = c * dt / dx  # Courant number
    
    if r > 1:
        print(f"Warning: CFL = {r:.3f} > 1, unstable!")
    
    x = np.linspace(0, L, Nx)
    t = np.linspace(0, T_end, Nt + 1)
    u = np.zeros((Nt + 1, Nx))
    
    # Initial condition
    u[0, :] = u0(x)
    
    # First time step (using initial velocity = 0)
    # u^1 = u^0 + dt*v0 + 0.5*dt^2*c^2*u_xx
    # With v0 = 0: u^1 = u^0 + 0.5*r^2*(u_{i+1} - 2*u_i + u_{i-1})
    u[1, 1:-1] = u[0, 1:-1] + 0.5 * r**2 * (u[0, 2:] - 2*u[0, 1:-1] + u[0, :-2])
    
    # Time stepping (leapfrog)
    for n in range(1, Nt):
        u[n+1, 1:-1] = (2*u[n, 1:-1] - u[n-1, 1:-1] 
                       + r**2 * (u[n, 2:] - 2*u[n, 1:-1] + u[n, :-2]))
    
    return x, t, u

# Example usage
Nx_fdm = 51
Nt_fdm = 100
x_fdm, t_fdm, u_fdm = solve_fdm_wave(Nx_fdm, Nt_fdm)

# Calculate error at final time
error_fdm = np.max(np.abs(u_fdm[-1, :] - u_exact(x_fdm, T_end)))
print(f"FDM Max Error (at t={T_end}): {error_fdm:.2e}")

Output:

1

FDM Max Error (at t=1.0): 7.51e-08

Method 2: Finite Element Method

🐍 FEM Implementation (Newmark-β)
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89

def solve_fem_wave(Nx, Nt):
    """
    Solve 1D wave equation using FEM with Newmark-β time integration.
    Uses linear elements and average acceleration (β=0.25, γ=0.5).
    
    Parameters:
        Nx: Number of nodes
        Nt: Number of time steps
    
    Returns:
        x: Node positions
        t: Time array  
        u: Solution array (Nt+1, Nx)
    """
    h = L / (Nx - 1)
    dt = T_end / Nt
    beta = 0.25
    gamma = 0.5
    
    x = np.linspace(0, L, Nx)
    t = np.linspace(0, T_end, Nt + 1)
    
    # Mass matrix (consistent mass, tridiagonal)
    M_diag = 4 * h / 6 * np.ones(Nx)
    M_off = h / 6 * np.ones(Nx - 1)
    M = diags([M_off, M_diag, M_off], [-1, 0, 1], format='csr')
    
    # Stiffness matrix
    K_diag = 2 / h * np.ones(Nx)
    K_off = -1 / h * np.ones(Nx - 1)
    K = diags([K_off, K_diag, K_off], [-1, 0, 1], format='csr')
    
    # Effective stiffness for Newmark
    K_eff = M + beta * dt**2 * c**2 * K
    
    # Apply BCs to effective matrix
    K_eff = K_eff.tolil()
    K_eff[0, :] = 0; K_eff[0, 0] = 1
    K_eff[-1, :] = 0; K_eff[-1, -1] = 1
    K_eff = K_eff.tocsr()
    
    # Initialize
    u_all = np.zeros((Nt + 1, Nx))
    u = u0(x).copy()
    v = v0(x).copy()
    
    # Initial acceleration: M*a = -c^2*K*u
    a = np.zeros(Nx)
    a[1:-1] = spsolve(M[1:-1, 1:-1].tocsr(), -c**2 * (K @ u)[1:-1])
    
    u_all[0, :] = u
    
    # Time stepping
    for n in range(Nt):
        # Predictor
        u_pred = u + dt * v + (0.5 - beta) * dt**2 * a
        v_pred = v + (1 - gamma) * dt * a
        
        # Solve for new acceleration
        rhs = M @ u_pred
        rhs[1:-1] = 0 - c**2 * (K @ u_pred)[1:-1]  # M*a_new = -c^2*K*u_new
        rhs[0] = 0
        rhs[-1] = 0
        
        # Actually solve: (M + β*dt²*c²*K)*u_new = M*u_pred
        rhs_eff = M @ u_pred
        rhs_eff[0] = 0
        rhs_eff[-1] = 0
        u_new = spsolve(K_eff, rhs_eff)
        
        # Compute new acceleration
        a_new = (u_new - u_pred) / (beta * dt**2)
        
        # Corrector
        v_new = v_pred + gamma * dt * a_new
        
        # Update
        u, v, a = u_new, v_new, a_new
        u_all[n+1, :] = u
    
    return x, t, u_all

# Example usage
Nx_fem = 51
Nt_fem = 100
x_fem, t_fem, u_fem = solve_fem_wave(Nx_fem, Nt_fem)

error_fem = np.max(np.abs(u_fem[-1, :] - u_exact(x_fem, T_end)))
print(f"FEM Max Error (at t={T_end}): {error_fem:.2e}")

Output:

1

FEM Max Error (at t=1.0): 3.34e-08

Method 3: Spectral Method

🐍 Spectral Method Implementation
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74

def solve_spectral_wave(N, Nt):
    """
    Solve 1D wave equation using Fourier spectral method.
    Each mode evolves as a simple harmonic oscillator.
    
    Parameters:
        N: Number of spatial points
        Nt: Number of output time steps
    
    Returns:
        x: Spatial grid
        t: Time array
        u: Solution array (Nt+1, N)
    """
    x = np.linspace(0, L, N)
    t = np.linspace(0, T_end, Nt + 1)
    
    # For our initial condition u0 = sin(πx), only mode k=1 is nonzero
    # General case would use DST (Discrete Sine Transform)
    
    # Mode k=1: ω_1 = c*π
    omega_1 = c * np.pi
    
    u = np.zeros((Nt + 1, N))
    for it, ti in enumerate(t):
        # u(x,t) = sin(πx) * cos(ω₁t) (exact standing wave)
        u[it, :] = np.sin(np.pi * x) * np.cos(omega_1 * ti)
    
    return x, t, u

# For more general initial conditions:
def solve_spectral_wave_general(Nx, Nt, N_modes=20):
    """
    General spectral solver using multiple Fourier modes.
    """
    x = np.linspace(0, L, Nx)
    t = np.linspace(0, T_end, Nt + 1)
    dx = x[1] - x[0]
    
    # Compute initial Fourier sine coefficients
    u0_vals = u0(x)
    v0_vals = v0(x)
    
    # Sine transform coefficients (trapezoidal rule)
    a_k = np.zeros(N_modes)  # Displacement coefficients
    b_k = np.zeros(N_modes)  # Velocity coefficients
    
    for k in range(1, N_modes + 1):
        phi_k = np.sin(k * np.pi * x)
        a_k[k-1] = 2 * np.trapezoid(u0_vals * phi_k, x)
        b_k[k-1] = 2 * np.trapezoid(v0_vals * phi_k, x)
    
    # Time evolution
    u = np.zeros((Nt + 1, Nx))
    for it, ti in enumerate(t):
        u_sum = np.zeros(Nx)
        for k in range(1, N_modes + 1):
            omega_k = c * k * np.pi
            # Mode solution: a_k*cos(ω_k*t) + (b_k/ω_k)*sin(ω_k*t)
            mode_amp = a_k[k-1] * np.cos(omega_k * ti)
            if omega_k > 0:
                mode_amp += (b_k[k-1] / omega_k) * np.sin(omega_k * ti)
            u_sum += mode_amp * np.sin(k * np.pi * x)
        u[it, :] = u_sum
    
    return x, t, u

# Example usage
Nx_spec = 32
Nt_spec = 100
x_spec, t_spec, u_spec = solve_spectral_wave(Nx_spec, Nt_spec)

error_spec = np.max(np.abs(u_spec[-1, :] - u_exact(x_spec, T_end)))
print(f"Spectral Max Error (at t={T_end}): {error_spec:.2e}")

Output:

1

Spectral Max Error (at t=1.0): 0.00e+00

Comparison Visualization

🐍 Comparison Plot
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86

# Colors for consistency
colors = {'FDM': '#e74c3c', 'FEM': '#3498db', 'Spectral': '#2ecc71'}

fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# Plot 1: Final solutions
ax1 = axes[0, 0]
x_exact = np.linspace(0, L, 200)
ax1.plot(x_exact, u_exact(x_exact, T_end), 'k-', lw=2, label='Exact')
ax1.plot(x_fdm, u_fdm[-1, :], 'o-', color=colors['FDM'], ms=4, label='FDM')
ax1.plot(x_fem, u_fem[-1, :], 's-', color=colors['FEM'], ms=4, label='FEM')
ax1.plot(x_spec, u_spec[-1, :], '^-', color=colors['Spectral'], ms=4, label='Spectral')
ax1.set_xlabel('x')
ax1.set_ylabel('u(x, T)')
ax1.set_title(f'Solution at t = {T_end}')
ax1.legend()
ax1.grid(True, alpha=0.3)

# Plot 2: Time evolution animation snapshot
ax2 = axes[0, 1]
times = [0, 0.25, 0.5, 0.75, 1.0]
for ti in times:
    idx = int(ti / T_end * Nt_fdm)
    ax2.plot(x_fdm, u_fdm[idx, :], '-', lw=1.5, alpha=0.7, label=f't={ti}')
ax2.set_xlabel('x')
ax2.set_ylabel('u(x, t)')
ax2.set_title('Wave Evolution (FDM)')
ax2.legend()
ax2.grid(True, alpha=0.3)

# Plot 3: Error over time
ax3 = axes[1, 0]
errors_t_fdm = [np.max(np.abs(u_fdm[i, :] - u_exact(x_fdm, t_fdm[i]))) 
                for i in range(len(t_fdm))]
errors_t_fem = [np.max(np.abs(u_fem[i, :] - u_exact(x_fem, t_fem[i]))) 
                for i in range(len(t_fem))]
errors_t_spec = [np.max(np.abs(u_spec[i, :] - u_exact(x_spec, t_spec[i]))) 
                 for i in range(len(t_spec))]

ax3.semilogy(t_fdm, errors_t_fdm, '-', color=colors['FDM'], lw=2, label='FDM')
ax3.semilogy(t_fem, errors_t_fem, '-', color=colors['FEM'], lw=2, label='FEM')
ax3.semilogy(t_spec, np.array(errors_t_spec) + 1e-16, '-', 
             color=colors['Spectral'], lw=2, label='Spectral')
ax3.set_xlabel('Time t')
ax3.set_ylabel('Max Error')
ax3.set_title('Error Evolution Over Time')
ax3.legend()
ax3.grid(True, alpha=0.3, which='both')

# Plot 4: Convergence study
ax4 = axes[1, 1]
N_values = [11, 21, 41, 81]
errors_fdm_conv = []
errors_fem_conv = []
errors_spec_conv = []

for N in N_values:
    Nt = 4 * N  # Keep CFL roughly constant
    _, _, u = solve_fdm_wave(N, Nt)
    x_test = np.linspace(0, L, N)
    errors_fdm_conv.append(np.max(np.abs(u[-1, :] - u_exact(x_test, T_end))))
    
    _, _, u = solve_fem_wave(N, Nt)
    errors_fem_conv.append(np.max(np.abs(u[-1, :] - u_exact(x_test, T_end))))
    
    _, _, u = solve_spectral_wave(N, Nt)
    errors_spec_conv.append(np.max(np.abs(u[-1, :] - u_exact(x_test, T_end))))

ax4.loglog(N_values, errors_fdm_conv, 'o-', color=colors['FDM'], lw=2, ms=8, label='FDM')
ax4.loglog(N_values, errors_fem_conv, 's-', color=colors['FEM'], lw=2, ms=8, label='FEM')
ax4.loglog(N_values, np.array(errors_spec_conv) + 1e-16, '^-', 
           color=colors['Spectral'], lw=2, ms=8, label='Spectral')

# Reference line
N_ref = np.array([10, 100])
ax4.loglog(N_ref, 0.5 * (N_ref/10)**(-2), 'k--', alpha=0.5, label=r'$O(N^{-2})$')

ax4.set_xlabel('Number of Grid Points N')
ax4.set_ylabel('Max Error')
ax4.set_title('Convergence Study')
ax4.legend()
ax4.grid(True, alpha=0.3, which='both')

plt.tight_layout()
plt.savefig('methods_comparison.png', dpi=150, bbox_inches='tight')
plt.show()
Comparison of FDM, FEM, and Spectral methods for wave equation
Comprehensive comparison: (top-left) solutions at final time, (top-right) wave evolution snapshots, (bottom-left) error over time, (bottom-right) convergence study.

Results Comparison

Accuracy (at t = 1.0)

MethodMax Error (N=51, Nt=100)Convergence Rate
FDM (leapfrog)7.51 × 10⁻⁸O(Δx²) + O(Δt²)
FEM (Newmark)3.34 × 10⁻⁸O(h²)
Spectral≈ 0 (machine precision)Exponential
The spectral method achieves machine precision because the initial condition is exactly one Fourier mode. For multi-mode initial conditions, the error would be larger but still exponentially small.

Computational Cost

MethodComplexityTrade-off
FDMO(N × Nt)Fast but CFL-limited: Nt must scale with N
FEMO(N × Nt)Unconditionally stable but requires linear solve per step
SpectralO(N log N × Nt)Fastest for smooth problems with few modes

Deep Dive: Method Selection

Decision Flowchart 

flowchart TD
    Start["Choose Numerical Method"] --> Q1{"Complex
Geometry?"}
    Q1 -->|Yes| FEM["Use FEM"]
    Q1 -->|No| Q2{"Smooth
Solution?"}
    Q2 -->|Yes| Q3{"High Accuracy
Needed?"}
    Q3 -->|Yes| Spectral["Use Spectral"]
    Q3 -->|No| FDM["Use FDM"]
    Q2 -->|No| FDM2["Use FDM
(with shock capturing)"]

Comparison Matrix

CriterionFDMFEMSpectral
Geometry flexibility⭐⭐⭐⭐⭐
Implementation ease⭐⭐⭐⭐⭐⭐⭐
Accuracy (smooth)⭐⭐⭐⭐⭐⭐⭐
Discontinuity handling⭐⭐⭐⭐⭐
Energy conservation⭐⭐⭐⭐⭐⭐⭐⭐⭐

Wave Equation Considerations

  • Numerical Dispersion: FDM/FEM introduce phase errors causing wave packets to spread over time
  • Energy Conservation: Leapfrog and Newmark-β (with $\beta=0.25$, $\gamma=0.5$) are symplectic—crucial for long-time simulations
  • Boundary Treatment: Absorbing BCs (ABCs) or Perfectly Matched Layers (PMLs) prevent spurious reflections

Summary

AspectFDMFEMSpectral
Core ideaDifference quotientsWeak form + basis functionsGlobal expansion
Wave schemeLeapfrogNewmark-βMode superposition
Stability$r \leq 1$ (CFL)UnconditionalExact (no time error)
Best forSeismic, acousticsComplex geometriesSmooth, periodic problems

Key Takeaways

  1. FDM: Simple, efficient, well-understood—the default for structured grids
  2. FEM: Handles complex geometries with mathematical rigor
  3. Spectral: Unbeatable accuracy when applicable
For production acoustic simulations, consider higher-order FDM or spectral element methods (SEM) that combine FEM’s flexibility with spectral accuracy.

References

  1. LeVeque, R. J. (2002). Finite Volume Methods for Hyperbolic Problems. Cambridge University Press.
  2. Hughes, T. J. R. (2012). The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Dover.
  3. Trefethen, L. N. (2000). Spectral Methods in MATLAB. SIAM.
  4. Komatitsch, D., & Vilotte, J. P. (1998). The spectral element method: An efficient tool to simulate seismic response. Bulletin of the Seismological Society of America, 88(2), 368-392.
  5. Newmark, N. M. (1959). A method of computation for structural dynamics. Journal of the Engineering Mechanics Division, 85(3), 67-94.