企业官网建设流程全解析

简介：Timothy Sauer《数值分析》第三版的教师解答手册（Instructor's Solutions Manual），面向使用该教材修读数值计算方法课程的高校学生、考研复习者与授课教师，用于核对课后习题答案、理解算法推导并排查计算偏差。内容按教材章节编排，从多项式求值、二进制与浮点表示、有效数字损失，到二分法、不动点迭代、Newton 法、精度极限，再到高斯消元与 LU 分解、插值误差、Chebyshev 插值与三次样条、最小二乘与 QR 分解、数值微分与 Newton-Cotes、Romberg 与自适应积分，直至常微分方程初值问题、Runge-Kutta 与边值问题的打靶法和有限元，均给出对应习题的求解过程与结果。压缩包共 1 个 PDF 文件，约 6.75MB，单文件即含全章节答案，便于按目录跳转检索。已有 3381 人学习，适合作为教材配套练习与考前复盘的对照材料。

1. 对着答案抄数字没用，真正要抄的是崩在第几位

数值分析第三版的答案手册，我一般不是拿来对数字的。同一道 Horner 嵌套乘法题，实数域上两种展开顺序完全等价，落到双精度里能差出 4.76×10⁻¹²；同一道线性方程组，残差已经小到 10⁻¹⁶，解的误差还在 10⁻⁸ 量级晃。这份配套解答手册覆盖第 0 章到第 13 章，从多项式求值、浮点表示一路铺到最小二乘、常微分方程、FFT、DCT 和 SVD，每一节都给出手算步骤、计算机题的输出和误差分析。它适合三类人：正在跟课、需要核对推导的人；写完算法实现、结果对不上却不知道看哪一行的人；以及准备把书里算法抄成代码、想提前知道哪里会翻车的人。把答案当“标准值”用是最浪费的用法，把它当“误差参照物”用才有价值。

2. 从第 0 章开始对账：Horner 嵌套乘法与浮点舍入的三个实验

第 0 章看着像预备知识，实际上是全书的误差基线。这一章的答案密度最高，也最容易被跳过。跳过它的代价在后面每一章都会反复出现：你以为在对算法答案，其实在核对一个已经被舍入污染过的中间量。

2.1 嵌套乘法：同一道题的两种展开顺序

答案 0.1 节第一题给的是 $P(x)=1+x(1+x(5+x(1+x(6))))$，求 $P(1/3)$。手工从最内层往外推：$6$、$1+\frac13\cdot6=3$、$5+\frac13\cdot3=6$、$1+\frac13\cdot6=3$、$1+\frac13\cdot3=2$，结果是 2。换成直接按幂次展开 $6x^4+x^3+5x^2+x+1$ 再代入，实数域上结果一样，但计算路径长了一倍。

计算机题 0.1 节第 1 题把这件事放大了。系数取 51 个 1，$x=1.00001$，答案手册给的 MATLAB 调用是nest(50,ones(51,1),1.00001)，输出 51.01275208274999，与闭式 $(x^{51}-1)/(x-1)$ 相差 $4.76\times10^{-12}$。用 Python 复现同一件事：

def horner(coeffs, x): """coeffs 按降幂排列，acc 从最高次系数开始累积""" acc = 0.0 for c in coeffs: acc = acc * x + c return acc def power_sum(coeffs, x): """朴素做法：逐项算 x 的幂再乘系数，仅用于对照""" n = len(coeffs) - 1 return sum(c * x ** (n - i) for i, c in enumerate(coeffs)) x = 1.00001 coeffs = [1.0] * 51 v_horner = horner(coeffs, x) v_naive = power_sum(coeffs, x) v_exact = (x**51 - 1.0) / (x - 1.0) print(v_horner, v_naive, v_exact) print(abs(v_horner - v_exact), abs(v_naive - v_exact))

