数值分析上机题实战：迭代法、误差分析与避坑指南

简介：内容包括东南大学出版《数值分析》教材的MATLAB上机题解答，涵盖三次样条插值函数、重积分复化梯形公式与外推计算、常微分方程初值问题的RK4、AB4及AB4-AM4预测校正方法。文档给出每道题的算法推导、完整MATLAB实现思路、计算结果与误差对比，并附有作者的调试心得和反思，适合数值分析课程的学生作为实验参考或复习资料。资源共1个docx文件，压缩包约566KB，结构紧凑，便于下载后直接查看。目前已有904人学习下载。文中不仅呈现最终结果，还记录了从公式输入到程序调试的常见错误与解决办法，有助于提升读者对数值算法的理解和编程能力。

1. 数值分析上机题到底是什么：从一道不收敛的迭代式说起

某次做结构分析模拟，需要解一个带状稀疏线性方程组。我图省事，直接上了最普通的雅可比迭代，结果跑了三千多步数值还在震荡，最后换高斯-赛德尔加松弛因子才勉强压下去。这个坑，我最早其实是在数值分析上机题里撞过的——当时只当是作业，没沉淀成经验，结果工作里又翻了一次车。

数值分析上机题，对外行看是“用程序实现数学算法”，对做过的人来说，它真正练的是三件事：第一，把收敛性、误差估计这些理论结论落到具体算例上验证；第二，在有限精度、有限初值、有限计算量的约束下，判断一个算法到底能不能用；第三，学会用程序去发现“理论很完美、实践掉链子”的边界。无论是学生在赶作业，还是工程师想补数值计算的地基，这篇文章都按“题型判断 → 算法实现 → 验证交付 → 踩坑排查”展开，新手能照着一套步骤从零跑通，熟手可以重点看参数选择和边界条件。

2. 数值分析上机题的高频题型与算法选型：先把判断框架立住

2.1 五类经典题型与核心算法谱系

数值分析上机题绕不开五个主题：非线性方程求根、线性方程组求解、插值与逼近、数值积分、常微分方程初值问题。配套的算法每一类都有好几套，但上机题真正考的不是“会不会调库”，而是能不能把一个算法的迭代过程、误差变化、稳定性边界用程序讲清楚。

上机题一般要求至少覆盖其中两类。线性方程组迭代法是出现频率最高的，因为它既能卡你算法实现，又能卡你矩阵性质判断。记住一点：选算法先看问题结构，再看理论收敛保障。比如非线性方程，二分法稳定但线性收敛，牛顿法二阶收敛但初值敏感，两边权衡要看函数是否光滑、导数是否好求、初值能不能给准。

2.2 上机题真正难点不在写代码：误差、稳定性与计算量

很多新手以为上机题是“把公式翻译成代码”，跑完就完事。实际上，上机题的判分点和学习重点都在误差分析。离散误差、舍入误差和收敛性怎么权衡，才是这类题目真正要练的东西。

离散误差来自用有限步代替无限过程。比如复化辛普森公式，把积分区间切成 n 段，步长减半，误差大约降为原来的 1/16，这就是四阶精度。拿 sin(x) 在 [0, pi] 上的积分做实验，解析值是 2，n=10 时误差在 1e-4 量级，n=20 时降到 1e-5 量级，亲眼看到这个下降趋势，才算真的理解了“收敛阶”是什么意思。

舍入误差则来自浮点数有限精度。实际做项目时你会发现，步长不能无限缩小。把辛普森积分步长从 1e-3 继续往 1e-5 减，误差反而会回升，因为每一步的计算误差在累积，这就是精度拐点。上机题恰好提供一个可控环境让你观察这个现象：套路是先验算解析解，再逐步加密步长画误差曲线，看到曲线先降后升才叫真正理解了数值方法。

