ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

数值分析上机题通关指南:从算法原理到代码实现

数值分析上机题通关指南:从算法原理到代码实现

数值分析这门课,很多人是上课听得懂、作业不会写,到了上机题更是对着编辑器半天憋不出几行。原因是这类题目表面上考数学定理,实际上考的是“把数学变成能跑的程序”的能力,而这中间缺了不少过渡。这篇内容我按数值分析的几道上机编程题为线索,把从读题、选算法、写代码,到验证误差的完整流程拆开讲一遍。我不是来复述教材定理的,我要讲的是实际操作中你会怎么想、怎么做、以及在哪儿容易翻车。如果你正准备对付这类作业,或者已经写得头疼,这篇文章应该能帮你少走很多弯路。

1. 先看透题目设计:上机题真正在检验什么

1.1 上机题和我们以为的“写代码”不一样

很多学生的第一反应是:上机题不就是把公式翻译成代码吗?理论上好像确实是这样,但实际做起来你会发现,每一步都在做决定,而且每个决定都在挑战你对算法的理解。

比如教材里写二分法,收敛条件是区间两端函数值异号,迭代格式就是把区间一分为二。你觉得这还不简单?可一旦写代码,你要回答的问题就变成了:循环什么时候停?是看函数值绝对值足够小,还是看区间长度足够小?迭代次数超过多少算失败?如果一开始给的区间两端同号,程序是直接报错还是自动处理?这些细节教材不会替你回答,但上机题会考。

数值分析的上机题,本质上是把“算法的可靠性和效率”从抽象概念变成具体代码。老师给你一个方程、一个矩阵、一组离散数据点,你不仅要算出结果,还要在有限次迭代内得到满足精度要求的结果。这要求你对算法的适用条件、收敛性、误差来源都有清晰认识。所以我个人觉得,做上机题之前,最该做的一件事不是马上打开编译器,而是先问自己:这道题考的是哪个算法?这个算法在什么条件下能用?它的误差从哪里来?

1.2 上机题无非就是这几类,逃不掉的

把常见的数值分析上机题翻来覆去地看,你会发现出题范围其实相当稳定,基本集中在五个方向:非线性方程求根、线性方程组求解、插值与拟合、数值积分、常微分方程初值问题。

我通常会给学生一张对照表,让他们拿到题目先归类:

题目特征关键词对应的算法簇典型考核点
方程求根、迭代收敛、初值敏感二分法、牛顿法、割线法收敛条件、收敛阶、停机准则
线性方程组、矩阵分解、条件数高斯消元、LU分解、迭代法主元选取、存储方式、误差分析
已知数据点、插值、多项式逼近拉格朗日插值、牛顿插值、最小二乘拟合插值余项、龙格现象、系数计算
定积分、数值积分、区间分割梯形公式、辛普森公式、复合求积公式截断误差、步长选取、代数精度
微分方程初值问题、步进求解欧拉法、改进欧拉法、龙格库塔法局部截断误差、稳定性、步长控制

这张表看起来简单,但很有用。分类完成以后,很多问题已经被解决了一大半,因为你已经知道这道题在《数值分析》教材里属于哪一章,可以去翻对应的算法流程了。

2. 四类必考题型怎么拆解,代码怎么落地

2.1 非线性方程求根:别只知道套公式,要会处理迭代失败

求根类题目是上机题里最常见的,典型问法有“用二分法求方程 x³ + x − 1 = 0 在区间 [0,1] 内的根,误差不超过 10⁻⁵”或者“用牛顿法求 f(x) = x² − 2 的正根,比较不同初值的迭代情况”。

二分法的代码实现并不复杂,核心就是区间折半。但我会建议你在写之前把停机准则想清楚。一般有两个判断:一是区间长度小于给定容差,二是函数值绝对值小于容差。如果只判断函数值,存在风险,因为函数在某些区域很“平”,函数值已经很小,但离真实根还很远。如果只判断区间长度,那没问题,但要注意,如果题目要求精度是 10⁻⁵,区间初始长度为 1,二分法每轮减半,大约需要 log₂(1/10⁻⁵) ≈ 17 次迭代。这个预估很有用,算出来以后你就知道自己的循环次数设置成 100 是足够安全的,设置成 10 则很可能跑不到精度。