horner的循环体只有一次乘法和一次加法，acc始终是当前部分和，不会出现大数相减；power_sum的每一项都要单独求幂，误差逐项累积。这里的关键参数是x与 1 的距离：当 $x$ 接近 1 时，$(x^{51}-1)/(x-1)$ 的分子分母同时趋近于 0 和 0，闭式本身就处在相消边缘，任何一点舍入都会被放大。答案手册在这一题后面专门标注了差值 $4.76\times10^{-12}$，这个数字不是装饰，它告诉你“嵌套乘法比朴素展开少一次数量级的误差”这句话在具体输入下的量级。

答案 0.1 节第 6 题还给了运算量对比：$P(x)=a_0+x^5(a_5+x^5(a_{10}+x^5a_{15}))$ 这类稀疏结构，嵌套写法需要 6 次乘法和 3 次加法；$x^7$ 版本需要 9 次乘法和 4 次加法。这个数字要在写代码时记住——Horner 省的不只是精度，还有乘法次数，在泰勒展开、多项式拟合的循环里这是实打实的性能差。

2.2 二进制展开与浮点编码：把答案里的位串还原出来

答案 0.2 节给了大量十进制转二进制的中间过程，比如 $(79)_{10}=(1001111)2$、$(227){10}=(11100011)2$、$(31/64){10}=(0.011111)2$、$(0.1){10}=(0.00011)_2$ 后面跟着省略号。这些省略号是整章的重点：有限小数在二进制里未必有限。

def int_to_bin(n): if n == 0: return "0" bits = [] while n: bits.append(str(n % 2)) n //= 2 return "".join(reversed(bits)) def frac_to_bin(f, maxbits=24): """f 在 [0,1) 内，反复乘 2 取整数位""" bits = [] for _ in range(maxbits): f *= 2 if f >= 1: bits.append("1") f -= 1 else: bits.append("0") if f == 0: # 有限二进制表示，提前退出 break return "".join(bits) print(int_to_bin(79)) # 1001111 print(int_to_bin(227)) # 11100011 print(frac_to_bin(0.125)) # 001 print(frac_to_bin(0.1)) # 000110011001100110011001，进入 0011 循环

int_to_bin用的是除 2 取余再逆序，这是最不容易出错的手写方式；frac_to_bin用乘 2 取整，maxbits控制最多输出多少位，0.1那一行会稳定输出0001之后接1001的循环节，永远退不出去。答案 0.2 节第 3 题第 (f) 小问写的是(0.1)10 = (0.00011)2，后面那串省略号对应的就是0011循环，不是排版偷懒。

再往下一层是浮点编码。答案 0.3 节第 1 题把 $\mathrm{fl}(1/3)$ 展开成 $+1.0101\cdots01\times2^{-2}$，第 1 题第 (d) 小问把 $\mathrm{fl}(0.9)$ 因为第 53 位非零而向上舍入。这些描述可以直接用位操作验证：

import struct def ieee754_parts(x): bits = struct.unpack(">Q", struct.pack(">d", x))[0] sign = (bits >> 63) & 0x1 exp = (bits >> 52) & 0x7FF frac = bits & ((1 << 52) - 1) return sign, exp, frac sign, exp, frac = ieee754_parts(1/3) print(sign, exp, exp - 1023, hex(frac)) # 0 1021 -2 0x5555555555555

exp - 1023得到真实指数 $-2$，和答案里的 $2^{-2}$ 对上；frac是0x5555555555555，二进制正好是 52 位的0101…01，对应答案里那串交替位。参数说明：struct.pack(">d", x)按大端打包 8 字节双精度，>Q再按无符号 64 位整数解出来，这样符号位、指数位、尾数位的位置关系和教材里的定义完全一致，不需要记偏移量。

2.3 有效位丢失：两种可以当场复现的翻车

