企业官网建设流程全解析

日拱一卒这个系列，我给自己定的规矩很简单：每天解决一个具体问题，把这个问题的推导、验证过程和踩过的坑写下来。今天这卒，拱在亥姆霍兹方程验证上。

亥姆霍兹方程，只要是沾过电磁、声学、量子力学或者结构振动的人，肯定都见过这个形式：∇²u + k²u = 0。它本质上是时域波动方程在做时间简谐假设之后得到的频域方程，用处极广。但公式会写，和真正动手算出解、再把数值结果和解析结果对上，完全是两码事。这篇文章记录的就是这么一次完整的验证——从波动方程推导出亥姆霍兹方程；在一维区间和二维矩形域中给出解析解；用有限差分法从零实现数值求解；最后用多个误差指标做对比。

整篇适合三类人：一类是刚学完数理方程、想动手验证理论的学生，一类是做声学或者微波仿真、需要确认求解器正确性的工程师，还有一类是想学有限差分法但不想一上来就碰商业软件的朋友。

1. 从波动方程到亥姆霍兹方程：验证思路的拆解

1.1 时间简谐假设：为什么一个偏微分方程能变成椭圆型方程

亥姆霍兹方程最标准的出身是标量波动方程：

∂²u/∂t² = c²∇²u

这里的 u 可以是声压、电磁场的某个分量，也可以是弹性介质里的位移分量；c 是介质中的波速。实际工程问题里，绝大多数稳态情况都是单频激励。比如扬声器发出的正弦声波、通信天线的连续波，都是单一频率。这种条件下，场的时变规律可以用一个复指数 e^(-iωt) 统一描写，也就是做时间简谐假设：

u(x, y, z, t) = φ(x, y, z) · e^(-iωt)

代入原波动方程后，时间二阶导数直接变成代数因子 -ω²，于是：

∇²φ + (ω²/c²)φ = 0

把 ω/c 记作 k，就是标准亥姆霍兹方程：

∇²φ + k²φ = 0

这个推导本身不难。但为什么要强调这一层？因为很多人面对亥姆霍兹方程时会忽略它的物理语境，直接把它当成一个纯偏微分方程来解。实际上，k 的意义就是波数，它和频率、波速、波长直接相关。在验证过程中，k 的取值决定了内部场的振荡次数，也就决定了数值网格到底需要多密。如果把 k 取错了，后面误差分析的结果会非常难看，而且很难排查。

1.2 验证方案的整体设计：解析路线、数值路线和误差校核

我这次验证没有直接拿现成的商业软件，而是把路线拆成三条：

• 解析线：选择区域形状规则、边界条件简单的场景，把解用显式公式写出来。一维用三角函数解，二维矩形域用正弦级数解。
• 数值线：用有限差分法把偏微分方程离散成代数方程组，编写代码求解。视图是先用最简单的格式把流程跑通，不追求高阶精度。
• 校核线：在相同的采样点上对比解析解和数值解，分别计算 L2 误差和最大绝对误差，并观察网格加密时误差的下降速度。

之所以不直接上有限元软件，是因为这个问题的核心是“验证”，不是“解决复杂工程问题”。有限差分法的离散过程完全透明，误差来源清晰，适合把整个知识链条打通。用有限元软件当然能算，但那是黑盒，得到的是一堆云图和数值，中间的数学结构反而被掩盖了。对于日拱一卒型的学习任务，我宁可每一步都见底。

2. 解析解推导：一维本征解和二维矩形域模态

2.1 一维亥姆霍兹方程的边值问题与本征解

一维亥姆霍兹方程非常直接：

d²u/dx² + k²u = 0

这是一个线性常系数二阶微分方程，通解是：

u(x) = A·sin(kx) + B·cos(kx)

如果考虑一个长度为 L 的区间 [0, L]，两端都是固支边界，即 u(0) = 0，u(L) = 0，那么代入边界条件后，B 必须等于零，同时要求 sin(kL) = 0。于是 k 不能任意取，它只能取离散的一系列值：

kₙ = nπ/L，n = 1, 2, 3, ...

对应每个 kₙ，解是：

uₙ(x) = sin(nπx/L)

这就是这个边值问题的本征解。换个角度说，k 的取值不是指定的，而是由边界条件和区域形状决定的。数值验证里如果想要非平凡解，就必须在这个本征值框架下理解问题，而不是随便给一个 k 就去线性方程组里找非零解——那样只会得到零解，这一点后面会专门展开。

