数值方法引论

该节由 Alice_Rain 编写

Wir müssen wissen. Wir werden wissen. —— David Hilbert（希尔伯特）

我们必须知道，我们必将知道。

1. 数值方法的诞生

各位同学一定对集合论并不陌生。19世纪末，康托尔的集合论一度被认为可以为整个数学提供统一基础。但 1901 年，伯特兰·罗素提出了著名的罗素悖论：

设集合 $S = \{x \mid x \notin x\}$ （即所有不属于自身的集合组成的集合）。

那么 $S \in S$ 是否成立？

如果 $S \in S$ ，根据定义， $S$ 不应属于 $S$ 。
如果 $S \notin S$ ，根据定义， $S$ 应该属于 $S$ 。

这个悖论的重要性不在于它本身有多“猎奇”，而在于它暴露出一个更深的问题：数学究竟建立在什么基础上？什么样的对象才算“存在”？一个结论如果不能被明确构造出来，是否仍然算真正的数学结论？

20 世纪上半叶，数学界围绕这些问题展开了长期讨论，其中影响最大的三种立场通常被概括为逻辑主义、直觉主义和形式主义。它们并不是数值方法的直接“发明者”，但它们重新定义了人们对证明、构造与算法的理解，而这些观念后来都汇入了数值计算的思想中。

I. 逻辑主义（Logicism）

逻辑主义试图说明：数学本质上是逻辑的延伸，数学概念可以被还原到一套严格的符号体系中。

这种追求形式化的努力虽然没有彻底统一数学，却留下了一个极其重要的遗产：数学过程可以被写成规则、符号和步骤。对计算机科学来说，这种思想后来发展为形式系统、程序语义和自动推理；对数值方法来说，它提醒我们，算法并不是“灵感”，而是一套可以被精确定义和重复执行的操作流程。

II. 直觉主义（Intuitionism）

直觉主义强调：数学对象不是凭空“存在”的，只有当你能把它构造出来时，它才真正有意义。换句话说，存在即被构造。

这与数值方法的关系非常直接。对于许多实际问题，我们并不满足于“某个解存在”的抽象证明，而更关心：能不能设计一个有限步骤的算法，把这个解逐步逼近出来？例如，当我们说某个方程有解时，数值方法真正追问的是：能否通过迭代求出它的小数近似，能否把误差控制在可接受范围内？

III. 形式主义（Formalism）

形式主义则更进一步地把数学看作一套按规则操作符号的系统。只要规则清晰、推演合法，数学就可以在纸面上被机械地执行。

这种思想直接推动了“可计算性”问题的提出，并最终通向图灵机与现代计算机的理论基础。也正是在这个意义上，数值方法才真正拥有了现实载体：连续的函数、导数和积分，不再只停留在黑板上，而是被转化成计算机能够逐步执行的离散运算。

从这个角度看，数值方法并不是凭空出现的一套技巧，而是现代数学与计算思想共同发展的结果。逻辑主义强调形式化，直觉主义强调可构造性，形式主义强调规则化与可计算性；三者汇合之后，便自然导向了数值方法最核心的立场：面对复杂的连续问题，我们关心的不是“是否存在一个完美的解析答案”，而是“能否设计一个可执行、可逼近、可检验误差的算法”。

2. 数值方法是什么，为什么需要它？

数值方法（Numerical Methods）是利用计算机求解数学问题近似解的算法与理论，是数学和计算机领域的结合体，是求解现实问题的桥梁，它属于计算数学核心领域，又称数值分析。它将复杂的连续数学模型（如微分方程）离散化为有限的算术运算，用于解决科学工程中无法求出解析解的难题，核心在于平衡精度、收敛性和计算效率。

如果你有志于使用计算机进行物理模拟，准备数学建模比赛，或者对图形着色器（shader）的编写感兴趣，对人工智能的底层算法的实现好奇，那么学习数值方法对你是很有用的。

首先我们需要在数学上达成一个共识：许多看似简单的积分、微分方程和偏微分方程，都无法写成我们熟悉的闭式解析解。