答案 0.4 节的标题就是“Loss of Significance”，这一节最有价值的不是结论，而是它挑的例子。第一个是经典的二次方程求根：$x^2-2\times10^8x+1=0$。求根公式里 $-b-\sqrt{b^2-4ac}$ 这一步，$\sqrt{b^2-4ac}$ 和 $b$ 几乎相等，相减直接把有效位砍掉。

import math a, b, c = 1.0, -2e8, 1.0 disc = math.sqrt(b*b - 4*a*c) x1 = (-b + disc) / (2*a) x2_naive = (-b - disc) / (2*a) # 相消发生在这一步 x2_stable = c / (a * x1) # 用韦达定理 x1*x2 = c/a 绕开相减 print(x1, x2_naive, x2_stable) print(abs(x2_naive - x2_stable) / abs(x2_stable))

x2_naive和x2_stable的相对差会落在 $10^{-8}$ 量级，而双精度的机器精度是 $2.2\times10^{-16}$，中间差了整整八个数量级。参数上唯一的开关是a、b、c的量级比：只要 $b^2 \gg 4ac$，这个坑就成立，和具体系数无关。答案手册在这一节给的修法就是韦达定理，代码里x2_stable = c / (a * x1)一行解决，代价是多一次除法，换来八个数量级的精度。

第二个实验更极端，验证机器精度的吞噬效应：

assert 1.0 + 2.0**-53 == 1.0 # 半个 ulp 被吞掉 assert 1.0 + 2.0**-52 > 1.0 # 一个 ulp 才留下 assert (1.0 + (2**-51 + 2**-52 + 2**-53)) - 1.0 == 2.0**-50

第三行对应答案 0.3 节第 5 题第 (b) 小问：$2^{-51}+2^{-52}+2^{-53}=3.5\epsilon_{\text{mach}}$，加 1 之后落在 $3\epsilon$ 和 $4\epsilon$ 的正中间，按“四舍六入五成双”的规则向偶数尾数靠，取 $4\epsilon$，减回去就是 $2^{-50}$。这三行断言写进测试里，比手算一遍管用，因为它把“舍入规则”这件事钉死成了可执行条款。

这三组实验做完，第 0 章就不再是预备知识，而是一把尺子。后面每一章的答案，凡是出现小数位对不上、迭代提前停住、残差和误差不同步，回头用这把尺子量一下，基本能定位到是舍入、相消还是迭代终止条件的问题。

3. 解方程与线性方程组：把收敛判据写成可执行断言

第 1 章和第 2 章是全书最容易被低估的部分。二分法、不动点迭代、牛顿法看起来都是老面孔，答案手册在这一块真正给的是“迭代到哪一步该停”的判断依据，以及残差和误差什么时候会分家。

3.1 二分法的迭代次数不是一个公式就完事

答案 1.1 节的标准结论是：区间 $[a,b]$ 上二分到容差 $\text{tol}$，需要 $n\ge\log_2\frac{b-a}{\text{tol}}$ 步。这句话在纸面上没问题，落到代码里有两个额外约束：一是tol不能小于区间端点的 ulp，否则中点会塌缩到端点，循环永远不收敛；二是函数在端点同号时，公式再对也没有根。

import math def bisect(f, a, b, tol=1e-12, maxiter=200): fa, fb = f(a), f(b) if fa * fb > 0: raise ValueError("端点同号，区间内不一定有根") n_need = math.ceil(math.log2((b - a) / tol)) for k in range(min(n_need, maxiter)): mid = 0.5 * (a + b) if mid == a or mid == b: # 中点无法再细分，提前退出 return mid, k fm = f(mid) if fm == 0.0: return mid, k if fa * fm < 0: b, fb = mid, fm else: a, fa = mid, fm return 0.5 * (a + b), k

