Appearance
数值分析（可选）
概念
math/10-calculus.md 里已经用过一堆数值方法——中心差分求导、梯形法、辛普森法、欧拉法。 但那里只说了"误差阶是多少"，没说"误差从哪来"。
这一篇补的就是那一句：
它面对三重限制，缺一不可地解释"为什么不能天真地照公式写代码"：
| 限制 | 后果 |
|---|---|
| ① 字长有限（双精度 53 位有效位） | 舍入误差；"两个几乎相等的数相减"会吃掉全部有效位 |
| ② 速度有限（不能取无限小步长、不能做无限项求和） | 截断误差；必须"取到第几项/多大步长" |
| ③ 问题本身可能是病态的 | 再好的算法也救不回来——输入微扰会让输出剧烈变化 |
这三条对应的三个核心概念：
| 概念 | 回答的问题 |
|---|---|
| 误差 | "我这个结果离真值有多远，误差的可信位数是多少？" |
| 收敛阶 | "多迭代一步，误差缩小多少倍？" |
| 条件数与稳定性 | "是问题本身难算，还是我的算法不好？" |
主线落点：
| 落点 | 数值分析在做什么 |
|---|---|
arch 浮点运算 | eps、舍入模式、误差累积；"为什么浮点加法不满足结合律" |
math/10 的全部数值方法 | 本篇给出它们的误差来源与最优步长 |
soft 数值计算/机器学习 | 梯度下降是迭代法；矩阵求逆的稳定性 |
| 性能测量的可信度 | "这个数到底有几位可信" |
⚠️ 本篇是可选篇——主线不硬依赖它，但读完它，
math/10那几段"为什么步长不能太小""为什么辛普森更准"就不再是结论，而是可推导的东西。
原理
一、两个误差来源：截断与舍入
任何数值方法的误差都由两部分组成：
| 误差 | 与步长 h 的关系 | 减小 h 会怎样 |
|---|---|---|
| 截断误差 | O(h^p)，随 h 减小而减小 | 变小 ✔ |
| 舍入误差 | 约 eps/h，随 h 减小而增大 | 变大 ✘ |
两者一高一低，于是总误差有个最优步长**——这就是 math/10-calculus.md 例 2 里"最优 h ≈ 1e-5"的完整解释。
机器 epsilon（双精度）：
含义：1 + ε 与 1 在双精度下就已经分不开了（当 ε < 1.11e-16，即半个 ulp 时）。它决定了"最多能有约 15~16 位十进制有效数字"。
有效位数的粗略规则：双精度给约 16 位十进制有效位；每做一次可能损失若干位，做 n 次运算可能损失约 √n 位（随机游走模型）。
★ 灾难性抵消（catastrophic cancellation）：
两个非常接近的数相减，结果的绝对误差不变，但相对误差爆炸——因为有效位被"抵消掉了"。
经典例子：1 − cos x（x 很小时）
| 写法 | 问题 |
|---|---|
直接 1 - cos(x) | cos x → 1，两个接近 1 的数相减 → 有效位全丢；x < 1e-8 时直接得 0 |
改用 2·sin²(x/2) | sin(x/2) 本身是小量，平方后仍是精确的小量 → 结果准确 |
同样的坑还有：√(x+1) − √x（应改写成 1/(√(x+1)+√x)）、e^x − 1（应用 expm1）、ln(1+x)（应用 log1p）。
一条通用原则："两个大数相减得小数"就是危险信号，试着把式子改写成"直接算那个小数"的形式。 例 1、例 5 都会把这件事量化。
二、条件数与算法稳定性：是问题难，还是算法差
条件数衡量"问题本身对输入扰动的敏感度"（函数 f 在 x 处的相对条件数）：
| 条件数 | 含义 |
|---|---|
≈ 1 | 良态——输入有多少位精度，输出就有多少位 |
≫ 1 | 病态——输入有 k 位有效位，输出可能只剩 k − log₁₀ cond 位 |
关键区分（务必记住）：
| 概念 | 说的是谁 |
|---|---|
| 条件数 | 问题本身的性质（与算法无关，选不出更好的） |
| 稳定性 | 算法对舍入误差的放大程度（这个可以选） |
"病态问题"与"不稳定算法"是两件事：病态问题用稳定算法也只能保住"输入给的那么多位"；而稳定算法在良态问题上能把精度保到几平
eps的量级。 反过来，一个不稳定的算法会让良态问题也变成灾难——例 5 的二次方程就是一个精确控制的例子。
三、方程求根：二分法与牛顿法
问题：求 f(x) = 0 的根。
| 方法 | 思路 | 收敛阶 | 优缺点 |
|---|---|---|---|
| 二分法 | 每步把区间折半，保留有根的那半 | 线性（每步误差 ×1/2，即每步多得 1 个二进制位） | 只要有符号变化就必然收敛；但慢 |
| 牛顿法 | x_{n+1} = x_n − f(x_n)/f'(x_n)（切线法） | 二次（误差平方） | 快得惊人；但可能不收敛（初值不好、导数为 0） |
| 不动点迭代 | x = g(x)，反复代入 | 线性，收敛条件是 ∣g'∣ < 1 | 最简单，判据也最明确 |
牛顿法的二次收敛是什么意思：
"上一步有 1 位小数错，下一步就只剩 2 位；再下一步 4 位、8 位、16 位"——实际效果是"有效位数每步翻倍"。 例 2 会把它算出来。
准确起见，收敛阶的定义：若 lim |e_{n+1}| / |e_n|^p = C ≠ 0，则称 p 阶收敛。 p = 1 线性、p = 2 二次。
四、插值与拟合
插值：找一个多项式穿过给定的 n+1 个点。 n+1 个点唯一确定一个 n 次多项式（拉格朗日/牛顿差商形式）。
⚠️ 龙格现象：在等距节点上，高次插值会在区间两端剧烈振荡，节点越密、振荡越厉害（反直觉！）。
| 节点选择 | 现象 |
|---|---|
| 等距节点 | 端点附近振荡发散（龙格现象） |
切比雪夫节点（x_k = cos((2k+1)π/(2n+2)) 投到区间上） | 振荡被压到最小（因为节点向两端聚集） |
工程结论：别用高次全局插值。 要细分就"分段低次"——分段线性（折线）、分段三次（样条）——这才是图形学与数据插值的实际做法。
拟合与插值的区别：
| 插值 | 拟合 | |
|---|---|---|
| 目标 | 严格穿过每个点 | 整体误差最小（如最小二乘） |
| 数据 | 点少且被认为准确 | 点多且带噪声 |
| 解 | 解线性方程组 | AᵀA x = Aᵀb（正规方程） |
"最小二乘"的几何意义就是 math/11-linalg.md 里的"把 b 投影到 A 的列空间上"——残差向量与列空间正交，这就是"最小"的来源。
五、数值积分与微分的误差项（回链 math/10）
math/10-calculus.md 给出了"误差阶"，这里给出"误差的系数"：
三个直接结论：
| 结论 | 说明 |
|---|---|
| 辛普森对三次以下多项式精确 | 因为 f^{(4)} = 0——这解释了 math/10 例 3 里"辛普森算 x² 误差是 0" |
| 误差与导数的阶数挂钩 | 被积函数越"弯"，梯形法误差越大 |
h → h/2 时误差按 2^{-p} 缩小 | 梯形 1/4、辛普森 1/16 |
理查德森外推（把"已知阶数"变成"更高阶"）：如果用步长 h 和 h/2 各算一次，可以组合出误差 O(h^{p+1}) 甚至 O(h^{p+2}) 的估计——辛普森法本身就是对梯形法做了一次理查德森外推的产物。
数值微分的最优步长（把第一节的两条误差合起来）：
推导思路：把"截断 h^p"与"舍入 eps/h"相加，对 h 求极值——前向差分 p = 1 得 √eps，中心差分 p = 2 得 eps^{1/3}。 实测值就在 1e-5 附近——与 math/10 例 2 的扫描结果一致（那张表里误差最小的正是 h = 1e-5）。
六、线性方程组：直接法与迭代法
直接法：高斯消元，代价约 n³/3 次乘法（比完整矩阵乘法的 n³ 少 2/3，因为只用消一半）；必须选主元（math/11-linalg.md 已讲）。
迭代法（适合稀疏大矩阵）：
| 方法 | 一句话 |
|---|---|
| 雅可比迭代 | 用上一轮的全部新值算本轮 |
| 高斯-赛德尔迭代 | 本轮就算、立刻用（用更新的值） |
收敛判据（回链 math/11 的特征值）：
"谱半径 = 最大特征值的绝对值"——迭代法是"收敛速度由最大特征值决定"，与 math/12-probability.md 里"马尔可夫链收敛速度由第二大特征值决定"是同一类结论。
对角线占优（|a_ii| > Σ_{j≠i}|a_ij|）是"一定收敛"的充分条件——这也是"选主元"在迭代法里的对应物。
示例
例 1：灾难性抵消——同一个数学式，两种写法的精度差多少
python
import math, unicodedata
def w(s):
return sum(2 if unicodedata.east_asian_width(c) in "WF" else 1 for c in s)
def pad(s, n):
return s + " " * max(0, n - w(s))
print("=== ① 1 − cos x（真值 ≈ x²/2） ===")
print(" " + pad("x", 8) + pad("直接算 1-cos(x)", 24) + pad("改写 2sin²(x/2)", 24)
+ pad("真值 x²/2", 24) + "直接的相对误差")
for k in (1, 4, 6, 8, 9):
x = 10.0 ** (-k)
a = 1 - math.cos(x)
b = 2 * math.sin(x / 2) ** 2
t = x * x / 2
rel = abs(a - t) / t if a != 0 else float("inf")
print(" " + pad("1e-%d" % k, 8) + pad("%.18e" % a, 26) + pad("%.18e" % b, 26)
+ pad("%.18e" % t, 26)
+ ("100%（结果直接是 0）" if a == 0 else "%.3e" % rel))
print(" ★ x = 1e-8 时「直接算」得到 0，而真值是 5e-17 —— 有效位被完全抵消")
print(" ★ 改写后的 2sin²(x/2) 一直保持机器精度：因为它「先算出那个小数，再平方」")
print()
print("=== ② √(x+1) − √x（真值 ≈ 1/(2√x)） ===")
print(" " + pad("x", 14) + pad("直接相减", 24) + pad("改写 1/(√(x+1)+√x)", 24) + "直接的相对误差")
for x in (1e6, 1e8, 1e10, 1e12, 1e14):
a = math.sqrt(x + 1) - math.sqrt(x)
b = 1 / (math.sqrt(x + 1) + math.sqrt(x))
print(" " + pad("%g" % x, 14) + pad("%.18e" % a, 26) + pad("%.18e" % b, 26)
+ "%.3e" % (abs(a - b) / b))
print(" ★ 分子有理化一下，误差从「大数相减」变成「直接算小量」—— 一行改写换来 9 位精度")预期输出：
=== ① 1 − cos x（真值 ≈ x²/2） ===
x 直接算 1-cos(x) 改写 2sin²(x/2) 真值 x²/2 直接的相对误差
1e-1 4.995834721974179438e-03 4.995834721974234081e-03 5.000000000000000971e-03 8.331e-04
1e-4 4.999999969612645145e-09 4.999999995833334195e-09 5.000000000000000105e-09 6.077e-09
1e-6 5.000444502911705058e-13 4.999999999999582876e-13 4.999999999999999899e-13 8.890e-05
1e-8 0.000000000000000000e+00 5.000000000000000512e-17 5.000000000000000512e-17 100%（结果直接是 0）
1e-9 0.000000000000000000e+00 5.000000000000000358e-19 5.000000000000000358e-19 100%（结果直接是 0）
★ x = 1e-8 时「直接算」得到 0，而真值是 5e-17 —— 有效位被完全抵消
★ 改写后的 2sin²(x/2) 一直保持机器精度：因为它「先算出那个小数，再平方」
=== ② √(x+1) − √x（真值 ≈ 1/(2√x)） ===
x 直接相减 改写 1/(√(x+1)+√x) 直接的相对误差
1e+06 4.999998750463419128e-04 4.999998750000625262e-04 9.256e-11
1e+08 5.000000055588316172e-05 4.999999987499999612e-05 1.362e-08
1e+10 4.999994416721165180e-06 4.999999999875000369e-06 1.117e-06
1e+12 5.000038072466850281e-07 4.999999999998749341e-07 7.614e-06
1e+14 5.029141902923583984e-08 4.999999999999987201e-08 5.828e-03
★ 分子有理化一下，误差从「大数相减」变成「直接算小量」—— 一行改写换来 9 位精度例 2：牛顿法求 ∛2，看清"二次收敛"
python
import math, unicodedata
def w(s):
return sum(2 if unicodedata.east_asian_width(c) in "WF" else 1 for c in s)
def pad(s, n):
return s + " " * max(0, n - w(s))
target = 2 ** (1 / 3)
print("=== 牛顿法求 ∛2：解 f(x) = x³ − 2 = 0 ===")
print(" 迭代式：x ← x − (x³ − 2)/(3x²) = (2x + 2/x²)/3")
print(" 初值 x0 = 1.0")
print()
print(" " + pad("迭代 n", 8) + pad("x_n", 22) + pad("误差 |x_n − ∛2|", 18) + pad("误差比 e_n/e_{n-1}²", 20) + "有效位数")
x = 1.0
prev_err = None
for n in range(1, 7):
x = (2 * x + 2 / (x * x)) / 3
err = abs(x - target)
ratio = "—" if prev_err in (None, 0) else "%.4f" % (err / prev_err ** 2)
digits = "%.1f" % (-math.log10(err)) if err > 0 else "15+（机器精度）"
print(" " + pad(str(n), 8) + pad("%.20f" % x, 24) + pad("%.6e" % err, 18)
+ pad(ratio, 20) + digits)
prev_err = err
print()
print(" ∛2 = %.20f（Python 用 2**(1/3) 算的）" % target)
print(" ★ 关键看「误差比 e_n/e_{n-1}²」这一列：它趋于一个常数 → 这就是 p = 2 的定义")
print(" ★ 有效位数每步翻倍：1 → 2 → 4 → 8 → 16 位，4 步就到机器精度")
print()
print("=== 收敛阶的定义（把它写清楚） ===")
print(" 若 lim |e_{n+1}| / |e_n|^p = C ≠ 0，则称方法 p 阶收敛")
print(" p = 1（线性）：误差每步乘 C（二分法 C = 1/2，每步多 1 个二进制位）")
print(" p = 2（二次）：误差每步被平方（牛顿法；有效位数翻倍）")预期输出：
=== 牛顿法求 ∛2：解 f(x) = x³ − 2 = 0 ===
迭代式：x ← x − (x³ − 2)/(3x²) = (2x + 2/x²)/3
初值 x0 = 1.0
迭代 n x_n 误差 |x_n − ∛2| 误差比 e_n/e_{n-1}² 有效位数
1 1.33333333333333325932 7.341228e-02 — 1.1
2 1.26388888888888883955 3.967839e-03 0.7362 2.4
3 1.25993349344997707107 1.244356e-05 0.7904 4.9
4 1.25992105001776977247 1.228966e-10 0.7937 9.9
5 1.25992104989487319067 0.000000e+00 0.0000 15+（机器精度）
6 1.25992104989487319067 0.000000e+00 — 15+（机器精度）
∛2 = 1.25992104989487319067（Python 用 2**(1/3) 算的）
★ 关键看「误差比 e_n/e_{n-1}²」这一列：它趋于一个常数 → 这就是 p = 2 的定义
★ 有效位数每步翻倍：1 → 2 → 4 → 8 → 16 位，4 步就到机器精度
=== 收敛阶的定义（把它写清楚） ===
若 lim |e_{n+1}| / |e_n|^p = C ≠ 0，则称方法 p 阶收敛
p = 1（线性）：误差每步乘 C（二分法 C = 1/2，每步多 1 个二进制位）
p = 2（二次）：误差每步被平方（牛顿法；有效位数翻倍）例 3：插值与龙格现象——等距节点为什么会振荡
python
import math, unicodedata
def w(s):
return sum(2 if unicodedata.east_asian_width(c) in "WF" else 1 for c in s)
def pad(s, n):
return s + " " * max(0, n - w(s))
def lagrange_eval(xs, ys, x):
"""拉格朗日插值：λ_i(x) = Π_{j≠i} (x−x_j)/(x_i−x_j)"""
s = 0.0
n = len(xs)
for i in range(n):
term = ys[i]
for j in range(n):
if j != i:
term *= (x - xs[j]) / (xs[i] - xs[j])
s += term
return s
RUNGE = lambda x: 1 / (1 + 25 * x * x)
N = 11 # 12 个节点 → 11 次多项式
xs_eq = [-1 + 2 * k / (N - 1) for k in range(N)]
xs_ch = [math.cos((2 * k + 1) * math.pi / (2 * N + 2)) for k in range(N - 1, -1, -1)]
GRID = [-1 + 2 * k / 400 for k in range(401)]
print("=== 龙格函数 f(x) = 1/(1+25x²)，插值多项式取 11 次 ===")
print(" " + pad("节点选择", 22) + pad("最大误差", 16) + "最大误差出现在")
for name, xs in (("等距节点（12 个）", xs_eq), ("切比雪夫节点（12 个）", xs_ch)):
ys = [RUNGE(x) for x in xs]
errs = [(abs(lagrange_eval(xs, ys, x) - RUNGE(x)), x) for x in GRID]
mx, at = max(errs)
print(" " + pad(name, 22) + pad("%.6f" % mx, 16) + "x = %+.4f" % at)
ys_eq = [RUNGE(x) for x in xs_eq]
print()
print(" 等距节点插值在若干点上的表现（真值 vs 插值）：")
print(" " + pad("x", 10) + pad("真值 f(x)", 16) + pad("插值 p(x)", 18) + "误差")
for x in (0.0, 0.5, 0.7, 0.9, 0.95, 1.0):
p = lagrange_eval(xs_eq, ys_eq, x)
print(" " + pad("%+.2f" % x, 10) + pad("%.6f" % RUNGE(x), 16) + pad("%.6f" % p, 18)
+ "%.6f" % abs(p - RUNGE(x)))
print(" ★ 表里取的都是「非节点」位置（0.5、0.7、0.9、0.95）；节点上插值必然精确，没有信息量")
print(" ★ 两端振荡越靠近端点越剧烈 —— 这就是龙格现象：等距节点越密，端点振荡反而越大")
print(" ★ 换切比雪夫节点（向两端聚集）能把最大误差压下去一个量级以上")
print(" ★ 工程做法：与其升次数，不如分段低次（折线/样条）")预期输出：
=== 龙格函数 f(x) = 1/(1+25x²)，插值多项式取 11 次 ===
节点选择 最大误差 最大误差出现在
等距节点（12 个） 1.915643 x = +0.9400
切比雪夫节点（12 个） 0.182758 x = +0.0000
等距节点插值在若干点上的表现（真值 vs 插值）：
x 真值 f(x) 插值 p(x) 误差
+0.00 1.000000 1.000000 0.000000
+0.50 0.137931 0.253755 0.115824
+0.70 0.075472 -0.226196 0.301668
+0.90 0.047059 1.578721 1.531662
+0.95 0.042440 1.923631 1.881191
+1.00 0.038462 0.038462 0.000000
★ 表里取的都是「非节点」位置（0.5、0.7、0.9、0.95）；节点上插值必然精确，没有信息量
★ 两端振荡越靠近端点越剧烈 —— 这就是龙格现象：等距节点越密，端点振荡反而越大
★ 换切比雪夫节点（向两端聚集）能把最大误差压下去一个量级以上
★ 工程做法：与其升次数，不如分段低次（折线/样条）例 4：数值积分的误差系数与理查德森外推
python
import math, unicodedata
def w(s):
return sum(2 if unicodedata.east_asian_width(c) in "WF" else 1 for c in s)
def pad(s, n):
return s + " " * max(0, n - w(s))
def trap(f, a, b, n):
h = (b - a) / n
s = 0.5 * (f(a) + f(b))
for i in range(1, n):
s += f(a + i * h)
return h * s
def simp(f, a, b, n):
h = (b - a) / n
s = f(a) + f(b)
for i in range(1, n):
s += (4 if i % 2 else 2) * f(a + i * h)
return h * s / 3
f = math.exp
EXACT = math.e - 1
print("=== ∫_0^1 e^x dx = e − 1 = %.15f ===" % EXACT)
print(" " + pad("n", 8) + pad("梯形法误差", 16) + pad("误差比", 10) + pad("辛普森误差", 16) + "误差比")
pt = ps = None
for n in (4, 8, 16, 32, 64):
et = abs(trap(f, 0, 1, n) - EXACT)
es = abs(simp(f, 0, 1, n) - EXACT)
rt = "—" if pt is None else "%.4f" % (et / pt)
rs = "—" if ps is None else "%.4f" % (es / ps)
print(" " + pad(str(n), 8) + pad("%.6e" % et, 16) + pad(rt, 10)
+ pad("%.6e" % es, 16) + rs)
pt, ps = et, es
print(" ★ 梯形法误差比 → 1/4（O(h²)）；辛普森 → 1/16（O(h⁴)）—— 与误差项公式一致")
print()
print("=== 理查德森外推：用两个梯形值组合出辛普森 ===")
print(" R(h) = (4·T(h/2) − T(h)) / 3")
print(" " + pad("h", 10) + pad("T(h)", 20) + pad("T(h/2)", 20) + pad("R(h)", 20) + pad("S(h/2)", 20) + "R 与 S 之差")
for n in (4, 8, 16, 32):
h = 1 / n
T1 = trap(f, 0, 1, n)
T2 = trap(f, 0, 1, 2 * n)
R = (4 * T2 - T1) / 3
S = simp(f, 0, 1, 2 * n)
print(" " + pad("%g" % h, 10) + pad("%.15f" % T1, 20) + pad("%.15f" % T2, 20)
+ pad("%.15f" % R, 20) + pad("%.15f" % S, 20) + "%.3e" % abs(R - S))
print(" ★ 结论：辛普森法就是「对梯形法做一次理查德森外推」的产物（精度从 h² 提到 h⁴）")
print(" ★ 只要有「已知阶数 + 两个步长的结果」，就能白拿一阶更高的精度")预期输出：
=== ∫_0^1 e^x dx = e − 1 = 1.718281828459045 ===
n 梯形法误差 误差比 辛普森误差 误差比
4 8.940076e-03 — 3.701346e-05 —
8 2.236764e-03 0.2502 2.326241e-06 0.0628
16 5.593001e-04 0.2500 1.455928e-07 0.0626
32 1.398319e-04 0.2500 9.102726e-09 0.0625
64 3.495839e-05 0.2500 5.689702e-10 0.0625
★ 梯形法误差比 → 1/4（O(h²)）；辛普森 → 1/16（O(h⁴)）—— 与误差项公式一致
=== 理查德森外推：用两个梯形值组合出辛普森 ===
R(h) = (4·T(h/2) − T(h)) / 3
h T(h) T(h/2) R(h) S(h/2) R 与 S 之差
0.25 1.727221904557517 1.720518592164302 1.718284154699897 1.718284154699897 2.220e-16
0.125 1.720518592164302 1.718841128579995 1.718281974051892 1.718281974051892 6.661e-16
0.0625 1.718841128579995 1.718421660316327 1.718281837561771 1.718281837561771 4.441e-16
0.03125 1.718421660316327 1.718316786850094 1.718281829028016 1.718281829028015 6.661e-16
★ 结论：辛普森法就是「对梯形法做一次理查德森外推」的产物（精度从 h² 提到 h⁴）
★ 只要有「已知阶数 + 两个步长的结果」，就能白拿一阶更高的精度例 5：二次方程求根的稳定化改写（不稳定的算法会毁掉良态问题）
python
import math, unicodedata
def w(s):
return sum(2 if unicodedata.east_asian_width(c) in "WF" else 1 for c in s)
def pad(s, n):
return s + " " * max(0, n - w(s))
a, b, c = 1.0, 1e8, 1.0 # x² + 1e8·x + 1 = 0
print("=== 解 x² + %.0e·x + 1 = 0 ===" % b)
print(" 真值（解析）：x1 = −1e8，x2 = −1e-8（两数之积 = c/a = 1）")
print()
disc = b * b - 4 * a * c
sq = math.sqrt(disc)
print(" 判别式 b² − 4ac = %.0f，√判别式 = %.8f，b = %.1f" % (disc, sq, b))
print(" 注意：√判别式 与 b 相差约 %.3e —— 两个大数几乎相等" % (b - sq))
print()
print(" 【朴素公式】x = (−b ± √(b²−4ac)) / (2a)")
naive_big = (-b - sq) / (2 * a)
naive_small = (-b + sq) / (2 * a)
print(" x(大) = %.20e" % naive_big)
print(" x(小) = %.20e ← 真值 −1e-8" % naive_small)
rel_naive = abs(naive_small - (-1e-8)) / 1e-8
print(" 小根的相对误差 = %.4f（约 %.1f%%）—— 有效位几乎全丢" % (rel_naive, rel_naive * 100))
print()
print(" 【稳定写法】小根不用减法，改用两根之积：x小 = c / (a·x大)")
stable_big = (-b - sq) / (2 * a)
stable_small = c / (a * stable_big)
print(" x(大) = %.20e" % stable_big)
print(" x(小) = %.20e" % stable_small)
rel_stable = abs(stable_small - (-1e-8)) / 1e-8
print(" 小根的相对误差 = %.3e" % rel_stable)
print()
print(" ★ 同一道题、同一个判别式：只是把小根换成「积除以大根」，误差从 %.2f%% 降到 %.1e"
% (rel_naive * 100, rel_stable))
print(" ★ 这就是「算法稳定性」：问题本身良态，但朴素公式把它算成了灾难")预期输出：
=== 解 x² + 1e+08·x + 1 = 0 ===
真值（解析）：x1 = −1e8，x2 = −1e-8（两数之积 = c/a = 1）
判别式 b² − 4ac = 9999999999999996，√判别式 = 99999999.99999999，b = 100000000.0
注意：√判别式 与 b 相差约 1.490e-08 —— 两个大数几乎相等
【朴素公式】x = (−b ± √(b²−4ac)) / (2a)
x(大) = -1.00000000000000000000e+08
x(小) = -7.45058059692382812500e-09 ← 真值 −1e-8
小根的相对误差 = 0.2549（约 25.5%）—— 有效位几乎全丢
【稳定写法】小根不用减法，改用两根之积：x小 = c / (a·x大)
x(大) = -1.00000000000000000000e+08
x(小) = -1.00000000000000002092e-08
小根的相对误差 = 0.000e+00
★ 同一道题、同一个判别式：只是把小根换成「积除以大根」，误差从 25.49% 降到 0.0e+00
★ 这就是「算法稳定性」：问题本身良态，但朴素公式把它算成了灾难考点
考点
1. 两类误差与最优步长
- 截断误差
O(h^p)：随h减小而减小； - 舍入误差
≈ eps/h：随h减小而增大； - 两者相加取极值 → 最优步长：前向差分
√eps ≈ 1.5e-8；中心差分eps^{1/3} ≈ 6e-6； - 双精度
ε_mach ≈ 2.22e-16（约 15~16 位十进制有效数字）。
2. 灾难性抵消与三种标准改写
| 危险写法 | 稳定改写 |
|---|---|
1 − cos x | 2 sin²(x/2) |
√(x+1) − √x | 1/(√(x+1) + √x) |
e^x − 1 | expm1(x) |
ln(1+x) | log1p(x) |
二次方程小根 (−b+√disc)/(2a) | c/(a·x_大) |
判据一句话：凡是"两个大数相减得小数"，就是要改写的地方。
3. 条件数 vs 稳定性
- 条件数
cond = |x f'(x)/f(x)|：问题本身的敏感度，与算法无关； - 稳定性：算法对舍入误差的放大程度，可以选；
- 病态问题 + 稳定算法 = 只能保住输入带来的位数；良态问题 + 不稳定算法 = 灾难。
4. 方程求根的三条
| 方法 | 收敛阶 | 备注 |
|---|---|---|
| 二分法 | 线性（1/2） | 必收敛；每步多 1 个二进制位 |
| 牛顿法 | 二次 | 有效位数每步翻倍；可能不收敛 |
不动点 x=g(x) | 线性 | 收敛条件 ∣g'∣ < 1 |
- 收敛阶定义：
lim |e_{n+1}|/|e_n|^p = C ≠ 0； - 牛顿法对
∛2从x₀=1起：误差7.34e-2 → 3.97e-3 → 1.24e-5 → 1.23e-10 → 0。
5. 插值
n+1个点唯一确定一个n次多项式；- 龙格现象：等距节点 + 高次 → 两端剧烈振荡；
- 切比雪夫节点（向两端聚集）能压住振荡；
- 工程做法：分段低次（折线、样条），别升次数；
- 插值 vs 拟合：插值必过点，拟合求整体最小（最小二乘 = 投影到列空间）。
6. 数值积分的误差项
- 辛普森对三次以下多项式精确（
f^{(4)} = 0）； - 理查德森外推
R(h) = (4T(h/2) − T(h))/3——这就是辛普森法的来历； - 误差比：梯形
h→h/2时 ×1/4，辛普森 ×1/16。
7. 线性方程组的两个代价与一个判据
- 高斯消元约
n³/3次乘法（比矩阵乘n³少 2/3），必须选主元； - 雅可比 / 高斯-赛德尔迭代适合稀疏矩阵；
- 收敛判据：迭代矩阵谱半径
ρ < 1；对角线占优是充分条件。
8. 本节五个易错点
eps = 2.22e-16是"1+eps能分辨"的量级，不是"绝对误差下限"——实际误差还要乘上数值本身的大小；- 数值求导的
h不是越小越好（太小被舍入吃掉）； - 牛顿法不保证收敛：初值很差、
f'≈0、根附近有拐点都可能失败； - "更高次插值更准"是错的（龙格现象）；
- 条件数大 ≠ 算法差：病态是问题的属性，先分清是哪一个再改。
小结
- 数值分析 = 有限精度怎么算连续数学：三重限制——字长、速度、问题本身的病态。
- 两类误差：截断
O(h^p)（大h时主导）+ 舍入≈eps/h（小h时主导），两者交汇处就是最优步长。 ε_mach ≈ 2.22e-16：决定了双精度约 15~16 位有效数字。- 灾难性抵消：
1−cos x在x=1e-8时直接算出 0（真值5e-17）；改写2sin²(x/2)、1/(√(x+1)+√x)、c/(a·x_大)能拿回全部精度。 - 条件数（问题）与稳定性（算法）要分开：不稳定算法会把良态问题算成灾难（例 5 的小根误差 25%）。
- 求根：二分法线性必收敛；牛顿法二次、有效位翻倍；不动点需
|g'|<1。 - 插值：等距高次会龙格振荡；切比雪夫节点能压住；工程上用分段低次。
- 积分误差项：梯形
h²、辛普森h⁴；辛普森 = 梯形 + 一次理查德森外推。 - 迭代法收敛看谱半径
ρ < 1——与math/11的特征值、math/12的马尔可夫收敛速度同源。
回到主线：这一篇是对 math/10-calculus.md 的"补刀"。
math/10给出的都是"结论"：中心差分O(h²)、辛普森O(h⁴)、欧拉法O(h)、最优h ≈ 1e-5。本篇给出的是"为什么"：误差分成截断与舍入两半，一个随h变小、一个随h变大——最优步长就是这两条曲线的交点；辛普森的高阶来自理查德森外推；而"什么时候不能照公式写"的答案，全在"灾难性抵消"这四个字里。它在主线的另一处落点是
arch：浮点运算的舍入模式、eps的量级、"浮点加法不满足结合律"——那不是硬件 bug，而是本篇第一节的直接后果。一句话记住它的地位："数学给出的是实数域上的答案，计算机给出的是浮点数域上的答案；两者之差就是这门课。"
下一篇是复变函数与积分变换：它是 math 的最后一篇，也是 elec 信号与系统的直接前置——为什么"正弦信号在系统里形状不变"？为什么拉普拉斯变换能把微分方程变成代数方程？答案都在复数的指数形式与极点分布里。
下一篇：复变函数与积分变换（可选）
评论（0）
当前浏览器不允许本地存储，评论无法保存。
还没有评论，来说两句。