MMS：二维稳态导热
思路：先决定答案
概述
我想确认自己写的导热求解器是否正确，但复杂形状和边界条件没有解析解。怎么办？
自己构造答案就行了，这就是构造解方法（MMS）。比如选定 T = cos(πx/3)·sin(πy) 这样光滑的任意函数作为“答案”，代入导热方程，就能反推出得到这个答案所需的热源 q 和边界条件。把它们交给求解器求解，再把得到的温度与预先定好的答案比较。网格加密时误差按理论速度减小，就可以说代码实现正确。
方程与构造解
导热系数 $k = 1$ 的单位正方形。$x = 0$、$y = 0$、$y = 1$ 处以构造解的值作为温度给定，$x = 1$ 处给定由构造解求出的热流（诺伊曼条件）。$a$ 为改变构造解形状的参数。
收敛阶
$e_h$ 为网格间距 $h$ 时的误差。中心差分五点格式为二阶精度，$p$ 应趋近2。边界条件或源项实现有误时，$p$ 会低于理论值。
边界条件的两种写法
左边是在边界外设虚拟点的二阶精度写法，右边是单侧差分的一阶精度写法。把右边作为“错误”植入，测试MMS能否发现。
“由答案造问题”的思路
构造解的思想从20世纪80年代起就在数值计算研究者中使用，但把它系统整理为代码验证方法的是Roache等人，通过1998年的著作和21世纪初的论文推广开来。其划时代之处在于用逆向思维打破了“没有解析解就无法验证”的壁垒。如今许多商业和研究代码都公开MMS验证结果，ASME V&V 10和20也把它作为代码验证的核心方法。
算例
例1：正确实现的收敛
| 分割数 | 最大误差（a = 2） | 阶数 | 最大误差（a = 3） | 阶数 |
|---|---|---|---|---|
| 8 | 5.93×10⁻³ | — | 8.27×10⁻³ | — |
| 16 | 1.49×10⁻³ | 1.99 | 2.07×10⁻³ | 2.00 |
| 32 | 3.74×10⁻⁴ | 2.00 | 5.16×10⁻⁴ | 2.00 |
| 64 | 9.34×10⁻⁵ | 2.00 | 1.29×10⁻⁴ | 2.00 |
两种构造解下，网格每减半误差都降到1/4，确认了理论上的二阶收敛。
例2：边界条件错误
| 分割数 | a = 2（边界曲率为0） | 阶数 | a = 3（边界有曲率） | 阶数 |
|---|---|---|---|---|
| 8 | 5.97×10⁻³ | — | 5.77×10⁻³ | — |
| 16 | 1.50×10⁻³ | 2.00 | 4.13×10⁻³ | 0.48 |
| 32 | 3.74×10⁻⁴ | 2.00 | 2.39×10⁻³ | 0.79 |
| 64 | 9.35×10⁻⁵ | 2.00 | 1.28×10⁻³ | 0.90 |
a = 2时，明明有错误却看起来是二阶收敛！
这就是MMS的陷阱。单侧差分的误差与边界处温度的二阶导数（曲率）成正比。cos(πx/2)在 x = 1 处曲率恰好为0，错误的误差也为0而被掩盖。取 a = 3 时边界处保留曲率，阶数0.48→0.79→0.90趋近1，一阶精度的错误就清楚显现。构造解要选得让所有项都不为零地起作用，且在边界上不取特殊值。惯用做法是把三角函数与指数函数、多项式混合，使用略显“不规则”的函数。
MMS的实施
步骤
- 选择光滑、且所有项（扩散、对流、非线性项、边界导数）都不为零的构造解。
- 用符号计算软件推导源项和边界条件，避免手算错误。
- 在三到四级网格上求解，求出误差和阶数。
- 确认阶数趋近理论值；不符时逐一隔离各项和边界条件查找原因。
“误差很小，阶数却不对”
有个例子：MMS误差低于0.1%，足够小就判为“合格”，后来查出阶数其实是1。误差大小可以靠加密网格任意减小，单凭这一点不能证明正确。即使有一阶精度的错误，网格加密误差也会变小，只是比理论慢。实际问题中往往保证不了那么细的网格，错误就会影响结果。MMS要看的不是误差大小，而是网格加密时误差减小的方式，即阶数。
常见错误
错误与对策
| 错误 | 后果 | 对策 |
|---|---|---|
| 构造解过于简单 | 错误被掩盖 | 选让所有项起作用的函数 |
| 按误差大小判定 | 漏掉低阶错误 | 按阶数判定 |
| 手工推导源项 | 把推导错误误认为代码错误 | 用符号计算推导 |
| 只用粗网格看阶数 | 未进入渐近区 | 用三到四级网格看趋势 |
| 不验证边界条件 | 边界错误残留 | 对每种边界分别测试 |
我想进一步了解MMS和收敛阶。
整体思路见MMS概述，阶数评估见收敛阶评估，源项推导见构造解源项，代码验证的整体定位见代码验证，实际问题的数值误差估计见Roache的GCI。还可以用Richardson外推工具算阶数，用二维导热工具试试二维导热。