ARTICLE · INTELLIGENCE

战地情报 · 详情页

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

牛顿迭代法实验深度解析:从原理推导、收敛阶分析到发散陷阱与实战代码

牛顿迭代法实验深度解析:从原理推导、收敛阶分析到发散陷阱与实战代码

做计算方法课程实验做到第四次,主题是牛顿迭代法。说实话,一开始我多少有点轻视这个题目——课堂上公式推得明明白白,x_{n+1} = x_n - f(x_n)/f'(x_n)一行就能写完,打开编辑器可能十分钟就出结果了。可真正动手做实验,把初值换着给、把不同形态的函数拿来试、把收敛过程逐条打印出来看之后,我才反应过来,这个看似简单的迭代格式里藏着的门道远比想象的多:初值稍微偏一点,序列可能直接飞到十万八千里外;碰到重根,原本引以为傲的二阶收敛速度瞬间退化成一阶;停止准则选不好,打印出来的“精确解”根本就不是你想要的根。这篇文章把整套实验从原理到实现、从异常现象到误差分析的过程完整记录一遍,给正在做计算方法或数值分析实验的同学一份能直接参考的实操记录,也适合只想把牛顿迭代法真正搞懂,而不是停留在背公式层面的读者。

1. 先别急着敲代码:牛顿迭代法要解决的是哪一类问题

1.1 现实中的方程很少长成 ax²+bx+c=0 的样子

我们从小接受的方程训练基本都长一个样:移项、合并同类项、套求根公式。二次方程有现成的判别式和公式,三次四次虽然公式长得很劝退,但至少理论上存在解析解。可是到了五次及以上的一般多项式方程,数学家已经证明不存在通用的根式解;更别提实际工程里冒出来的那些混着指数、三角函数和隐式关系的方程,比如 x·e^x = 1,或者天体轨道计算里常见的 E - e·sin(E) = M 这种方程,它们压根没有“精确表达式”这回事,你只能在给定的精度要求下求一个数值近似解。

举个例子,电路分析里二极管方程 I_s·(e^(V/V_T) - 1) - I = 0 就是一个典型的超越方程,V 在指数上,没法通过初等变换解出来。再比如化工计算里求气体压缩因子 z 的状态方程,往往是含 z 的三次方程,系数还带着温度和压力,手推公式根本不现实。这些场景的共同点是:方程本身没有解析根,或者有但代价高到不可接受,这时候迭代法就是默认选择。

这个“数值近似”的过程,思路都差不多:先给一个猜测值,然后按照某种规则把猜测值不断修正,直到前后两次的结果差得足够小,就认为收敛到了方程的根。牛顿迭代法之所以是计算方法课程里绕不开的经典,就是因为它是迭代法中收敛速度最快的那一档,而且几何意义非常直观,背后的原理还能很自然地推广到非线性方程组和最小二乘问题,学会了它就是给后续课程打底子。

1.2 “局部收敛”四个字决定了实验怎么做

不过牛顿法有一个和二分法截然不同的脾气:它是局部收敛方法。所谓局部收敛,意思是只有当初值落在某个根的“引力范围”内,迭代序列才会乖乖收敛到那个根;初值给得不好,序列可能收敛到另一个根,也可能压根不收敛。二分法只需要你给一个两端函数值异号的区间,就能保证缩出根来,属于全局收敛,代价是线性速度,慢;牛顿法快是快,但非常依赖初值。

更扎心的事实是,一个非线性方程往往有不止一个根。以 f(x) = x³ - 3x + 1 为例,它有三个实根,你从不同的初值出发,最终会落到不同的根上。如果你把牛顿法放到复平面上观察,不同根的吸引域边界甚至可能是分形结构——这不是玄学,而是牛顿迭代作为离散动力系统的真实性质。所以你会发现,几乎所有计算方法实验都会要求你对同一个方程换若干个不同的初值试跑一遍。这不是课程设计故意刁难你,而是想让你亲手体会“局部收敛”到底是什么感觉。

第一次做这个实验的时候,我就在初值选择上吃了不少苦头,后面专门用一节来写这些翻车现场。这里先记住一个结论:牛顿法的实验报告里,初值敏感性测试和收敛结果一样重要,它比任何一个漂亮的输出值都更能体现你对这个方法的理解程度。