我在验证时就取 L = 1，n = 1，这样解析解就是 u(x) = sin(πx)，对应的 k = π。这个场景干净到极致，理论上最大误差应该只来自离散近似，很容易做定量评估。

2.2 二维矩形域分离变量法与模态组合

二维的路线和一维完全一致。考虑矩形区域 [0, a] × [0, b]，四条边都是零边界，即 u = 0。方程为：

∂²u/∂x² + ∂²u/∂y² + k²u = 0

用分离变量法，设 u(x, y) = X(x)Y(y)。代入后可以得到两个独立的一维问题，最终一组标准的本征解是：

uₘₙ(x, y) = sin(mπx/a)·sin(nπy/b)

对应的波数：

kₘₙ = π√(m²/a² + n²/b²)

这里 m、n 都是正整数。验证时取 a = b = 1，m = n = 1，于是：

k₁₁ = π√2

解析解就是 u(x, y) = sin(πx)·sin(πy)。这个函数在正方形区域中心取最大值1，四条边是零。无论是形状还是边界条件，都很适合作为数值验证的基准。

2.3 把解析解和物理意义挂上钩：声腔模态

如果觉得上面的推导只是纸面数学，我可以举一个常用的物理对应：空气填充的矩形硬壁腔体。四壁刚性边界在声学上近似为声压的法向导数为零，但这个例子里如果用声压为零的软边界，解的形式完全是同样套路。更常见的情况是硬壁腔体的速度势模态，本征函数同样是三角函数乘积，对应的本征频率由 k 换算：

fₘₙ = c·kₘₙ/(2π)

比如波速 c = 343 m/s，a = b = 1 m 时，第一阶模态频率大约是 77.2 Hz。数值求解亥姆霍兹方程的本质，就是在求这个腔体有哪些共振频率、每个频率对应的模态形状长什么样。这也解释了为什么这类方程在室内声学、扬声器箱体设计、微波谐振腔分析里永远绕不开。验证解析解时如果加入了物理视角，你会更容易判断算出来的模态到底是合理还是异常。

3. 数值验证实操：用有限差分法求解亥姆霍兹方程

3.1 为什么选有限差分法而不是有限元

有限差分法的思路一句话就能说清楚：把微分用差商代替。二阶导数在均匀网格上最常用的中心差分格式是：

u''(xᵢ) ≈ (uᵢ₊₁ - 2uᵢ + uᵢ₋₁) / h²

这个格式的截断误差是 O(h²)，也就是网格间距 h 减半时，误差大约变成原来的四分之一。对验证型工作来说，这个精度完全够用，而且处理流程非常透明。

相比之下，有限元法边界条件处理灵活、能适配复杂几何，但需要引入形函数、单元组装、数值积分等一大堆概念。对这个验证课题来说，属于杀鸡用牛刀。更重要的是，有限差分法用 Python 加 NumPy 就能实现，不需要接触大型依赖库，这在学习阶段能节省大量时间。

3.2 一维本征问题：矩阵组装与 Python 实现

一维离散后的方程可以写成矩阵形式。内部节点编号从 1 到 M，边界节点是 0 和 N。离散方程是：

(uᵢ₊₁ - 2uᵢ + uᵢ₋₁) / h² + k²uᵢ = 0

移项以后，它等价于一个矩阵特征值问题：

A·v = -k²v

矩阵 A 是二阶导数算子的离散矩阵，是三对角结构：

A[i, i] = -2/h² A[i, i-1] = 1/h² A[i, i+1] = 1/h²

求解这个特征值问题，得到的特征值就是 -k²，特征向量就是本征函数。我实际使用的代码如下：

import numpy as np L = 1.0 N = 200 h = L / N M = N - 1 # 构造二阶导数中心差分矩阵 A = np.zeros((M, M)) for i in range(M): A[i, i] = -2.0 / h**2 if i > 0: A[i, i - 1] = 1.0 / h**2 if i < M - 1: A[i, i + 1] = 1.0 / h**2 # 求解特征值问题 A v = -k^2 v eigvals, eigvecs = np.linalg.eigh(A) # 最小正波数对应最接近零的负特征值，也就是最大特征值 k2_num = -eigvals[-1] k_num = np.sqrt(k2_num) # 取出对应特征向量，补上两端边界零值 v = eigvecs[:, -1] u_num = np.concatenate((np.zeros(1), v, np.zeros(1))) # 归一化到解析解的量纲（解析解最大值为 1） u_num = u_num / np.max(np.abs(u_num)) # 解析解与误差 x = np.linspace(0, L, N + 1) u_exact = np.sin(np.pi * x) err = u_num - u_exact print("解析 k =", np.pi) print("数值 k =", k_num) print("L2 误差 =", np.linalg.norm(err) * np.sqrt(h)) print("最大绝对误差 =", np.max(np.abs(err)))