牛顿法的实现要复杂一点,因为多了一个导数的计算。我写过一个标准流程,大概长这样:

def newton(f, df, x0, tol=1e-8, max_iter=50): x = x0 for k in range(max_iter): fx = f(x) if abs(fx) < tol: return x, k + 1 d = df(x) if abs(d) < 1e-14: raise ValueError("导数接近零,牛顿法不适用") x = x - fx / d return x, max_iter

这里特别提醒一句:牛顿法对初值非常敏感。同样的方程,初值选得不好,可能会迭代到另外的根,甚至直接发散。所以我交作业时一般会准备两组初值,一组靠近真实根,一组故意选差一点的,然后在报告里写出观察结果。老师看到这种对比,会觉得你是真的理解了这个算法的性质,而不是照着书抄了一遍代码。

割线法就是不需要显式导数的替代方案,用差商代替导数。它的收敛阶低于牛顿法,但胜在每次迭代只需要算一个函数值。很多时候上机题会故意让你把牛顿法和割线法都实现一遍,这时候写一个统一的迭代框架会更方便对比。

2.2 线性方程组求解:主元选不好,结果差到离谱

线性方程组的上机题,一开始通常都很直接:给定一个系数矩阵和右端项,让你用高斯消元或 LU 分解求解。看起来就是三层循环的事,但里面有一个关键细节会直接决定你的结果是可信还是废掉——主元。

我不止一次见过学生直接按教材的“高斯消元”流程写,不处理主元为零的情况。打个比方,如果有矩阵

[ A = \begin{pmatrix} 0.01 & 1 & 1 \ 1 & 1 & 0 \ 1 & 0 & 1 \end{pmatrix} ]

第一行第一个元素 0.01 很小但不为零,如果直接把它当主元消去下方元素,中间会产生很大的中间数值,计算时会出现严重的舍入误差,甚至把真实解淹没掉。解决办法就是列主元高斯消元:每一步从当前列下方选择绝对值最大的元素,交换到主元位置,再做消元。

这个操作表面上只是多了几行交换代码,但效果差别很大。我在模拟项目X里测过同一个三阶矩阵,不选主元时解的误差可能到 10⁻³ 量级,选主元之后能到 10⁻¹⁴ 左右,差了十几个数量级。所以如果你发现自己的计算结果和手算结果对不上,第一反应应该回去检查主元处理。

LU 分解一般是在高斯消元的基础上顺手完成的。你现在消元过程中用到的乘数 (m_{ik}) 就是 L 矩阵下三角部分的元素,消元结束后的上三角矩阵就是 U。如果你的课程要求你实现“Doolittle 分解”或“Crout 分解”,那也逃不开这个思路。代码实现时要注意存储顺序,很多作业会用紧凑格式,把 L 和 U 放在同一个矩阵里,省存储空间。这种做法在实际工程里很常见,但刚开始写的时候容易把 l 和 u 搞混,给自己留点注释会省很多调试时间。

2.3 插值问题:别把高次多项式神话,也要防龙格现象

插值类的上机题,常见的任务是“已知某函数在若干点上的值,构造拉格朗日插值多项式或牛顿插值多项式,并计算指定点的函数值”。这类题目的代码难度不大,但有一个概念陷阱。

拉格朗日插值的公式长得很整齐:

[ L_n(x)=\sum_{i=0}^{n} y_i \prod_{j\neq i}\frac{x-x_j}{x_i-x_j} ]

但它的缺点是,每增加一个节点,所有基函数都要重新计算,相当于推倒重来。牛顿插值用差商构造,好处是增加节点时可以利用前面算过的差商表,计算效率高。我在作业里一般会两种都实现,然后用同一组数据对比结果,确保两个方法得到的插值多项式在相同节点上的函数值完全一致。这其实是很好的自检方式,因为两种方法理论上是等价的,如果你写出来的结果不一样,那一定是某处代码有 bug。