三体问题： 庞加莱证明了三体问题没有一般的解析解（无法写成简单的公式）。
物理上的微分方程： 现实中的许多微分方程都没有通用的解析解，例如 N-S 方程，它描述了流体运动的基本规律。我们不仅难以写出它的通解，甚至连“解是否始终存在且光滑”都是著名的千禧年数学难题。但借助数值方法，我们仍然可以通过有限体积法（FVM）让飞机在虚拟风洞里跑起来。
微积分： 运算在初等函数环上不封闭，很容易产生非初等解。例如 $e^{-x^2}$ 这样极其简单的函数，其原函数（高斯误差函数）就无法表示为初等函数的组合。这在概率论（正态分布）中随处可见。

解析数学研究的是连续对象，但计算机只能处理有限步、有限精度的离散数据。我们不可能在机器里取无穷多个点，也不可能无限精确地表示实数，因此必须把连续问题转化为计算机能够处理的离散问题。在数值方法中，我们常见的处理方式包括：

采样 (Sampling)： 我们只选取 $x_0, x_1, \dots, x_n$ 处的值。
插值： 在两个采样点之间，我们假设它是直线（线性插值）或某种简单的多项式。
网格 (Mesh)： 在 2D/3D 模拟中，这意味着将连续空间切分成一个个三角形或四边形。

离散化处理必然会带来误差。减小步长 $h$ 往往可以降低截断误差，但 $h$ 并不是越小越好，因为计算机的内存和算力是有限的，而且 float 或 double 遵循 IEEE 754 标准，只能存储有限位数的有效数字。

Tips

在计算机中，1.0 是可以用二进制完美表示的数，但是通过计算得到的 1.0 往往是不准确的。你可以尝试运行下方的代码，看看结果吧！

main.c

double x = 0.0;
for(int i = 0; i < 10; i++) x += 0.1;
if (x == 1.0) printf("Equal");
else printf("%.20f", x);

这是因为 $0.1$ 在二进制下是无限循环小数，计算机存储的是它的近似值。因此，浮点数比较通常不应直接使用 == 1.0，而应使用带容差的判断，例如 $fabs(x - 1.0) < \varepsilon$ ，其中容差需要结合具体问题的尺度来选择。

在实数系统中，不存在“比 1 大的下一个实数”；但在计算机的浮点数系统中，1.0 附近可表示的数是离散分布的，因此 1.0 与下一个可表示数之间存在一个最小间隔。这个量通常与 机器 epsilon ( $\epsilon_{mach}$ ) 有关。

对于 double 来说， $\epsilon \approx 2.22 \times 10^{-16}$ 。

当你尝试计算 $1.0 + 10^{-16}$ 时，结果可能因为有效位数不足而仍被舍入为 1.0。也就是说，当 $h$ 小到一定程度后，截断误差虽然继续减小，但舍入误差未必同步减小，甚至可能在大量运算中不断累积。

当然，如果 $h$ 取得过大，数值解也可能明显偏离真实解，甚至出现不稳定现象。这并不意味着“高阶算法天然更容易发散”，而是说明任何算法的精度结论都有适用范围；当步长过大时，原本依赖小步长成立的误差分析与稳定性结论都可能失效。在解释这一点前，先需要介绍数值方法里算法的概念，以及衡量它们准确性的误差。

3. 初值问题及其典型模型

在数值方法中，一个非常核心的任务是处理初值问题（IVP）：已知系统在 $t_0$ 时刻的状态，利用变化律（微分方程）推演它在未来时刻的状态。以下是几个经典的动力学模型：

天体轨道方程（万有引力）

根据牛顿第二定律和万有引力定律，天体的运动可以表示为二阶常微分方程：

\frac{d^2 \mathbf{r}}{dt^2} = -\frac{GM}{r^3} \mathbf{r}

在数值模拟中，我们通常通过引入速度变量 $\mathbf{v}$ ，将二阶方程改写为关于位置与速度的一阶 ODE 方程组。这也是三体问题难以写出通用解析解、必须依靠数值迭代的原因之一。