n_need是理论步数，maxiter是兜底；mid == a or mid == b这一行是实战补丁，对应答案 1.3 节“Limits of Accuracy”里说的精度天花板——双精度下区间宽度小于约 $10^{-16}$ 乘以端点量级时，二分法已经无法再提供新信息。用 $f(x)=x^3+x-1$ 在 $[0,1]$ 上跑，tol=1e-12时n_need是 40，实际返回的迭代次数通常也是 40 左右；把tol压到1e-18，理论步数变成 60，但实际会在 53 步左右被mid == a截断，多出来的 7 步是空转。这个差异就是答案手册里没有明写、但做计算机题时一定会遇到的边界。

3.2 不动点迭代与牛顿法：收敛阶怎么现场核

答案 1.2 节用 $g(x)=\cos x$ 做不动点迭代的例子，不动点是 0.7390851332151607。收敛速度由 $|g'(x^)|=|\sin x^|\approx0.6736$ 决定，每一步误差乘 0.67，要把误差从 $10^{-1}$ 压到 $10^{-10}$ 大约需要 $\log(10^{-9})/\log(0.6736)\approx52$ 步。

import math def fixed_point(g, x0, tol=1e-12, maxiter=500): x = x0 for k in range(maxiter): x_new = g(x) if abs(x_new - x) < tol: # 用步长做停止条件 return x_new, k x = x_new return x, maxiter root, steps = fixed_point(math.cos, 0.5) print(root, steps) # 约 0.7390851332，50 步上下

tol控制的是相邻两次迭代的差，不是真实误差；对压缩映射来说两者同阶，所以够用。步数用来和答案手册对照：如果你跑出来是 20 步，说明收敛因子估计错了，或者g写成了 $g(x)=\cos(x)+c$ 之类改变了导数。

牛顿法的验证方式不一样，看的是误差比值的稳定性。$f(x)=x^3+x-1$，$f'(x)=3x^2+1$，根在 0.6823278038 附近。

def newton(f, df, x0, tol=1e-14, maxiter=50): x = x0 hist = [] for k in range(maxiter): fx = f(x) x_new = x - fx / df(x) hist.append(abs(x_new - x)) if abs(x_new - x) < tol: x = x_new break x = x_new return x, hist f = lambda t: t**3 + t - 1 df = lambda t: 3*t*t + 1 root, hist = newton(f, df, 1.0) print(root) print([round(h, 18) for h in hist[:6]])

hist里相邻两项的比值 $e_{k+1}/e_k^2$ 会稳定在一个常数附近（这一题大约 0.3 到 0.5），这就是二次收敛的可执行判据。如果这个比值随着迭代发散，要么导数是手推错的，要么初值离根太远进入了另一个吸引域。答案 1.4 节的习题里专门有几道初值敏感的例子，用这段代码跑一遍就能看到“牛顿法并不是无条件收敛”这句话落在哪一步。

3.3 PA=LU 与残差-误差分离：线性方程组最贵的一课

答案 2.4 节的 PA=LU 分解，验证方式很直接：把置换矩阵、下三角、上三角乘回去，看是否等于原矩阵。SciPy 的lu返回的P满足 $PA=LU$，注意这里P已经在左边，不需要再转置。

import numpy as np from scipy.linalg import lu A = np.array([[0.0, 2.0, 1.0], [1.0, 1.0, 1.0], [2.0, 1.0, 3.0]]) P, L, U = lu(A) print(np.allclose(P @ A, L @ U)) # True print(np.allclose(np.tril(L), L)) # L 是单位下三角

P @ A和L @ U的差值应该落在机器精度量级；np.tril(L)检查L对角线以上是否为零、对角线是否为 1。答案 2.4 节给的手算过程里，行交换的顺序直接决定P长什么样，不同教材的编号习惯不一样，用P @ A == L @ U这个等式卡一遍最省事，不用去猜索引从 0 还是 1 开始。

真正值钱的是答案 2.3 节“Sources of Error”和 2.5 节迭代法之间的关系：残差小，不等于误差小。用 Hilbert 矩阵演示最直观。