计算量是第三维约束。高斯消去法对稠密 n 阶矩阵是 O(n^3)，迭代法每轮 O(n^2) 但可能要迭代很多轮。同样的 n=1000 矩阵，直接法和迭代法的耗时差距能到几十倍。我一般建议在手算推导时先估算一下计算量级，再决定是上迭代法还是直接法。

2.3 算法选型的三个判断依据：收敛阶、条件数与停机准则

判断一个上机题该用哪套算法，我建议从三个维度入手。

第一是收敛阶。理论课会告诉你牛顿法二阶收敛、二分法线性收敛，但上机题里真正重要的是“初值敏感性”。牛顿法二阶收敛的前提是初值落在收敛域内，很多题目给的是随意初值，一跑就飞。稳妥做法是先用二分法或画图缩小根的区间，再用牛顿法加速。这种“粗定位+精收敛”的组合，是我在大量上机题里见过最稳的套路。

第二是条件数。线性方程组里，残差小不等于误差小。同一套迭代法，在良态矩阵上 20 轮收敛，换到条件数 1e10 的病态矩阵上，可能迭代 1000 轮误差还在原地打转。遇到题目给了一个看起来“普通”的矩阵，先算一下 cond(A)，超过 1e6 就要小心，考虑换直接法或先做预处理。

第三是停机准则。很多上机题要求“迭代到误差小于 1e-6”，但误差是拿什么定义的？是用相邻两步解的插值，还是残差范数，还是相对误差？这直接影响结果可信度。我一般同时看“迭代步增量”和“残差范数”两个指标，只有两个都低于阈值才判定收敛，避免出现“数值不动但离真解很远”的假收敛。

3. 用最小可复现代码跑通线性方程组：雅可比迭代与高斯-赛德尔对照

3.1 从一道经典上机题切入：构造一个收敛条件可控的算例

上机题最容易翻车的地方是随便编一个矩阵就开始迭代，结果不收敛也不知道是自己代码错还是矩阵本身问题。我建议先从严格对角占优的三对角矩阵入手，这类矩阵理论上保证雅可比和高斯-赛德尔都收敛，方便你把注意力放在算法实现上。

构造方式很直接：主对角线放 4，相邻对角线放 -1，得到一个对称正定的三对角矩阵。取精确解为全 1 向量，右端项按 A 乘全 1 向量算出来，这样你就有了一个“已知答案”的对照，后面所有误差分析都有了锚点。

import numpy as np def make_tridiag(n): # 构造 n 阶三对角矩阵：主对角 4，副对角 -1 A = np.zeros((n, n)) for i in range(n): A[i, i] = 4.0 if i > 0: A[i, i - 1] = -1.0 if i < n - 1: A[i, i + 1] = -1.0 return A n = 5 A = make_tridiag(n) b = A @ np.ones(n) # 让精确解为全 1 向量 print("A =", A) print("b =", b)

这段代码的重点在最后一行：用矩阵乘全 1 向量生成右端项，而不是手写一组看起来差不多的数。这样一来，你的验证目标非常明确——程序跑出来的解必须逼近 1，否则一定是算法实现出了问题，而不是数据本身的问题。参数方面，n 从 5 开始最合适，矩阵规模小，迭代次数少，方便打印中间结果；等小规模跑通了再把 n 改成 50 或 100，验证代码在更大规模下是否依然稳定。

3.2 雅可比迭代与高斯-赛德尔的最小实现

雅可比迭代的核心是把系数矩阵拆成对角部分和非对角部分，每次用上一轮的全部分量来计算这一轮的新值。它的最大优势是天然适合向量化，缺点是需要额外保存一整份旧向量。

def jacobi(A, b, x0, tol=1e-8, max_iter=1000): # 雅可比迭代：x_new = D^{-1}(b - (L+U) x_old) n = A.shape[0] D = np.diag(np.diag(A)) R = A - D D_inv = np.linalg.inv(D) x = x0.copy() for k in range(max_iter): x_new = D_inv @ (b - R @ x) err = np.linalg.norm(x_new - x, ord=np.inf) x = x_new if err < tol: return x, k + 1, err return x, max_iter, err

