四组难度递进:A 概念判断 → B 推导与计算 → C 编程实现 → D 综合与开放;每题答案默认折叠(::: sol),先闭卷写,写不出来再点开。
建议用法:A 组限时 20 分钟、每题必须写出理由(只写「对 / 错」不给分);B 组手写完整步骤,不许写「略」;C 组按验收标准跑通并对照「参考实现要点」;D 组写成半页方案,四件事缺一不可。
分工:本页只管练与测;公式汇总、易混对照与复习路线看复习页(那页是「查」),课程地图与讲次导航看课堂精编导览。
题源限第 1–7 讲与第 9 讲。第 8 讲(10-12)的材料已归档,占位页已建但内容待补,本页凡涉及处一律按复习页的口径处理:以第 9 讲开头的复盘为据(唯一性、压缩条件下的收敛、收敛阶、牛顿法二阶收敛)。
三个约定
- 可复算:B 组每题都给出中间步骤与数值答案,凡是「数值」都可用四则运算复现;文末说明关键数值的复核方式。
- 标注纪律:凡不是这 8 讲课堂结论的延伸(补充篇/原书内容、需要另行推导的推论),一律标「需核实」,不冒充课堂结论。
- 术语:算法名首次出现给英文全称;「上凸 / 下凸」按课堂口径(f′′<0 / f′′>0),与英文 convex/concave 相反,答题时建议直接写 f′′ 的符号。
判定要求:判断 + 一句话理由;判为「错」的要指出错在哪个条件上。题源见每题括号内的小节。
A1. 「误差四分类中,舍入误差是本课程要控制、要研究的核心对象;截断误差由机器字长决定,只能控制其传播。」(第 1 讲 §四)
答案
错——两句恰好对调。 四类误差里:模型 / 观测误差无法消除;截断误差是数值方法本身带来的近似,是本课程能控制、要研究的核心对象;舍入误差由机器字长决定,可控制的只是它的传播。把「截断」当成机器决定的,整条误差红线就搭错了。
A2. 「绝对误差限变小,有效数字一定增加——例如 a=3.1579 的 ε(a) 从 0.0004 减到 0.0002 时,有效位数从 4 位变成 5 位。」(第 2 讲 §一)
答案
错。 位数不看误差限变没变小,只看「误差限是否不超过某一位的半个单位」:ε=0.0004 时区间 [3.1575,3.1583] 落在千分位的半单位区间 [3.1575,3.1585] 内 → 4 位;ε=0.0002 时区间 [3.1577,3.1581] 仍然落在同一半单位区间内 → 还是 4 位。要想多一位,误差限必须跨过万分位的半单位 0.00005。
A3. 「次可乘性 ∥AB∥≤∥A∥∥B∥ 对任何矩阵范数都可以证明,所以不必列进公理。」(第 2 讲 §四·§五)
答案
错。 对从属(算子)范数,次可乘性确实可证(x 取单位向量,用两次相容关系);但对一般矩阵范数它是公理要求(矩阵范数四公理的第 4 条)。把「对某类范数成立」读成「对所有范数成立」是常见失分点——这条区分正是第 2 讲结尾证明乘积不等式的意义所在。
A4. 「∥G∥<1 只保证『存在某个初值』使迭代收敛;要保证对任意初值收敛,还得再验证 ρ(G)<1。」(第 4 讲 §五、第 5 讲 §三)
答案
错——两条判据的强弱被颠倒了。 ∥G∥<1(某个范数)⇒ 解存在唯一且对任意初值收敛(充分,但结论已经足够强);ρ(G)<1 才是充要条件(对任意初值收敛 ⟺),只是大规模问题算不出特征根。正确表述:范数条件「好算但偏严」,谱半径条件「最准但难算」。
A5. 「SOR 里 ω 越大越好:ω→2 松弛越充分、收敛越快,取到 ω=2 是极限情形,仍有 ρ(G)<1。」(第 6 讲 §五)
答案
错。 ω∈(0,2) 只是必要条件;本页 B5 对具体矩阵算出 ω=2 时 ρ(G)=1(特征值 (−5±i11)/6,模恰为 1),充要条件不成立。C2 的实测也显示 ω=1.5 的步数(25)反而多于 ω≈1.02(5)——参数不是越大越好,得看 ρ(G(ω))。
A6. 「幂法改用 ∞-范数单位化后,ykTAyk 仍收敛到 λ1,只是慢一点。」(第 7 讲 §五)
答案
错。 ∞-范数下 ∥yk∥2=1,ykTyk=1,公式必须换成商 ykTAyk/ykTyk(板书在「∥⋅∥∞」那一行特意提示「另作调整」)。数值见 B6(3):∞-范数下 yTAy=5.0、yTy=1.381966,商才等于 3.618034=λ1。所以不是「慢一点」,而是公式本身要换。
A7. 「位移反幂法里,p 只要『接近』目标特征值 λk 就行,不必是最近的那个;A−pI 是否奇异也无所谓。」(第 7 讲 §八)
答案
错。 前提是 ∣λk−p∣ 严格小于其余所有 ∣λi−p∣——只有这时 A−pI 的目标特征值 λk−p 才是模最小的,反幂法才瞄得准;否则求到的是别的特征值。第二个说法也错:反幂法要求每步解方程组,A−pI 必须可逆。(p 恰等于某个特征值时 A−pI 奇异,这一点课堂未讨论,属需核实的延伸。)
A8. 「牛顿法的二阶收敛是四条全局条件保证的:没有 (i)–(iv) 就只能得到一阶收敛。」(第 9 讲 §二·§三·§四)
答案
错——把「阶」与「范围」混成了一件事。 二阶收敛来自 φ′(s)=0(φ=x−f/f′,φ′=ff′′/(f′)2,在 f(s)=0、f′(s)=0 处为零),这是收敛阶定理给的,与四条件无关;四条条件管的是收敛范围——把「初值充分靠近」升级成「初值落在 f(x0)f′′(x0)>0 指定的阴影区」。所以:没有四条件,牛顿法在合适的初值下仍是二阶,只是你不知道多近才算「合适」。
A9. 「单点割线法把 xk−1 钉死成 x0 后,收敛阶从 1.618 降到 1,代价是每步仍要算两次函数值。」(第 9 讲 §八)
答案
错——后半句错。 阶降为一阶是对的(实测误差比恒为 0.424,见 B7(5));但 f(x0) 只算一次并全程复用,每步只新增一次求值。它的处境正是「省一次求值,赔一档阶」:牛顿 4 步 / 割线 8 步 / 单点割线 31 步(阈值 10−12,见 C4)。
A10. 「插值多项式的存在唯一性依赖基的选取:只有换成拉格朗日基之后,存在唯一性才成立。」(第 9 讲 §十一·§十二)
答案
错。 存在唯一由范德蒙德(Vandermonde)行列式 detA=∏i<j(xj−xi)=0 保证,与基的选取无关(同一个多项式空间的不同基而已)。基决定的是系数好不好求:{1,x,x2,…} 的系数矩阵病态(C5 实测 n=12 时 cond2=6.8×109),拉格朗日基在节点上的取值矩阵是单位阵,系数直接读出。这正对应复习页的口诀「存在唯一 = 能算出来」。
要求:写出每一步中间量;数值答案保留 4 位有效数字(有精确值时给精确值)。全部可手算,必要时用计算器核第 4 位。
B1. 有效数字、相对误差限与四则误差限(第 1 讲 §六、第 2 讲 §一·§二)
(1) a=3.1579:分别取 ε(a)=0.0004 与 0.0008,判断有效位数、给出规范写法并算相对误差限;
(2) A=2.14、B=3.15 都是「末位取整」的近似值,求 A+B 的误差限与有效位数;
(3) 取 a=2.00(εa=0.005)、b=0.100(εb=0.0005),算 ε(a/b) 与 εr(a/b),并与「相对误差限相加」对照;再把 b 缩到 0.010(误差限不变)重算一次;
(4) 递推 In=n1−5In−1:判断正推是否稳定,写出稳定方向,并估计 10 步后初始误差被放大多少倍。
答案
(1) ε=0.0004:区间 [3.1575,3.1583]⊂[3.1575,3.1585](千分位半单位)→ 4 位,写作 3.158,εr=0.0004/3.1579=1.27×10−4。ε=0.0008:区间 [3.1571,3.1587] 越出该半单位区间,退到百分位 [3.155,3.165] 才包住 → 3 位,写作 3.16,εr=0.0008/3.1579=2.53×10−4(误差限翻倍,相对误差限也翻倍,但位数只退一位)。
(2) εA+B=0.005+0.005=0.01 超过了百分位的半单位 0.005,而十分位的半单位 0.05≥0.01 → A+B=5.29≈5.3,只剩 2 位有效数字(相加会把位数「吃掉」,这是误差限公式的直接后果)。
(3) ε(a/b)=b2∣b∣εa+∣a∣εb=0.010.1×0.005+2×0.0005=0.010.00075=0.15,a/b=20,εr=0.15/20=0.75%;而 εr(a)+εr(b)=0.25%+0.5%=0.75% ✓ 两者一致(乘除法的相对误差限近似相加)。b=0.010 时:ε(a/b)=10−40.01×0.005+2×0.0005=10.5,a/b=200,εr=5.25%——从 0.75% 涨到 5.25%,全部来自 εr(b)=5% 与分母里的 b2。小数不放分母的定量形态。
(4) 误差传播 εn=−5εn−1,放大因子 ∣−5∣=5>1 → 正向不稳定;反解 In−1=51/n−In,因子 −1/5,∣−1/5∣<1 → 反向稳定(从较大的 n 出发往回推)。10 步放大 510=9765625≈9.8×106 倍:初始 10−6 的误差会被放大到量级 10。判据只看放大因子,与递推式长什么样无关。
B2. 条件数估计与病态诊断(第 3 讲 §一、第 9 讲 §十一)
设 A=(1111.001),b=(2, 2.001)T。
(1) 手算 A−1、∥A∥∞、∥A−1∥∞ 与 cond∞(A)=∥A∥∞∥A−1∥∞;
(2) 把 b2 改成 2.002(右端相对扰动约 0.05%)后重解,比较解的相对变化,并与扰动界 cond×δ 对照;
(3) 问:改用列主元高斯消去法能不能把条件数降下来?为什么?
答案
(1) detA=1×1.001−1×1=0.001,
A−1=0.0011(1.001−1−11)=(1001−1000−10001000).
∥A∥∞=max(1+1, 1+1.001)=2.001;∥A−1∥∞=max(1001+1000, 1000+1000)=2001;
cond∞(A)=2.001×2001=4004.001≈4.0×103.
(与复习页 §2.3登记的量级一致。)log104004=3.60,数量级上对应 3~4 位有效数字的损失——注意「损失位数 ≈log10cond」这条结论属补充篇 2.4,课堂未讲,需核实。
(2) 原方程两式相减得 0.001x2=0.001,故原解 x=(1,1)T(代入可验)。新右端 (2, 2.002)T:相减得 0.001x2=0.002⇒x2=2,再回代 x1=0,新解 x′=(0,2)T。相对扰动 ∥Δb∥∞/∥b∥∞=0.001/2.001=5.0×10−4(0.05%);解的相对变化 ∥Δx∥∞/∥x∥∞=1/1=100%——放大 2×103 倍,与 cond∞ 同量级。而界给出 cond×δ=4004×5.0×10−4≈2.0(即 200%):实测 100%≤200% ✓,界成立且几乎被取到。
(3) 不能。 条件数是问题(矩阵与右端)的属性:cond=∥A∥∥A−1∥ 与用什么算法无关。列主元只把消去乘数控制到 ∣mik∣≤1,防止误差在消去过程中被继续放大;它救不了「两行几乎相同」这个几何事实——第二行只比第一行多在第二分量上加了一个 0.001,两个方程几乎重合,交点对右端极敏感。处置要换提法(换基、换模型),见 D1 与第 9 讲 §十一。
B3. 消去一步、主元与 LU 分解(第 3 讲 §四·§五、第 4 讲 §三)
设 A=0.512121143,b=(3,7,6)T。
(1) 算顺序主子式 Δ1,Δ2,判断不换行的高斯消去能走几步、卡在哪一步;detA=3=0 说明了什么?
(2) 不换行做第一步消去:写出 m21,m31 与消去后的第二、三行;再改用列主元做第一步,比较乘数;
(3) 用列主元把方程组解到底(含回代),给出解;
(4) 对 A′=24−2274375 做 Doolittle 分解(逐圈),写出 L 与 U。
答案
(1) Δ1=0.5=0,Δ2=0.5×2−1×1=0。由充要条件(前 n−1=2 个顺序主子式均不为零 ⟺ 能不换行消成上三角)知:第一步能走,第二步卡住(实测第一步后第二行为 (0,0,2),主元 a22(1)=0)。detA=3=0 说明 Ax=b 有唯一解——「能不能消」看前 n−1 个顺序主子式,「有没有唯一解」看第 n 个,这是两件事。
(2) m21=−0.51=−2,m31=−0.52=−4(∣m∣>1,误差会被放大);第二行 (1,2,4)−2(0.5,1,1)=(0,0,2),第三行 (2,1,3)−4(0.5,1,1)=(0,−3,−1)。列主元:第一列绝对值最大元为 2(第 3 行),交换第 1、3 行得 A=210.5121341、b=(6,7,3)T,此时 m21=−21、m31=−41,都满足 ∣m∣≤1——这就是列主元的全部目的。
(3) 第一步消去后 20011.50.7532.50.25,b=(6,4,1.5)T;第二步主元 1.5>0.75 不必交换,m32=−1.50.75=−21,得 20011.5032.5−1,b=(6,4,−0.5)T。回代:x3=0.5;x2=1.54−2.5×0.5=611=1.8333;x1=26−1×611−3×0.5=34=1.3333。即 x=(34,611,21)T,代回三行逐行验证 ✓。
(4) 逐圈(第 1 行 U → 第 1 列 L → 第 2 行 U → 第 2 列 L):u11=2, u12=2, u13=3;l21=4/2=2, l31=−2/2=−1;u22=7−2×2=3, u23=7−2×3=1;l32=34−(−1)×2=2。回代得最后一元 u33=5−(−1)×3−2×1=6,于是
L=12−1012001,U=200230316,
乘回去逐元素等于 A′ ✓(顺序主子式 Δ1=2,Δ2=6,Δ3=36 均非零,分解存在)。关键在安排顺序:每一圈只含一个未知数。
B4. 迭代矩阵、谱半径与收敛判定(第 5 讲 §四·§五)
设 A=(3−1−14),b=(1,1)T,三分 A=L+D+U。
(1) 写出 D、L、U 与 GJ=D−1(−(L+U)),算 ∥GJ∥∞ 与 ρ(GJ)(精确值 + 4 位小数);
(2) 用行对角占优判断 Jacobi 收敛性,并说明这一步用的是哪条判据;
(3) 写出 GGS=−(L+D)−1U 并算 ρ(GGS);
(4) 后验估计的常数 1−∥GJ∥∞∥GJ∥∞ 是多少?它有什么用?
答案
(1) D=diag(3,4),L=(0−100),U=(00−10),
GJ=D−1(−(L+U))=(041310).
∥GJ∥∞=max(31,41)=31=0.3333。形如 (0ba0) 的矩阵特征值为 ±ab,故 ρ(GJ)=121=231=0.2887。
(2) 行对角占优(strict diagonal dominance):∣3∣>∣−1∣ ✓、∣4∣>∣−1∣ ✓,于是
∥GJ∥∞=imax∣aii∣∑j=i∣aij∣=max(31,41)=31<1 ⇒ Jacobi 收敛.
这一步走的是充分条件(∥G∥∞<1,手算即可);(1) 的 ρ(GJ)<1 是充要条件,本例太小、两者都能算,大规模问题只能走前者。
(3) L+D=(3−104),(L+D)−1=(31121041),
GGS=−(L+D)−1U=(0031121) (下三角) ⇒ ρ(GGS)=121=0.0833.
(数值上 0.0833=0.28872;这种「平方关系」需要额外的一致次序(consistently ordered)条件,本 8 讲未讲,需核实。)
(4) 常数 =2/31/3=21。用途:后验误差估计 ∥x∗−x(k)∥∞≤21∥x(k)−x(k−1)∥∞——这就是「相邻差小就停机」的定量依据。C2 会用实测比值(0.36∼0.45)验证这条界成立。
B5. SOR:构造、ω=2 与最优参数(第 6 讲 §二·§三·§五)
仍用 A=(3−1−14)(n=2)。
(1) 写出 N=L+ω−1D、P=(ω−1−1)D−U,验证 N−P=A;
(2) 取 ω=1.2,用 ∣detG∣=∣ω−1∣n 算出 ∣detG∣,并用 detG=detP/detN 直接核对;
(3) 取 ω=2,写出 G=N−1P 并算出 ρ(G);为什么这说明 ω=2 不能用?
(4) 在 (0,2) 上扫描 ρ(G(ω)),找出使 ρ 最小的 ω;并说明「必要不充分」在这里意味着什么。
答案
(1) N=L+ω−1D=(3/ω−104/ω)(下三角),P=(ω−1−1)D−U=(3(ω−1−1)014(ω−1−1))(上三角)。相减:对角元 3/ω−3(ω−1−1)=3、4/ω−4(ω−1−1)=4,下三角 −1,上三角 −1 → N−P=A ✓(构造的检验式)。
(2) ω=1.2:∣detG∣=∣ω−1∣2=0.22=0.04。直接核对:detN=(3/ω)(4/ω)=12/1.44=8.3333,∣detP∣=ω−1−12×3×4=0.166672×12=0.33333,故 ∣detG∣=0.33333/8.3333=0.04 ✓。
(3) ω=2:N=(1.5−102),P=(−1.501−2),N−1=(3231021),
G=N−1P=(−1−2132−32),trG=−35,detG=1.
特征方程 λ2+35λ+1=0,根 λ=6−5±i11,∣λ∣2=3625+11=1 → ρ(G)=1。充要条件是 ρ(G)<1,而 ∣detG∣=∣ω−1∣n=1 恰好落在必要条件的边界上:特征值模之积为 1,就不可能「每个模都小于 1」。所以 ω=2 不保证收敛,本题更是 ρ=1(误差按模 1 的复特征值旋转、不衰减)。
(4) 以步长 10−3 扫描得 ω≈1.022 处最小:ρmin≈0.0218(ω=1.02 时 0.0354、ω=1.05 时 0.0500——极小点很窄)。经典教材的最优参数公式 ωopt=1+1−ρ(GJ)22=1.0217 与扫描值吻合,但此公式不在本 8 讲范围,需核实。「必要不充分」的意思:定理只保证 ω∈/(0,2) 一定不收敛;区间内部收不收敛要看矩阵——本题 A 对称正定,由 Ostrowski 定理 0<ω<2 是充要(区间内都收敛);换成一般矩阵,同样 ω∈(0,2) 也可能 ρ(G)≥1。
B6. 幂法一步、Rayleigh 商与反幂法(第 7 讲 §四·§五·§七)
设 A=(2113),特征值 λ1=25+5=3.6180339887、λ2=25−5=1.3819660113,取 u0=(1,0)T。
(1) 按 yk=∥uk∥2uk、uk+1=Ayk 迭代 4 步,列出 uk(可整型化)与 yk;
(2) 每步算 Rayleigh 商 Rk=ykTAyk=ukTukukTAuk,列出 ek=∣Rk−λ1∣ 与比值 ek+1/ek,说明它趋于什么;
(3) 若改用 ∞-范数单位化,yTAy 等于多少?该怎么读特征值?
(4) 用反幂法求 λ2(模最小):写出 A−1,按「每步解 Auk+1=yk」迭代 3 步,给出 Rk 序列。
答案
(1) uk 同乘正数不改变 yk,故可整型化:u0=(1,0)→u1=(2,1)→u2=(5,5)→u3=(15,20)→u4=(50,75)。对应
y1=5(2,1)≈(0.8944,0.4472),y2=2(1,1)≈(0.7071,0.7071),y3=(0.6,0.8),y4=13(2,3)≈(0.5547,0.8321).
(顺序不能反:先单位化,再乘 A。)
(2) R0=2(y0=(1,0))、R1=3、R2=3.5、R3=3.6、R4=1347=3.6153846;误差 e=0.618034, 0.118034, 0.018034, 0.002650,…,比值 0.191→0.153→0.147→0.146,趋于 0.1459=(λ1λ2)2=0.3819662。注意:特征向量的误差按 ∣λ2/λ1∣=0.382 衰减,而由 R 读出的特征值误差按它的平方 0.146 衰减——这条「平方加速」是本页的数值观察(课堂未展开,需核实)。
(3) ∞-范数极限单位向量为 y=(0.618034, 1),此时 yTAy=5.0、yTy=1.381966,两者不是同一个数;但
yTyyTAy=1.3819665.0=3.618034=λ1 ✓.
所以换范数后必须补上分母 ykTyk。
(4) A−1=51(3−1−12);迭代 uk+1=A−1yk(实际每步解 Auk+1=yk,不构造 A−1),整型化序列 (1,0)→(3,−1)→(2,−1)→(7,−4):
R0=2,R1=1.5,R2=1.4,R3=1318=1.384615,
继续两步得 1.382353、1.382022 → 报 λ2=1.3820,与 25−5=1.3819660 一致到 4 位以上小数(误差比同样趋于 0.146)。
B7. 二分法、不动点迭代、牛顿法与割线法族(第 7 讲 §九·§十、第 9 讲 §四·§七·§八)
统一取 f(x)=x3−x−1(根 s=1.3247179572447460)。
(1) 验证 f(1)f(2)<0,做三步二分(写出中点与保留区间),并算「区间长 <10−6」所需步数;
(2) 改写成 x=(x+1)1/3:在 [1,1.5] 上验证自映射与压缩条件(给出 L),迭代 4 步并算误差比;
(3) 另一变形 x=x+f(x)=x3−1:算 φ′(s) 并判断能否用;
(4) 牛顿法从 x0=1.5 出发迭代 4 步,用实测阶公式 p≈ln(ek/ek−1)ln(ek+1/ek) 算两个 p,并与理论因子 2f′(s)f′′(s) 对照;
(5) 割线法(x0=1,x1=2)与单点割线法(x0=2 钉住、x1=1.5)各迭代到误差 <10−12,列出步数与实测阶,并解释两者与牛顿法的差距。
答案
(1) f(1)=−1<0,f(2)=5>0 ✓(介值定理保证有根)。c1=1.5、f=0.875>0→[1,1.5];c2=1.25、f=−0.296875→[1.25,1.5];c3=1.375、f=0.224609→[1.25,1.375]。区间长 =2n1,需 2n1<10−6⇒n>log2106=19.93,即 20 步。
(2) φ(x)=(x+1)1/3 单增:φ(1)=1.2599∈[1,1.5],φ(1.5)=1.3572∈[1,1.5] → 自映射 ✓;φ′(x)=31(x+1)−2/3,φ′(1)=0.2100、φ′(1.5)=0.1810,取 L=0.21<1 → 压缩 ✓(两条件各管什么:第一条管迭代值不越界,第二条管误差收缩)。迭代:1→1.2599210→1.3122938→1.3223538→1.3242687;误差 6.48×10−2,1.242×10−2,2.364×10−3,4.492×10−4,比值 0.1917,0.1903,0.1900→φ′(s)=0.18995:线性(一阶)收敛,误差每步约乘 0.19。
(3) φ2′(x)=3x2,φ2′(s)=3×1.7548=5.2646>1 → 在根附近误差被放大,迭代发散(实测 1.5→2.375→12.40→1904→6.9×109)。变形不唯一,选错就发散——这正是收敛定理存在的理由。
(4) 牛顿迭代(xk+1=xk−f/f′,f′=3x2−1):
| k |
xk |
ek=∣xk−s∣ |
| 0 |
1.5 |
1.7528×10−1 |
| 1 |
1.3478260870 |
2.3108×10−2 |
| 2 |
1.3252003990 |
4.8244×10−4 |
| 3 |
1.3247181740 |
2.1675×10−7 |
| 4 |
1.3247179572 |
4.4×10−14 |
实测阶:k=2 处 p=1.9922,k=3 处 p=2.0000 → 二阶 ✓。理论因子 2f′(s)f′′(s)=2(3s2−1)6s=8.52937.9483=0.9319,而实测 e3/e22=0.931 ✓(二阶收敛的常数正是这个因子)。
(5) 割线法(差商代导数)序列:1,2,1.1666667,1.2531120,1.3372064,1.3238501,1.3247079,1.3247179654,1.3247179572,误差 3.25×10−1,6.75×10−1,1.58×10−1,7.16×10−2,1.25×10−2,8.7×10−4,1.0×10−5,8.1×10−9,7.6×10−14——到 10−12 共 8 步;稳定段实测阶 1.5270,1.6731,1.5958,1.6267≈1.618(黄金比,r2−r−1=0)。注意前两步误差反而增大(0.32→0.68):超线性只在根附近显现。
单点割线法(x0=2 钉住)序列:1.5,1.3939394,1.3532486,1.3366746,1.3297631,1.3268528,…,误差比恒为 0.424(常数!)→ 一阶;到 10−12 需 31 步。
差距的来源一句话:牛顿每步用 f 与 f′(二阶)、割线用两次函数值近似导数(阶 1.618)、单点割线把远端钉死(信息陈旧 → 一阶)。步数 4 / 8 / 31,而二分法要 40 步(同在 10−12)。
B8. 拉格朗日插值与误差界(第 9 讲 §十·§十二)
(1) 求过 (0,1),(1,2),(2,5) 的插值多项式:写出三个基函数 li(x) 与 L2(x),化简后对照 li 的筛选性质;
(2) 同一组节点上取 f(x)=ex(系数即 (1,e,e2)),算 x=0.5 与 x=1.5 的插值值与真实误差;
(3) 用余项式 R2(x)=f(x)−L2(x)=3!f(3)(ξ)∏i=02(x−xi) 给出这两点的误差上界,并给出 [0,2] 上的整体上界(提示:求 max∣x(x−1)(x−2)∣);
(4) 这个界「紧」吗?
答案
(1) l0(x)=(0−1)(0−2)(x−1)(x−2)=2(x−1)(x−2),l1(x)=(1)(1−2)x(x−2)=−x(x−2),l2(x)=(2)(2−1)x(x−1)=2x(x−1)。于是
L2(x)=1⋅l0+2⋅l1+5⋅l2=2x2−3x+2−2x2+4x+25x2−5x=x2+1.
筛选性质:li(xi)=1、li(xj)=0 (j=i)——代入节点时只有自己那项存活,系数直接就是函数值 (1,2,5),L2(xi)=xi2+1 与数据一致 ✓(一步方程都不用解,这就是换基的收益)。
(2) 同一组基、系数换成 (1,e,e2):(l0,l1,l2)∣x=0.5=(0.375,0.75,−0.125),
L2(0.5)=0.375+0.75e−0.125e2=1.49007936,e0.5=1.64872127 ⇒ R2=−0.15864191.
x=1.5:L2(1.5)=4.68460741,e1.5=4.48168907⇒R2=+0.20291834。
(3) f(3)=ex 在 [0,2] 上最大为 e2=7.389056,故 3!f(3)(ξ)≤67.389056=1.231509;而 x=0.5 与 1.5 处 ∣∏(x−xi)∣=0.5×0.5×1.5=0.375,于是
∣R2(0.5)∣, ∣R2(1.5)∣≤1.231509×0.375=0.461816.
整体上界:[0,2] 上 max∣x(x−1)(x−2)∣=0.384900(在 x=1±31=0.4226, 1.5774 取到),故 ∣R2(x)∣≤1.231509×0.384900=0.474008。
(4) 不紧。 实测误差只占界的一部分:x=0.5 处 34%、x=1.5 处 44%;[0,2] 上最大实测误差 0.2135 也只到全局界的 45%。余项界给的是有保证的上界,不是误差估计——要预测误差量级、判断「节点加密是否收敛」,得等下一讲(Runge 现象一类问题,需核实)。余项式本身出现在第 9 讲 §十二 页末预告,本页只做误差界。
依赖仅限 numpy。每题给可运行的验收标准——不通过就说明对应环节没接上;「参考实现要点」默认折叠(点开看),只给关键行。
C1. Doolittle 分解与三对角追赶法(Thomas algorithm)(第 3 讲 §四、第 4 讲 §三·§四)
实现 doolittle(A) -> (L, U) 与 thomas(lo, diag, up, d) -> x(前者前代 + 回代解 Ax=b)。
验收标准:① 对 A′=24−2274375:max|L@U - A| < 1e-12,且 np.diag(L) == 1 全部成立(B3(4) 的答案用来对表);② 用前代 + 回代解 A′x=(1,1,1)T 与 np.linalg.solve 的解差 <10−10;③ 三对角 T=tridiag(−1,2,−1)(n=3):L、U 的非零元数都等于 2n−1=5(无填充),并断言 np.count_nonzero(np.triu(L, 2)) == 0 且 np.count_nonzero(np.tril(U, -2)) == 0;④ 追赶法解 2−10−12−10−12x=(1,0,1)T 得 (1,1,1)T,误差 <10−12;⑤ 打印 n=1000 时追赶法(约 3n)与稠密 LU(约 2n3/3)的运算量估计。
参考实现要点
逐圈顺序不能乱——先 U[k, k:] = A[k, k:] - L[k, :k] @ U[:k, k:],再 L[k+1:, k] = (A[k+1:, k] - L[k+1:, :k] @ U[:k, k]) / U[k, k];追赶法的关键行是 den = diag[i] - lo[i-1]*c_prev(这一行就是「不产生填充」的全部内容),随后回代 x[i] = d_[i] - c_[i]*x[i+1]。
C2. Jacobi / Gauss–Seidel / SOR 的实测对比(第 5 讲 §四·§五、第 6 讲 §三)
对 A=(3−1−14)、b=(1,1)T、x(0)=(0,0)T:构造三种迭代矩阵,跑迭代并画误差-步数曲线。
验收标准:① 实测 ρ(GJ)=0.2886751(与 231 差 <10−9)、∥GJ∥∞=0.3333333(差 <10−12)、ρ(GGS)=0.0833333(差 <10−9);② 到 ∥x(k)−x∗∥∞<10−8 的步数:Jacobi 14、GS 7、SOR(ω=1.02175) 5、SOR(ω=1.5) 25(初值与阈值固定后,步数应完全一致);③ 后验估计:对每个 k 断言 err[k] <= 0.5*diff[k](Jacobi、∞-范数),并打印实测比值(0.36 与 0.45 交替,界是 0.5);④ SOR 用分量式更新,且 ω=1 时与 GS 逐位相同(差 <10−15);⑤ 扫描 ω∈(0,2) 画 ρ(G(ω)),标出 ω=2 处 ρ=1。
参考实现要点
SOR 的关键行 x[i] = (1-w)*x[i] + w*(b[i] - A[i,:]@x + A[i,i]*x[i])/A[i,i]——右端一定要用本轮已更新的 x(「形式隐式、实际显式」);写成「整向量同时更新」就退化成 JOR,ω=2 的实验会得出错误结论。
C3. 幂法、反幂法与位移反幂法(第 7 讲 §四·§五·§七·§八)
对 A=(2113):① 幂法(2-范数单位化)迭代到 ∣Rk−λ1∣<10−10;② 记录每步误差比;③ 反幂法:禁止 np.linalg.inv,每步 u = np.linalg.solve(A, u/np.linalg.norm(u));④ 位移反幂法:取 p=3.6,对 A−3.6I 重复 ③。
验收标准:① 幂法步数 ≤15(实测 13),且 abs(R - 3.618033988749895) < 1e-10;② 特征向量残差 ∥Ay−Ry∥2<10−4(实测 1.3×10−5;残差比特征值误差落后一步,别把两个阈值设反);③ 第 5 步之后误差比落在 0.1459±0.005;④ 反幂法给出的 λ2 与 1.38196601125 差 <10−10;⑤ 源码里 linalg.inv 出现次数为 0(grep -c 断言);⑥ 位移反幂法在 3 步内把 ∣R−λ1∣ 压到 10−10 以下(实测第 2 步 2.6×10−8、第 3 步 1.7×10−12)。
参考实现要点
幂法顺序不能反(u = A @ (u/np.linalg.norm(u)));每步 R 用 y @ A @ y;位移反幂法只是把 solve(A, ...) 换成 solve(A - p*np.eye(n), ...),最后一步别忘 +p。
C4. 四种方法的实测收敛阶与步数(第 7 讲 §九、第 9 讲 §六·§七·§八)
对 f=x3−x−1(s=1.3247179572447460)实现二分法、牛顿法(x0=1.5)、割线法(x0=1,x1=2)、单点割线法(x0=2 钉住,x1=1.5)。
验收标准:① 牛顿实测阶(取最后两个稳定 k)∈[1.9,2.1](实测 1.99、2.00);② 割线实测阶 ∈[1.5,1.75](实测 1.53∼1.67);③ 单点割线法不要套阶数公式(数值噪声会把 p 算飞),改为断言相邻误差比最后 5 个都落在 0.42±0.01;④ 到 ∣xk−s∣<10−12 的迭代次数严格满足 牛顿(4)< 割线(8)< 单点割线(31)< 二分法(40);⑤ 二分法步数与预测 log2εb−a 相差 ≤1;⑥ 断言割线法前两步误差上升(3.2×10−1→6.8×10−1),从第 4 步才进入超线性段。
参考实现要点
牛顿的理论常数 2f′(s)f′′(s)=0.9319,可顺手打印 e[k+1]/e[k]**2 与它对照;单点割线法的分母 xk−x0f(xk)−f(x0) 要防 x_k == x_0;收敛阶估计式 p≈ln(ek/ek−1)ln(ek+1/ek) 是本页给的重测方法(不是课堂结论)。
C5. 范德蒙德条件数、拉格朗日插值与误差界复验(第 9 讲 §十一·§十二)
取节点 xi=i/n(i=0,…,n),构造 V=np.vander(nodes, n+1, increasing=True):① 打印 n=2,4,6,8,10,12 的 cond2(V);② 随机取系数 a,算 b=Va 后用 np.linalg.solve 反解 a^,报告相对误差 ∥a^−a∥∞/∥a∥∞;③ 实现 lagrange(nodes, vals, x),验证 li(xj)=δij;④ 复验 B8 的误差界:在 [0,2] 上取 20001 个点,断言 ∣L2(x)−ex∣≤6e2max∣x(x−1)(x−2)∣。
验收标准:① cond2 随 n 指数上升——n=8 时 >106、n=12 时 >109(实测 2.0×106、6.8×109);② n=8 时反解系数的相对误差比 eps 大至少 3 个数量级(实测 6.0×10−11≈2.7×105eps),n=12 时 >10−9(实测 3.7×10−8)——与 cond2×eps 同量级;③ 拉格朗日路线不解任何方程,「节点处精确」到 10−12,n=12 时也不碰范德蒙德矩阵;④ 误差界不等式在全部 20001 个点上成立,且实测最大误差不超过界的 0.5 倍(实测 0.45)。
参考实现要点
分母 ∏j=i(xi−xj) 只算一次并缓存(节点固定),分子每次重算——这就是「按点加权求和」的代价与收益;比较两路线时注意它们追求的是同一个多项式,差别只在求解系数的方式。
每题必须给方案 + 理由 + 失败判据 + 验证计划(四项缺一不可),篇幅半页到一页。
D1. 病态方程组的诊断与处置流程
背景:B2 的 A=(1111.001)(cond∞≈4.0×103,两行近似平行)。给出从「发现病态」到「交一份可信结果」的完整流程。
参考方案与评分要点
方案(四步):① 诊断——算(或估)cond=∥A∥∥A−1∥:小规模用范数手算(第 3 讲练的正是 ∥A∥1、∥A∥∞),大规模用已知解算例做扰动实验(把 b 扰动 δ,看解跑多远);② 定风险——把 cond 与「要求保留的有效位数」对照,预留 cond 倍的前向误差(补充篇 2.4 的界,需核实);③ 处置(按代价从低到高)——换提法/换基(插值问题从范德蒙德基换到拉格朗日基,第 9 讲 §十二)→ 正则化/增广(把病态问题换成良态问题的近似)→ 提高浮点精度(float64 → float128 / 任意精度库)→ 若只关心残差,改用残差判据停止;④ 绝不指望列主元或迭代改进「治好」条件数。
理由:条件数是问题的属性(第 9 讲的红线:存在唯一 = 能算出来);列主元只把消去乘数控制到 ∣m∣≤1,避免误差在消去过程中被继续放大,它不改变 ∥A∥∥A−1∥。
失败判据:① 残差小但解误差大——B2 的扰动实验里残差可以做到机器精度量级而解的相对误差仍是 100%,所以「残差小」不能作为收工标准;② 把 b 换成更高精度后解仍大幅变化(说明问题不在右端精度,而在条件数);③ 换基后条件数没有下降(换错了基,或问题根本不是「基」的问题)。
验证计划:小规模时用已知解构造算例,扫 cond 从小到大,记录「实测前向误差 / 理论界」的比值(B2 实测 100%/200%=0.5);大规模时用「残差 + 与高精度解对比」双轨;把每步的 cond 估计、采取的措施、剩余不确定度写进报告,形成可复现的检查表。
评分要点:四步流程具体(1 分);理由指到「问题属性 vs 算法改良」(0.5 分);失败判据可测量(0.5 分);验证计划含高精度对照(1 分)。
D2. 为什么 SOR 的 ω 不能取 2
不只回答「定理说不行」,而要给出可验证的完整论证。
参考方案与评分要点
方案:三步走——(a) 用必要条件算边界:∣detG∣=∣ω−1∣n,ω=2 时恰为 1;(b) 对具体矩阵直接把 G(2) 写出来求谱(B5(3):G=(−1−2132−32),trG=−35、detG=1,特征值 6−5±i11,模为 1);(c) 扫描 ρ(G(ω)) 并检查两个端点与内部的行为。
理由:ρ(G)=1 使充要条件 ρ(G)<1 不成立——特征值模为 1 时误差既不衰减也不放大,而是按复特征值旋转;且 ∣detG∣=1 意味着「模之积为 1」,不可能推出「每个模都小于 1」。这正是第 6 讲 §五「区间外一定不收敛」的端点形态:ω=2 卡在必要条件的边界上,既不能保证收敛,具体到本题更是 ρ=1。
失败判据:① ω→0+ 时 N=L+ω−1D 的元素发散(1−ω−1→∞),迭代失去意义——注意区间是开的(第 6 讲考注 1);② ω∈(0,2) 内部、但矩阵非对称正定时 ρ(G)≥1——只查区间不查谱就会漏(Ostrowski 定理只对对称正定给充要);③ 「迭代若干步后误差不再下降」时,先怀疑 ρ≥1,而不是「精度不够」。
验证计划:对同一个 A 扫 ω∈(0.01,3.99)(步长 10−3)画 ρ(G(ω)),标出 ρ=1 的两个点;对 ω=2 单独跑一次迭代,打印前 20 步误差(应不衰减);用 ω≈1.022 对照(5 步进入 10−8)。脚本留在仓库里,作为「SOR 参数不能凭直觉」的证据。
评分要点:给出 ρ(2)=1 的推导而不只是结论(1 分);理由说到「必要条件的边界」(0.5 分);失败判据含「内部也可能不收敛」(0.5 分);验证计划可执行(1 分)。
D3. 从这 8 讲里挑一个算法:它在现代库里的替代品与理由
在「幂法 / 反幂法」「Jacobi–GS–SOR」「三对角追赶法」中挑一个(下面以第一个为例)。
参考方案与评分要点
方案:以「幂法 / 反幂法」为例:课堂算法只求一个特征对(模最大 / 模最小 / 最接近 p),收敛率被 ∣λ2/λ1∣ 支配;现代做法是 Lanczos / Arnoldi + 隐式重启(稀疏:scipy.sparse.linalg.eigsh / eigs,即 ARPACK),稠密对称全谱问题交 LAPACK 的对称特征求解器(dsyevr / dsyevd);此外还常用「位移-反迭代精修」把近似特征对打到机器精度(课堂的位移反幂法就是它的种子)。若挑「三对角追赶法」:现代对应物是结构保持的带状求解器(LAPACK dgtsv),课堂强调的 O(n)、无填充正是它的设计目标;若挑「Jacobi–GS–SOR」:它们在现代退到预条件子 / 平滑子的位置,主干换成 Krylov 子空间方法(CG / GMRES,第 4 讲 §五已点名 GMRES)。
理由:① 幂法的收敛率由谱间隙决定,间隙很小时慢到不可用;重启动 Krylov 每次用整个子空间的信息(等价于把多个幂法步叠在一起),对谱间隙不敏感。② 幂法 / 反幂法每步都要矩阵-向量乘或解方程组,Lanczos 只做乘,天然适合稀疏与矩阵免乘法(matrix-free)场景。③ SOR 的最优 ω 事先不知道(B5 里靠扫描;本 8 讲没有最优 ω 的公式),而 GMRES 无参数、CG 只需对称正定——少一个需要猜的参数在工程上就是少一类失败模式。
失败判据:① eigsh 的残差 ∥Ay^−λ^y^∥ 不随迭代下降 → 先核对目标语义(which='LM' 不等于「按模最大」的直觉,必要时用 sigma 做位移);② 稀疏矩阵却调了稠密 np.linalg.eig → 内存爆(对比 n=106 的存储需求);③ SOR 的 ω 扫描在目标矩阵上找不到明显极小 → 该矩阵不适合 SOR,换 Krylov + 预条件;④ 追赶法的 O(n) 只在结构保持时成立——一旦为选主元换行破坏了带宽,填充立刻出现,O(n) 失效。
验证计划:同一稀疏对称正定算例上比较自写幂法与 eigsh 的迭代数、残差、耗时(双对数图);用已知谱的小矩阵核对正确性(例如 B4 的 A,谱为 27±5);把「三种目标(模最大 / 模最小 / 最接近 p)」与库参数的对照表写下来(本 8 讲的「一个算法解决三个问题」在库里对应不同参数,这张表就是迁移的桥)。
评分要点:替代品具体到库与函数(1 分);理由指到收敛率 / 存储 / 参数敏感性(0.5 分);失败判据可测量(0.5 分);验证计划含已知谱核对(1 分)。
D4. 线性方程组求解器的选型流程(综合第 3–6 讲)
三个待求解问题:(a) n=100 的稠密非奇异方程组(可能病态);(b) n=106 的一维两点边值问题差分离散(三对角、行对角占优);(c) n=107 的三维离散问题(稀疏、非三对角)。给出选型流程与判据。
参考方案与评分要点
方案:一张按顺序问四个问题的流程表——
| 顺序 |
问题 |
结论 |
| 1 |
有结构吗? |
三对角 / 带状 → 追赶法类 O(n)(第 4 讲 §四);稠密 → 下一步 |
| 2 |
要反复解同一个 A 吗? |
要 → 先做 LU 分解,一次分解反复使用(第 4 讲 §二);只解一次 → 列主元消去 |
| 3 |
规模大且稀疏吗? |
大 → 迭代法:有行对角占优或对称正定才敢用 Jacobi / GS / SOR(第 5 讲 §四、第 6 讲 §八);否则上 Krylov(GMRES/CG,第 4 讲 §五) |
| 4 |
病态吗? |
cond 大 → 按 D1 处置(换提法/换基、加精度),别指望算法改良 |
| 三个问题的落点:(a) n=100 稠密 → LU(多右端复用时)或列主元消去,重点查 cond;(b) 三对角 → 追赶法,约 3n 次运算;(c) 大而稀疏 → 迭代法 / Krylov + 预条件。 |
|
|
| 理由:复杂度与稳定性是两条独立的轴。三对角 O(n)(约 3n)对比朴素消去 O(n3):n=106 时前者 3×106 次、后者约 32×1018 次(第 3 讲 §四的 1018 量级),差十几个数量级。稳定性一侧:行对角占优给的是充分条件(好算),ρ(G)<1 是充要但算不出(第 5 讲 §二·§三),SOR 的 ω 还额外依赖矩阵类型。 |
|
|
| 失败判据:① 用 ∥G∥∞<1 停机,却忽略「∥G∥ 接近 1 时相邻差小不能断定接近真解」的警告(第 5 讲 §三)→ 停机判据失效;② 对 n=106 稀疏矩阵用稠密 LU → 填充 + 内存爆;③ 对非对角占优的三对角矩阵直接追赶(不选主元)→ 主元可能为零,重演 B3 的卡壳——这时要么先做选主元的带状消去(可能破坏带宽),要么换迭代法;④ 病态问题用迭代法「磨」到收敛 → 得到的仍是病态问题的近似解。 |
|
|
| 验证计划:① 已知解算例:造 xtrue、算 b=Axtrue,比较解相对误差与归一化残差 ∥A∥∞∥x^∥∞∥b−Ax^∥∞;② 结构测试:断言追赶法在 n=1000 时 L、U 的非零元数为 2n−1;③ 规模测试:n 从 102 扫到 104,双对数图上斜率应为 1(三对角)与 3(稠密);④ 病态测试:把 B2 的 A 扩成 n×n 的「相邻行近似成比例」版本,观察 cond 增长与解质量退化。 |
|
|
| 评分要点:流程表四问清晰(1 分);理由含复杂度与稳定性两条轴(0.5 分);失败判据可测量(0.5 分);验证计划含斜率与残差归一化(1 分)。 |
|
|
| 组 | 题量 | 每题 | 小计 | 建议用时 | 判分要点 |
| A 概念判断 | 10 | 2 | 20 | 20 min | 判定对错 1 分 + 理由 1 分;只写「对/错」不给分 |
| B 推导与计算 | 8 | 6 | 48 | 90 min | 关键中间步骤 4 分 + 数值答案 2 分;缺步骤按比例扣 |
| C 编程实现 | 5 | 4 | 20 | 120 min(含上机) | 每条验收标准都要跑到;少跑一条扣 0.5 分 |
| D 综合与开放 | 4 | 3 | 12 | 60 min | 方案 1 分 + 理由 0.5 分 + 失败判据 0.5 分 + 验证计划 1 分;四件缺一不可 |
| 合计 | 27 | — | 100 | 约 290 min | — |
档位:≥90 扎实;75∼89 基本掌握(按下面的漏点定位表补两处);60∼74 需回读讲次页对应小节;<60 从 A 组错题对应的讲次重来一遍,再重做 B 组。
计分注意:B 组「数值答对但理由错」按半分处理——这 8 讲失分最多的地方是条件(ρ 与 ∥G∥、必要与充分、局部与全局),不是算术。
关键数值的复核方式
本页所有可复算的数值——B2 的条件数与扰动实验、B3 的消元与 LU、B4 / B5 的谱半径与 SOR 扫描、B6 的幂法 / 反幂法序列、B7 的四种迭代步数与实测阶、B8 的插值误差界、C5 的范德蒙德条件数——均已用 NumPy(v2.4.6,float64)复核,脚本片段如下(与本文数值一一对应):
import numpy as np
A = np.array([[1.,1],[1,1.001]]); print(np.linalg.cond(A, np.inf)) # 4004.001 —— B2
Ab = np.array([[3.,-1],[-1,4]]) # B4:ρ(G_J)=1/(2√3)
GJ = -np.linalg.solve(np.diag(np.diag(Ab)), np.tril(Ab,-1)+np.triu(Ab,1))
print(np.linalg.eigvals(GJ), np.linalg.norm(GJ, np.inf)) # ±0.2887, 0.3333
def G_sor(A, w): # B5:ω=2 的谱
D = np.diag(np.diag(A)); L = np.tril(A,-1); U = np.triu(A,1)
return np.linalg.solve(L + D/w, (1/w - 1)*D - U)
print(np.linalg.eigvals(G_sor(Ab, 2.0))) # (-5±i√11)/6,|λ|=1
f = lambda x: x**3 - x - 1; s = 1.3247179572447460 # B7:实测阶
e = [abs(1.5 - s)]; x = 1.5
for _ in range(4): x = x - f(x)/(3*x**2 - 1); e.append(abs(x - s))
print([np.log(e[k+1]/e[k]) / np.log(e[k]/e[k-1]) for k in (2, 3)]) # [1.9922, 2.0000]
复核记录:cond∞=4004.001;ρ(GJ)=0.2886751、∥GJ∥∞=0.3333333;ω=2 时特征值 6−5±i11、∣λ∣=1;幂法第 4 步 R4=47/13=3.615385;步数(阈值 10−12)牛顿 4 / 割线 8 / 单点割线 31 / 二分 40;范德蒙德(节点 0,…,1,n=12)cond2=6.8×109。
凡标「需核实」的结论——log10cond 与有效数字损失的对应、Rayleigh 商的「平方加速」、SOR 最优 ω 的显式公式、ρ(GGS)=ρ(GJ)2 的「平方关系」——都不是这 8 讲的课堂结论,属补充篇 / 原书或需另行推导的延伸;答题时若要引用,请按需核实处理。