import numpy as np from scipy.linalg import hilbert n = 8 A = hilbert(n) x_true = np.ones(n) b = A @ x_true x_hat = np.linalg.solve(A, b) residual = np.linalg.norm(A @ x_hat - b) error = np.linalg.norm(x_hat - x_true) print(f"cond(A) = {np.linalg.cond(A):.3e}") print(f"residual = {residual:.3e}, error = {error:.3e}")

cond(A)在 $n=8$ 时大约是 $1.5\times10^{10}$，残差会落在 $10^{-17}$ 量级（已经是舍入噪声），而误差会落在 $10^{-8}$ 量级。这两行输出的比值大致等于条件数，正好对应答案 2.3 节里的误差上界不等式。参数上唯一要调的是n：$n$ 每加 1，条件数涨一个数量级以上，n=10时误差会直接吃掉全部有效位，解出来的x_hat连符号都不可信。

把这段代码和答案 2.5 节的 Jacobi、Gauss-Seidel 放在一起看，会得到一个很实用的判断：迭代法的收敛速度看谱半径，迭代解的可靠性看条件数，这两件事互相独立。谱半径小于 1 只保证收敛，不保证收敛到的东西有有效位。

4. 插值、最小二乘与积分：误差表的定量读法

第 3 到第 5 章是答案手册里表格最密集的部分。Chebyshev 插值、三次样条、QR 分解、Romberg 积分，每一节的习题答案都是一张误差表。这些表不是让你抄的，是让你对照自己的实现，看误差下降的阶数对不对。

4.1 Chebyshev 节点与三次样条的边界条件

答案 3.3 节的 Chebyshev 节点定义是 $x_k=\cos\frac{(2k+1)\pi}{2n+2}$，$k=0,\dots,n$，落在 $[-1,1]$ 上。它和等距节点的差别，用 Runge 函数 $f(x)=1/(1+25x^2)$ 一测就出来。

等距节点在 $n$ 增大时最大误差反而上升，这就是 Runge 现象；Chebyshev 节点把误差压着往下走。答案 3.3 节的习题答案给的误差值和我这里跑出来的在同一个量级，差异来自节点公式的端点约定（有的书写 $\cos\frac{2k+1}{2n+2}\pi$，有的写 $\cos\frac{2k-1}{2n}\pi$，两者相差一个平移）。

import numpy as np from numpy.polynomial.chebyshev import chebpts1, chebfit, chebval def runge(x): return 1.0 / (1.0 + 25.0 * x * x) for n in (6, 10, 16): x_eq = np.linspace(-1, 1, n + 1) x_ch = chebpts1(n + 1) # n+1 个 Chebyshev 点 xg = np.linspace(-1, 1, 2001) p_eq = np.polyfit(x_eq, runge(x_eq), n) p_ch = chebfit(x_ch, runge(x_ch), n) err_eq = np.max(np.abs(np.polyval(p_eq, xg) - runge(xg))) err_ch = np.max(np.abs(chebval(xg, p_ch) - runge(xg))) print(n, f"{err_eq:.2e}", f"{err_ch:.2e}")

chebpts1(n+1)返回 $n+1$ 个第一类 Chebyshev 点，chebfit和chebval在 Chebyshev 基下拟合和求值，避免把 Chebyshev 点硬塞进单项式基导致的病态。这个细节对应答案 3.3 节后半段关于“基的选择影响条件数”的讨论，很多人抄了节点公式却仍然用np.polyfit，误差还是炸，原因就在基。

三次样条这边，答案 3.4 节的重点是边界条件。自然样条要求两端二阶导为零，clamped 样条要求两端一阶导给定，not-a-knot 要求前两个区间用同一个三次多项式。SciPy 的CubicSpline用bc_type参数切换：

import numpy as np from scipy.interpolate import CubicSpline x = np.array([0.0, 1.0, 2.0, 3.0]) y = np.array([0.0, 1.0, 0.0, 1.0]) cs_nat = CubicSpline(x, y, bc_type="natural") cs_nak = CubicSpline(x, y, bc_type="not-a-knot") print(cs_nat(0.5), cs_nak(0.5))