这里有一个细节值得停留：为什么取特征值数组的最后一个，而不是第一个？因为np.linalg.eigh返回的特征值是升序排列。一维 Dirichlet 问题里所有特征值都是负的，数值最大的那个负特征值最接近零，对应最小的 k²，也就是基模 n = 1。如果取第一个特征值，得到的是最高阶模态，那就完全对不上了。

3.3 二维五点差分：从一维推广到稀疏矩阵

二维的离散格式是一维的自然扩展。在均匀网格上用五点差分格式：

(uᵢ₊₁,ⱼ - 2uᵢ,ⱼ + uᵢ₋₁,ⱼ) / hₓ² + (uᵢ,ⱼ₊₁ - 2uᵢ,ⱼ + uᵢ,ⱼ₋₁) / hᵧ² + k²uᵢ,ⱼ = 0

如果直接用二维数组下标硬解，代码会很啰嗦。更干净的做法是构造二维拉普拉斯矩阵，它等于两个一维矩阵的 Kronecker 和：

L₂ = Lₓ ⊗ I + I ⊗ Lᵧ

然后用稀疏矩阵求解特征值问题。这里我给出核心片段：

import numpy as np a = b = 1.0 Nx = Ny = 40 hx = a / Nx hy = b / Ny Mx = Nx - 1 My = Ny - 1 def lap_1d(M, h): A = np.zeros((M, M)) for i in range(M): A[i, i] = -2.0 / h**2 if i > 0: A[i, i - 1] = 1.0 / h**2 if i < M - 1: A[i, i + 1] = 1.0 / h**2 return A A1x = lap_1d(Mx, hx) A1y = lap_1d(My, hy) Ix = np.eye(Mx) Iy = np.eye(My) A2d = np.kron(Iy, A1x) + np.kron(A1y, Ix) eigvals, eigvecs = np.linalg.eigh(A2d) k2_num = -eigvals[-1] k_num = np.sqrt(k2_num) print("解析 k11 =", np.pi * np.sqrt(2)) print("数值 k11 =", k_num)

注意这里 grid 是 Nx × Ny 的矩形网格，内部自由度是 (Nx-1) × (Ny-1)。当 Nx=Ny=40 时，自由度是 1521，用稠密矩阵np.linalg.eigh仍然很快。如果想要更高分辨率，比如 100×100 网格，自由度接近一万，稠密矩阵还能勉强，再大就应该切到scipy.sparse.linalg.eigsh了。这也是为什么工程上大网格都倾向用稀疏迭代求解器的原因。

3.4 收敛性验算：二阶精度是怎么来的

误差分析不能只看一组数据，至少要跑三组网格，看误差是不是按预期下降。下面是我在一维问题上实际跑出来的数据，N 分别取 40、80、160：

k 的解析值是 π ≈ 3.1416。可以看到，网格从 40 翻到 80 时，误差大约是原来的四分之一，从 80 翻到 160 时又是四分之一。这就是中心差分二阶精度的典型特征。特征函数的最大绝对误差也按同样规律下降，说明不仅特征值收敛，形状也收敛。

二维问题收敛规律类似，但因为同时存在 x 和 y 两个方向，总自由度增加很快，我的实测经验是 Nx=Ny=20、40、80 三组数据足够看到二阶趋势。跑更高分辨率时内存和时间成本上升明显，需要权衡。

4. 常见问题与排查技巧实录

4.1 零解陷阱：为什么线性方程组解出来全为零

这是新手做亥姆霍兹方程最容易踩的坑。如果直接给定一个 k，再配合零边界条件，去解离散方程 A·u + k²u = 0，得到的解必然是零。这不是代码写错了，而是数学上这个边值问题在非本征频率下只有零解。只有 k 恰好等于本征值时，零边界条件下才存在非平凡解，但那时矩阵又变成奇异的，直接求逆也得不到稳定结果。

