Float16 舍入误差如何影响热方程的 von Neumann 稳定性分析？
方法二适合作为起点：按照实际执行顺序逐步建立舍入模型，将每一步的舍入误差整理为附加扰动，再分析差分格式如何传播这些扰动。
对于常系数、均匀网格上的经典分析，精确算术下的条件仍是 0 ≤ r ≤ 1/2。Float16 会引入额外误差，误差也会累积，但仅凭精度格式无法得到通用的“修正 CFL 条件”。还需要明确运算顺序、边界条件、系数计算方式，以及希望保证哪一种稳定性。
采用什么浮点误差模型？
为避免与温度变量 u 混淆，下面用 μ 表示 unit roundoff。
假设 Float16 指 IEEE 754 binary16，并采用 round-to-nearest。binary16 的有效精度为 11 个二进制位，因此：
[ \mu=2^{-11}\approx4.8828125\times10^{-4}. ]
它与 1 到下一个浮点数之间的间隔不同，后者是 2⁻¹⁰。NVIDIA 的 FP16 格式说明给出了 binary16 的位数布局。
在没有溢出且下溢不破坏相对误差界的前提下，四则运算均可写为：
[ \operatorname{fl}(x\mathbin{\mathrm{op}}y) =(x\mathbin{\mathrm{op}}y)(1+\delta), \qquad |\delta|\le\mu, \quad \mathrm{op}\in{+,-,*,/}. ]
原问题给出的加减法模型也可以使用。在上述前提下，取 α = β = δ 就能得到其中一种有效表示。逐步推导时，为每个运算结果对应一个 δ 通常更方便。Cornell 的浮点误差分析讲义说明了这种模型及其适用限制。
相减发生消去时，输入已有的误差相对于较小的差值可能变得很大，但不能因此直接认定标准舍入模型失效。若出现 subnormal 或 flush-to-zero，应另外处理绝对误差，不能继续无条件套用纯相对误差界。
如何按方法二建立模型？
用 ρ 表示程序实际使用的浮点系数，它可能与数学定义的 r = aΔt/Δx² 不完全相同。暂时固定一个网格点和时间层，记：
[ v_-=\tilde u_{i-1}^n,\qquad v_0=\tilde u_i^n,\qquad v_+=\tilde u_{i+1}^n. ]
方法二对应的伪代码为：
# fl16 表示将该次运算的精确结果舍入到 binary16
# 所有输入均为已经存储的浮点值
t1 = fl16(2 * v0)
t2 = fl16(vplus - t1)
t3 = fl16(t2 + vminus)
t4 = fl16(rho * t3)
vnew = fl16(v0 + t4)
对每个舍入点分别建模：
[ \begin{aligned} t_1&=2v_0(1+\delta_1),\ t_2&=(v_+-t_1)(1+\delta_2),\ t_3&=(t_2+v_-)(1+\delta_3),\ t_4&=\rho t_3(1+\delta_4),\ \tilde u_i^{n+1}&=(v_0+t_4)(1+\delta_5), \end{aligned} \qquad |\delta_j|\le\mu. ]
这些 δ 随网格点、时间步和数据变化，不需要假定它们相互独立或均值为零。
对于 binary16，浮点数乘以 2 在结果不溢出时可以精确表示，因此这里可取 δ₁ = 0。代回可得：
[ \tilde u_i^{n+1} =(1+\delta_5)v_0 +\rho(1+\delta_5)(1+\delta_4)(1+\delta_3) \left[(1+\delta_2)(v_+-2v_0)+v_-\right]. ]
这个表达式保留了各舍入因子的来源。方法一遗漏了最终加法等运算的舍入影响，不能直接作为该执行顺序的完整模型。若要使用整体扰动形式，应先从逐步模型推导其系数和误差界。
将舍入误差整理为每步扰动
定义精确算术下、使用实际系数 ρ 的差分算子：
[ (A_\rho v)i =(1-2\rho)v_i+\rho v{i+1}+\rho v_{i-1}. ]
浮点计算可以写成：
[ \tilde u^{n+1}=A_\rho\tilde u^n+\eta^n, ]
其中 ηⁿ 是这一时间步全部舍入操作产生的残差。该等式可以精确定义 ηⁿ，不要求将舍入误差视为随机噪声。
例如，在上述逐步模型成立时，令：
[ \gamma_k=\frac{k\mu}{1-k\mu},\qquad k\mu<1. ]
利用至多四个舍入因子的乘积界，可以得到一个保守的局部绝对误差界：
[ |\eta_i^n| \le\mu|v_0| +|\rho|\gamma_4 \left(|v_+|+2|v_0|+|v_-|\right). ]
因此，在解幅值有界时，每步残差可用 O(μ) 乘以局部解幅值来控制。这个界没有利用误差之间的抵消，通常会高估实际误差。
这里的 ηⁿ 只包括使用 ρ 更新时的算术舍入误差。若要与使用数学系数 r 的格式比较，还需加入系数误差：
[ \tilde u^{n+1} =A_r\tilde u^n +(\rho-r)L\tilde u^n+\eta^n, \qquad (Lv)i=v{i+1}-2v_i+v_{i-1}. ]
如何纳入 von Neumann 稳定性分析？
先假设采用周期边界、均匀网格和常系数。在这些条件下，Fourier 模态能够对角化差分算子。其他边界条件需要分析相应的边界离散算子。
对模态 exp(i iθ)，精确算子的放大因子为：
[ G_\rho(\theta) =1+\rho(e^{i\theta}-2+e^{-i\theta}) =1-4\rho\sin^2(\theta/2). ]
要求所有波数满足 |Gρ(θ)| ≤ 1，就得到经典条件：
[ 0\le\rho\le\frac12. ]
这一推导可参见 MIT 的 von Neumann 稳定性分析讲义。
接着令 Uⁿ 表示使用同一系数 ρ、但在精确算术下计算的离散解，并定义：
[ e^n=\tilde u^n-U^n. ]
两种更新相减后得到：
[ e^{n+1}=A_\rho e^n+\eta^n. ]
对它作离散 Fourier 变换：
[ \widehat e^{,n+1}(\theta) =G_\rho(\theta)\widehat e^{,n}(\theta) +\widehat\eta^{,n}(\theta). ]
于是：
[ \widehat e^{,n}(\theta) =G_\rho(\theta)^n\widehat e^{,0}(\theta) +\sum_{m=0}^{n-1} G_\rho(\theta)^{n-1-m}\widehat\eta^{,m}(\theta). ]
其中，Gρ 描述已有误差的传播，η 描述每一步新产生的误差。
舍入操作依赖数据，通常也会破坏空间平移不变性，因此完整的浮点更新一般不能用一个固定的“带舍入误差放大因子”来描述。上面的受迫递推仍然成立，但误差源本身可能包含不同 Fourier 模态。
满足 CFL 条件后，误差会怎样累积？
在周期边界下，若 0 ≤ ρ ≤ 1/2，则 Aρ 在离散 2-范数下满足：
[ |A_\rho^k|_2\le1. ]
因此：
[ |e^n|_2 \le|e^0|2 +\sum{m=0}^{n-1}|\eta^m|_2. ]
如果另有每步残差上界 ‖ηᵐ‖₂ ≤ E，就得到：
[ |e^n|_2\le|e^0|_2+nE. ]
这是最坏情形界。由于局部残差界依赖计算解的幅值，用它估计 E 时，还需要控制相关时间区间内的解幅值。没有随机性假设时，不能直接将累计误差写成 O(√n μ)。
对于 |Gρ(θ)| < 1 的模态，历史误差会衰减。但周期边界下 Gρ(0) = 1，零频误差没有这种衰减；当 ρ = 1/2 且网格包含 θ = π 时，最高频模态还满足 Gρ(π) = −1，同样没有幅值衰减。
舍入误差会改变 CFL 条件吗？
需要区分三个结论：
- 经典差分算子的稳定区间仍是 0 ≤ ρ ≤ 1/2。这是精确算子传播误差的性质，浮点精度不会改变上述 Gρ 的表达式。
- 程序应核查实际使用的 ρ。如果 r 通过多次浮点运算计算，还应单独分析这些运算的误差，不能只检查纸面上的 r。
- 满足 CFL 条件不等于舍入误差永远有界。它控制齐次误差传播；持续注入的舍入误差仍可能累积。
因此，不能仅凭 Float16 就给出通用的 r ≤ 1/2 − Cμ，并声称这保证了完整浮点算法在任意步数下稳定。选择严格小于 1/2 的实际系数，可以为高频模态提供衰减余量，但周期问题的零频模态仍没有衰减。
如果研究目标是在固定 μ 下不断细化网格，还应注意：时间步数增加后，累计舍入误差界可能变大。精确算术下的一致性与稳定性，并不自动保证固定低精度计算的误差随网格细化持续下降。
实现时应确认哪些细节？
- 确认 Float16 是存储格式还是实际运算格式。中间结果若使用更高精度，舍入点和误差界会变化。
- 确认编译器是否采用 FMA 或重排表达式。融合运算只在融合结果处舍入，模型应随实际执行方式调整。
- 使用上一时间层的数据计算全部新值，避免原地覆盖意外改变差分格式。
- 单独处理溢出、subnormal 和 flush-to-zero。即使最终更新值有限，中间量也可能超出表示范围。
- 留意小增量被最终加法舍去的现象。时间步过小时，计算可能出现数值停滞，因此减小 Δt 并不总能提高 Float16 结果的准确度。
这里分析的是浮点解相对于精确离散解的误差。若要评价相对于热方程真解的总误差，还需加入离散截断误差以及初始、边界数据的相关误差。
备注：内容仅供参考。