Chapter 9: Numerical Methods and Precision

图形程序里的数值问题通常先表现为画面问题：一条运动轨迹开始抖动，一段曲线速度忽快忽慢，两个远处表面闪烁，粒子弹簧越算越发散，或者同一段动画在不同帧率下结果不同。本章要建立的能力是：看到这些现象时，能够把它们追踪到离散化、插值、浮点精度、时间步长、误差累计和修正策略，而不把所有异常都归因于 shader 写错或 GPU 性能不足。

贯穿本章的材料是一帧“沿曲线飞行的发光粒子”。粒子沿 Catmull-Rom 曲线移动，位置传入 vertex shader，fragment shader 输出发光颜色，物理侧用一个弹簧把相机或拖尾点拉向目标，最终深度缓冲决定前后遮挡。这个例子很小，却同时经过 CPU 更新、动画插值、GPU 浮点计算、深度测试和视觉调试。数值误差会在这些位置进入系统，并在屏幕上放大成可见异常。

数值方法（Numerical Methods）指用有限步骤、有限精度和近似规则求解连续问题的方法。在图形学里，连续问题包括曲线求值、运动积分、光照计算、深度比较、物理约束和图像重建。精度（Precision）指一个数值表示和一次计算能保留多少有效信息。二者连接起来后，渲染工程中的判断顺序变成：先确认输入范围，再确认近似模型，再确认时间步长和迭代规则，最后检查浮点表示、epsilon 和可视化证据。

本章不把数学写成独立课程。每个公式和判断都落回图形管线中的对象：曲线参数、顶点属性、uniform、buffer、depth value、动画状态、solver 迭代和屏幕像素。

9.1 Approximation and Iteration Error Sources

近似误差的入口是“连续目标被有限步骤表达”。沿曲线飞行的粒子在数学上可以处在曲线参数 $t$ 的任意位置，渲染时却只能在每一帧取一个有限的 $t$ 值，再把结果写入 buffer 或 uniform。帧之间的真实连续运动被一串离散状态替代，这个替代过程带来离散化误差。

离散化（discretization）把连续空间、连续时间或连续函数切成有限样本。曲线被采样成若干点、物理时间被切成多个 $\Delta t$、屏幕被切成像素网格、深度被写成有限位数。离散化解决的是“GPU 和 CPU 只能处理有限数据”的问题；它的代价是中间信息被近似表示。采样密度越低，曲线越容易出现折线感；时间步长越大，物理状态越容易跨过碰撞面或越过稳定区间。

迭代误差来自“用多次近似更新逼近目标解”。弹簧相机、布料约束、粒子碰撞和后处理中的某些滤波都可能使用迭代。以弹簧为例，目标是让当前位置 $x$ 逐步接近目标位置 $target$，速度 $v$ 根据力和阻尼更新。一次更新很难直接得到全局精确解，程序会在每帧执行一步或多步近似积分。迭代次数、时间步长和收敛条件共同决定误差大小。

下面这段简化代码展示了一个常见的半隐式 Euler 更新。它只保留误差进入状态的关键路径：加速度先改速度，速度再改位置，$\Delta t$ 同时影响两次更新。

struct SpringState {
    float x;
    float v;
};

void stepSpring(SpringState& state, float target, float stiffness, float damping, float dt) {
    float force = stiffness * (target - state.x) - damping * state.v;
    float acceleration = force;
    state.v += acceleration * dt;
    state.x += state.v * dt;
}

这段代码里有三类误差。第一类是模型误差：代码把复杂运动简化成一维弹簧，并假设质量为 1。第二类是截断误差：连续微分方程被一阶步进公式替代，时间步长越大，近似越粗。第三类是舍入误差：每次乘法和加法都在 float 范围内舍入，长时间累计后会改变最终位置。

Solver tolerance 是迭代求解器的停止阈值。它回答的问题是“误差小到什么程度就接受当前结果”。在碰撞修正或约束求解中，阈值过大时物体会停在可见穿透状态；阈值过小时迭代次数增加，结果仍受浮点下限限制。工程上需要把 tolerance 绑定到视觉尺度，例如角色半径、屏幕像素、碰撞壳厚度或世界单位；缺少尺度语境的小常数很难支撑稳定判断。

