构造解方法（MMS）概述与步骤
理论：先确定解
概述
我想验证自己写的热分析代码，但热导率随温度变化的问题找不到精确解。怎么办？
用构造解方法（MMS）。把思路反过来，先决定“希望得到的解”。把这个解代入方程，左右两边就会有不相等的部分，把它作为源项加到方程中，决定的解就正好是该方程的解。然后在代码中加入该源项求解，确认网格加密时误差是否按格式理论速度减小。关键是用误差的减小方式判断，而不是看解是否看起来对。
例题方程
热导率依赖于解的非线性扩散方程。构造解光滑，在两端分别为0和1。
源项与收敛阶
用符号计算（如sympy）求源项，可以避免手算错误。由网格间距减半时的误差之比求表观收敛阶。
逆向思维：构造解
有解析解的问题数量有限，几乎不适用于实际软件所处理的非线性、多维、复杂边界方程。20世纪80年代Steinberg和Roache提出先确定解、用源项使方程成立的方法，21世纪初Roache、Salari、Knupp等人将其整理为程序验证的标准步骤。如今广泛用于包括湍流模型的流体代码以及固体、热、电磁软件的验证，ASME验证标准也把它作为代码验证的核心方法。
算例：发现小bug
用有限差分比较两种实现：正确实现在面上用两侧平均值计算热导率，错误实现只用左侧值（一个字符的差别）。用Picard迭代求解，求节点最大误差：
| 分段数 | 正确实现误差 | 阶数 | 错误实现误差 | 阶数 |
|---|---|---|---|---|
| 10 | 1.13×10⁻² | — | 4.62×10⁻² | — |
| 20 | 2.85×10⁻³ | 1.98 | 2.13×10⁻² | 1.11 |
| 40 | 7.10×10⁻⁴ | 2.00 | 1.04×10⁻² | 1.03 |
| 80 | 1.78×10⁻⁴ | 2.00 | 5.15×10⁻³ | 1.02 |
| 160 | 4.44×10⁻⁵ | 2.00 | 2.56×10⁻³ | 1.01 |
错误实现误差也只有1%左右，而且随网格加密减小。这样恐怕察觉不到。
这正是MMS的长处。只看误差大小，错误实现也像是“大致对、加密就会变好”，画成图几乎与构造解重合。但误差减小方式明显不同：正确实现为二阶，错误实现为一阶。在面上用单侧值计算系数，使本应二阶的格式混入了一阶误差。实际问题中，这种差别会作为细网格下也不消失的百分之几误差残留。MMS通过阶数这一数字，让看不出的错误浮现出来。
步骤与技巧
- 选择使所有项（对流、扩散、源、非线性系数）都起作用的构造解；过于简单的解会掩盖错误。
- 用符号计算求源项和边界条件，自动输出到代码中。
- 在四级以上网格上求解，求误差范数和阶数。
- 改变边界条件、网格类型和非线性迭代设置，逐一确认。
- 把结果做成自动测试，每次修改代码都确认阶数。
“只有边界是一阶”
用MMS验证某流体代码时，内部网格点误差按二阶减小，整体最大误差却只按一阶减小。看误差分布，只有紧邻壁面边界的点误差大，边界条件实现为一阶精度——壁面梯度用了一阶单侧差分近似。改为二阶单侧差分后整体变为二阶。MMS中不仅看误差大小，还看分布，就能定位错误所在。
常见错误
错误与对策
| 错误 | 后果 | 对策 |
|---|---|---|
| 只凭误差大小判断 | 漏掉阶数下降 | 务必求阶数 |
| 构造解过于简单 | 部分项不起作用，错误被掩盖 | 使所有项起作用 |
| 手算源项 | 源项本身出错 | 用符号计算自动化 |
| 迭代收敛不充分 | 阶数紊乱 | 严格收敛判据 |
| 只用两套网格判断 | 不知道是否处于渐近范围 | 四级以上 |
我想了解相关内容。
可以用Richardson外推工具试算收敛阶。相关页面参考MMS收敛阶、MMS源项、二维导热MMS、Navier-Stokes方程MMS和代码验证。