逻辑说明：每一步的 x_new 完全由上一轮的 x 决定，所以雅可比天然可以矩阵运算。D_inv 是对角矩阵的逆，这里直接用了 numpy 的 inv，实际工程里可以按元素相除省掉一次矩阵求逆。参数说明：tol=1e-8 是相邻两步解的最大绝对差，ord=np.inf 取无穷范数，也就是向量中变化最大的那个分量；max_iter=1000 是安全阀，防止不收敛时程序无限跑下去。

高斯-赛德尔的不同之处是：计算第 i 个分量时，前面 i-1 个分量已经更新了，要立刻用新值。这个“边算边用”的策略通常能加快收敛，但它天然是串行的，难以直接向量化。

def gauss_seidel(A, b, x0, tol=1e-8, max_iter=1000): # 高斯-赛德尔：第 i 个分量用本轮已更新的 x_new[:i] 和上轮的 x[i+1:] n = A.shape[0] x = x0.copy() for k in range(max_iter): x_new = x.copy() for i in range(n): s = b[i] - A[i, :i] @ x_new[:i] - A[i, i+1:] @ x[i+1:] x_new[i] = s / A[i, i] err = np.linalg.norm(x_new - x, ord=np.inf) x = x_new if err < tol: return x, k + 1, err return x, max_iter, err

这段代码有个细节值得注意：s 的计算里，前半部分用的是 x_new，后半部分用的是 x，这正是高斯-赛德尔与雅可比的本质区别。如果这里抄错，代码不会报错，但收敛速度会退化成雅可比甚至不收敛。参数上，内层循环逐行做点积，n 很大时会慢，后续可以用 scipy.sparse 的矩阵格式优化，但作为上机题演示已经足够。

3.3 主脚本设计：记录迭代次数与误差曲线

有了两个基础算法，下一步是写主脚本把它们放在同一个起跑线上对比。我用同一个初值、同一个停机准则、同一个矩阵，记录各自的迭代次数和最终误差。

x0 = np.zeros(n) x_j, it_j, err_j = jacobi(A, b, x0) x_gs, it_gs, err_gs = gauss_seidel(A, b, x0) print(f"Jacobi: {it_j} iters, err = {err_j:.2e}") print(f"Gauss-Seidel: {it_gs} iters, err = {err_gs:.2e}") print(f"Jacobi solution: {x_j}") print(f"GS solution: {x_gs}")

运行后你会发现，高斯-赛德尔的迭代次数通常只有雅可比的一半左右。这个现象符合理论：对这类三对角矩阵，高斯-赛德尔迭代矩阵的谱半径是雅可比谱半径的平方，收敛速度自然快一档。这里可以顺手做个试验，把 tol 改成 1e-12，观察两个算法的迭代次数如何增长，这比任何理论推导都直观。需要强调的是，这个谱半径平方关系只对特定矩阵类成立，不要当成普适结论。

如果还想看误差下降的曲线，就在算法里加一句 history.append(err)，最后用 matplotlib 画对数坐标图。观察重点不是谁更低，而是收敛曲线的斜率——斜率越陡，说明每一步误差衰减越快，这才是算法收敛性的直观体现。

3.4 参数怎么改：停机准则、松弛因子与初值的影响

上机题做到这里，很多人会误以为“跑通就等于会了”，其实参数实验才是拿高分和真正理解算法的关键。第一个参数是停机准则，把 tol 从 1e-6 收紧到 1e-12，迭代次数会明显增加，但解的精度的提升幅度未必成正比。要记录下“tol 每收一个量级，迭代次数增加多少”，这才是分析报告里能写出深度的地方。

第二个参数是 SOR 松弛因子。在高斯-赛德尔基础上引入一个 omega，公式变成加权平均的形式：

