MMS收敛阶的计算方法
理论:由误差之比读出阶数
概述
用MMS算出误差后,阶数怎么计算?只是简单相除吗?
公式很简单,但陷阱不少:弄错网格比、加密过头碰到舍入误差底限、用单精度计算、解不光滑、只在边界用了低阶公式——都会使阶数偏离理论值。要把它们与代码缺陷区分开,就得知道各自表现如何。我们用二阶中心差分求解一维泊松方程逐一确认。
观测阶数公式
$e_1, e_2$ 为粗、细网格误差,$r$ 为网格间距之比;只有 $r = 2$ 时才能写成 $\log_2(e_1/e_2)$。
舍入误差底限
二阶导数差分要除以 $h^2$,所以 $h$ 越小舍入误差越大。比离散误差与舍入误差之和最小的网格更细时,误差反而增大。
舍入误差与网格细度
数值微分中步长太小会因舍入误差导致精度下降,这从计算机诞生初期就已知道。浮点数表示在1985年统一为IEEE 754,双精度约16位有效数字,单精度约7位。随着GPU普及,单精度和半精度计算重新被使用,加密网格精度却上不去的问题如今又成了身边的注意事项。
算例:一维泊松方程
−u″ = f,制造解 u = sin(πx) + x³,二阶中心差分,两端狄利克雷条件,误差取最大值范数。
双精度与单精度
| 等分数 | 双精度误差 | 阶数 | 等分数 | 单精度误差 | 阶数 |
|---|---|---|---|---|---|
| 100 | 8.23e-5 | 2.00 | 30 | 9.14e-4 | 2.00 |
| 1,000 | 8.22e-7 | 2.00 | 100 | 8.27e-5 | 2.00 |
| 10,000 | 8.24e-9 | 2.00 | 300 | 3.03e-5 | 0.91 |
| 100,000 | 3.26e-10 | 1.40 | 1,000 | 4.01e-4 | −2.15 |
| 1,000,000 | 6.85e-7 | −3.32 | 3,000 | 2.90e-3 | −1.80 |
其他陷阱
| 情况 | 得到的阶数 | 应有阶数 |
|---|---|---|
| 解为x^1.5(最大值范数,512等分) | 1.44 | 2(光滑解时) |
| 解为x^1.5(L2范数,512等分) | 1.47 | 2 |
| 网格比1.5用log₂计算 | 1.17 | 2.00 |
| 诺伊曼端用单侧差分(一阶) | 0.99 | 2 |
| 诺伊曼端用虚点(二阶) | 2.00 | 2 |
即使双精度,100万等分时误差也会增大。并不是越细越好。
二阶导数差分把相邻值的微小差除以h²。h变小时,这个差接近舍入误差大小,再被除法放大。双精度约在10万等分碰到底限,单精度仅300等分就失效。用MMS确认阶数时,只用阶数稳定在理论值的范围,碰到底限的点要剔除。另一件要紧的事是阶数下降时要分辨原因:解不光滑时稳定在1.44这样不上不下的值;边界公式是一阶时阶数干脆降到1;舍入误差则会突然崩溃变成负值。这些表现的差异是区分缺陷与其他原因的线索。
求阶数的步骤
- 以恒定比构造四级以上网格,记录比r。
- 按 p = ln(e₁/e₂)/ln r 对每对相邻网格求阶数并列表。
- 把阶数稳定在理论值的范围作为“渐近范围”使用。
- 细网格上误差增大时怀疑舍入误差底限,提高精度确认。
- 制造解用光滑函数,同时确认边界条件公式的阶数。
“移到GPU后阶数崩了”
某团队为在GPU上加速,把导热代码切换为单精度,结果MMS验证中细网格的阶数变成负值。他们怀疑移植出错查了好几天,原因却是单精度舍入误差;改回双精度后阶数稳定在2.00。于是改为计算主体用单精度、只有残差计算和部分求和用双精度,在保持速度的同时获得了所需精度。
常见错误
错误与对策
| 错误 | 后果 | 对策 |
|---|---|---|
| 网格比不是2却用log₂ | 阶数算错 | 除以ln r |
| 加密过头 | 舍入误差底限 | 只用渐近范围 |
| 用单精度验证 | 过早失效 | 用双精度确认 |
| 制造解不光滑 | 阶数不上不下 | 用光滑函数 |
| 不看边界公式阶数 | 整体变为一阶 | 边界也用同阶 |
我想了解相关内容。
相关页面参考MMS概述、MMS:二维弹性体、MMS:二维导热、代码验证和Richardson外推。