截断误差（truncation error）来自公式阶数。用一阶方法更新运动时，每一步只使用当前斜率近似一段曲线；真实曲线在这段时间内可能已经弯曲。二阶、四阶方法能降低局部截断误差，但会增加采样次数和状态读取成本。实时图形程序常在稳定性、可预测性、成本和视觉容忍度之间选择，不会单纯追求数学阶数最高。

在渲染帧里，误差通常沿下面的路径传播：CPU 用 $\Delta t$ 更新粒子参数，曲线求值得到世界坐标，坐标进入 shader 后经过矩阵变换，clip space 中的结果做透视除法，深度值写入 depth buffer。任何一步的误差都可能被后续变换放大。大坐标下的小位移丢失，会让顶点抖动；near/far 设置不合理，会让深度精度集中在相机附近，远处表面出现闪烁。

定位这类问题时，先把误差来源分层：输入数据是否已经离散过粗，时间步长是否随帧率剧烈变化，迭代次数是否不足，停止阈值是否与视觉尺度匹配，最后再检查 float 表示能否容纳当前数值范围。这个顺序能把“看起来随机”的画面异常拆回可验证对象。

9.2 插值法与样条曲线

插值的目标是在已有样本之间构造中间值。沿曲线飞行的粒子不会每帧手工存储一个位置，程序通常只保存关键点、关键帧或控制点，然后在运行时计算中间状态。插值质量直接决定运动是否平滑、速度是否稳定、方向是否突变，以及顶点属性是否在三角形内部连续变化。

线性插值（linear interpolation, lerp）在两个值之间按比例混合。给定 $a$、$b$ 和参数 $t$，公式是 $a(1-t)+bt$。它的优点是成本低、行为可预测、适合颜色、位置、标量权重等线性量。它的边界也清楚：当对象处在曲面、球面方向或旋转空间时，直接 lerp 可能改变长度或旋转速度。

下面的 GLSL 片段展示了 lerp 在 shader 中的典型位置。这里的 mix 对应线性插值，用于在两个颜色之间生成当前粒子的发光颜色。

vec3 baseColor = vec3(0.1, 0.4, 1.0);
vec3 hotColor = vec3(1.0, 0.8, 0.2);
float glow = clamp(particleEnergy, 0.0, 1.0);
vec3 finalColor = mix(baseColor, hotColor, glow);

这段代码的数值判断点在 glow。如果 particleEnergy 的范围没有提前固定，颜色插值会超出预期范围。clamp 把输入限制到 $0$ 到 $1$，使颜色混合保持在两个端点之间。它解决的是值域稳定性问题，并不解决能量计算本身是否合理。

球面线性插值（spherical linear interpolation, slerp）用于单位方向或四元数旋转。它沿球面最短弧移动，适合相机朝向、骨骼旋转和法线方向过渡。slerp 的成本高于 lerp，但能保持单位长度和角速度更稳定。对方向向量直接 lerp 后再 normalize 也常见，它成本较低，适合角度变化较小或视觉要求较宽的情况。

重心插值（barycentric interpolation）是光栅化阶段的重要插值方式。三角形三个顶点的属性进入 rasterizer 后，硬件会根据片元在三角形内的位置计算权重，再插值得到 fragment shader 的输入。位置、颜色、UV、法线都可能通过这条路径进入片元。透视投影下的属性插值还需要透视校正，否则纹理会在远近变化时发生扭曲。

Hermite 曲线使用端点和端点切线控制曲线形状。它适合动画曲线、速度可控的路径和需要指定进入/离开方向的运动。Catmull-Rom 曲线可以看成一种由相邻控制点自动推导切线的三次样条，适合让路径穿过给定点。沿曲线飞行的粒子使用 Catmull-Rom 时，控制点少，路径连续，编辑成本低。

Catmull-Rom 的关键边界是参数化方式。统一参数化会让相邻控制点距离不同的曲线出现速度变化和局部回折风险；centripetal 参数化能改善这类问题。工程中需要先问“控制点是否均匀分布”和“速度是否需要稳定”。如果只追求穿点和平滑，Catmull-Rom 很方便；如果需要严格控制速度、加速度或弧长，通常要做弧长表或改用更明确的动画曲线。

