1. 波动方程与有限差分法概述

波动方程是描述波在介质中传播规律的核心偏微分方程。在一维均匀介质中，若记位移为 u(x,t)，波速为 c，则齐次波动方程写作：

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

这个式子表达了一个非常直观的物理事实：介质中某一点的加速度（对时间的二阶导数）正比于该点附近位移分布的弯曲程度（对空间的二阶导数）。弯曲越剧烈，回复力越大，质点被拉回平衡位置的趋势越强，波便以速度 c 向两侧传播。波速 c 由介质本身决定，例如弦上横波 c = √(T/ρ)，其中 T 为张力、ρ 为线密度；声波在流体中 c = √(K/ρ)，K 为体积模量。方程本身不含耗散项，因此理想波动方程描述的是无衰减、无频散的传播过程。

要得到唯一解，仅有方程还不够，必须配合初始条件与边界条件。初始条件通常给出初始位移 u(x,0) 与初始速度 ∂u/∂t|_t=0；边界条件则描述端点行为，常见的有固定端（Dirichlet，u=0）、自由端（Neumann，∂u/∂x=0）以及吸收边界（如一阶旁轴近似）。这三者共同构成一个适定的初边值问题，也是数值仿真必须完整设定的输入。

有限差分法的基本思想

解析解只在少数简单几何与均匀介质下可得，实际工程中的非均匀介质、复杂边界、变系数问题几乎无法手算。有限差分法的思路是把连续的求解区域离散成网格，用差商代替微商，把偏微分方程转化为代数方程组，再逐步推进求解。

以一维情形为例，将空间区间划分为步长 Δx 的网格点 x_i = iΔx，时间划分为步长 Δt 的时刻 t_n = nΔt，记 u_iⁿ ≈ u(x_i,t_n)。利用中心差分：

∂²u/∂t² ≈ (u_i^{n+1} − 2u_i^n + u_i^{n−1}) / Δt²
∂²u/∂x² ≈ (u_{i+1}^n − 2u_i^n + u_{i−1}^n) / Δx²

代入波动方程并整理，得到显式递推格式：

u_i^{n+1} = 2u_i^n − u_i^{n−1} + r² (u_{i+1}^n − 2u_i^n + u_{i−1}^n)
其中 r = c·Δt/Δx 称为库朗数（Courant number）

该格式的物理含义清晰：下一时刻的位移由当前时刻位移、上一时刻位移以及左右邻点的空间差分共同决定，相当于在离散网格上重现“加速度由曲率驱动”的机制。启动时需要两个时间层，因此第一层通常由初始速度条件通过半步或一阶差分构造。

除稳定性外，还需关注两个精度问题。其一是数值频散：离散网格对不同波长的波传播速度略有差异，短波（接近网格尺度）传播偏慢，导致波形逐渐“拖尾”。减小 Δx（即每波长更多网格点，一般建议不少于 10～20 点）可显著抑制频散。其二是截断误差，中心差分在空间与时间上均为二阶精度，误差随步长平方减小。

数值仿真在波动问题中的作用

数值仿真的价值在于把“不可解”的问题变成“可算”的问题。它至少承担三类作用：第一，预测。在给定介质参数、震源或激励、边界条件后，仿真可以给出任意时刻、任意位置的波场分布，用于评估结构响应、噪声传播或地震动强度。第二，解释。当实测波形与理论预期不符时，通过调整介质速度模型、边界反射系数或衰减参数进行正演对比，可反推实际传播机制。第三，设计。在换能器布设、隔振方案、超声检测参数选择等工程任务中，仿真可在物理实验之前快速筛选方案，降低试错成本。

典型应用场景包括：地震波场正演与偏移成像（石油勘探、地震危险性分析）、结构健康监测中的导波传播（管道、板壳缺陷定位）、超声无损检测（缺陷回波时程预测）、声学与水下声传播建模、以及弦、膜、杆等基础力学振动教学演示。这些场景的共同点是：介质非均匀、边界复杂或激励任意，解析解难以覆盖，而有限差分格式结构简单、易于并行、对复杂速度模型适应性强，因此成为最常用的入门与工程方法之一。

一个可动手的最小示例

为把上述概念落到可执行层面，考虑一根两端固定的弦，长度 L = 1 m，波速 c = 100 m/s，初始位移为半正弦 u(x,0) = sin(πx/L)，初始速度为 0。取 Δx = 0.01 m（100 个网格），由稳定条件选 r = 0.9，则 Δt = rΔx/c = 9×10⁻⁵ s。第一层用 u_i¹ = u_i⁰ + ½r²(u_i+1⁰ − 2u_i⁰ + u_i−1⁰) 构造（对应初始速度为零）。随后按递推式逐层推进，两端点始终置零。运行若干周期后，可观察波形在两端反射、叠加，若 r 取 1.05 则会看到振幅指数增长——这正是库朗条件被违反的直接表现，也是初学者最应亲手验证的一次“失败实验”。