4. 一维波动方程显式差分格式

一维波动方程是描述弦振动、杆中纵波、声波在管道中传播等物理过程的基本模型，其标准形式为

∂²u/∂t² = c² ∂²u/∂x²

其中 u(x,t) 为位移或声压等物理量，c 为波速。要在计算机上求解它，需要把连续的时间和空间离散成网格，用差分商代替偏导数。本节聚焦最经典、最常用的显式格式：空间二阶中心差分加时间二阶中心差分，并讨论它的推导、稳定性条件（CFL 条件）与数值色散特性。

一、网格划分与差分近似

先约定离散方式：空间方向取步长 Δx，节点编号 i = 0,1,…,N；时间方向取步长 Δt，时间层编号 n = 0,1,2,…。记 u_iⁿ ≈ u(iΔx, nΔt)。

对时间二阶导数，在 t = nΔt 处用中心差分：

∂²u/∂t² ≈ (u_iⁿ⁺¹ − 2u_iⁿ + u_iⁿ⁻¹) / Δt²

对空间二阶导数，在 x = iΔx 处同样用中心差分：

∂²u/∂x² ≈ (u_i+1ⁿ − 2u_iⁿ + u_i−1ⁿ) / Δx²

两者截断误差均为 O(Δt²) 和 O(Δx²)，因此整体是二阶精度格式。

二、显式格式的推导

把上述两个近似代入波动方程，得到

(u_iⁿ⁺¹ − 2u_iⁿ + u_iⁿ⁻¹)/Δt² = c² (u_i+1ⁿ − 2u_iⁿ + u_i−1ⁿ)/Δx²

引入网格比（库朗数）r = cΔt/Δx，整理后把未知的 n+1 层单独解出：

u_iⁿ⁺¹ = 2u_iⁿ − u_iⁿ⁻¹ + r² (u_i+1ⁿ − 2u_iⁿ + u_i−1ⁿ)

这就是一维波动方程的显式三点格式。它之所以叫“显式”，是因为 u_iⁿ⁺¹ 可以直接由前两个时间层上的已知值算出，不需要解线性方程组。计算时只需保存 n 层和 n−1 层两行数据，逐点推进即可，内存占用小、实现简单。

启动计算需要两个初始层：n=0 层由初始位移 u(x,0) 给出；n=1 层通常由初始速度 u_t(x,0) 结合一次向前差分构造，例如 u_i¹ = u_i⁰ + Δt·v_i⁰ + (r²/2)(u_i+1⁰ − 2u_i⁰ + u_i−1⁰)，以保证整体仍为二阶精度。边界条件则按物理问题设定：固定端取 u=0，自由端可用虚拟节点镜像处理。

三、稳定性条件（CFL 条件）

显式格式并非无条件可用。用冯·诺依曼稳定性分析，设误差分量为 u_iⁿ = λⁿ e^{jk iΔx}，代入格式得到增长因子 λ 满足的二次方程：

λ² − 2(1 − 2r² sin²(θ/2))λ + 1 = 0，θ = kΔx

要使 |λ| ≤ 1 对所有波数成立，必须满足

r = cΔt/Δx ≤ 1，即 Δt ≤ Δx/c

这就是著名的 CFL（Courant–Friedrichs–Lewy）条件。它的物理含义很直观：在一个时间步内，波传播的距离 cΔt 不能超过一个空间网格 Δx，否则数值依赖域无法覆盖真实依赖域，计算必然发散。实际编程时通常取 r = 0.8～0.9 留出安全裕量，而不是贴着 r=1 运行。

四、数值色散分析

即使格式稳定，数值解与真实解之间仍存在相位误差，这就是数值色散。把离散格式的色散关系写出来：

sin(ωΔt/2) = r sin(kΔx/2)

真实连续介质满足 ω = ck，即数值相速度应为 c。但由上式可得数值相速度

c_num = ω/k = (2/(kΔt)) arcsin(r sin(kΔx/2))

可见 c_num 依赖波数 k，不同频率的波以不同速度传播，短波（kΔx 接近 π）的相速度明显偏低，长波则接近真实波速。这意味着一个初始脉冲在传播中会逐渐“散开”，出现伪振荡和波形畸变。网格越粗（Δx 越大），色散越严重；r 越接近 1，色散误差越小。因此工程上常要求每个波长至少布置 10～20 个网格点，以把色散控制在可接受范围。

五、一个可执行的算例

设弦长 L = 1 m，波速 c = 100 m/s，两端固定，初始位移为 u(x,0) = sin(πx/L)，初始速度为 0。取 N = 100 个空间区间，则 Δx = 0.01 m。按 CFL 条件 Δt ≤ Δx/c = 1×10⁻⁴ s，取 r = 0.9，则 Δt = 0.9×10⁻⁴ s。核心推进循环可写成：

r = 0.9; r2 = r*r
for n in range(1, nsteps):
    for i in range(1, N):
        u_new[i] = 2*u[i] - u_old[i] + r2*(u[i+1] - 2*u[i] + u[i-1])
    u_new[0] = 0; u_new[N] = 0      # 固定端边界
    u_old, u, u_new = u, u_new, u_old

运行后应观察到：波形以周期 T = 2L/c = 0.02 s 往复振荡，振幅基本守恒；若把 r 改为 1.05 再运行，几十步内位移幅值迅速增大直至溢出，这正是违反 CFL 条件的典型表现。若把 N 减到 20（Δx = 0.05 m），虽然仍满足 CFL，但波形在传播数个周期后明显变形、出现高频毛刺，这就是色散误差的直接体现。

六、实践要点小结

• 先由 Δx 和 c 确定 Δt 上限，再取 r ≈ 0.8～0.9，不要贴着 1 运行。
• 初始速度层要用二阶精度公式构造，否则会引入 O(Δt) 的启动误差。
• 每个波长至少 10～20 个网格点，以控制色散；对长时程传播问题应做网格收敛性检验。
• 用能量（离散形式的动能加势能）随时间的变化监控稳定性：稳定格式下能量应在小范围内振荡而非单调增长。

掌握这套显式格式后，读者已具备求解一维波动问题的最小完整工具链：网格、差分、推进、边界与稳定性判据。后续更复杂的格式与多维问题，都可以在此基础上扩展。