一条曲线路径的计算可以拆成三层。第一层是参数 $u$ 如何从时间得到，例如 u += speed * dt。第二层是曲线如何从 $u$ 求位置，例如 Catmull-Rom 在四个控制点之间求值。第三层是方向、速度和相机姿态如何从位置变化得到。抖动可能出在任意一层：时间步长跳变会让 $u$ 不稳定，控制点分布会让速度不均匀，方向求导会放大相邻帧的位置误差。

float segmentU = localTime / segmentDuration;
segmentU = std::clamp(segmentU, 0.0f, 1.0f);
Vec3 position = catmullRom(p0, p1, p2, p3, segmentU);
Vec3 tangent = normalize(catmullRomTangent(p0, p1, p2, p3, segmentU));

这段代码有两个检查点。segmentU 的范围必须稳定，因为曲线公式通常假设参数处在 $0$ 到 $1$。tangent 在控制点重合或曲线局部极短时可能接近零向量，normalize 会放大噪声。稳定实现会给切线长度设下限，并在异常时沿用上一帧方向或使用相邻段方向。

插值方法的选择可以按对象类型判断：标量和颜色优先 lerp，方向和旋转优先 slerp 或归一化后的向量插值，三角形内部属性依赖重心插值，路径和动画曲线使用 Hermite、Catmull-Rom 或更完整的样条系统。这个分类的目的在于把“平滑”拆成可检查约束：值域、长度、角速度、穿点、连续性和计算成本。

9.3 浮点运算精度与舍入误差分析

浮点数把实数近似存进有限位宽。图形程序大量使用 float，因为 GPU 对 32 位浮点通常有高吞吐，buffer、shader、矩阵和纹理格式也围绕它建立了成熟路径。float 的代价是有效位数有限，同样大小的位移在小坐标和大坐标中可见程度不同。

常见的 32 位浮点可理解为符号、指数和有效数字三部分。指数决定数值范围，有效数字决定相对精度。数值越大，相邻可表示数之间的间距越大。这意味着在世界坐标接近十万、百万时，小于某个间距的位移可能在表示时丢失。粒子在原点附近平滑移动，搬到大世界坐标后开始抖动，经常和这个规律有关。

舍入误差（rounding error）发生在每次计算结果写回有限精度表示时。单次误差很小，循环累计后会改变轨迹、能量或约束状态。粒子系统、积分器、骨骼动画和摄像机平滑都可能受到累计误差影响。累计求和尤其敏感，因为小增量加到大总量上时，小增量可能无法改变结果。

下面的例子展示了大坐标下的小位移问题。这个问题来自有限精度表示，和具体 API 无直接绑定。

float worldX = 1000000.0f;
float smallMove = 0.01f;
float nextX = worldX + smallMove;

如果当前数量级下 float 已经无法表达 0.01 级别的差异，nextX 可能和 worldX 表示成同一个值。渲染侧的表现是对象移动输入存在，屏幕位置却停住或跳格。大世界渲染常用 floating origin、camera-relative rendering 或分块局部坐标，把进入 shader 的坐标范围压回相机附近，从而保留局部位移。

深度缓冲也是精度问题的集中位置。透视投影会把 view space 的深度映射到 depth buffer。常见透视设置下，深度精度更多集中在 near plane 附近，远处可分辨间隔变粗。near 设得过小、far 设得过大时，远处两个接近的表面会争夺同一批深度值，表现为 z-fighting。具体 NDC 深度范围和剪裁约定受图形 API 影响，OpenGL、Vulkan、Direct3D 和 Metal 的细节需要按当前管线配置核对。

z-fighting 的排查顺序很明确。先检查两个表面是否真的共面或距离过小，再检查 near/far 比例，再检查 depth format，例如 24 位深度、32 位浮点深度或反向 Z 设置。最后检查是否存在多 pass 重画同一几何、polygon offset 使用不稳定、或者模型导入阶段产生重复面。每一步都对应可观察证据：线框重合、相机远近变化、depth visualization、draw call 顺序和 depth state。

矩阵乘法中的精度误差也会影响画面。模型矩阵、视图矩阵和投影矩阵连续相乘后，顶点从 model space 进入 clip space。矩阵元素范围差距过大时，小尺度局部变化可能被大平移项淹没。相机相对坐标能减轻这类问题：CPU 侧先用双精度或高精度场景坐标管理世界，提交到 GPU 前减去相机原点，把 shader 中的 position 保持在较小范围。