def sor(A, b, x0, omega=1.2, tol=1e-8, max_iter=1000): # SOR 超松弛：omega=1 时退化为高斯-赛德尔 n = A.shape[0] x = x0.copy() for k in range(max_iter): x_new = x.copy() for i in range(n): s = b[i] - A[i, :i] @ x_new[:i] - A[i, i+1:] @ x[i+1:] x_new[i] = (1 - omega) * x[i] + omega * s / A[i, i] err = np.linalg.norm(x_new - x, ord=np.inf) x = x_new if err < tol: return x, k + 1, err return x, max_iter, err

参数含义：omega > 1 是超松弛，通常能加速收敛；omega < 1 是低松弛，用于一些特殊问题。对上面的三对角矩阵，可以写一个循环从 1.0 试到 1.8，打印每种 omega 下的迭代次数，画一条“omega-迭代次数”曲线，找到最优松弛因子。这个实验本身就是极好的上机题素材。

第三个参数是初值。把 x0 从全 0 改成全 100，对比迭代次数变化。对收敛性有保障的矩阵，初值只会影响迭代轮数，不会改变是否收敛；但对条件数大的矩阵，初值可能直接决定发散方向。这个对比实验能帮你建立对“初值敏感性”的直觉，比背理论结论有用得多。

4. 上机题交付与验证套路：把实验做严谨再写报告

4.1 上机报告的五个标准模块

上机题最终交付物一般是程序加报告。我批改过不少作业，发现一个好报告和普通报告的分界线特别明显。普通报告是“贴代码+贴运行结果”，好报告则按下面五个模块组织，每一块都有明确目的。

第一个模块是问题陈述与数学模型。用数学语言把题目重新表达一遍，列出方程、边界条件、已知量。这个步骤看起来多余，实际上能逼你确认自己没有理解偏题意。很多翻车都是从一开始题意就理解错了，后边所有实现都是在错误的方向上狂奔。

第二个模块是算法选取与理由。回答一个问题：为什么选这个算法而不选另一个？要结合矩阵性质、收敛条件来说。比如选雅可比迭代，是因为矩阵严格对角占优、收敛理论完备；不选高斯消去，是因为题目规模大且稀疏，迭代法单轮成本更低。理由可以简单，但不能不写。

第三个模块是程序实现要点。这里写清楚迭代格式、停机准则、初值设置，以及代码里最关键的几行做了什么。不要贴全部源码，贴核心片段并解释。重点说明你如何保证实现与理论公式一致。

第四个模块是数值结果与误差分析。这是最见功力的部分。需要包含：不同参数下的迭代次数表、误差下降曲线、与解析解或参考解的对比。分析不能停在“结果对了”，要指出误差趋势是否符合理论收敛阶、在哪里出现偏离、为什么出现。

第五个模块是结论与讨论。写一两个实验中的观察，比如“松弛因子在 1.2 附近最优，但最优位置随矩阵规模的增大而偏移”。不必长篇大论，三到五条结论即可，但每条必须由前面实验数据支撑。

4.2 误差对比表与拐点分析：一组数据该怎么呈现

上机题报告里的数据呈现，我强烈建议用表格加拐点分析，而不是贴一堆 print 输出。下面是我常用的对比表格式，把矩阵规模、方法、迭代次数、最终残差、相对误差放在一行，一眼就能看出趋势。

以上数字是示意值，实际以你运行环境为准，但趋势是稳定的：n 增大时迭代次数大致线性增长，高斯-赛德尔始终约为雅可比的一半。写报告时不要只贴这张表，要指出拐点——n 从 5 增长到 50，迭代次数从 24 涨到 320，增长了大约 13 倍，这符合谱半径随矩阵规模逼近 1 的理论预期。