还有一类题目特别容易拿龙格现象说事:对函数 f(x) = 1/(1+25x²) 在区间 [-1,1] 上取等距节点做高次插值,次数越高,插值结果在边界附近震荡越厉害。我在第一次做这个实验的时候也觉得很反直觉,明明插值节点加密了,逼近效果怎么反而更差?这本质上是因为高次多项式的灵活性太强,在边界处为了满足所有节点的取值,会产生剧烈的振荡。实际教学中如果有这个题目,你最好把不同次数插值结果的图像画出来,在报告中写清楚“高次牛顿插值并不一定优于低次”,这比写一堆理论分析更直观,也更有说服力。

如果要提高逼近效果,一般会用到分段低次插值或样条插值。有些作业会进阶要求你用三次样条插值,这时候关键在于边界条件的选择:自然边界条件是两端二阶导为零,固定边界条件是两端一阶导已知。代码实现起来比拉格朗日麻烦不少,但对数据的平滑性改善很显著。这类题做一次,你对插值算法的理解就能上一个台阶。

2.4 数值积分与常微分方程:误差项才是真正的主角

数值积分的上机题有一个非常固定的套路:给一个定积分,让你分别用梯形公式和辛普森公式计算,然后比较不同等分数 n 下的误差变化。

复合梯形公式的误差阶是 (O(h^2)),复合辛普森公式的误差阶是 (O(h^4))。这意味着什么?如果你把步长 h 缩小一半,梯形公式的误差大约变成原来的 1/4,而辛普森公式的误差大约变成原来的 1/16。我每次做这类题目都会把不同 n 下的误差列成一张表,然后用这个规律去验证自己的实现。如果你发现辛普森公式的误差没有按照预期的速率下降,那很可能是在实现时把偶数节点的系数搞错了。

辛普森公式的代码逻辑很清晰,但一定要注意 n 必须是偶数。因为复合辛普森公式把区间分成 n 个小子区间,每两个区间对应一组 (4, 2, 4) 的系数结构,如果 n 是奇数,最后会剩下一段处理不了。很多同学会在这里踩坑,写着写着报错或者算出来的数明显不对。我一般会在代码开头加一个判断:

if n % 2 != 0: raise ValueError("辛普森公式要求等分数n为偶数")

至于常微分方程初值问题,最基础的是欧拉法和改进欧拉法,再进阶就是经典的四阶龙格库塔法。上机题问法经常是“用步长 h=0.1 和 h=0.05 求解初值问题,比较结果与精确解的误差”。这里你需要抓住一个核心概念:全局误差和步长的关系。欧拉法全局误差是 (O(h)),经典龙格库塔法是 (O(h^4))。如果你用 h=0.1 和 h=0.05 各算一遍,然后比较误差的比值,会发现 RK4 的误差激烈地缩小了大约 16 倍,而欧拉法只缩小了大约 2 倍。这种数值实验比你空口说“RK4精度更高”要有力得多。

写 RK4 的时候,我建议不要一次性把四步全都写在同一个表达式里,那样调试会非常痛苦。拆开写:

k1 = f(tn, yn) k2 = f(tn + h/2, yn + h*k1/2) k3 = f(tn + h/2, yn + h*k2/2) k4 = f(tn + h, yn + h*k3) yn += h*(k1 + 2*k2 + 2*k3 + k4)/6

这样每一步都有对应含义,出了问题一眼就能看出来是哪一步的系数或函数值不对。

3. 完整实操:一道典型题从拆题到验收的全过程

3.1 拆解数学模型,明确输入输出

这里我用一道比较综合的题目做例子,比如:“用二分法和牛顿法分别求方程 (x^3 - x - 1 = 0) 在区间 [1,2] 内的根,要求精确到 10⁻⁶,打印迭代过程中每个节点的值。”

