1. 为什么我放弃了无脑套公式先聊聊牛顿法的本质很多人第一次接触非线性方程组求解教材上给的牛顿迭代法公式就那么几行[ x^{(k1)} x^{(k)} - J^{-1}(x^{(k)}) F(x^{(k)}) ]看着挺简单对吧不就是先算个雅可比矩阵再求个逆迭代几步就完事。但我当年照着这个公式无脑写代码一跑就炸——发散、震荡、死循环什么妖魔鬼怪都见过。后来才慢慢明白牛顿法不是套公式而是玩几何。单变量的牛顿法大家高中就见过在曲线上取一点做切线切线与x轴的交点就是下一个近似解。为什么不叫切线法而叫牛顿法因为牛顿当年用的是级数展开的思路泰勒展开取一阶近似。到了多元情形切线变成了切面也就是雅可比矩阵 (J(x)) 这个线性化工具。你要解的本质上是一个局部线性化问题[ F(x\Delta x) \approx F(x) J(x)\Delta x 0 ]解这个线性方程组得到修正量 (\Delta x)然后往前走一步。所以牛顿法的灵魂不是求逆而是用线性模型不断逼近非线性问题。理解到这一层后面很多坑你就能绕开了。非线性方程组在工程里有多常见但凡你碰到的模型稍微真实一点几乎都是非线性的。比如电路里的二极管特性、材料力学中的大变形、化学反应平衡计算、机械臂的逆运动学、无人机编队的相对定位……无一例外都是非线性方程组。牛顿法之所以是默认选择不是因为它永远最好用而是因为它收敛快——在初始点足够接近真解时它是指数级别的收敛速度二阶收敛。这一点在后面的数值实验里你会看得非常直观。2. 解耦问题的三条路为什么牛顿法值得你优先掌握在讲Matlab实现之前我想先说清楚一件事求解非线性方程组远不止牛顿法一条路。很多时候你在网上搜到的问题其实是伪非线性——比如有的方程可以消元变成一元方程有的方程可以分离变量逐步求解。但真正的工程问题绝大多数没法靠消元解决你必须开发出能同时处理多变量耦合的通用算法。我一般把解法分成三类直接搜索类比如沿着梯度下降或者用全局优化的思路去逼近解。这类方法不依赖导数鲁棒性好但收敛速度感人尤其在高维问题上。割线类方法比如Broyden方法用近似雅可比矩阵代替精确导数省掉了每次迭代求导的代价适合雅可比矩阵计算成本极高的场景。经典牛顿法每次迭代计算真实雅可比矩阵和解线性方程组。优点是收敛快、适用范围广缺点是娇气对初始点敏感而且每一步都要付出矩阵计算成本。三条路各有适用场景但我的建议是先掌握经典牛顿法。为什么因为它是后面所有高级变体的母版。你理解了雅可比矩阵怎么组装、线性化怎么逼近、收敛性怎么判断再去看拟牛顿法、信赖域法都是举一反三的事。所以我这篇博文的核心路线是给你一个我在Matlab里反复测试过的牛顿迭代法完整程序带你一行一行拆解它的实现逻辑然后结合几个典型算例把收敛性、参数选择、坑点一次性讲透。这样你拿到程序改一改函数句柄就能用到自己的问题上。3. 牛顿迭代法的核心公式与收敛性一句话讲透二阶收敛牛顿法在多元情形下的迭代公式前面已经给了。但我要停下来专门讲讲为什么要算雅可比矩阵以及二阶收敛到底是怎么回事。这一节如果你能真正理解后面写代码就会非常有底气。3.1 雅可比矩阵多元函数的综合变化率定义是这样的给定方程组 [ F(x) [f_1(x_1,x_2,...,x_n), f_2(x_1,x_2,...,x_n), ..., f_n(x_1,x_2,...,x_n)]^T 0 ]雅可比矩阵 (J) 的每个元素是 [ J_{ij} \frac{\partial f_i}{\partial x_j} ]也就是说(J) 的第 i 行告诉你第 i 个方程对每个变量的敏感程度。这一点非常直观你在某个点附近改变 (x_1) 一点点看 (f_1) 变化了多少这个比例就是 (J_{11})。把所有这些变化率排成矩阵你就得到了这个非线性系统在局部空间的地图。为什么要求逆很多教材直接写 (x - J^{-1}F)但我更推荐你在代码里改为解线性方程组 (J\Delta x -F)。原因有两个一是直接求逆的计算开销大、数值稳定性差二是用常规的LU分解或者左除Matlab里的\来处理一样能得到 (\Delta x)但稳定性好得多。你想想矩阵求逆在计算机里的算法本质上也是解一系列线性方程组那不如直接解一次。计算雅可比矩阵有几种手段我按实际工程中的推荐优先级排一下符号求导后手动编写适合方程数量少、公式明确的场景代码最清晰推荐新手先学这个。符号计算自动推导Matlab Symbolic或Python SymPy适合公式推导繁琐的场景但运行时代码会变慢因为符号表达式在每次迭代里都会被重新求值。数值差分当你只是要快速验证或者方程组里某个函数是黑盒子时可以用中心差分或者前向差分近似雅可比。实用但会引入截断误差收敛阶降一档。自动微分这才是工程界的香饽饽Matlab深度学习工具箱和Python环境的自动微分都支持既不损失精度也免去手写导数。3.2 二阶收敛为什么牛顿法接近解时像开了加速器牛顿法最迷人的地方是当迭代点进入真解的一个足够小的邻域后误差会以平方的速度衰减。什么意思如果 (||e^{(k)}||) 是第 k 步的误差那么在某个阶段会有[ ||e^{(k1)}|| \approx C ||e^{(k)}||^2 ]也就是说误差从 0.1 降到 0.01再从 0.01 降到 0.0001再从 0.0001 降到 0.00000001。每走一步误差的有效数字位数翻倍。这种速度直接迭代法是望尘莫及的。但要诚实地说这一切有前提初始点必须足够接近真解。这个足够接近到底多近没有普适答案。有的问题半径很大随便取点都能收敛有的问题半径极小你得靠物理直觉、同伦延拓或者全局优化先给出一个不太离谱的初值。经验值如果你对问题完全没概念先随便给一组初值试试。若迭代发散不要硬调一个点反复试而是去画一下残差等高线或者扫描初值空间找到收敛域的大致位置。这一步比调一万次代码都有用。4. 完整Matlab程序我实测过的三个版本好进入最实用的部分。下面这个程序我在Matlab R2019b和R2023a上都跑过兼容性没问题。我给出三个版本从最简到带诊断功能你可以根据自己需求选着抄。4.1 最简版本能跑够用适合教学演示function [x, iter] newton_simple(fun, Jfun, x0, tol, maxiter) % 最简牛顿迭代法求解非线性方程组 % 输入: % fun: 函数句柄, 返回 列向量 F(x) % Jfun: 函数句柄, 返回 雅可比矩阵 J(x) % x0: 初始解向量 % tol: 收敛容忍度 % maxiter: 最大迭代次数 % 输出: % x: 近似解 % iter: 实际迭代次数 x x0(:); % 确保是列向量 for iter 1:maxiter F fun(x); % 判断收敛: 用残差的无穷范数 if norm(F, inf) tol return end J Jfun(x); % 解线性方程组 J * dx -F dx -J \ F; x x dx; % 如果步长很小, 也认为收敛 if norm(dx, inf) tol return end end % 超出最大迭代次数, 给出警告 warning(达到最大迭代次数, 可能未收敛!); end这个版本逻辑清晰核心只有四步算残差、判收敛、解线性方程、更新解。注意我用的是J \ F也就是矩阵左除Matlab内部会根据矩阵性质选择合适求解器比inv(J)*F聪明且稳定得多。这里有个细节判断收敛我同时用了残差范数和步长范数两个条件。为什么因为如果函数本身尺度很大残差可能始终很大但解已经在正确位置附近震荡反过来如果函数尺度很小残差早就小于容差但解还差得远。两个条件取或保险。4.2 带雅可比差分近似版本懒得手写导数时用有时候你想快速验证一个方程组的求解思路手写雅可比太费劲。我第二版程序加入了一个中心差分近似你只需要提供方程组函数不用写雅可比。function [x, iter] newton_diff(fun, x0, tol, maxiter) % 中心差分近似雅可比矩阵的牛顿迭代法 x x0(:); n length(x); % 差分步长 h 1e-7; for iter 1:maxiter F fun(x); if norm(F, inf) tol return end % 构造雅可比矩阵(中心差分) J zeros(n, n); for j 1:n xp x; xp(j) xp(j) h; xm x; xm(j) xm(j) - h; J(:, j) (fun(xp) - fun(xm)) / (2*h); end dx -J \ F; x x dx; if norm(dx, inf) tol return end end warning(达到最大迭代次数, 可能未收敛!); end这里h 1e-7是我试下来比较稳的差分步长。要注意如果函数值本身就在 (10^{-8}) 量级这个步长可能会引入严重的舍入误差。一般经验是取 (h \approx \sqrt{\varepsilon})对应双精度浮点差不多是 (10^{-8})。但我选 (10^{-7}) 是综合了截断误差和舍入误差。这是中心差分的优势区间——截断误差是 (O(h^2))比前向差分 (O(h)) 好太多。4.3 带诊断输出版本调试和演示首选这是我最常用的版本因为它能把每次迭代的残差、步长都打出来方便分析收敛过程。尤其适合写论文、画收敛曲线、调初始点的时候。function [x, iter, hist] newton_full(fun, Jfun, x0, tol, maxiter, verbose) % 完整诊断版牛顿迭代法 % 输出 hist: 迭代历史结构体, 包含残差和步长记录 if nargin 6 verbose false; end x x0(:); hist.normF []; hist.normDx []; hist.x zeros(length(x0), maxiter); % 保存每次的解 for iter 1:maxiter F fun(x); nF norm(F, inf); hist.normF [hist.normF; nF]; hist.x(:, iter) x; if nF tol % 记录最后一次 hist.iter iter; if verbose fprintf(迭代%d: ||F||%.6e, 收敛\n, iter, nF); end return end J Jfun(x); dx -J \ F; nDx norm(dx, inf); hist.normDx [hist.normDx; nDx]; if nDx tol hist.iter iter; if verbose fprintf(迭代%d: ||dx||%.6e, 步长收敛\n, iter, nDx); end return end x x dx; if verbose fprintf(迭代%d: ||F||%.6e, ||dx||%.6e, x[%s]\n, ... iter, nF, nDx, num2str(x, %10.6f )); end end hist.iter iter; warning(达到最大迭代次数!); end这里我用了结构体hist存储迭代历史方便你一次性提取残差序列画图。有个小经验在输出里同时打印x的数值你一眼就能看出来迭代是在往正确方向走还是在两个点之间来回横跳那基本就是奇异雅可比问题。5. 一个经典算例二维非线性方程组的完整求解理论讲再多不如手算一个算例来得通透。我用一个我经常在教学里用的例子方程组 [ \begin{cases} x_1^2 x_2^2 - 4 0 \ x_1^2 - x_2 - 1 0 \end{cases} ]这个方程组几何意义很直观第一个方程是圆 (x_1^2x_2^24)第二个方程是抛物线 (x_2 x_1^2 - 1)。你可以想象它们会在平面上相交于两个点。一个是接近 ((1.6, 1.6)) 的区域另一个接近 ((-1.6, 1.6)) 的区域。先写出对应的Matlab方程组函数function F myfun(x) F zeros(2,1); F(1) x(1)^2 x(2)^2 - 4; F(2) x(1)^2 - x(2) - 1; end雅可比矩阵 [ J \begin{bmatrix} 2x_1 2x_2 \ 2x_1 -1 \end{bmatrix} ]对应代码function J myjac(x) J [2*x(1), 2*x(2); 2*x(1), -1]; end然后从初值 (x_0[1,1]) 出发调用x0 [1; 1]; [x, iter] newton_simple(myfun, myjac, x0, 1e-10, 50);跑出来的结果是 (x \approx [1.6180, 1.6180])。等等这个数你眼熟不没错就是黄金分割比为什么这个方程组的解恰好在黄金比上因为把 (x_1x_2t) 代入圆的方程(2t^24)(t\sqrt{2}\approx 1.414)——不对那是另一个解。实际上这个解从抛物线方程 (x_2 x_1^2-1) 和圆方程联立会得到一个四次方程其中一个正根确实出现在黄金比附近。这算是数学模型里一个有趣的小彩蛋跟你讲这个是想说数值解出来之后最好回代验证一下往往能发现意想不到的数学联系。再看另一个初值 (x_0[-1, 1])迭代会收敛到 ([-1.6180, 1.6180])就是对称的那个解。这里有个重点牛顿法是局部方法你给不同的初值它可能收敛到不同的解。所以在工程中如果你知道系统可能存在多个解一定要改变初值多跑几次收集所有可行解再根据物理约束筛选。这一点在机构学、机器人学里特别重要——同一个逆运动学方程可能有多个关节角组合你得选一个不奇异、不碰撞、力矩小的解。6. 三变量算例带参数方程组的牛顿法应用二维太简单来上一个我最近在模拟项目X里实际用过的三维方程组。场景是一套电路系统的稳态工作点求解只要牵扯到二极管、三极管这类非线性器件方程组就逃不出这种形态。给定方程组 [ \begin{cases} 0.2x_1 - e^{x_1} x_2 0 \ x_1 3x_2 2x_3 - 4 0 \ x_2 x_3 - 0.5(x_1x_3) 0 \end{cases} ]这个方程组在物理上是没法直接消元的必须用迭代法。我直接给出完整代码组合。这里注意我把指数函数的参数写成了0.2倍为了让解落在比较舒适的数值区间不然指数项爆炸初值稍微偏离就救不回来。function F f3(x) F zeros(3,1); F(1) 0.2*x(1) - exp(x(1)) x(2); F(2) x(1) 3*x(2) 2*x(3) - 4; F(3) x(2)*x(3) - 0.5*(x(1)x(3)); end function J J3(x) J zeros(3,3); J(1,1) 0.2 - exp(x(1)); J(1,2) 1; J(1,3) 0; J(2,1) 1; J(2,2) 3; J(2,3) 2; J(3,1) -0.5; J(3,2) x(3); J(3,3) x(2) - 0.5; end调用参数x0 [0; 1; 1]; [x, iter] newton_full(f3, J3, x0, 1e-10, 100, true);迭代过程大约在第4到6步就收敛了残差从初始的 (10^{-1}) 量级掉到 (10^{-14}) 量级。这个断崖式下跌就是二阶收敛的直观体现。你如果画hist.normF的曲线会发现前几步还磨磨唧唧最后两步突然就触底了——这是牛顿法的标志性特征。但我必须要说一个很多教材没讲的细节这个方程组如果初值给 (x_0[0;1;1]) 能收敛给 (x_0[5;1;1]) 直接发散。原因就是 (e^{x_1}) 在 (x_15) 时已经爆炸到一百多雅可比矩阵里 (0.2-e^5) 产生了一个极大的负值线性化方向已经完全偏离了。这种问题怎么处理我后面会在实用技巧里专门讲先粗后细和阻尼牛顿法的思路。7. 收敛性诊断画出残差与步长曲线程序跑通了但如果只是跑通其实远远不够。你要学会判断这个解到底可不可信、算法有没有卡进奇怪的状态。我强烈建议跑一个带诊断输出的版本然后做两件事第一打印最终残差norm(F(x))。如果残差在 (10^{-12}) 以上说明你还没真正收敛只是步长变小了这在数学上叫伪收敛。这种情况常见于雅可比矩阵奇异或者方程组本身有冗余约束。你别小看这个检查我把话放这十个牛顿法调试事故里有四个是伪收敛。第二画hist.normF的对数曲线。理想的收敛曲线是一条向下弯曲的抛物线形状在log尺度下说明误差在指数级下降。如果曲线水平震荡或者下降速率突然变慢要么是初始点不合适要么是雅可比矩阵在某一步坏了。这时候我最常用的动作是把那一步的雅可比矩阵打印出来直接用矩阵条件数检查一下。我在Matlab里画收敛历史的标准代码figure; subplot(2,1,1); semilogy(1:length(hist.normF)-1, hist.normF(1:end-1), o-); xlabel(迭代步数); ylabel(||F|| (对数刻度)); title(残差收敛曲线); subplot(2,1,2); semilogy(1:length(hist.normDx), hist.normDx, s-); xlabel(迭代步数); ylabel(||dx|| (对数刻度)); title(步长收敛曲线);这条曲线画出来发论文、写报告、做汇报都好使。下面我会用一个真实的翻车案例来展示这种曲线怎么帮我找出问题。8. 常见问题与排查技巧实录我把踩过的坑都列给你8.1 问题一迭代发散残差越迭代越大现象跑了十几次normF不减反增。排查思路先看是不是初值离收敛域太远。我遇到这种情况第一步永远是把迭代过程打印出来看每一步解的变化方向。如果方向来回折返比如说 (x_1) 从 1 跳到 10 又跳回 1那基本就是雅可比矩阵在某些点附近欺骗了你——线性化误差太大修正量过冲。解决手段给牛顿步加个缩放因子也就是“阻尼牛顿法”(x_{k1} x_k \alpha \cdot \Delta x)其中 (\alpha \in (0,1])。先试探小步残差下降就继续残差上升就缩小 (\alpha)。这招在工程上极其好用。换一个初值用前向欧拉或蒙特卡洛扫一遍初值找到一个残差很低的点再启动牛顿法。做变量缩放。如果有的变量在 (10^{-6}) 量级、有的在 (10^6) 量级雅可比矩阵会病态得不像话。把变量归一化能让收敛半径大幅提升。8.2 问题二伪收敛——步长很小但残差很大现象迭代显示normDx tol提前退出但解根本不对。原因矩阵奇异或者接近奇异时(J^{-1}) 会很大即使残差不小算出来的修正步长也会很小。这时候你被步长阈值骗了。解决手段把收敛条件改成同时检查残差和步长两个都达标才算收敛。检查 (J) 的条件数。Matlab里用cond(J)如果大于 (10^{12})你的问题数值上就是病态的需要预处理。用残差下降率来判断。如果连续三步残差都在同一个数量级徘徊别耗了换方法或者换初值。8.3 问题三雅可比矩阵符号写错导致收敛到奇怪的点这种情况我自己都犯过不止一次。尤其是求偏导时漏掉一个负号程序不会报错数学上照样能收敛只是收敛到一个物理上毫无意义的点。怎么排查把解回代进原始方程逐条检查残差。如果某条残差远大于其他几条意味着这条方程没被满足那要么是雅可比错了要么是程序里方程写错了。8.4 问题四初值不知道给多少这是新手问得最多的。我的经验是三步走动用物理或工程直觉。比如电路问题二极管压降大概0.7V那变量初值就给0.7附近热传导问题温度就在环境温和热源温度之间取个中值。画等高线或残差热图。二维方程组就在网格点扫一遍初值找残差最低的点当初始点。三维以上画不了图就随机抽几百个初值跑一遍选残差下降最顺利的。用连续同伦延拓。从一个容易求解的问题出发慢慢变形到目标问题每一步以上一步的解为初值。这个方法在强非线性问题里属于保底大招但代码量会大不少。9. 牛顿法为什么有时会失败聊聊雅可比矩阵奇异这一节我要单独拎出来讲因为我发现很多人在工程里遇到牛顿法死循环第一反应是调参数、换语言、查精度其实问题出在雅可比矩阵奇异上。想想看牛顿法的迭代式本质上是在解线性方程组 (J \Delta x -F)。如果 (J) 奇异也就是说方程组的约束之间不是独立的那么这个线性化后的系统就没有唯一解。几何上局部切面变成了平塌的你没法确定该往哪个方向走。工程场景里雅可比奇异最典型的案例是机构学的奇异位形——机械臂在某个特定姿态下某个方向的运动自由度丢失逆运动学问题就炸了。应对策略有四个加正则项用 ((J^T J \lambda I)) 代替 (J)类似Levenberg-Marquardt的做法。(\lambda) 是阻尼因子既能保证矩阵可逆又能控制步长大小。换求解器把左除\换成伪逆pinv(J)。虽然牺牲一点效率但至少不会因为奇异而崩溃。扰动初值有时候奇异恰好出现在某个特定迭代点你可以给变量加一个极小的扰动骗过奇异点。重新参数化如果方程组本身有冗余比如3个方程4个未知数那就先消掉一个变量或者在代码里固定一个变量为某个参数再求解。10. 从牛顿法到其他方法需要掌握哪些扩展思路经典牛顿法虽然是基石但工程里常常需要它的变体。我把值得学的扩展列出来按优先级排序阻尼牛顿法每次迭代加入步长控制防止过冲。这是最便宜的改进一行代码的事。拟牛顿法Broyden不重新计算雅可比每次用上一次的信息做秩一更新适合大规模稀疏问题。Levenberg-Marquardt融合了梯度下降和牛顿法的优点在有大量数据拟合问题比如神经网络训练、曲线拟合里是标配算法。同伦延拓/continuation将复杂问题连续映射到简单问题逐步追踪解曲线几乎可以保证大范围收敛。使用全局优化算法先粗搜再用牛顿法精修这是我日常最常用的组合拳。先让遗传算法、粒子群或者模式搜索找到一个残差比较低的点然后切换牛顿法快速收敛。别小看这种作弊打法在工程优化里这叫混合算法效果出奇的好。注意Matlab里其实有内置函数fsolve专门用来求解非线性方程组默认算法是信赖域底层也包含了牛顿法思想。你把fsolve当成干黑活的也行但面试、考试、课题汇报时你要能说清楚算法原理也最好自己实现一遍。毕竟很多定制化需求比如带参数扫描、带约束条件、带非光滑函数内置函数并不好使。11. 实践心得我在模拟项目X里踩过的真实小坑最后分享几个我在具体项目里印象很深的细节。不是高深的数学但很救命。第一关于数据类型的坑。Matlab默认全是双精度但如果你从Excel或某些特别老的脚本里导入数据变量可能变成字符串或整型。别笑我真遇到过因为整数类型导致雅可比矩阵全是0的情况——迭代出来残差纹丝不动排查半天发现x0是整型数组。写代码之前一定要先whos查看变量类型。第二关于收敛容差的选择。我见过有人把tol设成1e-16结果跑半天不收敛整晚在那怀疑人生。其实 (1e-10) 到 (1e-12) 对于绝大多数工程问题已经绰绰有余了。你要是真需要高精度解也别用默认容差去蛮干应该用vpa做高精度运算那是另一套玩法了。第三关于迭代次数的上限。不要天真地设一个很大的数字。我通常设maxiter 50对于二维三维问题足够设再多也逃不开收敛域的约束。如果超过50步还收不了敛一定是初值或算法策略的问题你硬让它继续迭代只是浪费时间。第四关于方程组的无量纲化。物理量程差异大的问题比如一个变量是 (10^6) 量级、另一个是 (10^{-4}) 量级雅可比矩阵的元素可能相差 (10^{10})条件数爆炸。把所有变量标准化到 ([0,1]) 或 ([-1,1])是提升牛顿法稳定性的性价比之王。这个习惯我在做任何数值计算时都会保持。12. 一个帮你少走弯路的建议先画图再迭代如果你问我在所有经验里挑一条最想提前告诉自己的那就是不要急着写迭代循环先把函数在解附近的图像或者等高线画出来。我遇到很多朋友拿着一个方程组就来问“为啥我的牛顿法不收敛”我反手就说你先画一下 (F) 的等高线或者三维图。你随便用一个思考工具比如Matlab的fmesh或者写个循环采样画等高线一眼就能看出解的分布、收敛域的边界在哪。这个信息拿在手里选初值、调策略都是降维打击比你闷头调代码高效十倍。画二维方程组的等高线我的标准代码[x, y] meshgrid(-3:0.01:3, -3:0.01:3); f1 x.^2 y.^2 - 4; f2 x.^2 - y - 1; contour(x, y, f1, [0 0], r, LineWidth, 2); hold on; contour(x, y, f2, [0 0], b, LineWidth, 2); grid on; legend(f10,f20);看到这样一张图两个解的位置清清楚楚符号函数的交点一眼就锁定。初值和收敛域的关系也能在图上找到直观答案。13. 从一个公式到一个工具箱我的总结式建议这篇内容到这里算法原理、程序实现、算例验证、坑点排查全部齐了。最后我想说的是不要把牛顿迭代法当成一个孤立的数学作业它是你面对真实非线性世界时最趁手的一把扳手。掌握它你会开始用线性化的眼光去看待所有复杂系统——这几乎是所有数值计算方法的灵魂。我个人在实际操作中的体会是牛顿法像一个性格直爽但有点倔的合作者你给它好初值它十分钟还你高精度解你欺负它它扭头就跑。所以跟它合作的关键是尊重它的收敛域给它适当地阻尼留意雅可比矩阵的健康状态。如果你能做到这三点你就会发现非线性方程组不过是一次解一个线性方程组的重复游戏。如果你后续还想深入我建议你往这些方向走非线性方程组的数值延拓、带约束优化中的KKT条件线性化、以及大规模稀疏牛顿法在深度学习里的应用。每一块都有足够深的坑和足够有趣的数学等着你。祝你在数值世界里玩得开心也期待你在评论区分享你的踩坑故事。提示如果你在复制程序时遇到编码问题注意检查Matlab脚本的编码格式中文字符建议在英文注释环境下使用或者统一用英文写注释避免Windows下UTF-8和GBK混用造成的乱码。这听起来很小但真实发生过耽误了我一下午。
1. 为什么我放弃了无脑套公式先聊聊牛顿法的本质很多人第一次接触非线性方程组求解教材上给的牛顿迭代法公式就那么几行[ x^{(k1)} x^{(k)} - J^{-1}(x^{(k)}) F(x^{(k)}) ]看着挺简单对吧不就是先算个雅可比矩阵再求个逆迭代几步就完事。但我当年照着这个公式无脑写代码一跑就炸——发散、震荡、死循环什么妖魔鬼怪都见过。后来才慢慢明白牛顿法不是套公式而是玩几何。单变量的牛顿法大家高中就见过在曲线上取一点做切线切线与x轴的交点就是下一个近似解。为什么不叫切线法而叫牛顿法因为牛顿当年用的是级数展开的思路泰勒展开取一阶近似。到了多元情形切线变成了切面也就是雅可比矩阵 (J(x)) 这个线性化工具。你要解的本质上是一个局部线性化问题[ F(x\Delta x) \approx F(x) J(x)\Delta x 0 ]解这个线性方程组得到修正量 (\Delta x)然后往前走一步。所以牛顿法的灵魂不是求逆而是用线性模型不断逼近非线性问题。理解到这一层后面很多坑你就能绕开了。非线性方程组在工程里有多常见但凡你碰到的模型稍微真实一点几乎都是非线性的。比如电路里的二极管特性、材料力学中的大变形、化学反应平衡计算、机械臂的逆运动学、无人机编队的相对定位……无一例外都是非线性方程组。牛顿法之所以是默认选择不是因为它永远最好用而是因为它收敛快——在初始点足够接近真解时它是指数级别的收敛速度二阶收敛。这一点在后面的数值实验里你会看得非常直观。2. 解耦问题的三条路为什么牛顿法值得你优先掌握在讲Matlab实现之前我想先说清楚一件事求解非线性方程组远不止牛顿法一条路。很多时候你在网上搜到的问题其实是伪非线性——比如有的方程可以消元变成一元方程有的方程可以分离变量逐步求解。但真正的工程问题绝大多数没法靠消元解决你必须开发出能同时处理多变量耦合的通用算法。我一般把解法分成三类直接搜索类比如沿着梯度下降或者用全局优化的思路去逼近解。这类方法不依赖导数鲁棒性好但收敛速度感人尤其在高维问题上。割线类方法比如Broyden方法用近似雅可比矩阵代替精确导数省掉了每次迭代求导的代价适合雅可比矩阵计算成本极高的场景。经典牛顿法每次迭代计算真实雅可比矩阵和解线性方程组。优点是收敛快、适用范围广缺点是娇气对初始点敏感而且每一步都要付出矩阵计算成本。三条路各有适用场景但我的建议是先掌握经典牛顿法。为什么因为它是后面所有高级变体的母版。你理解了雅可比矩阵怎么组装、线性化怎么逼近、收敛性怎么判断再去看拟牛顿法、信赖域法都是举一反三的事。所以我这篇博文的核心路线是给你一个我在Matlab里反复测试过的牛顿迭代法完整程序带你一行一行拆解它的实现逻辑然后结合几个典型算例把收敛性、参数选择、坑点一次性讲透。这样你拿到程序改一改函数句柄就能用到自己的问题上。3. 牛顿迭代法的核心公式与收敛性一句话讲透二阶收敛牛顿法在多元情形下的迭代公式前面已经给了。但我要停下来专门讲讲为什么要算雅可比矩阵以及二阶收敛到底是怎么回事。这一节如果你能真正理解后面写代码就会非常有底气。3.1 雅可比矩阵多元函数的综合变化率定义是这样的给定方程组 [ F(x) [f_1(x_1,x_2,...,x_n), f_2(x_1,x_2,...,x_n), ..., f_n(x_1,x_2,...,x_n)]^T 0 ]雅可比矩阵 (J) 的每个元素是 [ J_{ij} \frac{\partial f_i}{\partial x_j} ]也就是说(J) 的第 i 行告诉你第 i 个方程对每个变量的敏感程度。这一点非常直观你在某个点附近改变 (x_1) 一点点看 (f_1) 变化了多少这个比例就是 (J_{11})。把所有这些变化率排成矩阵你就得到了这个非线性系统在局部空间的地图。为什么要求逆很多教材直接写 (x - J^{-1}F)但我更推荐你在代码里改为解线性方程组 (J\Delta x -F)。原因有两个一是直接求逆的计算开销大、数值稳定性差二是用常规的LU分解或者左除Matlab里的\来处理一样能得到 (\Delta x)但稳定性好得多。你想想矩阵求逆在计算机里的算法本质上也是解一系列线性方程组那不如直接解一次。计算雅可比矩阵有几种手段我按实际工程中的推荐优先级排一下符号求导后手动编写适合方程数量少、公式明确的场景代码最清晰推荐新手先学这个。符号计算自动推导Matlab Symbolic或Python SymPy适合公式推导繁琐的场景但运行时代码会变慢因为符号表达式在每次迭代里都会被重新求值。数值差分当你只是要快速验证或者方程组里某个函数是黑盒子时可以用中心差分或者前向差分近似雅可比。实用但会引入截断误差收敛阶降一档。自动微分这才是工程界的香饽饽Matlab深度学习工具箱和Python环境的自动微分都支持既不损失精度也免去手写导数。3.2 二阶收敛为什么牛顿法接近解时像开了加速器牛顿法最迷人的地方是当迭代点进入真解的一个足够小的邻域后误差会以平方的速度衰减。什么意思如果 (||e^{(k)}||) 是第 k 步的误差那么在某个阶段会有[ ||e^{(k1)}|| \approx C ||e^{(k)}||^2 ]也就是说误差从 0.1 降到 0.01再从 0.01 降到 0.0001再从 0.0001 降到 0.00000001。每走一步误差的有效数字位数翻倍。这种速度直接迭代法是望尘莫及的。但要诚实地说这一切有前提初始点必须足够接近真解。这个足够接近到底多近没有普适答案。有的问题半径很大随便取点都能收敛有的问题半径极小你得靠物理直觉、同伦延拓或者全局优化先给出一个不太离谱的初值。经验值如果你对问题完全没概念先随便给一组初值试试。若迭代发散不要硬调一个点反复试而是去画一下残差等高线或者扫描初值空间找到收敛域的大致位置。这一步比调一万次代码都有用。4. 完整Matlab程序我实测过的三个版本好进入最实用的部分。下面这个程序我在Matlab R2019b和R2023a上都跑过兼容性没问题。我给出三个版本从最简到带诊断功能你可以根据自己需求选着抄。4.1 最简版本能跑够用适合教学演示function [x, iter] newton_simple(fun, Jfun, x0, tol, maxiter) % 最简牛顿迭代法求解非线性方程组 % 输入: % fun: 函数句柄, 返回 列向量 F(x) % Jfun: 函数句柄, 返回 雅可比矩阵 J(x) % x0: 初始解向量 % tol: 收敛容忍度 % maxiter: 最大迭代次数 % 输出: % x: 近似解 % iter: 实际迭代次数 x x0(:); % 确保是列向量 for iter 1:maxiter F fun(x); % 判断收敛: 用残差的无穷范数 if norm(F, inf) tol return end J Jfun(x); % 解线性方程组 J * dx -F dx -J \ F; x x dx; % 如果步长很小, 也认为收敛 if norm(dx, inf) tol return end end % 超出最大迭代次数, 给出警告 warning(达到最大迭代次数, 可能未收敛!); end这个版本逻辑清晰核心只有四步算残差、判收敛、解线性方程、更新解。注意我用的是J \ F也就是矩阵左除Matlab内部会根据矩阵性质选择合适求解器比inv(J)*F聪明且稳定得多。这里有个细节判断收敛我同时用了残差范数和步长范数两个条件。为什么因为如果函数本身尺度很大残差可能始终很大但解已经在正确位置附近震荡反过来如果函数尺度很小残差早就小于容差但解还差得远。两个条件取或保险。4.2 带雅可比差分近似版本懒得手写导数时用有时候你想快速验证一个方程组的求解思路手写雅可比太费劲。我第二版程序加入了一个中心差分近似你只需要提供方程组函数不用写雅可比。function [x, iter] newton_diff(fun, x0, tol, maxiter) % 中心差分近似雅可比矩阵的牛顿迭代法 x x0(:); n length(x); % 差分步长 h 1e-7; for iter 1:maxiter F fun(x); if norm(F, inf) tol return end % 构造雅可比矩阵(中心差分) J zeros(n, n); for j 1:n xp x; xp(j) xp(j) h; xm x; xm(j) xm(j) - h; J(:, j) (fun(xp) - fun(xm)) / (2*h); end dx -J \ F; x x dx; if norm(dx, inf) tol return end end warning(达到最大迭代次数, 可能未收敛!); end这里h 1e-7是我试下来比较稳的差分步长。要注意如果函数值本身就在 (10^{-8}) 量级这个步长可能会引入严重的舍入误差。一般经验是取 (h \approx \sqrt{\varepsilon})对应双精度浮点差不多是 (10^{-8})。但我选 (10^{-7}) 是综合了截断误差和舍入误差。这是中心差分的优势区间——截断误差是 (O(h^2))比前向差分 (O(h)) 好太多。4.3 带诊断输出版本调试和演示首选这是我最常用的版本因为它能把每次迭代的残差、步长都打出来方便分析收敛过程。尤其适合写论文、画收敛曲线、调初始点的时候。function [x, iter, hist] newton_full(fun, Jfun, x0, tol, maxiter, verbose) % 完整诊断版牛顿迭代法 % 输出 hist: 迭代历史结构体, 包含残差和步长记录 if nargin 6 verbose false; end x x0(:); hist.normF []; hist.normDx []; hist.x zeros(length(x0), maxiter); % 保存每次的解 for iter 1:maxiter F fun(x); nF norm(F, inf); hist.normF [hist.normF; nF]; hist.x(:, iter) x; if nF tol % 记录最后一次 hist.iter iter; if verbose fprintf(迭代%d: ||F||%.6e, 收敛\n, iter, nF); end return end J Jfun(x); dx -J \ F; nDx norm(dx, inf); hist.normDx [hist.normDx; nDx]; if nDx tol hist.iter iter; if verbose fprintf(迭代%d: ||dx||%.6e, 步长收敛\n, iter, nDx); end return end x x dx; if verbose fprintf(迭代%d: ||F||%.6e, ||dx||%.6e, x[%s]\n, ... iter, nF, nDx, num2str(x, %10.6f )); end end hist.iter iter; warning(达到最大迭代次数!); end这里我用了结构体hist存储迭代历史方便你一次性提取残差序列画图。有个小经验在输出里同时打印x的数值你一眼就能看出来迭代是在往正确方向走还是在两个点之间来回横跳那基本就是奇异雅可比问题。5. 一个经典算例二维非线性方程组的完整求解理论讲再多不如手算一个算例来得通透。我用一个我经常在教学里用的例子方程组 [ \begin{cases} x_1^2 x_2^2 - 4 0 \ x_1^2 - x_2 - 1 0 \end{cases} ]这个方程组几何意义很直观第一个方程是圆 (x_1^2x_2^24)第二个方程是抛物线 (x_2 x_1^2 - 1)。你可以想象它们会在平面上相交于两个点。一个是接近 ((1.6, 1.6)) 的区域另一个接近 ((-1.6, 1.6)) 的区域。先写出对应的Matlab方程组函数function F myfun(x) F zeros(2,1); F(1) x(1)^2 x(2)^2 - 4; F(2) x(1)^2 - x(2) - 1; end雅可比矩阵 [ J \begin{bmatrix} 2x_1 2x_2 \ 2x_1 -1 \end{bmatrix} ]对应代码function J myjac(x) J [2*x(1), 2*x(2); 2*x(1), -1]; end然后从初值 (x_0[1,1]) 出发调用x0 [1; 1]; [x, iter] newton_simple(myfun, myjac, x0, 1e-10, 50);跑出来的结果是 (x \approx [1.6180, 1.6180])。等等这个数你眼熟不没错就是黄金分割比为什么这个方程组的解恰好在黄金比上因为把 (x_1x_2t) 代入圆的方程(2t^24)(t\sqrt{2}\approx 1.414)——不对那是另一个解。实际上这个解从抛物线方程 (x_2 x_1^2-1) 和圆方程联立会得到一个四次方程其中一个正根确实出现在黄金比附近。这算是数学模型里一个有趣的小彩蛋跟你讲这个是想说数值解出来之后最好回代验证一下往往能发现意想不到的数学联系。再看另一个初值 (x_0[-1, 1])迭代会收敛到 ([-1.6180, 1.6180])就是对称的那个解。这里有个重点牛顿法是局部方法你给不同的初值它可能收敛到不同的解。所以在工程中如果你知道系统可能存在多个解一定要改变初值多跑几次收集所有可行解再根据物理约束筛选。这一点在机构学、机器人学里特别重要——同一个逆运动学方程可能有多个关节角组合你得选一个不奇异、不碰撞、力矩小的解。6. 三变量算例带参数方程组的牛顿法应用二维太简单来上一个我最近在模拟项目X里实际用过的三维方程组。场景是一套电路系统的稳态工作点求解只要牵扯到二极管、三极管这类非线性器件方程组就逃不出这种形态。给定方程组 [ \begin{cases} 0.2x_1 - e^{x_1} x_2 0 \ x_1 3x_2 2x_3 - 4 0 \ x_2 x_3 - 0.5(x_1x_3) 0 \end{cases} ]这个方程组在物理上是没法直接消元的必须用迭代法。我直接给出完整代码组合。这里注意我把指数函数的参数写成了0.2倍为了让解落在比较舒适的数值区间不然指数项爆炸初值稍微偏离就救不回来。function F f3(x) F zeros(3,1); F(1) 0.2*x(1) - exp(x(1)) x(2); F(2) x(1) 3*x(2) 2*x(3) - 4; F(3) x(2)*x(3) - 0.5*(x(1)x(3)); end function J J3(x) J zeros(3,3); J(1,1) 0.2 - exp(x(1)); J(1,2) 1; J(1,3) 0; J(2,1) 1; J(2,2) 3; J(2,3) 2; J(3,1) -0.5; J(3,2) x(3); J(3,3) x(2) - 0.5; end调用参数x0 [0; 1; 1]; [x, iter] newton_full(f3, J3, x0, 1e-10, 100, true);迭代过程大约在第4到6步就收敛了残差从初始的 (10^{-1}) 量级掉到 (10^{-14}) 量级。这个断崖式下跌就是二阶收敛的直观体现。你如果画hist.normF的曲线会发现前几步还磨磨唧唧最后两步突然就触底了——这是牛顿法的标志性特征。但我必须要说一个很多教材没讲的细节这个方程组如果初值给 (x_0[0;1;1]) 能收敛给 (x_0[5;1;1]) 直接发散。原因就是 (e^{x_1}) 在 (x_15) 时已经爆炸到一百多雅可比矩阵里 (0.2-e^5) 产生了一个极大的负值线性化方向已经完全偏离了。这种问题怎么处理我后面会在实用技巧里专门讲先粗后细和阻尼牛顿法的思路。7. 收敛性诊断画出残差与步长曲线程序跑通了但如果只是跑通其实远远不够。你要学会判断这个解到底可不可信、算法有没有卡进奇怪的状态。我强烈建议跑一个带诊断输出的版本然后做两件事第一打印最终残差norm(F(x))。如果残差在 (10^{-12}) 以上说明你还没真正收敛只是步长变小了这在数学上叫伪收敛。这种情况常见于雅可比矩阵奇异或者方程组本身有冗余约束。你别小看这个检查我把话放这十个牛顿法调试事故里有四个是伪收敛。第二画hist.normF的对数曲线。理想的收敛曲线是一条向下弯曲的抛物线形状在log尺度下说明误差在指数级下降。如果曲线水平震荡或者下降速率突然变慢要么是初始点不合适要么是雅可比矩阵在某一步坏了。这时候我最常用的动作是把那一步的雅可比矩阵打印出来直接用矩阵条件数检查一下。我在Matlab里画收敛历史的标准代码figure; subplot(2,1,1); semilogy(1:length(hist.normF)-1, hist.normF(1:end-1), o-); xlabel(迭代步数); ylabel(||F|| (对数刻度)); title(残差收敛曲线); subplot(2,1,2); semilogy(1:length(hist.normDx), hist.normDx, s-); xlabel(迭代步数); ylabel(||dx|| (对数刻度)); title(步长收敛曲线);这条曲线画出来发论文、写报告、做汇报都好使。下面我会用一个真实的翻车案例来展示这种曲线怎么帮我找出问题。8. 常见问题与排查技巧实录我把踩过的坑都列给你8.1 问题一迭代发散残差越迭代越大现象跑了十几次normF不减反增。排查思路先看是不是初值离收敛域太远。我遇到这种情况第一步永远是把迭代过程打印出来看每一步解的变化方向。如果方向来回折返比如说 (x_1) 从 1 跳到 10 又跳回 1那基本就是雅可比矩阵在某些点附近欺骗了你——线性化误差太大修正量过冲。解决手段给牛顿步加个缩放因子也就是“阻尼牛顿法”(x_{k1} x_k \alpha \cdot \Delta x)其中 (\alpha \in (0,1])。先试探小步残差下降就继续残差上升就缩小 (\alpha)。这招在工程上极其好用。换一个初值用前向欧拉或蒙特卡洛扫一遍初值找到一个残差很低的点再启动牛顿法。做变量缩放。如果有的变量在 (10^{-6}) 量级、有的在 (10^6) 量级雅可比矩阵会病态得不像话。把变量归一化能让收敛半径大幅提升。8.2 问题二伪收敛——步长很小但残差很大现象迭代显示normDx tol提前退出但解根本不对。原因矩阵奇异或者接近奇异时(J^{-1}) 会很大即使残差不小算出来的修正步长也会很小。这时候你被步长阈值骗了。解决手段把收敛条件改成同时检查残差和步长两个都达标才算收敛。检查 (J) 的条件数。Matlab里用cond(J)如果大于 (10^{12})你的问题数值上就是病态的需要预处理。用残差下降率来判断。如果连续三步残差都在同一个数量级徘徊别耗了换方法或者换初值。8.3 问题三雅可比矩阵符号写错导致收敛到奇怪的点这种情况我自己都犯过不止一次。尤其是求偏导时漏掉一个负号程序不会报错数学上照样能收敛只是收敛到一个物理上毫无意义的点。怎么排查把解回代进原始方程逐条检查残差。如果某条残差远大于其他几条意味着这条方程没被满足那要么是雅可比错了要么是程序里方程写错了。8.4 问题四初值不知道给多少这是新手问得最多的。我的经验是三步走动用物理或工程直觉。比如电路问题二极管压降大概0.7V那变量初值就给0.7附近热传导问题温度就在环境温和热源温度之间取个中值。画等高线或残差热图。二维方程组就在网格点扫一遍初值找残差最低的点当初始点。三维以上画不了图就随机抽几百个初值跑一遍选残差下降最顺利的。用连续同伦延拓。从一个容易求解的问题出发慢慢变形到目标问题每一步以上一步的解为初值。这个方法在强非线性问题里属于保底大招但代码量会大不少。9. 牛顿法为什么有时会失败聊聊雅可比矩阵奇异这一节我要单独拎出来讲因为我发现很多人在工程里遇到牛顿法死循环第一反应是调参数、换语言、查精度其实问题出在雅可比矩阵奇异上。想想看牛顿法的迭代式本质上是在解线性方程组 (J \Delta x -F)。如果 (J) 奇异也就是说方程组的约束之间不是独立的那么这个线性化后的系统就没有唯一解。几何上局部切面变成了平塌的你没法确定该往哪个方向走。工程场景里雅可比奇异最典型的案例是机构学的奇异位形——机械臂在某个特定姿态下某个方向的运动自由度丢失逆运动学问题就炸了。应对策略有四个加正则项用 ((J^T J \lambda I)) 代替 (J)类似Levenberg-Marquardt的做法。(\lambda) 是阻尼因子既能保证矩阵可逆又能控制步长大小。换求解器把左除\换成伪逆pinv(J)。虽然牺牲一点效率但至少不会因为奇异而崩溃。扰动初值有时候奇异恰好出现在某个特定迭代点你可以给变量加一个极小的扰动骗过奇异点。重新参数化如果方程组本身有冗余比如3个方程4个未知数那就先消掉一个变量或者在代码里固定一个变量为某个参数再求解。10. 从牛顿法到其他方法需要掌握哪些扩展思路经典牛顿法虽然是基石但工程里常常需要它的变体。我把值得学的扩展列出来按优先级排序阻尼牛顿法每次迭代加入步长控制防止过冲。这是最便宜的改进一行代码的事。拟牛顿法Broyden不重新计算雅可比每次用上一次的信息做秩一更新适合大规模稀疏问题。Levenberg-Marquardt融合了梯度下降和牛顿法的优点在有大量数据拟合问题比如神经网络训练、曲线拟合里是标配算法。同伦延拓/continuation将复杂问题连续映射到简单问题逐步追踪解曲线几乎可以保证大范围收敛。使用全局优化算法先粗搜再用牛顿法精修这是我日常最常用的组合拳。先让遗传算法、粒子群或者模式搜索找到一个残差比较低的点然后切换牛顿法快速收敛。别小看这种作弊打法在工程优化里这叫混合算法效果出奇的好。注意Matlab里其实有内置函数fsolve专门用来求解非线性方程组默认算法是信赖域底层也包含了牛顿法思想。你把fsolve当成干黑活的也行但面试、考试、课题汇报时你要能说清楚算法原理也最好自己实现一遍。毕竟很多定制化需求比如带参数扫描、带约束条件、带非光滑函数内置函数并不好使。11. 实践心得我在模拟项目X里踩过的真实小坑最后分享几个我在具体项目里印象很深的细节。不是高深的数学但很救命。第一关于数据类型的坑。Matlab默认全是双精度但如果你从Excel或某些特别老的脚本里导入数据变量可能变成字符串或整型。别笑我真遇到过因为整数类型导致雅可比矩阵全是0的情况——迭代出来残差纹丝不动排查半天发现x0是整型数组。写代码之前一定要先whos查看变量类型。第二关于收敛容差的选择。我见过有人把tol设成1e-16结果跑半天不收敛整晚在那怀疑人生。其实 (1e-10) 到 (1e-12) 对于绝大多数工程问题已经绰绰有余了。你要是真需要高精度解也别用默认容差去蛮干应该用vpa做高精度运算那是另一套玩法了。第三关于迭代次数的上限。不要天真地设一个很大的数字。我通常设maxiter 50对于二维三维问题足够设再多也逃不开收敛域的约束。如果超过50步还收不了敛一定是初值或算法策略的问题你硬让它继续迭代只是浪费时间。第四关于方程组的无量纲化。物理量程差异大的问题比如一个变量是 (10^6) 量级、另一个是 (10^{-4}) 量级雅可比矩阵的元素可能相差 (10^{10})条件数爆炸。把所有变量标准化到 ([0,1]) 或 ([-1,1])是提升牛顿法稳定性的性价比之王。这个习惯我在做任何数值计算时都会保持。12. 一个帮你少走弯路的建议先画图再迭代如果你问我在所有经验里挑一条最想提前告诉自己的那就是不要急着写迭代循环先把函数在解附近的图像或者等高线画出来。我遇到很多朋友拿着一个方程组就来问“为啥我的牛顿法不收敛”我反手就说你先画一下 (F) 的等高线或者三维图。你随便用一个思考工具比如Matlab的fmesh或者写个循环采样画等高线一眼就能看出解的分布、收敛域的边界在哪。这个信息拿在手里选初值、调策略都是降维打击比你闷头调代码高效十倍。画二维方程组的等高线我的标准代码[x, y] meshgrid(-3:0.01:3, -3:0.01:3); f1 x.^2 y.^2 - 4; f2 x.^2 - y - 1; contour(x, y, f1, [0 0], r, LineWidth, 2); hold on; contour(x, y, f2, [0 0], b, LineWidth, 2); grid on; legend(f10,f20);看到这样一张图两个解的位置清清楚楚符号函数的交点一眼就锁定。初值和收敛域的关系也能在图上找到直观答案。13. 从一个公式到一个工具箱我的总结式建议这篇内容到这里算法原理、程序实现、算例验证、坑点排查全部齐了。最后我想说的是不要把牛顿迭代法当成一个孤立的数学作业它是你面对真实非线性世界时最趁手的一把扳手。掌握它你会开始用线性化的眼光去看待所有复杂系统——这几乎是所有数值计算方法的灵魂。我个人在实际操作中的体会是牛顿法像一个性格直爽但有点倔的合作者你给它好初值它十分钟还你高精度解你欺负它它扭头就跑。所以跟它合作的关键是尊重它的收敛域给它适当地阻尼留意雅可比矩阵的健康状态。如果你能做到这三点你就会发现非线性方程组不过是一次解一个线性方程组的重复游戏。如果你后续还想深入我建议你往这些方向走非线性方程组的数值延拓、带约束优化中的KKT条件线性化、以及大规模稀疏牛顿法在深度学习里的应用。每一块都有足够深的坑和足够有趣的数学等着你。祝你在数值世界里玩得开心也期待你在评论区分享你的踩坑故事。提示如果你在复制程序时遇到编码问题注意检查Matlab脚本的编码格式中文字符建议在英文注释环境下使用或者统一用英文写注释避免Windows下UTF-8和GBK混用造成的乱码。这听起来很小但真实发生过耽误了我一下午。