bc_type="natural"对应答案里的自然样条，"not-a-knot"是 MATLABspline的默认行为。两端给出的数值会不一样，答案手册里的具体小数对应某一种边界条件，抄答案前先确认你用的是哪一种，否则差在小数点后第三位还以为是代码写错了。

4.2 正规方程与 QR：条件数的平方律

答案 4.1 节讲正规方程 $A^TAx=A^Tb$，4.3 节讲 QR 分解。两节的习题答案经常给出同一组数据的不同拟合系数，差别就在条件数。

import numpy as np from scipy.linalg import qr t = np.linspace(0, 1, 25) y = np.sin(2 * np.pi * t) + 0.05 * np.random.default_rng(0).standard_normal(t.size) A = np.vander(t, 11) # 10 次多项式，11 个基函数 print(f"cond(A) = {np.linalg.cond(A):.3e}") print(f"cond(A^T A) = {np.linalg.cond(A.T @ A):.3e}") c_normal = np.linalg.solve(A.T @ A, A.T @ y) Q, R = qr(A, mode="economic") c_qr = np.linalg.solve(R, Q.T @ y) print(np.max(np.abs(c_normal - c_qr)))

cond(A^T A)约等于cond(A)的平方，这就是正规方程在 10 次多项式上开始失准的原因。参数t.size和11是这一节唯二要调的量：基函数个数超过 12 时，正规方程解出来的高次项系数会开始乱跳，QR 还能撑住；超过 16 时两者都会崩，那时候要换 SVD 或者降低次数。答案 4.3 节给出的 QR 结果比 4.1 节的正规方程结果更接近真值，原因就在这个平方律。

4.3 Newton-Cotes、Romberg 与 Gauss 的精度阶验证

答案 5.2 到 5.5 节是积分四连：Newton-Cotes、Romberg、自适应、Gauss。验证精度阶最直接的办法是拿一个已知积分值的函数，把步长减半，看误差下降的倍数。

import numpy as np def romberg(f, a, b, levels=6): R = np.zeros((levels, levels)) for i in range(levels): n = 2 ** i h = (b - a) / n xs = np.linspace(a, b, n + 1) R[i, 0] = h * (0.5 * f(xs[0]) + np.sum(f(xs[1:-1])) + 0.5 * f(xs[-1])) for j in range(1, i + 1): R[i, j] = R[i, j-1] + (R[i, j-1] - R[i-1, j-1]) / (4**j - 1) return R R = romberg(np.exp, 0, 1, 6) print(np.round(R, 12)) print(f"exact = {np.e - 1:.12f}")

R[i, j]的第j列对应 $O(h^{2j+2})$ 的精度：第 0 列是梯形法 $O(h^2)$，第 1 列是 Simpson $O(h^4)$，第 2 列是 Boole 规则 $O(h^6)$。levels=6时最右下角的元素和np.e - 1的差会落在 $10^{-16}$ 量级，也就是机器精度封顶。看 Romberg 表有个技巧：如果某一列不再下降而是开始抖动，说明已经触到舍入噪声，再加层数没有意义，这个现象在答案 5.3 节的计算机题里有对应。

Gauss 积分的验证方式不同，看的是多项式的精确次数。$n$ 点 Gauss-Legendre 能精确积分 $2n-1$ 次多项式：

import numpy as np from numpy.polynomial.legendre import leggauss for n in (2, 3, 4): x, w = leggauss(n) # 在 [-1,1] 上的节点和权重 for deg in (2*n - 2 <p> <a href="https://download.csdn.net/download/sunyongjie1983/10942404" style="color:#ec7500;font-size:14px;"> 本文还有配套的精品资源，点击获取 </a> <img alt="menu-r.4af5f7ec.gif" src="https://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif" style="width:16px;margin-left:4px;vertical-align:text-bottom;cursor:text;"> </p>