数值分析上机编程题看着比理论考试轻松无非是求根、解线性方程组、插值拟合、数值积分、常微分方程这几大类。但真正坐到电脑前开始写代码才会发现教材公式和能跑出正确结果的程序之间隔着大量细节。我带实验课那几年每一轮都会见到类似的翻车现场牛顿法直接除零、高斯消元结果飘到没边、插值曲线末端翘得离谱、定积分越算越离谱。这篇文章就把几道典型上机题里最容易踩坑的地方整理出来给正在做这类作业、准备相关实践考试、或者自学数值计算的人做个参考。我不会只贴代码更多是解释每一步为什么要这么写以及检查结果时容易忽略什么。1. 非线性方程求根牛顿法既送分又送命的几个细节1.1 题目背后的考察点不只是“会写循环”牛顿迭代公式x_{k1} x_k - f(x_k)/f(x_k)确实一行就能写完。但上机题真正想考察的是三件事能不能选对初值、有没有把收敛条件写对、迭代终止条件设得合不合理。我见过不少同学把max_iter设成 100循环结束后直接打印结果根本不检查残差输出的点是离真实根很远的数值自己还以为算完了。这种问题在报告里几乎藏不住——老师只要运行一遍你的代码或者扫一眼误差表立刻就能发现问题。上机题里出现频率最高的例子是f(x)x^3 - x - 1 0正实根在 1.3247 附近要求精度到1e-6。看起来简单但如果你把初值随手取成 0 或者 5牛顿法的表现会完全不同。选初值不是“随便给一个数”而是要大致知道根在哪个区间。最稳妥的做法是先画个函数图或者用二分法夹逼出一个小范围再交给牛顿法精细逼近。1.2 牛顿法代码骨架导数、终止条件、除零检查这里给一个可以直接用的 Python 示例。核心是四点解析导数优先、终止条件要双保险、迭代次数要设上限、除零要显式处理。def newton(f, df, x0, tol1e-8, max_iter50): x x0 for k in range(max_iter): fx f(x) dfx df(x) if abs(dfx) 1e-14: raise ValueError(导数为零或过小请更换初始点) x_new x - fx / dfx if abs(x_new - x) tol and abs(f(x_new)) tol: return x_new, k 1 x x_new raise RuntimeError(迭代不收敛请调整初值或改用其他方法)为什么终止条件要同时看|x_new - x|和|f(x_new)|只判断前者有一个隐蔽问题遇到重根时牛顿法会退化为线性收敛相邻两步之间的差可能已经很小但函数值离 0 还很远。比如(x-1)^2 0初值取 1.2迭代几次后Δx看着没问题残差却还停在1e-3量级。双条件可以避免这种“看起来收敛了其实没到根”的假象。还有一点经验如果题目给的函数求导很麻烦不要自己去写一个差分近似然后祈祷步长选得合适。差分步长选大了精度差选小了浮点相减把你吞掉。这种时候直接用弦截法它不需要导数收敛阶约 1.618比二分法快得多代码也就多两行。1.3 牛顿法翻车现场与弦截法兜底牛顿法最怕两类情况一是f(x)在迭代过程中接近零二是一开始就飞出了收敛域。f(x)很小意味着切线几乎是水平的下一步会跳到很远的地方经常一步就飞到10^6量级。这时候不要急着调初值拿二分法先固定出一个区间或者直接换弦截法更省事。def secant(f, x0, x1, tol1e-8, max_iter50): for k in range(max_iter): if abs(f(x1) - f(x0)) 1e-14: raise ValueError(弦截法分母为零请换初始点) x2 x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0)) if abs(x2 - x1) tol and abs(f(x2)) tol: return x2, k 1 x0, x1 x1, x2 raise RuntimeError(弦截法迭代不收敛)弦截法的思想是用前两步的割线斜率代替切线斜率。它的优势是不求导数但要求两个初始点不能太离谱。如果题目明确要求比较方法报告里可以放一张简单的对比表说明你在什么情况下选什么方法。方法收敛阶是否需要导数主要风险二分法1否收敛慢但稳定牛顿法2是导数接近零时发散弦截法约1.618否分母接近零时失效1.4 用收敛阶检测程序有没有写对这个习惯我非常推荐算完根之后把每次迭代的误差序列存下来观察相邻误差的比值。对牛顿法二阶收敛意味着e_{k1} / e_k^2大致趋于一个常数对弦截法比值接近1.618对二分法比值接近0.5。如果你实现的是牛顿法却测出e_{k1} / e_k趋于常数那基本可以断定代码写错了或者这个根是重根。实际操作可以这样先自己用高精度算出一个参考根r然后记录每一步的e_k |x_k - r|在循环里把比值打印出来。这一步看似多余却是确认程序正确性的黄金标准。上机报告里加一列“收敛比”老师一看就知道你不是单纯把公式抄成代码而是真的理解了迭代法的行为。2. 线性方程组上机主元策略和迭代收敛是扣分高发区2.1 顺序消元的隐患比想象中严重线性方程组的经典上机题是“用高斯消元法求解Ax b”。很多同学直接按教材的教科书式子写顺序消元主元一路取A[k][k]。这种写法在课堂例题上表现很好因为例题的矩阵都是精心设计的主元不会太小。但一旦遇到现实中的病态矩阵比如希尔伯特矩阵H_ij 1/(ij-1)5 阶就能让人头疼。顺序消元的结果可能和小数点后几位就开始出错而列主元版本明显更稳。原因不复杂浮点数的相对误差会被除法放大。顺序消元遇到一个绝对值很小的主元相当于把后面所有误差都放大了好几倍。列主元做的事情很简单——每次消元前在当前列的下方找绝对值最大的元素把那一行换到主元位置。这个操作代价极小收益却巨大。很多上机题没明说“必须用列主元”但评分时参考实现就是默认带主元策略的版本。2.2 列主元高斯消元的实现细节下面这段代码可以直接抄进作业但更重要的是看懂每一个步骤import numpy as np def gauss_pivot(A, b): A np.array(A, dtypefloat).copy() b np.array(b, dtypefloat).copy() n len(b) for col in range(n): pivot col int(np.argmax(np.abs(A[col:, col]))) if np.abs(A[pivot, col]) 1e-12: raise ValueError(主元过小矩阵可能奇异) if pivot ! col: A[[col, pivot]] A[[pivot, col]] b[[col, pivot]] b[[pivot, col]] for r in range(col 1, n): factor A[r, col] / A[col, col] A[r, col:] - factor * A[col, col:] b[r] - factor * b[col] x np.zeros(n) for i in range(n - 1, -1, -1): x[i] (b[i] - A[i, i1:] x[i1:]) / A[i, i] return x需要注意几个点。第一判断主元太小时不是只用 0而是设一个阈值1e-12因为浮点下的“数值奇异”比“精确奇异”更常见。第二行交换时A和b要一起换漏掉b是高频失误。第三回代是从最后一个方程倒着来方向错了结果全错。如果你还想展示进阶理解可以在报告里提一句列主元消元本质上就是带排列矩阵的 LU 分解。消元过程中把factor存到矩阵下三角部分就得到了 LU 因子。这样同一个A对应多个不同的b时可以只消元一次反复回代。这个点在“数值分析”课程里绝对加分。2.3 Jacobi与Gauss-Seidel收敛判断别只看对角占优迭代法的上机题通常要求你“分别用 Jacobi 迭代和 Gauss-Seidel 迭代求解比较迭代次数”。很多教材都会提“严格对角占优是收敛的充分条件”于是有同学判断收敛时只看对角线是不是足够大。但老师设的题往往不会那么简单对角占优的条件可能不满足算法却仍然收敛反过来有些矩阵看起来符合对角占优但迭代矩阵谱半径大于 1照样发散。正确的做法是计算迭代矩阵的谱半径。设A D S_low S_up其中D是严格对角部分S_low、S_up分别是严格下、上三角部分不带对角。那么 Jacobi 迭代矩阵是B_j -D^{-1}(S_low S_up)Gauss-Seidel 迭代矩阵是B_gs -(D S_low)^{-1}S_up。谱半径ρ(B) 1才是充要的收敛条件。D np.diag(np.diag(A)) S_low np.tril(A, -1) S_up np.triu(A, 1) Bj -np.linalg.inv(D) (S_low S_up) Bgs -np.linalg.inv(D S_low) S_up rho_j max(abs(np.linalg.eigvals(Bj))) rho_gs max(abs(np.linalg.eigvals(Bgs))) print(fJacobi谱半径: {rho_j:.6f}, Gauss-Seidel谱半径: {rho_gs:.6f})经验是Gauss-Seidel 通常比 Jacobi 收敛快但“通常”不等于“总是”。也存在某些非对称矩阵Jacobi 收敛而 Gauss-Seidel 发散。不要在报告里写“Gauss-Seidel 一定更好”——这句话会被老师一眼看穿。老老实实列出两个谱半径和实际迭代次数才是科学的交代。2.4 病态方程组残差小不等于解准这是一道很经典的“坑”题。你辛辛苦苦解出一个x代回去算残差r A x - b发现残差只有1e-12于是认为答案很准。但拿去和真解一比可能差了好几个数量级。问题不在你的算法而在矩阵本身病态。量化病态程度的指标是条件数cond(A)。误差界可以写成||x - x_true|| / ||x_true|| cond(A) * ||r|| / ||b||也就是说残差小只说明“剩余小”但误差的放大倍率由条件数决定。希尔伯特矩阵是典型的病态例子5 阶条件数已经到1e5量级再大几阶就肉眼可见地离谱。用这种矩阵测算法列主元消元仍然会有误差但比顺序消元好得多。上机报告里如果能顺手输出np.linalg.cond(A)并解释“为什么残差小但解不准”这比做十道题都更能体现你理解了数值分析的灵魂。3. 插值与拟合差商表递推比拉格朗日公式更值得写3.1 拉格朗日表达式与两个坑插值题里最容易被写出来的就是拉格朗日形式因为它和教材公式一模一样一堆连乘再求和。但直接用有一个明显的问题每增加一个节点所有基函数都要重新计算一遍。如果只是求三个点的插值函数值无伤大雅如果上机题要求你“在密集网格上画出插值曲线”拉格朗日的计算量就会让人很难受。还有一个隐藏细节直接照抄公式时分母的x_i - x_j在i j时需要跳过初学者经常忘记结果除零。此外当节点数量较多时乘积项出现极小数或极小数浮点累加会带来额外误差。不是说拉格朗日不能用而是对于上机题来说牛顿差商形式几乎总是更优的选择。3.2 牛顿差商表的实现过程牛顿插值的核心是差商表。上机题里最方便的实现是用一个数组滚动计算下面是标准写法def divided_diff(x, y): n len(x) c np.array(y, dtypefloat) for j in range(1, n): for i in range(n - 1, j - 1, -1): c[i] (c[i] - c[i - 1]) / (x[i] - x[i - j]) return c def newton_eval(c, x, t): n len(c) p c[n - 1] for i in range(n - 2, -1, -1): p p * (t - x[i]) c[i] return p这里的循环方向特别容易写错。内层循环是从后往前覆盖c[i]分母是x[i] - x[i - j]不是x[i] - x[i - 1]。j表示当前计算的是几阶差商分母必须跨j个节点。如果你把分母写成相邻节点的差结果完全不对。为什么用差商表更好因为它是递推结构增加一个新节点时前面的差商可以复用不需要从头重新计算。而且求值时用的是嵌套乘法的 Horner 形式计算效率高浮点累积误差也小。实验课上让一组用拉格朗日、另一组用牛顿差商跑同样的 20 个点插值时间差距肉眼可见。3.3 龙格现象等距节点高次插值的典型翻车插值题里有一个非常经典的“钓”人案例f(x) 1 / (1 25x^2)在区间[-1, 1]上等距取节点做高次插值。取 11 个点做 10 次插值曲线在两端剧烈振荡你以为加密节点到 21 个会变好结果振荡更严重。这不是编程 bug而是等距节点高次插值的固有缺陷。改进方案有两个方向。一是改用切比雪夫节点x_k cos((2k 1)π / (2n 2)), k 0, 1, ..., n切比雪夫节点在两端更密、中间更疏能显著抑制边缘振荡。二是放弃一个多项式穿过所有点改用分段低次插值比如三次样条。只要被插值函数在该区间上没有突变分段低次的稳定性远好于单条高次多项式。上机报告里我建议你画两张图一张是等距节点的插值结果一张是切比雪夫节点的插值结果叠加上原始函数曲线。这样一个直观的对比比你写十行文字都管用。3.4 观测数据带噪声时不要硬插值还有一种题型不是“插值”而是“拟合”给你一堆实验测量点让你求一条近似曲线。很多同学第一反应还是让曲线穿过每一个点于是用了高次插值。结果曲线在数据点之间剧烈抖动完全看不出趋势。原因在于测量数据本身带有噪声噪声没有规律硬让曲线穿过每个点等于把噪声也拟合进去了。这时候应该用最小二乘拟合。先尝试一次或二次多项式看残差分布是否随机如果残差有明显结构再考虑升高次数。但不要一上来就拟合到 10 次以上——高次多项式拟合两端同样会剧烈波动而且正规方程条件数极差。上机题常见的做法是算最小二乘的正规方程或者直接用现成工具。但用工具不等于可以做你得解释“为什么拟合次数选 3 而不是 15”这个判断依据才是评分重点。4. 数值积分与常微分方程复化公式只是及格线误差验证才是加分项4.1 复化梯形和辛普森代码简单但精度差异巨大数值积分题最基础的是复化梯形公式和复化辛普森公式代码确实不长。def trapz(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 s * h def simpson(f, a, b, n): if n % 2 1: n 1 h (b - a) / n s f(a) f(b) for i in range(1, n): s 4 * f(a i * h) if i % 2 1 else 2 * f(a i * h) return s * h / 3从数学性质上说复化梯形是二阶精度误差大致是O(h^2)复化辛普森是四阶精度误差大致是O(h^4)。对光滑函数同样的n辛普森的结果会漂亮得多。很多同学只实现了梯形公式就交了结果老师一问“你为什么不用辛普森”就答不上来。如果题目没指定方法优先做辛普森它比梯形多不了几行代码但精度完全不同。4.2 用减半步长法验证收敛阶上机题里最容易加分的动作是“验证收敛阶”。做法是比如用n10, 20, 40, 80分别计算积分值观察误差随步长减半的变化。对梯形公式误差比应该接近 4对辛普森公式误差比应该接近 16。方法n10n20误差比约复化梯形8.90e-32.20e-34复化辛普森4.20e-52.60e-616如果你算出来的误差比和理论值差很多意味着要么代码有 bug要么被积函数不够光滑。这个“自检”动作在报告里写出来比贴十行代码更有说服力。4.3 高斯-勒让德上机题里的隐藏进阶方向某些上机题不满足于复化求积公式会直接要求用两点或三点高斯-勒让德公式求∫_{-1}^{1} f(x) dx。两点公式的节点是±1/√3权重都是 1三点公式的节点是0和±√(3/5)权重分别是8/9和5/9。一般区间[a, b]要先做线性变换x (a b) / 2 (b - a) * t / 2积分结果再乘(b - a) / 2。高斯公式的厉害之处在于代数精度高两点公式对三次多项式精确成立三点公式对五次多项式精确成立。但要注意它的节点不是等距的所以你不能把它当成“加密版梯形公式”来理解。上机实验时用同一个函数比较复化辛普森和高斯公式的误差会非常直观。4.4 常微分方程初值问题欧拉、改进欧拉与RK4的步长烦恼常微分方程的上机题最常见的是这个模型y y, y(0) 1求y(1)精确值是e ≈ 2.718281828。用欧拉法步长取0.1算出来大概是2.5937误差明显用四点经典 RK4 方法同样步长算出来已经能对上2.7183的若干位。差距不是一点点而是数量级的。def rk4(f, a, b, y0, n): h (b - a) / n t a y y0 for _ in range(n): k1 f(t, y) k2 f(t h / 2, y h * k1 / 2) k3 f(t h / 2, y h * k2 / 2) k4 f(t h, y h * k3) y y h * (k1 2 * k2 2 * k3 k4) / 6 t t h return y如果步长减半欧拉法的误差大致也减半而 RK4 的误差大致降到原来的1/16。这就是“四阶方法”的含义。上机题要求比较不同方法时你要在报告里写清楚这个误差变化规律而不是只贴“谁算得更准”。还有一点要强调显式方法遇到刚性问题时会表现得很惨。如果你发现步长稍微放大数值解直接爆炸不要急着怀疑代码有 bug先想想是不是方程本身的性质导致稳定性限制。4.5 自适应步长的加分思路如果学有余力可以在 RK4 的基础上实现一个简单的自适应步长先以步长h算一步再以两个h/2算一步两者之差就是误差估计。如果估计误差大于容差就把步长减半如果误差远小于容差就尝试放大步长。这个思想会贯穿你以后遇到的所有数值仿真问题上机题里主动提一嘴老师会觉得你已经“跳出作业本身”了。5. 上机报告的隐形评分点可复现性、误差分析与查错路线5.1 报告里最容易被忽略的“环境与可复现”很多同学只交一份源代码里面连注释都没有更不要说运行环境说明。老师打开后要么缺库跑不起来要么版本不对报错第一印象就差了。上机报告至少应该包含语言和版本、依赖库、输入输出示例、运行一次的结果截图或文本。代码里不要写绝对路径不要依赖当前目录下别人没有的数据文件所有全局参数放在文件开头集中定义。这不是形式主义是“可复现性”的最低要求。如果你的代码用到了随机数比如某些蒙特卡洛验证记得固定随机种子。不然别人复现不出你的结果只会认为你的结果不可信。5.2 三个典型“看起来对、实际错”的案例第一个案例牛顿法只看Δx终止。前面提过重根会让相邻迭代差变小但残差仍然很大。你输出1e-8的Δx觉得精度很高实际函数值还有1e-3量级。修正方法是同时看残差。第二个案例高斯消元得到了一个天文数字般的解。打印主元会发现某个A[k][k]已经小到1e-15。这时候程序没有“报错”只是输出了一个无意义的结果。很多同学第一反应是“是不是回代公式写错了”实际上应该先检查主元。第三个案例被积函数在区间内震荡积分区间取得很大。复化梯形公式用n1000仍然误差巨大甚至数值结果出现负值。这种问题用再密的等距步长也未必解决需要用自适应积分或先对被积函数做变换。报告里能意识到“这个结果为什么不对”比结果本身更有价值。5.3 通用查错路线我发现很多同学调试数值程序时完全没有章法报错了就乱改参数最后改对了也不知道自己错在哪。下面这条路线基本能覆盖大多数上机题现象可能的根因检查方法牛顿法除零初值导致导数太小打印每次迭代的导数值线性方程解爆炸主元过小或矩阵奇异打印主元、计算条件数插值曲线末端翘起等距节点高次插值改用切比雪夫节点积分结果异常步长不够或被积函数震荡减半步长观察误差比ODE解发散步长超过稳定性限制缩小步长检查刚性先拿一个已知解析解的小例子验证算法线性方程组用 3×3 对角占优矩阵ODE 用y y积分用多项式函数。小例子跑对了再上复杂题目。过程中打印中间变量不要只打印最终结果。看到NaN或Inf时往前翻输出通常能找到源头。5.4 输出格式与浮点精度的几个提醒最后提几个琐碎但影响体验的细节。输出数值结果时尽量使用科学计数法比如%.6e这样误差对比一目了然而不是默认打印出十几位小数让人眼花。做浮点比较时不要用“完全相等”要留一个容差。循环里不要反复计算不依赖循环变量的常量比如区间长度先存到变量里。这些小习惯看起来不起眼但批改上机报告时整洁的输出格式能让你在第一眼就和其他人拉开差距。如果时间允许我建议每道题都补一小节“误差分析”哪怕只是两三行为什么选择这个方法、步长或节点数如何确定、结果和解析解的差距在什么量级。这一部分才是数值分析上机题真正想考察的东西——不是把公式翻译成代码而是知道你写的代码在什么条件下可靠、在什么条件下会失效。这个能力说实话比代码本身值钱得多。我自己后来做仿真项目遇到诡异结果时用的还是当年上机课总结出来的这套思路先查解析解、再打中间量、最后验证收敛阶。一次性能把事情想清楚省下来的调试时间远比多写几行代码多。