从解析解到数值解：为什么“算出来了”还不够
发布时间：
数值程序有一个容易被忽视的特点：只要方程组能够被求解，它通常总会返回某个结果。曲线平滑、参数变化规律合理，甚至不同模型之间也呈现出预期的趋势，但这些现象本身并不足以说明数值解已经可靠。
原因在于，数值方法从连续物理问题走到最终结果，中间经历了空间离散、边界截断、数值求导以及线性方程组求解等多个环节。每一步都可能引入误差，而这些误差最终混合在同一组结果中。仅凭最终曲线的“合理程度”，很难判断其中是否存在系统偏差。
如果一个问题恰好同时具有解析解和数值解，情况会简单许多。解析解可以作为参考基准，使误差分析从“结果是否看起来合理”转变为“误差究竟由哪个环节产生”。
一维层状介质中的大地电磁正演就是这样一个例子。本文以这一问题为背景，讨论一个更一般的数值计算问题：
当数值解与解析解存在偏差时，如何定位主要误差来源，并在精度与计算成本之间找到合理的改进方向。
为突出方法本身，文中的模型参数与误差数值只保留必要的量级和趋势，不对应任何特定数据集；实现代码也不在本文讨论范围内。
1. 为什么解析解适合作为数值基准
考虑一维水平层状介质。每一层具有恒定电阻率，最下方为无限延伸半空间。在忽略位移电流的条件下，对于角频率
\[\omega = 2\pi f,\]第 $j$ 层可以定义传播常数
\[k_j=\sqrt{\frac{i\omega\mu_0}{\rho_j}},\]以及本征阻抗
\[Z_{0j}=\sqrt{i\omega\mu_0\rho_j}.\]从最底部半空间出发，可以逐层向地表递推阻抗。一种常用形式为
\[Z_j= Z_{0j} \frac{ Z_{j+1}+Z_{0j}\tanh(k_jh_j) }{ Z_{0j}+Z_{j+1}\tanh(k_jh_j) }.\]得到地表阻抗后，可以计算视电阻率和相位：
\[\rho_a=\frac{|Z|^2}{\omega\mu_0},\] \[\phi=\arg(Z).\]对于理想的一维层状模型，这种递推关系能够直接给出高精度参考结果，并不需要显式离散地下空间。因此，它很适合作为有限差分结果的基准。
这一点很重要。若一个复杂模型本身没有已知标准答案，数值解出现偏差时往往难以判断问题究竟来自物理模型还是数值实现。相反，在具有解析解的简单模型上，两种方法求解的是同一个物理问题，因而可以将差异主要归因于数值近似。
从数值验证的角度看，解析解的价值不只是“提供一个更准确的答案”，更重要的是提供了一把可以逐项检查误差来源的标尺。
2. 从连续方程到最终响应，误差在哪里进入
有限差分法需要首先将连续地下空间离散为有限数量的节点。连续电场
\[E(z)\]在计算中被表示为
\[E_0,E_1,E_2,\ldots,E_n,\]空间导数则由相邻节点值近似。连续微分方程由此转化为代数方程组。
这一过程至少涉及四类常见误差：
| 环节 | 主要问题 | 常用诊断方式 |
|---|---|---|
| 空间离散 | 有限网格产生截断误差 | 系统加密网格并检查收敛 |
| 地层界面 | 界面与网格位置不匹配 | 对齐界面并进行控制变量比较 |
| 计算边界 | 有限区域替代无限介质 | 逐步扩大计算域 |
| 场量提取 | 用离散节点近似地表导数 | 比较不同阶数的导数公式 |
对于大地电磁问题，最后一项尤其值得注意。
阻抗需要地表电场与磁场，而磁场可由电场的空间导数得到。因此，即使有限差分已经较好地求出了地下电场，最终仍需要在地表进行一次数值求导。
整个计算链条可以写成
\[\text{连续模型} \rightarrow \text{空间离散} \rightarrow E(z)\text{ 的数值解} \rightarrow E'(0)\text{ 的近似} \rightarrow Z \rightarrow \rho_a,\phi.\]最终误差是这些步骤共同作用的结果。只看视电阻率和相位曲线，很难直接判断主要误差究竟来自哪里。
图 1 从连续模型到最终响应的计算链。解析解的意义之一，是让空间离散、边界与后处理误差可以被分别诊断。
3. 网格加密是必要的，但不应成为唯一手段
面对数值误差，最直接的改进通常是加密网格。
设特征网格尺度为 $h$。一阶有限差分的截断误差通常具有
\[O(h)\]的量级；二阶方法在进入渐近收敛区间以后，则可能表现为
\[O(h^2).\]因此，当网格逐渐加密时，数值结果应逐步向解析解靠近。这是验证离散方法的基本手段。
然而，“加密以后误差变小”只回答了一个问题：空间分辨率确实影响结果。它并不能说明原始误差主要来自哪个环节。
如果某个局部步骤本身只采用了一阶近似，那么通过整体加密网格来补偿这一误差，往往意味着用更多计算资源去弥补一个可以直接修改的低阶环节。此时更有效的做法，是在保持其他条件不变的情况下逐项检查误差来源。
这也是解析解真正开始发挥作用的地方。
图 2 数值解随网格加密向参考解收敛的示意。图中曲线仅表达数值关系，不使用原始实验数据。
4. 地表导数可能限制整体精度
最简单的地表导数可以由前两个节点近似：
\[E'(0)\approx\frac{E_1-E_0}{h}.\]这是标准的一阶前向差分，其截断误差为
\[O(h).\]如果地下内部的离散方案已经具有较高精度，而最终阻抗仍依赖这样一个一阶导数近似，那么整体计算精度就可能受到地表提取步骤的限制。
对于非均匀网格，若地表以下前两个单元宽度分别为 $a$ 和 $b$，可以利用前三个节点构造二阶近似：
\[E'(0)\approx -\frac{2a+b}{a(a+b)}E_0 +\frac{a+b}{ab}E_1 -\frac{a}{b(a+b)}E_2.\]这一修改的特点在于，它不改变地下有限差分方程，也不需要改变已经得到的电场数值解，只改变最后的地表导数提取方式。
在代表性测试中，仅提高这一局部步骤的近似阶数，就足以使视电阻率相关的整体误差下降接近两个数量级，相位误差也获得约一个数量级的改善。具体数值会随模型与网格变化，但趋势说明了一个更普遍的问题：
数值计算中，最终精度可能由某个局部的低阶步骤控制，而并非由整个求解器平均决定。
因此，在投入更高计算成本之前，先识别误差瓶颈通常更有效。
5. 如何把不同误差源分离出来
误差分析最有用的方法之一，是控制变量。
例如，在完全相同的地下电场数值解上，只改变地表导数公式。如果最终误差明显下降，那么可以较有把握地将改进归因于地表求导，而不必同时猜测网格、边界条件或其他参数的作用。
同样地，可以固定网格和求导方法，只改变底部计算边界的位置。若边界逐渐加深后误差先下降、随后趋于稳定，则表明边界误差已经不再是主要限制因素。
从概念上，可以将总误差写成
\[\varepsilon_{\mathrm{total}} \sim \varepsilon_{\mathrm{grid}} + \varepsilon_{\mathrm{surface}} + \varepsilon_{\mathrm{boundary}} + \varepsilon_{\mathrm{interface}} +\cdots\]这并不是严格的线性误差分解。不同误差项之间可能发生耦合，也可能相互抵消。这个表达式更适合作为一种分析框架：将“数值结果有误差”进一步拆成若干可以单独实验的问题。
对于视电阻率和相位，还可以分别定义带符号误差，例如
\[e_\rho(f)=100 \left( \frac{\rho_a^{\mathrm{num}}(f)} {\rho_a^{\mathrm{ref}}(f)}-1 \right),\]以及
\[e_\phi(f)= \phi^{\mathrm{num}}(f) - \phi^{\mathrm{ref}}(f).\]逐频率观察这两个量，比单独给出一个综合误差指标更容易识别某些局部频段的异常。
6. 为什么局部加密有时不会让总误差单调下降
局部网格加密的结果有时会出现一个看似反常的现象：某个区域被加密以后，总误差不但没有明显下降，甚至可能略有增加。
这并不意味着更细的网格降低了离散精度。
假设两个误差项分别为
\[\varepsilon_1>0, \qquad \varepsilon_2<0,\]那么总误差为
\[\varepsilon= \varepsilon_1+\varepsilon_2.\]| 如果两项原本恰好部分抵消，总误差可能表现得很小。局部加密只减小了其中一个误差项以后，原有的抵消关系会被改变，最终观察到的 ( | \varepsilon | ) 反而可能暂时增加。 |
因此，一次修改前后的误差大小并不足以证明方法是否真正改善。
更可靠的办法是采用一系列规律变化的网格尺度，例如
\[h,\quad \frac{h}{2},\quad \frac{h}{4},\quad \frac{h}{8},\]并观察误差是否呈现稳定的收敛规律。
7. 收敛阶比单独一个“小误差”更有意义
如果某一组计算的误差达到 $10^{-4}$，看起来已经相当小，但这个数字本身并不能充分证明数值方法可靠。
它可能来自真正高精度的离散，也可能只是若干误差项在这一组参数下发生了偶然抵消。
相反，如果随着网格尺度减半，误差近似满足
\[\varepsilon(h) \rightarrow \frac{\varepsilon(h)}{4} \rightarrow \frac{\varepsilon(h)}{16},\]那么方法正在表现出典型的二阶收敛特征。
这种收敛关系比单独报告“数值解和解析解很接近”更有说服力，因为它给出了误差随离散尺度变化的规律。换句话说，我们不仅知道当前结果比较准确，还知道继续加密时它会以怎样的方式趋近参考解。
在实际测试中，当地表求导与空间离散都采用较高阶处理以后，随着网格分辨率提高，误差逐渐呈现出接近二阶的收敛趋势。这个现象比某一个具体误差值更能说明数值方案进入了稳定的渐近区间。
图 3 一阶与二阶方法的典型收敛趋势。相比单独报告某一组误差，收敛阶能够更直接地检验离散方法是否按理论预期工作。
8. 计算边界需要“足够深”，而不是“尽可能深”
有限差分模型必须截断计算区域，而真实地下介质并不存在这样的人工边界。
常见的处理方法是将底边界设置到足够深的位置，使该处的场已经显著衰减，从而减弱边界条件对地表响应的影响。
但底边界并非越深越好。
逐步增加计算域深度时，误差通常会在开始阶段有所变化；达到一定深度以后，进一步扩大计算域对最终误差几乎没有影响。这说明边界误差已经被压低到次要水平，继续增加深度只会增加节点数量和求解成本。
因此，合理的数值设置并不是把所有参数取到尽可能大的值，而是找到误差已经不再对结果产生实质影响的尺度。
9. 趋肤深度提供了一个自然的网格尺度
大地电磁问题中，趋肤深度为网格设计提供了一个重要的物理尺度：
\[\delta= \sqrt{\frac{2\rho}{\omega\mu_0}} \approx 503\sqrt{\frac{\rho}{f}} \ \mathrm{m}.\]频率越低、电阻率越高，电磁场能够有效影响的深度尺度越大。
这意味着不同频率对应的空间尺度并不相同。若所有频率统一使用固定网格，高频部分可能因为网格过粗而出现较大离散误差，低频部分则可能使用远超需要的空间分辨率。
一种更自然的处理方式，是让网格尺度与趋肤深度相关：
\[\Delta z\sim\frac{\delta}{P},\]其中 $P$ 可以理解为每个趋肤深度内希望保留的离散单元数。
这样，网格分辨率就与电磁场自身的变化尺度建立了联系。相比单纯规定“每一层划分多少单元”，这种设计更容易兼顾精度与计算规模。
10. 精度优化和效率优化往往是同一个问题
误差来源被分离以后，优化方向会更加清楚。
整体网格加密能够降低离散误差，但计算量也随之增加；提高地表导数阶数则可能在几乎不增加节点数的情况下显著改善精度。不同方法的成本并不相同，因此不能只比较“谁的误差更小”，还应比较达到同一目标精度需要付出多少计算资源。
还有一些优化来自数值实现本身。
一维有限差分形成的系数矩阵具有明显的带状稀疏结构。若直接构造完整的 (N\times N) 矩阵，其存储量大致按
\[O(N^2)\]增长；若只存储有限数量的非零对角线，存储规模则可以接近
\[O(N).\]这种改变并不影响物理模型，也没有增加新的近似，只是避免存储和处理大量本来为零的矩阵元素。
从这个角度看，数值优化的核心并不是单纯让计算机进行更多运算，而是识别哪些计算真正决定精度，哪些开销可以在不损失结果质量的情况下被减少。
11. 一个更可靠的数值验证流程
将上述分析整理以后，可以得到一个相对通用的数值验证流程。
首先，在可能的情况下选择一个具有解析解、半解析解或高精度参考解的简化模型。它不需要和最终研究对象完全一致，只需要包含核心物理过程。
随后，为最终观测量定义明确的误差指标，并检查误差随频率、空间位置或其他自变量的分布，而不是只报告一个综合数字。
接着，通过控制变量分别检查网格、边界、界面位置和后处理步骤。若某个局部修改能够在较低成本下显著减小误差，应优先解决这个瓶颈。
在此基础上，再进行系统的网格收敛测试，并判断结果是否呈现理论预期的收敛阶。
最后，才是精度与效率之间的权衡：确定满足目标误差所需的最小计算规模，而不是无限制地提高分辨率。
这个流程可以概括为
\[\text{参考解} \rightarrow \text{误差度量} \rightarrow \text{误差定位} \rightarrow \text{收敛检验} \rightarrow \text{成本—精度权衡}.\]它并不限于大地电磁问题。有限元、有限体积、谱方法以及其他偏微分方程数值求解，都面对类似的问题。
12. 结语
解析解与数值解之间的比较，表面上只是一次精度验证，实际上能够提供关于数值方法本身的更多信息。
解析解给出理想模型的参考答案，数值方法则允许我们进一步处理复杂几何、非均匀介质以及难以获得解析表达式的问题。在一个两者都能求解的简化模型上进行比较，可以提前暴露离散、边界和后处理中的误差来源，并验证算法是否具有预期的收敛性质。
因此，在存在解析解的情况下额外实现数值解，并不是重复计算。更准确地说，这是在进入复杂模型之前，对数值方法进行一次可解释、可量化的校准。
数值程序最终总会给出某个结果。
真正需要建立的，是我们为什么可以相信这个结果。