正确的处理方式有两种：要么把问题定义成本征值问题，像上面代码里那样求解矩阵的特征值和特征向量；要么给系统加非零激励或者非零边界条件，比如某一边设 u = 1，此时 k 取任意非本征值都能得到唯一非零解，再和解析解对比。我之前一次就是卡在这，调了一下午矩阵，最后意识到是问题本身提法不对，而不是求解器问题。

4.2 边界条件处理错误导致模态形状变形

Dirichlet 边界，也就是 u = 0，处理起来最简单，直接把对应节点从未知量里删掉即可。真正容易出错的是 Neumann 边界，也就是法向导数为零的边界。这类边界在硬壁声场问题里非常常见。差分格式里需要边界外的虚拟节点，做法是：

∂u/∂n = 0 近似为 (u₊₁ - u₋₁)/(2h) = 0

也就是说边界外虚拟节点的值和边界内侧节点的值相同，然后把这个关系代回到离散方程里消掉虚拟节点。很多人在这一步会遗漏，导致边界附近的解明显偏离解析解，整体模态形状被拉歪。验证时我总是先画一维解和解析解叠在一起的图，边界附近误差突变就是边界条件写错了的最直观信号。

4.3 数值色散：网格不够密时波数会悄悄偏移

离散后的网格和连续介质之间其实有不可忽视的差异。中心差分格式的离散色散关系是：

sin²(k_h·h/2) / (h/2)² ≈ k²

也就是说，网格算出来的有效波数 k_h 和真实波数 k 并不完全相等，二者的偏差随 h 增大而增大。当网格间距接近波长的一半时，偏差会大到离谱，解完全失真。工程经验是每个波长最少要有 10 个网格点，也就是 h ≤ λ/10。如果波数 k = π，那 h ≤ 0.2；如果 k = 10π，那么 h ≤ 0.02。网格数太少时，即使本征值收敛趋势正确，绝对误差也会卡在一个偏高的平台上不去，这时候不是代码问题，是分辨率跟不上。

4.4 特征向量振幅不定：归一化是必做的一步

本征值问题的特征向量乘以任意非零常数仍然是特征向量。数值求解器给出的向量，振幅没有任何物理意义。因此拿数值特征函数和解析特征函数对比之前，必须先归一化。

我的做法是先找到数值解的最大绝对值，然后把整个向量除以这个最大值。解析解构造时也让最大值为 1，比如 sin(πx) 在 [0,1] 上的最大值恰好是 1。这样两者才有可比性。如果不做这一步，哪怕求解结果完全正确，误差指标也会因为缩放因子不同而变得很大，容易被误判为代码错误。

4.5 符号约定不统一和相位问题

亥姆霍兹方程推导里最容易让两个不同教材对不上号的，是时间因子用 e^(-iωt) 还是 e^(+iωt)。两个约定下，波数 k 的定义相差一个符号，本征函数里相位也会有差异。做理论推导时可能无所谓，但一旦要和实验数据的相位对齐，就必须统一约定。我习惯用 e^(-iωt)，因为电磁场领域大多数参考书的复向量记法都采用这个约定。建议你在自己的项目里固定一种约定，并在代码注释里写清楚，省得几个月后回看代码时自己都迷惑。

5. 再分享一点扩展思路

这次验证其实还留了不少可以继续拱的空间。比如把一维标准差分换成更高精度的四点差分格式，收敛阶会从二阶变成四阶，误差下降会更快，但矩阵带宽也会从三对角变成五对角，代码复杂度会上去。又比如把矩形区域换成圆形区域，底柱坐标下会出现贝塞尔函数，解析解变成 Jₙ(kr)cos(nθ) 的形式，数值验证的难度会上一个台阶，但和声波在圆形膜上的振动、光波在环形波导里的传播都能直接对应。这些方向我已经记录在下一批日拱一卒的选题里。

最后再分享一个实际工作中的小技巧：验证偏微分方程求解器时，不要只看一组网格的误差，一定要跑网格收敛性测试，并且把误差随网格数变化的趋势记录下来。如果误差斜率不符合预期，说明离散格式、边界条件或后处理里至少有一个环节存在问题。这个习惯能在后续处理更复杂的声学或电磁仿真问题时帮你省下大量排查时间。