拐点分析的价值在于，它把一个“运行结果”升级成了“对算法行为的解释”。如果只写“n=50 时迭代 320 次”，导师或读者看不出你理解了原理；补一句“因为迭代矩阵谱半径随 n 增大而趋向 1，收敛因子接近 1 导致轮数线性上涨”，这就说明你是带着理论在分析数据。

4.3 验证技巧：用解析解做地基，用高精度库做天花板

上机题最怕的问题是没有参照系，你不知道自己程序跑出来的解到底对不对。我有两个常用的验证技巧。

第一个是构造有解析解的算例。就像前面那组三对角矩阵，精确解取全 1 向量，右端项由矩阵乘解析解生成。这样做的好处是你的验证标准是绝对明确的，不依赖任何外部库。实际执行时，对比相对误差而不是绝对误差，因为病态矩阵的解可能数值很大，绝对误差会误导判断。

第二个是高精度库做天花板。对线性方程组，最直接的参照是 numpy.linalg.solve，它走的是 LAPACK 的高精度直接法。把你的迭代法结果和它对拍，如果两者差在 1e-6 以内，基本不会有大问题；如果差异大，优先检查迭代法实现而不是怀疑库有问题。

对积分和微分方程也有对应的天花板方案。数值积分用 sympy 的解析积分做对照，常微分方程用 scipy.integrate.solve_ivp 的高精度模式做参照。这些库本身可以当黑匣子用，但它们的作用不是替代你手写算法，而是提供一个“正确答案”来暴露你的实现问题。

5. 数值分析上机题避坑实录：四条踩坑记录与排查路径

5.1 现象：代码看着没问题，迭代就是不收敛

这是我见过频率最高的翻车现场。初学上机题，自己写了一个迭代法，矩阵和右端项都用随机数填充，结果怎么调初值都不收敛，最后怀疑是自己代码写错了。

原因拆解：随机生成的矩阵大概率不满足对角占优条件，迭代矩阵谱半径大于 1，理论就不保证收敛。这不是代码 bug，是算例构造错误。另一个隐蔽原因是迭代矩阵写错了，比如雅可比迭代里把 R 的符号搞反，相当于每一步都在往错误方向修正。

解决路径：第一步，检查谱半径。对雅可比迭代，迭代矩阵是 I - D^{-1}A，用下面的代码算一下最大特征值模：

def check_convergence(A): # 仅适用于雅可比迭代：检查迭代矩阵谱半径是否小于 1 D = np.diag(np.diag(A)) B = np.eye(A.shape[0]) - np.linalg.inv(D) @ A rho = np.max(np.abs(np.linalg.eigvals(B))) print("谱半径 =", rho) return rho < 1

第二步，如果谱半径小于 1 但仍不收敛，再怀疑迭代格式写错，打印前两三轮 x 的变化，手算一遍核对公式。第三步，换一个严格对角占优的矩阵重试，排除数据本身的问题。先把“算例没问题”证明出来，再谈算法优化。

5.2 现象：误差输出一直是 0.0，但手算结果对不上

第二个常见坑是误判“完美收敛”。有一次我调试一个非线性方程求根的二分法，程序跑完打印误差 0.0，我以为是精度极高，结果把根代回原方程发现差了 0.1。

原因拆解：这是典型的浮点比较陷阱。代码里如果用相邻两步的差做停机判断，当迭代进入平台期或步长被机器精度截断时，相邻两步的差会变成 0，但解的绝对值精度并未达到理论预期。另一些情况是初值恰好就是精确解，迭代一步就“收敛”了，这个误差为 0 实际是偶然，并不是算法表现。

解决路径：停机判断一律用范数差，而不是x_new == x。同时把初值设得离精确解远一点，比如np.full(n, 100.0)，强迫算法真正走完整条迭代路径。最后在报告里同时报告“相邻步差”和“残差范数”两个指标，残差才是物理上真正要消掉的量。

5.3 现象：收敛了，但精度远低于预期