捕食者-被捕食者模型（Lotka-Volterra）

在生物数学或游戏 AI 的生态模拟中非常常用。它描述了两个物种（如狼 $y$ 与羊 $x$ ）数量随时间的变化：

\frac{dx}{dt} = \alpha x - \beta xy

\frac{dy}{dt} = \delta xy - \gamma y

由于方程包含非线性项 ( $xy$ )，它没有简单的初等函数解，必须依靠积分器来观察物种数量波动的动态平衡。

热传导方程（Heat Equation）

如果你想在模拟物体表面的温度扩散或烟雾消散，这就是典型的偏微分方程（PDE）：

\frac{\partial u}{\partial t} = \alpha \nabla^2 u

在计算时，我们会将空间离散化为网格，将拉普拉斯算子 $\nabla^2$ 转换为邻居节点的差分运算，从而将 PDE 降级为计算机可解的代数运算。

在常微分方程初值问题中，完成这种时间步进过程的算法通常被称为积分器。下面列出几类常见积分器，在后续教学中我们会逐步介绍它们的思想与用法。

显式欧拉法 (Explicit Euler)： 数值计算的“开山鼻祖”。它的核心思想是局部线性化，即假设在极短的时间步长 $h$ 内，物体的运动速度（斜率）是不变的。
四阶龙格-库塔法 (RK4)： 在一个步长 $h$ 内进行四次“试探性”采样，通过加权平均得到一个更稳健的步进方向，常用且可靠。
辛积分器 (Symplectic Integrator)： 在长时间的物理模拟（如行星运行）中，普通积分器（如 RK4）会产生“能量漂移”，导致轨道慢慢向外螺旋或向内塌缩。辛算法能够保持相空间的几何性质，其最大的特点是长期能量守恒。
自适应步长积分器（Adaptive Step-size）： 固定的步长 $dt$ 往往是低效的：在物理状态剧烈变化（如碰撞瞬间）时需要极小步长，而在平稳滑行时可以用大步长。它会同时运行两套精度不同的算法（如 RK4 和 RK5），比较两者的差值。如果差值（估计误差）太大，就减小步长重算；如果差值很小，就增大步长加速。

几乎所有步进式积分器的都依赖泰勒展开 (Taylor Series)。为了计算方便，我们不得不砍掉无穷级数后面的项，这些被“截断”的部分就是截断误差。

f(x+h) = \underbrace{f(x) + f'(x)h}_{\text{线性近似}} + \underbrace{\frac{f''(x)}{2!}h^2 + \dots + \frac{f^{(n)}(x)}{n!}h^n}_{\text{高阶修正}} + \underbrace{R_n(x)}_{\text{余项/误差}}

显式欧拉法只保留到一阶项，其局部截断误差是 $O(h^2)$ 。
RK4 巧妙地保留到了四阶项，其局部截断误差是高阶无穷小 $O(h^5)$ 。

模拟总时间为 $T$ 的过程，需要走 $N = T/h$ 步。每一步的误差都会逐渐累积。对常见单步法而言，在适当的光滑性与稳定性条件下，全局误差通常会比局部截断误差低一阶。

欧拉法：局部 $O(h^2)$ $\to$ 全局 $O(h^1)$ （一阶算法）。
RK4：局部 $O(h^5)$ $\to$ 全局 $O(h^4)$ （四阶算法，所以叫 RK4）。

当为了追求性能而大幅增加步长时，高阶算法并不会自动保持原有的精度优势，因为它们的误差分析通常建立在 $h$ 足够小的前提下。若步长不再足够小，高阶截断项和稳定性问题都可能变得显著，数值解也可能迅速偏离真实轨道，甚至直接“飞出屏幕”。这说明高阶方法更强，并不等于它在任意步长下都更安全。

在接下来的教程中，我们会从插值与拟合讲起，再逐步介绍一些典型的积分器，帮助同学们建立对数值方法的基本认识。