法线和方向向量的精度问题表现为光照闪烁。normal map、切线空间、插值法线和变换后的法线都需要接近单位长度。经过插值后，法线长度通常会偏离 1，fragment shader 中需要 normalize。这里的 normalize 承担输入修正作用，它把光照公式 dot(N, L) 的输入恢复到预期范围。长度漂移会改变亮度，尤其在高光和法线贴图中更明显。

浮点误差的工程判断可以归纳成三句话。第一，检查数值范围，因为范围决定可用精度。第二，检查累计路径，因为循环更新会放大小误差。第三，检查比较规则，因为浮点值很少适合用完全相等判断。epsilon 应该来自对象尺度和视觉容忍度，例如碰撞厚度、像素大小、单位向量长度误差或归一化阈值。

9.4 提高数值稳定性的检查顺序

数值稳定性指输入存在小扰动、步进存在近似、计算存在舍入时，输出仍保持在可接受范围内。对沿曲线飞行的粒子来说，稳定意味着路径连续、速度可解释、深度遮挡不闪烁、弹簧不过度发散、shader 输出不会出现 NaN 或极端颜色。稳定性来自检查顺序，单个万能公式无法覆盖不同管线阶段的误差来源。

第一步检查取值范围。每个进入 shader 或 solver 的量都应有明确范围：曲线参数在 $0$ 到 $1$，颜色权重在 $0$ 到 $1$，法线长度接近 $1$，roughness 为正且处在材质模型允许范围，深度相关输入和当前投影矩阵一致。范围错误会让后续算法处在未设计区域内，任何高阶方法都很难修复。

第二步检查归一化。方向、法线、四元数、权重和概率分布都依赖归一化约束。归一化解决的是“量的语义保持一致”：方向向量表达方向，权重表达比例，四元数表达旋转。缺少归一化时，后续公式会把长度误差当作物理或光照信息使用。图形调试中常把法线编码成颜色输出，通过屏幕颜色是否跳变观察归一化问题。

第三步检查时间步长。动画和物理更新应区分渲染帧率和模拟步长。渲染帧率可能随负载变化，直接使用可变大步长更新弹簧或碰撞会导致结果依赖帧率。常见做法是使用固定模拟步长累计推进，渲染时对两个模拟状态做插值。这样能让物理结果更可复现，也让视觉平滑和模拟稳定分开处理。

accumulator += frameDt;
while (accumulator >= fixedDt) {
    previousState = currentState;
    stepSimulation(currentState, fixedDt);
    accumulator -= fixedDt;
}

float alpha = accumulator / fixedDt;
RenderState renderState = interpolate(previousState, currentState, alpha);

这段代码把 simulation update 和 render interpolation 分开。stepSimulation 使用固定 fixedDt，减少帧率变化对物理结果的影响；渲染状态通过 alpha 插值得到，减少画面跳动。检查这段逻辑时要关注两点：frameDt 需要上限，长时间卡顿后不应一次补过多步；插值只用于显示状态，不能反向写回模拟状态。

第四步检查 epsilon。epsilon 是“足够接近”的阈值，它必须带有单位和语境。判断向量长度是否接近零时，epsilon 对应世界单位；判断 UV 是否相等时，epsilon 对应纹理坐标范围；判断深度是否接近时，epsilon 对应当前 depth mapping 和 depth format。一个全局 1e-6 常量很难同时适配所有问题。

第五步检查守恒量或约束量。物理模拟里可以观察能量、速度上限、距离约束、碰撞穿透深度；动画系统里可以观察骨骼长度、权重总和、四元数长度；渲染里可以观察颜色范围、法线长度、roughness 下限和深度分布。守恒量不一定严格不变，但它们应该在可解释范围内变化。范围突然爆炸通常说明某次更新进入了错误分支。

第六步使用可视化输出。数值问题经常无法从最终美术效果直接定位，需要把中间量画出来。法线可以映射成 RGB，深度可以线性化后显示，速度可以用箭头或热力颜色表示，曲线参数可以用渐变色表示，NaN 可以输出成醒目颜色。可视化的价值在于把 buffer、shader 输入和 solver 状态转成屏幕证据。

下面的检查表可以作为图形工程中的默认顺序。