拿到题,我习惯先画一条线:输出什么、输入什么、精度怎么判断。输入是方程表达式、区间左端点、右端点、目标精度;输出是近似根、迭代次数、每轮迭代的区间或节点值。这个过程叫抽象,很多人省略了它直接写代码,结果写到一半开始混乱,不知道该从哪里返回。

然后我会做一次手算或心算预估:f(1) = -1,f(2) = 5,所以区间两端异号,二分法可用。真实根大约在 1.32 附近(因为 1.3³ ≈ 2.197,减去 1.3 再减 1 约等于 -0.103;1.32 迭代逼近),这不是为了绕过程序得到答案,而是方便后面判断程序输出是否合理。我一直觉得,写数值程序前手算一次近似值,能帮你过滤掉大量低级错误。

3.2 写出主循环,把停机条件放在最前面

二分法的主循环一般有两种写法:一种是固定的 for 循环,比如设定最多迭代 40 次;另一种是 while 循环,条件是区间长度大于容差且迭代次数没超上限。我更推荐第二种,因为当算法不收敛时你不会卡死在无限循环里。

def bisection_solve(f, a, b, tol): if f(a) * f(b) > 0: raise ValueError("区间两端函数值同号,不满足二分法条件") fa = f(a) k = 0 while (b - a) / 2 > tol and k < 100: c = (a + b) / 2 fc = f(c) print(f"k={k}, a={a:.8f}, b={b:.8f}, c={c:.8f}, f(c)={fc:.8e}") if fc == 0: return c, k if fa * fc < 0: b = c else: a = c fa = fc k += 1 return (a + b) / 2, k

这里一个容易被忽略的细节是,更新左端点时要把 fa 同时更新为 fc,而不是下一轮重新算一遍 f(a)。虽然重新算也能得到正确结果,但多了一次函数求值,如果函数本身比较昂贵,效率就低下来了。老师如果仔细看代码,这点小优化是很加分的。

3.3 用误差和收敛阶做自检,而不是只看“算出来了”

很多同学交上机题,判断正确的标准是“程序跑完没报错,输出一个数”。但这个数到底对不对,完全没有验证。我自己的习惯是至少做两层验证。

第一层,把数值解代回原方程,看残差绝对值是否小于题目要求。比如求根得到 x=1.324717957,代回去 f(x) 应该接近 10⁻⁷ 量级。如果残差很大,说明程序的停机准则或者算法逻辑有问题。

第二层,做收敛阶实验。对于一个真正收敛且实现正确的算法,随着迭代次数增加或步长缩小,误差应该按照理论速率下降。这种实验不需要精确的解析解也能做,比如牛顿法在单根附近是二阶收敛的,你可以记录相邻两次迭代的误差比值,如果这个比值趋近于一个常数区域,而不是震荡,那说明程序大概率是对的。这个方法,我后来发现比单纯对答案好用得多,因为它能验证实现逻辑是否真的符合算法理论。

3.4 报告里哪些内容容易拿分,哪些画蛇添足

上机题往往要求交报告。我自己批改这种报告时最看重的是三个东西:程序是否完整可运行、结果表是否有明确的迭代过程、是否有对算法行为的分析与观察。

结果表里,我建议至少包含迭代序号、近似解、函数值、相邻近似解的差这四列。没有迭代过程的报告,就像一道只有答案没有过程的数学题,即使结果正确也说服力不足。牛顿法尤其要打印相邻两步的差值,因为判断它是否真正进入收敛区域,要比看函数值绝对值变化更直观。

至于“画蛇添足”的内容,最常见的是把一大段教材原话抄进报告,然后不加任何自己的分析。老师不笨,一眼就知道哪些是你自己写的,哪些是复制粘贴的。报告的价值在于你的观察和踩坑记录,不在字数。

4. 调试数值代码时踩过的那些坑,我帮你提前踩一遍

4.1 判断条件写成“函数值等于零”导致死循环或漏解