2. 一条切线引出的迭代公式:原理推导与收敛速度分析

2.1 几何图像:切线穿过 x 轴的位置就是下一个猜测点

先把公式放一边,看图像。假设你在坐标平面上画了一条曲线 y = f(x),目标是找到曲线和 x 轴的交点。现在你手上有一个猜测点 x_n,对应曲线上一点 (x_n, f(x_n))。牛顿法的做法是:过这个点作曲线的切线,然后看这条切线与 x 轴相交在哪里,把交点横坐标当作下一次的猜测值 x_{n+1}。

这个过程非常像你用一把直尺去“卡”曲线:切线方向代表曲线在这个位置的变化趋势,切线越陡,说明函数在这一带变化越剧烈,猜测值的修正量越小;切线越平,说明函数接近水平,切线与 x 轴的交点就会跑到很远的地方去。这个几何直觉在后面理解“初值不好为什么会发散”时特别有用——一旦切线接近水平,一步就可能把你送到千里之外。

2.2 泰勒展开推导:从二阶截断得到迭代格式

几何上想通了,代数上一推就出来。设 x* 是方程的根,在 x_n 处把 f(x*) 做泰勒展开:

f(x*) = f(x_n) + f'(x_n)(x* - x_n) + (1/2)f''(ξ)(x* - x_n)²

左边 f(x*) = 0,如果丢掉最后那个二次余项,就得到一个线性近似:

0 ≈ f(x_n) + f'(x_n)(x* - x_n)

解出 x*,得到:

x* ≈ x_n - f(x_n) / f'(x_n)

右边这个近似值就是我们构造出来的下一个迭代点,于是牛顿迭代公式顺理成章:

x_{n+1} = x_n - f(x_n) / f'(x_n)

注意这里每一步都要重新计算 f(x_n) 和 f'(x_n),因为 x_n 一直在变,切线也一直在变。这跟最优化里每一步都要重新算梯度的道理是一样的:你不能拿着一个固定斜率走到底,那样很快就偏离目标了。

2.3 收敛阶的由来:为什么说它是二阶收敛

把上面丢掉的二次余项捡回来,就能分析误差。记 e_n = x_n - x* 为第 n 步的误差,重新整理泰勒展开式可以得到:

e_{n+1} ≈ ( f''(x*) / (2 f'(x*)) ) · e_n²

这个式子说明,牛顿法的误差是“平方级”缩小的。用具体数字感受一下:假如当前误差是 0.1,平方后变成 0.01,再下一步变成 0.0001,接着是 10^(-8)、10^(-16)。每一步的有效数字位数差不多翻一倍,这就是为什么对一个形态正常的方程,牛顿法往往四五步就能压到双精度浮点数的极限。

也正因为这个平方关系,牛顿法被称为二阶收敛方法。这里有个容易被忽略的细节:上面的结论成立的前提是 f'(x*) ≠ 0,也就是说根必须是单根。一旦碰到重根,分子的阶数变了,收敛阶会从 2 掉回 1,这个后文会专门展开讲,因为这是实验里最容易踩的坑之一。

3. 实验代码实现:Python 版牛顿迭代器与停止准则设计

3.1 导数怎么来:解析求导与数值微分的取舍

教科书里的例题总是默认 f'(x) 是现成的,手推一下就好。但真到实验和工程里,函数往往复杂到解析求导不现实,或者你拿到的就是一个只返回函数值的黑盒接口,这时候就得用数值微分。最常用的是中心差分公式:

f'(x) ≈ ( f(x + h) - f(x - h) ) / (2h)

这个公式的截断误差是 O(h²),比前向差分( f(x+h) - f(x) )/h的 O(h) 要好一档。步长 h 的选择有讲究:太大,截断误差占主导;太小,两个相近的大数相减会把舍入误差放大。在双精度下,经验上取 h ≈ 10^(-6) 左右比较平衡,理论上最优量级是机器精度的 1/3 次方,大约 6×10^(-6),所以取 1e-6 是合理的务实选择。

如果函数本身能手动求导,我还是建议用解析导数。原因很实际:实验的对照组要求验证收敛阶,解析导数能让你把“方法本身的误差”和“数值微分引入的额外误差”分开,报告写起来更有说服力。我在实验里两种都做了,先把解析导数版本跑通,再换成数值导数版本对比,这样既稳妥又多了一个可以讨论的维度。

3.2 停止准则:只看 |Δx| 或只看 |f(x)| 都会翻车

迭代总得有个停下来的判断。课程里常见两种准则:一种是看相邻两次迭代的差 |x_{n+1} - x_n| < tol,另一种是看残差 |f(x_{n+1})| < tol。很多同学偷懒只写其中一种,我在测试中遇到过两个典型的翻车案例。

只靠残差准则,会在函数很平坦时误判。比如 f(x) = x³,根在 0 附近,函数值本身小得非常快,x 还在 10^(-2) 量级时 |f(x)| 就已经到了 10^(-6),如果 tol 设成 1e-6,程序会提前宣布“收敛”,可实际上误差离机器精度还差着十万八千里。反过来,只靠 |Δx| 准则,在函数很陡峭的地方也可能出问题,因为函数陡意味着很小的 Δx 也能对应很大的残差。稳妥的做法是双条件:|Δx| 和 |f(x)| 都小于容差才停,同时把最大迭代次数加上,防止死循环。

另外还有一点必须提醒:容差不要设到 1e-18 这种低于机器精度的值。双精度浮点数的机器精度大约 2.2×10^(-16),当容差设得比它还小,循环条件永远无法满足,只会白白迭代到 max_iter 然后把异常抛出来。我第一次跑实验就犯了这么个想当然的错,调了半天才发现是容差写得比我以为的还要小。

3.3 一个可以直接跑的代码与运行输出

综合上面几点,我写了一个带历史记录的牛顿迭代函数。它返回每次迭代的近似根、步长和残差,方便后面画图和分析误差阶:

def newton(f, df, x0, tol=1e-14, max_iter=50): """牛顿迭代法求解 f(x)=0。 参数 ---- f : 目标函数 df : 导函数 x0 : 初始猜测 tol : 停止容差,同时检查 |dx| 与 |f(x)| max_iter : 最大迭代次数 返回 ---- (近似根, 每次迭代记录列表) """ x = x0 history = [] for k in range(1, max_iter + 1): fx = f(x) dfx = df(x) if abs(dfx) < 1e-14: print("导数接近零,算法提前终止") break x_new = x - fx / dfx dx = abs(x_new - x) residual = abs(f(x_new)) history.append((x_new, dx, residual)) print(f"k={k:2d} x={x_new:.15f} |dx|={dx:.3e} |f(x)|={residual:.3e}") if dx < tol and residual < tol: return x_new, history x = x_new return x, history

拿课程实验最经典的例子 f(x) = x³ - 2x - 5、初值 x0 = 2 来跑,输出大致是这样:

k= 1 x=2.100000000000000 |dx|=1.00e-01 |f(x)|=6.10e-02 k= 2 x=2.094568130000000 |dx|=5.43e-03 |f(x)|=1.86e-04 k= 3 x=2.094551500000000 |dx|=1.66e-05 |f(x)|=2.07e-07 k= 4 x=2.094551481542326 |dx|=1.85e-08 |f(x)|=1.78e-16

注意看第四步,残差直接从 10^(-7) 跳到了 10^(-16),这就是二阶收敛的直观体现:最后一步几乎是“一步到位”式地把误差压到了机器精度。这个例子常被用来做牛顿法的入门演示,因为它的导数在根附近既不接近零、函数形态也规整,收敛过程非常干净利落。

4. 实验里必须踩一遍的三个坑:发散、重根与震荡

做完上面那个漂亮的例子,如果就此收手写报告,说明实验做得还不够深。接下来这几个“反面教材”才是我觉得整个实验最有价值的部分,每一个都是实际跑出来的,不是书上看来的。

4.1 初值落在驻点附近:切线一平,迭代直接跑飞

牛顿迭代法翻车的第一种常见姿势,是初值落在驻点(导数接近零的极值点)附近。回到切线图像就很容易理解:切线越平,它与 x 轴的交点就被推到越远的地方,一步下来猜测值可能就跑到几十上百开外。

我用 f(x) = x³ - 3x + 1 试着从 x0 = 0.99 起步。这个函数在 x = 1 附近有一个极值点,f'(0.99) 只有 -0.0597 左右,而 f(0.99) 大约是 -0.9997,于是第一步修正量是 -0.9997 / -0.0597 ≈ 16.7,下一个点直接蹦到 -15.76。从一个接近极值点的地方出发,一眨眼就跑出去十几个单位,后续的迭代轨迹也变得不可控。更极端的例子是把初值精确取在极值点上,比如 x0 = 1,此时 f'(1) = 0,直接除零崩溃。

这个坑给实验带来的启示很直接:写代码时一定要判断导数绝对值是否接近零,否则程序会以除零或者巨大跳跃的方式失败。更重要的是,它解释了为什么实验总要你记录“初值取哪些值不收敛”——这不是为了凑字数,而是让你直观理解局部收敛的含义:初值不是随便给的,它直接决定了算法的生死。

4.2 重根附近的高阶收敛退化成线性

第二个坑比发散更隐蔽:函数有重根时,牛顿法并不是二阶收敛的。取 f(x) = (x - 1)³,根 x* = 1 是一个三重根,初值取 x0 = 0.9。按理说这个函数非常“好”,处处可导,根也明确,但从 0.9 出发时,误差序列是这样的:

  • e0 = -0.1000
  • e1 = -0.0667
  • e2 = -0.0444
  • e3 = -0.0296

每一步误差大约乘以 2/3,这是一个典型的线性收敛序列。算一算就明白为什么:对 m 重根,牛顿迭代式在根附近会退化成x_{n+1} ≈ x_n - (x_n - x*)/m,也就是误差乘以 (1 - 1/m)。m = 3 时就是每步缩小到原来的 2/3,误差从 10^(-1) 压到 10^(-15) 需要大约 90 步,和简单根场景下 4 步完成形成天壤之别。

知道原因之后,改进方法也顺手可得:把迭代格式改成x_{n+1} = x_n - m·f(x_n)/f'(x_n),也就是把修正量放大 m 倍,就能在重根附近恢复二阶收敛。前提是你得先知道重数 m,这在实际问题里往往又是一道坎。我在实验报告里把“三重根用原格式 vs 修正格式”的收敛过程画在了一起,一个直线下降、一个加速下坠,对比效果非常直观。

4.3 周期性震荡:算法陷入 0 和 1 之间反复横跳

第三个坑最反直觉:初值给得挺规矩,算法也能“稳定迭代”,但永远不收敛,因为序列陷入了周期循环。最经典的例子是 f(x) = x³ - 2x + 2,初值取 x0 = 0:

x0 = 0 时,f(0) = 2,f'(0) = -2,所以 x1 = 0 - 2/(-2) = 1;到了 x1 = 1 时,f(1) = 1,f'(1) = 1,于是 x2 = 1 - 1/1 = 0。序列在 0 和 1 之间无限循环,程序如果不设最大迭代次数,就会在这里空转到天荒地老。

这个例子我第一次跑的时候还以为是代码写错了,反复检查了半天才意识到算法本身就陷入了周期点。后来查了相关分析才知道,牛顿迭代法作为非线性动力学系统,存在周期轨道和分形边界,只是这个例子恰好给出了一个最简洁的二周期。对付这种情况没有万能解,但有两道防线一定能派上用场:一是 max_iter 必须设;二是当你发现近几步的取值在反复弹跳时,果断换初值,或者改用二分法兜底。这些防御性处理在实验报告里可以归到“算法鲁棒性讨论”,是实打实的加分内容。

5. 从“跑通”到“讲清楚”:误差阶验证、方法对比与报告写法

5.1 用实验数据反推收敛阶

实验不能只停留在“跑出了根”,还要用数据证明它确实二阶收敛。收敛阶的定义是:如果存在常数 C 和阶数 p,使得|e_{n+1}| / |e_n|^p → C,就称方法为 p 阶收敛。用实验数据验证 p 的经典做法是取相邻三步的误差,套对数公式:

p ≈ ln(|e_{n+1}| / |e_n|) / ln(|e_n| / |e_{n-1}|)

对前面 f(x) = x³ - 2x - 5 的例子,先用高精度结果 2.0945514815423265 作为近似真值,算出每步误差,再代入公式,得到的 p 会稳定在 2 附近。我跑出来的结果大约是 1.998、2.001 这样的数值,和理论完全吻合。这部分内容直接写进实验报告,比干巴巴地说一句“本方法二阶收敛”有说服力得多。

另外有个小技巧:把每步误差画成对数坐标图。二阶收敛的误差序列在对数图上表现为斜率逐级加倍、曲线急剧下坠;一阶收敛则是一条近似直线的均匀下降。两张图放在一起,收敛速度的差异一目了然,配图在实验报告里属于很提气的部分。

5.2 牛顿法、二分法、弦截法放在一张表里对比

计算方法课程通常会把牛顿法和二分法、弦截法放在同一轮实验里做对比。我的做法是拿同一个方程分别用三种方法求解,记录迭代次数和每步精度,然后整理成一张表:

方法收敛阶每次迭代开销初值要求特点
二分法1(线性)1 次函数求值给定变号区间全局收敛,稳但慢
牛顿法2(二阶)函数值 + 导数值初值靠近根最快,但可能发散
弦截法约 1.6181 次函数求值两个初始点次线性,不用求导

这张表值得展开理解。二分法是稳健型选手,只要区间两端异号就一定能找到根,代价是每步只缩小一半区间,精度要求高时迭代次数相当可观。牛顿法是速度型选手,一次迭代要花两次函数求值的开销(函数值和导数值),换来的是平方级收敛,但前提是初值陪跑到位。弦截法则是不想求导时的折中方案,用前两步的函数值做差商替代导数,收敛阶约 1.618,介于两者之间。

实验里有个很常见的现象值得留意:对同一个问题,二分法可能跑了五六十步,牛顿法四步就结束了,于是大家本能地觉得牛顿法全面碾压。但换个角度想,牛顿法每一步都要算导数,而且如果初值没给好,前面花费的迭代可能全部白费甚至反方向跑飞。所以“哪个方法更好”这种问题没有标准答案,取决于你对函数信息掌握多少、对初值有多大的把握、以及对计算开销的容忍度。实验报告里能写出这种取舍,就说明真的想明白了。

5.3 三个值得写进实验心得的经验

最后把我自己做完这轮实验沉淀下来的三条心得写在这里,都是可以直接用在后续课程和实际项目里的。

第一条,导数能用解析就给解析,给不了再用数值微分,并且用中心差分而不是前向差分,步长控制在 1e-6 量级。这样能最大程度降低“方法本身的误差”被“导数误差”污染的风险,方便你做收敛阶验证;如果你在报告里同时给出解析导数和数值导数的两组结果对比,还能额外体现你对误差来源的敏感度。

第二条,停止准则永远设双条件,|Δx| 和 |f(x)| 都判断,同时加上容差下界和最大迭代次数两个保护。我在测试里见到的绝大多数“不收敛”或“假收敛”,都不是算法本身的问题,而是停止条件写得太粗糙;把这条改好,很多莫名其妙的 bug 都会自动消失。

第三条,遇到重根和震荡,先别急着怀疑代码,回去看一眼函数形态和初值位置。重根时收敛慢不是 bug,震荡死循环是牛顿法作为非线性映射的正常动力学行为。把这些异常现象理解成方法本身的性质,实验报告的讨论深度就完全不一样了。

实验做完之后,我把这套带历史记录的newton函数存成了模板,后面写课程项目时遇到非线性方程,第一反应就是把它调出来改两行参数直接复用。换初值、换函数、换容差都只需要动一行代码,这种顺手的感觉让我觉得,这轮实验真正值回票价的地方,不是那个收敛漂亮的经典例题,而是把牛顿法从一个公式,变成了一个可以随时调用的趁手工具。

RELATED READING

延伸阅读

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