本研究由E. Haber和U. M. Ascher完成，作者分别来自加拿大不列颠哥伦比亚大学（University of British Columbia）的计算机科学系与地球和海洋科学系。该论文发表于《SIAM Journal on Scientific Computing》第22卷第6期，页码为1943至1961页，于2001年2月21日在线发表。
本研究的学术背景属于计算电磁学与地球物理数值模拟交叉领域。研究的核心问题是求解三维电磁问题中具有高度不连续系数（如电导率、磁导率和介电常数）的Maxwell方程组。在实际地球物理勘探和医学成像等应用中，地下介质的电导率和磁导率往往在空间上剧烈变化，跨越多个数量级，并且在介质界面处产生解的跳变不连续。传统的Yee氏方法虽然在电磁场时域模拟中广泛应用，但当频率较低、系数跳变剧烈时，由于旋度算子存在非平凡零空间，迭代求解离散代数方程组的收敛速度极慢，尤其是对于电场源问题。作者在先前工作中已经针对恒定磁导率情形提出过基于位函数和Helmholtz分解的方法，本研究将该方法推广至磁导率也可能剧烈变化的情形。研究目标是构造一种能够快速、准确求解该类三维电磁正问题的数值算法，以适应地球物理反演中对正演求解器需反复调用的需求。
研究的方法流程可以概括为以下若干步骤。第一步是对Maxwell方程组进行解析重构。作者首先将电场(e)进行Helmholtz分解，分解为无散部分(a)和梯度部分(\mathrm{grad}\,\varphi)，其中(\nabla \cdot a = 0)。将Helmholtz分解（Helmholtz decomposition）代入频域Maxwell方程，得到关于(a)与(\varphi)的耦合方程组。为了消除旋度算子的零空间对迭代求解的不利影响，作者在方程中减去一个梯度项(\mathrm{grad}(\mu^{-1}\nabla \cdot a))，该项在原方程中由于(\nabla \cdot a = 0)而为零，但在离散格式中起到稳定化作用，使所得微分算子成为强椭圆型算子。进一步地，作者引入广义电流密度(\hat{j} = \hat{\sigma}e = (\sigma - i\omega\varepsilon)e)和磁通量(b = \mu h)，并将系统写成一阶混合形式，包含五个方程：(\nabla \times h - \mathrm{grad}\,\psi - \hat{j} = j_s)，(i\omega\mu h - \nabla \times a = 0)，(\mu\psi - \nabla \cdot a = 0)，(-\nabla \cdot \hat{j} = \nabla \cdot j_s)，以及(\hat{j} - \hat{\sigma}(a + \mathrm{grad}\,\varphi) = 0)。其中(\psi)为辅助标量场。该系统中(\psi)通过边界条件和散度关系被证明恒为零，从而在离散层面保证(a)的无散性。
第二步是离散化方案的设计。作者采用交错网格有限体积法（staggered grid finite volume discretization），将(\hat{j})和(a)定义在单元面中心，(h)定义在单元棱边中点，(\varphi)和(\psi)定义在单元体心。这种变量排布方式与传统的Yee氏网格中(e)在棱边、(h)在面心的布置形成对偶。对于电导率(\hat{\sigma})在界面处的平均，作者使用调和平均（harmonic average），因为电流在界面法向呈串联关系；对于磁导率(\mu)在棱边处的平均，作者使用算术平均（arithmetic average），因为磁通在共享棱边的四个单元间呈并联关系。这种平均方式的选取对于保证界面物理通量的近似精度至关重要。离散算子满足离散的矢量恒等式，如(\nabla_h \times \mathrm{grad}_h = 0)和(\nabla_h \cdot \nabla_h \times = 0)，从而保持了与连续情形类似的微分算子结构。
第三步是线性代数系统的求解。通过消去(h)、(\psi)和(\hat{j})，最终得到关于(a)和(\varphi)的耦合大型稀疏复线性系统。该系统具有分块结构，主对角块分别对应于离散的矢量拉普拉斯型算子和标量拉普拉斯型算子，二者之间的耦合较弱。作者利用这一结构特点构造了基于块不完全LU分解（block ILU）的预处理子，同时将SSOR预处理子作为对比。由于在低频条件下矩阵的对角块占优，块ILU预处理子能够显著加速Krylov子空间方法（具体为BiCGSTAB）的收敛。
数值实验部分包含两个主要算例集。第一个算例集考察一个电导率和磁导率均匀的立方体嵌入均匀半空间中的模型，电导率对比度(\tilde{\sigma})从10¹到10⁶，磁导率对比度(\tilde{\mu})从10¹到10³，频率从0到10⁶ Hz。实验对比了本文提出的((a,\varphi))方法与Yee氏方法在电场源和磁场源下的迭代次数。结果以表格形式呈现：对于电场源，在频率为10⁰ Hz时，本文方法需要46次迭代而Yee方法需要9856次迭代；当频率升高至10⁶ Hz时，本文方法为97次而Yee方法为11212次；对于磁场源，两种方法的迭代次数差异很小。该结果表明，对于非无散源，传统Yee格式因旋度算子零空间中残差分量难以消除而收敛极慢，本文方法则显著更快，尤其在电场源和低频情形下优势超过两个数量级。进一步的实验考察了系数跳变大小对本文方法迭代次数的影响，发现(\tilde{\sigma})增大对迭代次数影响很小，而(\tilde{\mu})增大则显著增加迭代次数，从(\tilde{\mu}=10^1)时约30–60次迭代增加到(\tilde{\mu}=10^3)时约150–365次迭代。网格加密实验显示迭代次数随未知数个数以约三分之一次幂的速率增长，且块ILU预处理比SSOR所需迭代次数更少，但内存开销更大。
第二个算例集采用随机介质模型，其中每个网格单元以概率(p=0.5)取两组电导率和磁导率之一，用以模拟随机地球模型。在均匀44³网格上，对于电导率对比度从10²到10⁶、磁导率对比度从10¹到10³、频率从0到10⁶ Hz的多种组合进行了测试。结果表明本文方法对高度变化的电导率仍具有良好的求解效率。然而当频率较高且电导率较大，使得网格间距与趋肤深度不匹配时（(\omega\mu\sigma h)超过1），块ILU预处理子可能失效，表现为在800次迭代内无法达到10⁻⁷的残差。此时需要更细的网格以保持预处理子的有效性。
本研究的结论是：通过将Helmholtz分解、库仑规范稳定化处理、一阶混合形式转换以及交错网格有限体积离散相结合，能够有效求解具有显著电导率和磁导率不连续性的三维Maxwell方程。该方法对电场源问题相比传统Yee方法具有显著的速度优势，并且在随机介质等复杂模型中保持稳健效率。其科学价值在于为强不连续系数下的电磁正演问题提供了收敛性良好的离散化与预处理方案；其应用价值在于该算法已实际用于地球物理电磁反演中，在单处理器个人计算机上可在不到两分钟的时间内求解实际规模的正演问题，从而为反演算法的高效实现提供了可能。
本研究的亮点主要体现在以下几个方面：第一，将磁导率不连续性纳入位函数方法框架，扩展了先前仅适用于恒定磁导率的工作；第二，利用离散无散性质保证稳定化项在离散精确解处为零，从而避免了大对比度磁导率情形下罚函数方法可能引起的数值问题；第三，根据物理通量串联与并联关系分别对电导率和磁导率采用调和平均与算术平均，确保了界面通量的正确离散；第四，所构造的块ILU预处理子利用了解析重构后系统的弱耦合块结构，显著加速了Krylov子空间方法的收敛。此外，文中关于高频高电导率情形下预处理子失效的讨论也为实际应用中网格设计提供了重要参考。