第三种情况最容易被忽略：程序正常退出、残差达到设定阈值，但解和真解的相对误差在 1e-3 量级，远达不到题目的 1e-6 要求。

原因拆解：问题出在矩阵病态性和停机准则的错配上。残差范数小不等于误差小，它们是两个不同量纲的东西。对病态矩阵，一个微小的残差可能对应巨大的解误差，因为 A 的逆把舍入误差放大了 cond(A) 倍。这属于数值线代里的“条件数决定误差上限”问题，不是多迭代几轮能解决的。

解决路径：算一下 cond(A)，如果矩阵条件数大于 1e6，先考虑改换算法而非盲目加迭代。常见做法是换高斯消去或者对矩阵做对角缩放预处理，把条件数压下来再迭代。另一个实用操作是把停机阈值收紧到 1e-12 再跑一遍，如果误差有明显改善，说明之前只是停机太早；如果误差纹丝不动，那就是条件数在作祟，继续迭代也没用。

5.4 现象：程序跑得非常慢，怀疑是复杂度没控制好

最后一条是性能问题。上机题的数据规模通常不大，但如果用纯 Python 逐元素循环写高斯-赛德尔，n=1000 时耗时会到分钟级，实验做到一半人就没耐心了。

原因拆解：双层 for 循环里做数组切片和点积，每次运算都产生 Python 对象开销。雅可比可以向量化却用了循环，高斯-赛德尔虽难以完全向量化，但至少可以用 numpy 切片减少逐元素操作。

解决路径：把雅可比改成向量化写法，用矩阵乘一次算完所有分量，速度提升一个数量级。高斯-赛德尔如果逐元素循环慢，考虑用矩阵分裂形式做块更新，或者直接改用 scipy.sparse.linalg.gmres 这类成熟求解器。上机题本身不要求极致性能，但如果实验矩阵规模放大，先做好向量化再继续，否则一个普通的参数扫描实验要等半个小时。

6. 把上机题沉淀成自己的数值工具箱：一个能复用的验证框架

上机题做完就丢，是最亏的用法。我自己的习惯是，每做一道题，就把算法函数收敛成一个统一接口的版本，放进自己的工具库，下次遇到同类问题直接复用。统一接口的收益在于：换算法只改一行字典配置，不用改实验主体代码。

solvers = { "jacobi": jacobi, "gauss_seidel": gauss_seidel, "sor": lambda A, b, x0: sor(A, b, x0, omega=1.2), } def run_all_solvers(A, b, x0=None, n_list=None): if n_list is None: n_list = [5, 10, 20, 50] if x0 is None: x0 = np.zeros(A.shape[0]) for n in n_list: A_n = make_tridiag(n) b_n = A_n @ np.ones(n) for name, solver in solvers.items(): x, it, err = solver(A_n, b_n, np.zeros(n)) rel_err = np.linalg.norm(x - np.ones(n)) / np.sqrt(n) print(f"n={n:3d} {name:14s} iters={it:5d} rel_err={rel_err:.2e}")

这个框架的三个习惯，是我做上机题攒下来的最有用经验：第一，停机准则作为参数暴露，不要写死在函数里，实验时需要频繁调整；第二，算例构造函数和算法求解函数分离，换矩阵只改一处；第三，相对误差统一计算，不同规模的问题可比。

这些经验迁移到工程里，对应的是产线数值模块的调试习惯。比如在图像处理里做积分算子近似，或在交易系统里用数值方法求隐含参数，都需要先造一个小规模算例验证算法行为，再放开到真实数据。没有上机题训练出来的“先锚定解析解、再验误差曲线”的肌肉记忆，直接拿真实数据调参，很容易被表面收敛骗过去。

我自己在这类训练里最大的教训是：永远先验证最简情况，再往复杂里加。简单算例跑通了，后面出的每个问题都至少有个已知参照，不至于一出错就无从下手。希望这篇笔记帮到你。

本文还有配套的精品资源，点击获取