顺序	检查对象	需要回答的问题	常见视觉证据
1	取值范围	输入是否处在算法设计范围内	颜色溢出、曲线跳段、位置突变
2	归一化	方向、权重、旋转是否保持语义	光照闪烁、姿态缩放、方向翻转
3	时间步长	模拟是否被帧率变化影响	弹簧发散、碰撞穿透、动画忽快忽慢
4	epsilon	阈值是否匹配对象尺度	抖动、吸附失败、接触面闪烁
5	守恒量	状态是否超出可解释范围	能量增长、约束拉长、速度爆炸
6	可视化	中间量能否被屏幕验证	debug view 中出现断层、噪点或异常色

这个顺序的核心是从输入边界走向输出证据。先修正范围和约束，再讨论算法阶数和优化。很多画面异常在前两步就能定位，复杂 solver 和 GPU 工具分析应该放在基础数值条件确认之后。

9.5 动画插值与物理模拟中的误差控制

动画插值和物理模拟都会把误差带入连续画面。动画系统关注“状态之间如何过渡”，物理系统关注“状态如何随时间演化”。沿曲线飞行的粒子同时使用二者：路径位置由插值得到，拖尾或相机跟随由弹簧积分得到。二者的误差控制方式不同，混在一起处理会让问题更难定位。

姿态插值的主要风险是空间选错。位置可以用 lerp，旋转更适合四元数 slerp 或规范化后的四元数插值，缩放需要单独处理。骨骼动画中，一个关节的旋转误差会沿层级传播到子节点，末端手指、武器或相机挂点会放大微小误差。调试姿态时，应先看局部关节，再看父子组合矩阵，最后看 skinning 后的顶点位置。

动画曲线还要处理时间参数。关键帧时间分布不均匀时，按关键帧索引均匀推进会让速度异常。更稳定的做法是使用真实时间查找当前 keyframe 区间，再在区间内计算局部参数。路径动画需要进一步关注弧长，因为曲线参数均匀并不等于空间速度均匀。需要恒定速度时，可以预计算弧长表，用行进距离反查曲线参数。

物理模拟的误差控制从积分器开始。显式 Euler 成本低，但在刚性弹簧和大时间步长下容易增长能量；半隐式 Euler 常用于实时场景，稳定性更好；Verlet、RK 系列和约束投影适合不同模拟对象。选择积分器时要看对象：粒子拖尾、相机弹簧、布料、刚体碰撞和流体都对稳定性、能量表现和性能预算有不同要求。

碰撞修正会引入另一类误差：把已经穿透的状态推回合法区域。若每帧只做一次粗略修正，物体可能停在浅穿透状态；若修正过强，物体会抖动或弹飞。稳定做法通常包含接触容差、位置修正比例、速度投影和迭代次数控制。所有参数都应以世界尺度、对象半径和视觉容忍度解释。

沿曲线粒子的拖尾可以用弹簧点跟随主粒子。拖尾点追随目标时，如果 stiffness 过高、damping 过低、dt 过大，就会出现明显振荡。调参的顺序应先固定 dt，再调 damping 抑制振荡，再调 stiffness 控制响应速度，最后检查速度上限和位置误差。若先在可变帧率下调参数，得到的结果很可能随机器负载变化。

float maxFrameDt = 1.0f / 30.0f;
float frameDtClamped = std::min(frameDt, maxFrameDt);
stepSpring(trailPoint, targetX, stiffness, damping, frameDtClamped);

这段代码给帧时间设置上限，减少单次大步长把弹簧推入不稳定状态的概率。它适合简单视觉跟随，完整固定步长模拟还需要 accumulator、固定迭代和渲染插值。对于碰撞、约束或多人同步场景，固定步长和状态插值更可靠。对于纯视觉特效，上限裁剪加合理阻尼通常已经能满足需求。

误差还会从 CPU 和 GPU 的边界进入。CPU 可能用 double 保存编辑器和世界状态，GPU buffer 使用 float 渲染。提交数据时如果直接把大世界坐标写进 vertex buffer，shader 中的局部细节会丢失。更稳的路径是：CPU 侧保存高精度世界坐标；提交前转换到相机相对坐标；shader 中完成局部变换、投影和光照。这样可以把 float 的有效精度用在屏幕附近的几何细节上。