这是新手最容易犯的错误。很多人写二分法时喜欢用if fc == 0: return c,这个判断本身没错,但实际计算中,浮点数很少能精确等于零。更多时候,函数值在一个已经满足精度的区间里会像 10⁻¹⁶ 一样跳来跳去,永远不等于 0。所以停机准则应该用绝对值小于阈值,而不是严格等于零。这个阈值也不能随便拍脑袋,一般和题目要求的误差有关,我习惯设为 0.1 倍的目标误差甚至更小,避免提前停机导致结果不满足要求。

4.2 数组索引和循环边界错位,结果偏了却不知道

线性方程组的高斯消元,最容易出 bug 的地方是循环边界。三层循环里,外循环是消元的列 k,内循环是行 i 从 k+1 到 n-1,最内层是列 j 从 k+1 到 n-1。任何一层边界写错,都可能产生一个看似正常但实际错误的结果。更隐蔽的问题是当你实现列主元策略时,行交换以后右边端向量 b 也要跟着交换,很多人把矩阵 A 换了,却忘了换 b,结果算出来当然不对。这个错诊断起来特别难,因为输出结果看起来和正确答案有“一定的相似性”,但每一行都差一点。我后来干脆把参与运算的 A 和 b 打印出来检查,立刻就能发现。

4.3 步长和精度要求不同步,误差反而变大

做数值积分或常微分方程时,有人会觉得“步长取得越小越精确”,于是直接把步长设为 10⁻⁸。这个想法在理论上是成立的,但实际中会有两个问题:一是计算量成倍增加,程序可能跑很久;二是舍入误差会随计算次数累积,步长过小反而可能让误差上升。所以取值前需要做一个简单判断:题目给了精度要求,那步长应该选多大才够,但又不至于过度?比如复合梯形公式误差是 (O(h²)),要达到 10⁻⁶ 精度,h 大约需要 10⁻³ 量级,最多 10⁻⁴ 也就差不多了。设成 10⁻⁸ 只会让程序白白多算几百万次,还可能得到更差的结果。

4.4 用打印输出来调试太慢,学会在关键位置插桩

刚开始做上机题时,我会到处 print,把每个中间变量都打出来。这种方式排查简单问题够用,但程序复杂后,打印的信息量太多,反而让人眼花。后来我改用另一种方式:只打印抽象度更高的信息,比如每轮迭代的误差、每次消元的主元大小、每一步数值积分的结果。看到这些关键信息,我通常能快速定位是哪一步出了问题,然后再针对性打印那一步的中间变量。这种“先宏观定位,再微观排查”的思路,是我强烈推荐的调试习惯。

4.5 用库函数做交叉验证,但别直接抄答案

现在很多题目允许使用工具库,比如 Python 的 scipy、MATLAB 的内置函数。我的建议是:可以用内置函数来验证自己的答案,但报告里不要直接用内置函数的结果代替实现。上机题的目的本来就是让你手写一遍算法,理解其中的细节,直接调用封装好的函数等于把训练机会浪费掉了。我自己的做法是,手写算法之后,用库函数的结果作为“参考答案”做交叉验证。如果两者差别在合理误差范围内,说明实现正确;如果差别巨大,那大概率是我手写部分出了 bug,而不是库函数出问题。

最后说点实际的建议

如果你做上机题的时候觉得吃力,不要觉得是自己数学不行。更多情况下,是因为你从来没有把“算法流程”在脑子里走一遍,就开始写代码了。试着拿出一张白纸,先把算法的输入、输出、停机条件、每一步操作写成一个伪代码清单,再翻译成真正的代码。这个过程只需要十几分钟,但能帮你节省一个晚上的调试时间。

还有一个我个人的小经验:第一次写某一类算法时不要追求一次成功,先写一个能跑通小例子的版本,哪怕结果粗糙,然后逐步增加对精度、容错、效率的处理。客观地说,数值分析上机题是少数能让你同时感受到数学严谨性和工程味道的任务,做透了之后,你会发现自己以后看到其他数值仿真问题,心里也更有底。希望这篇内容能帮你顺利跑通作业,也真正把算法吃进肚子里。

RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进