动画和物理的最终检查点是视觉结果是否可复查。给路径输出曲线 debug 线，给每个控制点画编号，给速度画颜色梯度，给碰撞点画法线，给深度输出 debug view，给 NaN 设置专门颜色。数值稳定需要可观察证据确认：先把中间状态转成屏幕结果，再根据范围、时间步长、约束和精度逐项缩小问题。

本章建立的核心理解是：图形程序中的数值误差来源分布在多条路径上。它可能来自连续问题的离散化、插值空间选择、浮点位宽、时间步长、迭代停止条件、epsilon 和大世界坐标。稳定工程实践的关键，是把一帧中的视觉异常放回数据路径和管线阶段中，用固定顺序排查；直接替换算法或堆叠参数会掩盖真正的误差入口。

最小自检任务

你正在调试一个粒子拖尾效果：主粒子沿 Catmull-Rom 曲线飞行，拖尾点用弹簧跟随主粒子。画面中出现三个现象：曲线拐弯处速度突然变快，拖尾点在低帧率时振荡明显，远处地面和粒子阴影偶尔闪烁。请按本章的数值稳定性思路给出排查顺序，并说明每一步要观察的证据。

答案要点

先检查曲线参数和控制点分布。拐弯处速度变快通常说明参数均匀推进没有匹配空间弧长，或者 Catmull-Rom 控制点距离差异较大。证据是把曲线参数映射成颜色或在曲线上等时间采样画点，若点距明显变化，问题在路径参数化或弧长控制。

再检查弹簧更新的时间步长。低帧率振荡说明 dt 变大后积分进入不稳定区域。证据是记录每帧 frameDt、弹簧速度、位置误差和 overshoot 次数。稳定处理可以使用固定模拟步长，或至少给视觉弹簧设置 frameDt 上限，并按固定 dt → damping → stiffness → velocity limit 的顺序调参。

然后检查拖尾状态的约束和 epsilon。需要确认切线 normalize 前的长度下限、弹簧目标距离、速度上限和停止阈值是否与世界单位匹配。证据是 debug 输出中是否出现接近零向量、NaN、速度突增或误差在阈值附近来回跳动。

最后检查深度精度和几何关系。远处地面与阴影闪烁通常要观察表面是否共面、near/far 设置是否过宽、depth format 是否满足场景范围，以及是否多 pass 重画同一平面。证据包括线框重合、depth debug view、相机远近变化时闪烁强度变化、draw call 顺序和 depth state。

本章知识点总结

• 离散化误差：连续曲线、连续时间和连续空间进入图形程序后会变成有限样本，采样密度和时间步长决定近似质量。
• 迭代误差：弹簧、碰撞和约束求解通过多次近似更新逼近目标，迭代次数和停止阈值会影响最终状态。
• 截断误差：积分公式用有限阶近似替代连续演化，时间步长越大，单步近似和真实轨迹之间的差异越明显。
• Tolerance：停止阈值必须绑定对象尺度和视觉容忍度，随手设置的小常数很难适配不同图形问题。
• Lerp：线性插值适合标量、颜色和位置等线性量，关键检查点是参数范围和值域稳定。
• Slerp：球面线性插值适合方向和旋转，可以保持单位长度和角度变化更稳定。
• 重心插值：三角形内部属性由重心权重插值得到，透视投影下需要关注透视校正和属性连续性。
• 样条曲线：Hermite 和 Catmull-Rom 适合路径与动画曲线，参数化方式会直接影响速度稳定性。
• 浮点范围：float 的相对精度随数值量级变化，大世界坐标会降低局部小位移的可表示能力。
• 深度精度：near/far 比例、depth format、投影约定和共面几何共同决定 z-fighting 的风险。
• 归一化约束：方向、法线、四元数和权重需要维持语义范围，长度漂移会改变光照、姿态和插值结果。
• 固定步长：物理模拟使用固定步长能降低帧率变化对结果的影响，渲染插值负责补足视觉平滑。
• Epsilon：epsilon 应来自世界单位、纹理坐标、深度范围或像素尺度，每个比较场景都需要自己的阈值语境。
• 可视化证据：法线、深度、速度、曲线参数和 NaN 标记应转成 debug view，用屏幕证据定位数值问题。
• 排查顺序：稳定性问题应按取值范围、归一化、时间步长、epsilon、守恒量和可视化输出逐步检查。