简介：面向MATLAB数值计算学习者与科研人员，这份资源包系统整理了求解非线性方程组的多种数值方法。从基础的不动点迭代、经典牛顿法，到牛顿雅可比、牛顿-SOR、离散牛顿、牛顿下山、两点割线、拟牛顿、对称秩1、DFP/BFS、数值延拓、参数微分、最速下降、高斯牛顿、共轭梯度与阻尼最小二乘等，基本覆盖主流求解算法，便于对比不同方法的收敛性、适用场景与实现细节。压缩包共23个文件，含22个.m函数源代码和1个txt说明文档，整体仅11KB，轻量易用。源码函数命名规范，各算法独立封装，可直接调用或修改学习，适合课程作业、算法实验与工程参考。目前已有17047人学习下载，累计认可度较高。无论是数值分析初学者入门，还是需要快速选用合适算法的开发人员，都能从这套源码中获得清晰、可运行的MATLAB实现。

1. 十多种方法求解非线性方程组：MATLAB 源码包拆解与复现笔记

提到「MATLAB 非线性方程组」求解，多数人的第一反应是 fsolve 一把梭。但真正处理工程仿真里的强耦合方程时你会发现，fsolve 这个黑匣子对初值极其敏感：同一个方程组，换个初值就从收敛变成发散，报错信息还很难读。这份「十多种方法——求解非线性方程组 MATLAB」源码包把二分法、牛顿法、牛顿下山、简化牛顿、割线法、Broyden、SR1、DFP、BFGS、不动点迭代、梯度下降、高斯牛顿、Levenberg-Marquardt、信赖域、同伦延拓十五种算法全部实现成独立可跑的 MATLAB 函数，统一输入输出结构，自带三个测试方程组和一键对比脚本。

这套源码适合两类人：写论文需要可复现对比基线的研究者，以及被初值敏感性和 Jacobian 计算反复坑过的工程开发者。前者能直接拿收敛曲线和迭代次数当实验数据，后者把残差函数替换掉就能接入自己的仿真流程，内部实现基本不用动。

2. 项目定位与技术选型：为什么纯 MATLAB 能撑起一套方法库

2.1 三个工程痛点决定的方法集合

这套代码的初衷不是把教材算法机械翻译成 MATLAB，而是围绕我在仿真项目里反复遇到的三类问题来组织方法集合。看一遍方法列表就能读出它的设计逻辑：牛顿家族负责收敛速度，拟牛顿和无导数方法负责免去 Jacobian 的痛苦，同伦延拓和信赖域负责给初值发散兜底。

第一个痛点是初值选择。工程里给非线性方程组取初始猜测，大多数时候靠上一时步的解或者拍脑袋，系统状态一旦跳变，牛顿法的收敛域就接不住了。项目因此实现了三种带保护机制的算法：牛顿下山给每步迭代加上阻尼因子，让残差范数保持单调下降；信赖域把迭代步长限制在信任半径以内，模型预测不靠谱就缩小半径重来；同伦延拓更直接——把待解方程和另一个容易解的方程用参数 t 连续桥接，从 t=0 处的已知解一路追踪到 t=1 的目标解，初值给得再离谱也能沿着解流形慢慢摸过去。

第二个痛点是 Jacobian 求不出来。工程模型的残差函数往往对应一大段仿真逻辑，根本没有解析导数可用，手动推导又容易出错。项目默认用前向差分计算数值 Jacobian，差分步长、奇异矩阵兜底都封装在 utils 目录里；如果你手头有解析 Jacobian，也可以按约定传入，分发器会优先采用解析版本而不是覆盖它。

2.2 十五种方法的技术矩阵

选中 MATLAB 而不是 Python 或 C++ 是有原因的。矩阵左除J \ fx这种操作在 MATLAB 底层会自动选择 LU、Cholesky 或迭代求解路径，稠密稀疏都不用自己操心；核心求解路径不依赖任何工具箱，Symbooolic 工具箱只是可选加速手段，拿到裸 MATLAB 环境就能跑。再加上每个方法只有一两百行代码，改一个文件、跑一个脚本就能验证效果，这种迭代速度在工程项目里是实打实的优势。

设计上还留了一个容易被忽略的细节：所有方法文件都接收params结构体，而不是一长串位置参数。原因很简单，十五种方法各自需要的控制参数差异很大——信赖域要半径初值，同伦要 t 的推进步长，LM 要阻尼因子更新策略，如果全部塞进函数签名，调用方就得记十五套签名。统一成结构体之后，公共字段只有 method、tol、maxiter 三个，方法特有字段按需填充，没填就用默认值。这种「公共参数直通、私有参数可选」的协议，是这套代码能被快速二次开发的骨架。

2.3 测试口径：为什么对比实验要统一初值和统一容差

项目自带的三个测试方程组选得很有心机。变形 Rosenbrock 方程组制造了一个狭长弯曲的谷底，专门用来考验算法的「下山」能力；二极管电路方程组来自基尔霍夫定律，残差函数带上指数项，数值差分 Jacobian 在这种函数上很容易出现溢出；化学反应平衡浓度方程组则是典型的病态缩放问题，各变量量纲能差出四五个数量级。把这些方程放进同一个 run_all_demo 脚本里跑，等于把十五种方法拉到同一张考卷上。

对比口径也做了强制统一：所有方法从同一个 x0 出发，用同一个 tol 和 maxiter，返回结构体里都带residual_history字段。这样做的价值在于，论文里常见的「A 方法收敛 B 方法不收敛」的说法，在这套代码里是可以复现的——把两个方法的曲线画在同一张 semilogy 图上，谁在哪个迭代步卡住、谁的残差先穿过了阈值，一目了然。后面第五部分我会专门讲这里面的坑，对比实验的口径稍微偏一点，结论就可能翻车。

3. 目录结构与核心模块拆解：十五个求解器如何被一套接口串起来

3.1 先看目录结构再动手

解压源码包之后首先看到的是这样一个工程布局：

nonlinear_solver/ ├── run_all_demo.m # 一键对比脚本 ├── run_single_demo.m # 单方法演示脚本 ├── example_rosenbrock.m # 变形 Rosenbrock 方程组 ├── example_circuit.m # 二极管电路节点方程 ├── example_chemical.m # 化学反应平衡浓度 ├── core/ │ ├── solver_base.m # 统一入口/方法分发器 │ ├── bisection_method.m # 二分法 │ ├── newton_solver.m # 牛顿法 │ ├── newton_descent.m # 牛顿下山 │ ├── simplified_newton.m # 简化牛顿 │ ├── broyden_method.m # Broyden 拟牛顿 │ ├── quasi_newton_sr1.m # SR1 │ ├── quasi_newton_dfp.m # DFP │ ├── quasi_newton_bfgs.m # BFGS │ ├── secant_method.m # 割线法 │ ├── fixpoint_iter.m # 不动点迭代 │ ├── gradient_method.m # 梯度下降 │ ├── gauss_newton.m # 高斯牛顿 │ ├── levenberg_marquardt.m # Levenberg-Marquardt │ ├── trust_region.m # 信赖域 │ └── homotopy_method.m # 同伦延拓 ├── utils/ │ ├── jacobian_numeric.m # 数值差分 Jacobian │ ├── line_search.m # 线搜索 │ ├── check_converge.m # 收敛判据 │ └── plot_converge.m # 收敛曲线绘图 └── docs/ └── 方法对比表.md

core 目录下是十五个独立求解器，每个文件对应一种算法，不互相依赖。utils 目录放的是公共件：数值 Jacobian、线搜索、收敛判据和绘图函数。example 开头的三个文件是测试方程组，本质上是返回残差向量的函数句柄，可以被任何方法调用。run_all_demo 和 run_single_demo 是唯二需要你亲手运行的脚本，前者循环跑全部十五种方法并打印统计表，后者演示单个方法的完整调用方式。

把求解器和示例方程分开存放是刻意的。求解器只认「残差函数句柄 + 初值 + 参数」这个三元组，不知道也不关心残差背后是电路还是化学反应。你要接自己的工程模型，只需要新写一个返回残差的函数，其他什么都不用动。

3.2 统一入口 solver_base：一个分发器串起十五个方法

solver_base.m 是整个包的门面，它的核心是一个 switch 分发逻辑：

% core/solver_base.m % 统一分发入口：所有方法共用 F, x0, params 三个入参 function result = solver_base(F, x0, params) if nargin < 3, params = struct(); end if ~isfield(params, 'method'), params.method = 'newton'; end if ~isfield(params, 'tol'), params.tol = 1e-8; end if ~isfield(params, 'maxiter'), params.maxiter = 200; end switch params.method case 'bisection' % 仅适用于单变量方程 result = bisection_method(F, x0, params); case 'newton' result = newton_solver(F, x0, params); case 'newton_descent' result = newton_descent(F, x0, params); case 'simplified_newton' result = simplified_newton(F, x0, params); case 'secant' result = secant_method(F, x0, params); case 'broyden' result = broyden_method(F, x0, params); case 'sr1' result = quasi_newton_sr1(F, x0, params); case 'dfp' result = quasi_newton_dfp(F, x0, params); case 'bfgs' result = quasi_newton_bfgs(F, x0, params); case 'fixpoint' result = fixpoint_iter(F, x0, params); case 'gradient' result = gradient_method(F, x0, params); case 'gauss_newton' result = gauss_newton(F, x0, params); case 'lm' result = levenberg_marquardt(F, x0, params); case 'trust_region' result = trust_region(F, x0, params); case 'homotopy' result = homotopy_method(F, x0, params); otherwise error('solver_base: unknown method: %s', params.method); end result.method = params.method; % 回填方法名，便于后续统计 end

params 的默认值处理是这套协议的关键。用户只填 method 就能跑，tol 默认 1e-8、maxiter 默认 200，每个求解器内部都假设这两个字段一定存在。这种「先补默认再分发」的写法避免了十五个文件里重复写isfield判断，也让后续加新方法时只需要在 switch 里多一个 case，改动面被限制在一个文件里。

值得留意的是分发器最后回填了result.method。这样 run_all_demo 收集结果时不用自己记住循环到哪个方法，直接读每个 result 的 method 字段就能对齐做统计。你如果要加第十六种方法，复制一个现有求解器文件改算法，在 switch 里加一行 case 就算接入了。

3.3 公共件：数值 Jacobian 与双判据收敛检测

数值 Jacobian 是整套代码里被调用最频繁的函数，牛顿家族、高斯牛顿、LM、信赖域全都依赖它：

% utils/jacobian_numeric.m % 前向差分数值 Jacobian，返回 n x n 矩阵 function J = jacobian_numeric(F, x) n = numel(x); f0 = F(x); h = 1e-7; % 固定差分步长 J = zeros(n, n); for j = 1:n xp = x; xp(j) = xp(j) + h; % 只扰动第 j 个分量 J(:, j) = (F(xp) - f0) / h; end end

前向差分是工程里最不容易出幺蛾子的选择，实现只有七行，但代价是精度只有 O(h)。固定步长 1e-7 在大多数双精度问题里够用，如果你手头的方程各变量量纲差异很大，常见做法是改成 h = 1e-7 * max(abs(x(j)), 1) 的自适应形式，避免某个分量本身是 1e-10 却被一步 1e-7 的扰动淹没。另外要留意，如果残差函数里含exp、log这类强非线性，差分可能引入溢出，这也是为什么项目里二极管电路例子跑数值 Jacobian 时要格外小心。

收敛判据同样是个容易被低估的公共件：

% utils/check_converge.m % 双判据：残差范数与步长同时满足阈值才算收敛 function flag = check_converge(fx_new, x_new, x_old, tol) res_flag = norm(fx_new) < tol; % 残差足够小 step_flag = norm(x_new - x_old) < tol; % 步长足够小 flag = res_flag && step_flag; end

很多自写的求解器只盯着残差，结果遇到 Rosenbrock 这类谷底平坦的问题，残差已经降到 1e-10 但解还在缓慢漂移。这里用「残差 + 步长」双判据，两个条件同时满足才算收敛，能滤掉大部分假收敛。不过双判据也有副作用：有些方法在解附近步长天然极小，容易提前触发收敛，所以各求解器主循环里仍然保留了独立的残差判断作为快速路径，check_converge 主要在方法论验证场景下统一口径用。

4. 代码逐段精讲：一个方程组从输入到收敛的完整链路

4.1 测试方程怎么定义：以变形 Rosenbrock 为例

要理解整条调用链，先从残差函数的定义看起。关于这套代码的实际调用方式，我建议你直接在命令行敲run_single_demo跑一遍，然后按下面的顺序读代码：

% example_rosenbrock.m % 变形 Rosenbrock 方程组，真实解为 x* = (1, 1) function F = example_rosenbrock(x) F = zeros(2, 1); F(1) = 10 * (x(2) - x(1)^2); % 强耦合项，制造弯曲谷底 F(2) = 1 - x(1); % 线性约束把解钉在 x1 = 1 end

残差函数的协议非常简单：输入一个 n 维列向量 x，输出一个 n 维残差列向量 F，求解的目标就是让 F 的所有分量同时为 0。这里 F 的长度必须和 x 的长度一致——如果是欠定或超定方程组，这套代码并不直接支持，需要你先做变量代换或者改用最小二乘语义。example_circuit 和 example_chemical 两个文件遵循完全相同的协议，只是残差函数内部换成电路节点方程和化学平衡方程式而已。

调用侧同样简单：

% run_single_demo.m 中的核心调用片段 x0 = [-1.2; 1.0]; % 经典难初值，离解很远 params.method = 'newton_descent'; params.tol = 1e-9; params.maxiter = 300; r = solver_base(@example_rosenbrock, x0, params); fprintf('converged=%d, iter=%d, x=(%.6f, %.6f)\n', ... r.converged, r.iter, r.x(1), r.x(2));

注意到 solver_base 第一个参数传的是@example_rosenbrock函数句柄，而不是函数名字符串。这意味着你可以现场构造匿名函数，比如@(x) my_model(x, data)把外部数据通过闭包传进去，这也是后续做二次开发时把工程模型接入这套框架的标准姿势。初值 x0 必须是列向量，代码内部第一行x0(:)做了强制的形状规整，所以传行向量也不会报错，但返回的 result.x 一律是列向量。

4.2 牛顿法与牛顿下山：先保证能收敛，再谈收敛快

牛顿求解器是整个包的主干，其他大多数方法内部都会复用它：

% core/newton_solver.m % 标准牛顿法：x_{k+1} = x_k - J(x_k) \ F(x_k) function result = newton_solver(F, x0, params) x = x0(:); tol = params.tol; maxiter = params.maxiter; result.x = x; result.iter = 0; result.converged = false; result.residual_history = []; for k = 1:maxiter fx = F(x); res = norm(fx, 2); result.residual_history(k) = res; if res < tol result.converged = true; break; end J = jacobian_numeric(F, x); if rcond(J) < 1e-14 % Jacobian 接近奇异时改用伪逆，避免直接左除报错 delta = -pinv(J) * fx; else delta = -J \ fx; end x = x + delta; end result.x = x; result.iter = k; end

这个实现里有两个容易被忽略的工程细节。第一，rcond(J) < 1e-14的判断放在左除之前，rcond 返回的是矩阵条件数的倒数，越接近 0 说明矩阵越奇异。直接用J \ fx遇到奇异矩阵会抛出「Matrix is singular」的硬错误，整个脚本直接中断；换成pinv伪逆之后至少能继续迭代，虽然步长可能不准，但给了后续下山或者信赖域策略重新规划的机会。第二，residual_history 在每次迭代开始时记录残差，这样即使迭代中断，你仍然能把已经走过的收敛路径画出来，排查到底是哪一步开始发散的。

牛顿下山在牛顿法基础上加了一个线搜索保护：

% core/newton_descent.m % 牛顿下山：在完整牛顿方向上做阻尼，保证残差范数单调下降 function result = newton_descent(F, x0, params) x = x0(:); alpha = 1.0; for k = 1:params.maxiter fx = F(x); res = norm(fx); result.residual_history(k) = res; if res < params.tol result.converged = true; break; end J = jacobian_numeric(F, x); delta = -J \ fx; % 先求完整牛顿方向 % 从 alpha=1 开始，不满足下降条件就折半，最多折到 1e-6 while norm(F(x + alpha * delta)) >= res && alpha > 1e-6 alpha = alpha * 0.5; end x = x + alpha * delta; end result.x = x; result.iter = k; end

下山策略的哲学是：牛顿方向可能非常准，但当一个完整的牛顿步会把你推出收敛域时，沿着同一个方向走半步甚至四分之一步反而是更优的选择。while循环里的下降条件用的是残差范数的硬比较——新点的残差必须严格小于当前残差，否则就继续折半。alpha 从 1 开始意味着当牛顿步本身没问题时，下山法退化回纯牛顿法，不损失收敛速度。这也是初值问题排查时我第一个推荐的手段：如果纯牛顿法发散，换method='newton_descent'往往立竿见影。

4.3 BFGS 拟牛顿：不显式碰 Jacobian 也能保持超线性收敛

拟牛顿家族的思路和牛顿法完全不同：它不直接计算 Jacobian，而是用迭代过程中产生的 s 向量和 y 向量不断修正一个对 Hessian 逆的近似矩阵 H：

% core/quasi_newton_bfgs.m % BFGS 拟牛顿：每一步用 s 和 y 更新逆 Hessian 近似 H function result = quasi_newton_bfgs(F, x0, params) x = x0(:); n = numel(x); H = eye(n); % 初始逆 Hessian 近似取单位阵 fx = F(x); tol = params.tol; for k = 1:params.maxiter if norm(fx) < tol result.converged = true; break; end d = -H * fx; % 搜索方向，等价于 -Hessian^{-1} * 梯度 alpha = line_search(F, x, d, fx); s = alpha * d; % 实际走的步长向量 x_new = x + s; fx_new = F(x_new); y = fx_new - fx; % 残差变化量 % 曲率条件检查：s'*y > 0 才能保证 H 保持正定 if s' * y > 1e-12 rho = 1 / (y' * s); I = eye(n); H = (I - rho * s * y') * H * (I - rho * y * s') + rho * (s * s'); end x = x_new; fx = fx_new; end result.x = x; result.iter = k; end

这里最该读懂的是那个超长的 H 更新公式。BFGS 的核心思想是：我们希望更新后的 H 满足割线方程H * y = s，同时保持对称正定。上面这行更新公式就是从这两个约束里推出来的秩二修正——前半部分(I - rho*s*y') * H * (I - rho*y*s')负责保持对称性，后半部分rho * s * s'负责满足割线条件。s'*y > 1e-12这个曲率条件检查不能省：当步长太小或者残差变化异常时，割线方程不成立，强行更新会让 H 失去正定性，后续迭代直接崩。

line_search在这里的作用和下山法里的 alpha 折半很相似，但搜索准则更宽松——它只要求目标函数（这里是最小化 0.5*||F||^2）有足够下降，而不是残差范数严格单调。BFGS 对线搜索的精度要求不高，Armijo 型、回溯型都能配合，工程实现里用回溯折半最常见。这个文件里体现的设计取舍是：用矩阵-向量乘法和秩一修正替代了每次迭代都要做的 Jacobian 数值差分，对于 n 较大且 F 求值昂贵的工程模型，省下的成本非常可观。

4.4 同伦延拓：初值乱给时的最后防线

前面所有方法都有一个隐含前提：初值离真解不能太远。同伦延拓是唯一主动打破这个前提的方法，它的思路是把一个「不会解的方程」连续变形到一个「会解的方程」：

% core/homotopy_method.m % 同伦延拓：H(x,t) = F(x) - (1-t)*F(x0)，t 从 0 推进到 1 function result = homotopy_method(F, x0, params) x = x0(:); t = 0.0; dt = 0.1; % t 的推进步长 while t < 1 t = min(t + dt, 1); % 最后一步校准到 t=1 G = @(xx) F(xx) - (1 - t) * F(x0); % 当前同伦方程 inner = params; inner.maxiter = 50; % 内层牛顿迭代上限收紧 x = newton_solver(G, x, inner).x; if ~isfinite(norm(x)) % 数值爆炸就缩小 dt 回退重试 dt = dt * 0.5; t = t - dt; end end result.x = x; result.converged = norm(F(x)) < params.tol; end

理解这段代码的关键在于那个 G 的定义。当 t=0 时，G(x) = F(x) - F(x0)，显然 x0 就是 G 的精确解；当 t=1 时，G(x) = F(x)，G 的解就是我们真正要求的 F 的解。中间每个 t 值对应一个「比原方程简单一点」的方程，内层牛顿求解器从上一个 t 的解出发，通常几步就能收敛到当前 t 的解。这样一层层接力，初值再离谱也能被「带」到目标解流形上。

我个人的使用习惯是：先用同伦延拓跑一个粗糙解，再把这个解作为初值喂给 BFGS 或牛顿法精化。同伦延拓的收敛鲁棒性最好，但迭代次数多、计算量大，直接拿它做高精度求解不划算。dt 的默认值 0.1 意味着最多十段接力，如果中间某段内层牛顿失败，isfinite检查会捕获 NaN 并自动把 dt 减半重来，这个自我保护机制在处理强非线性方程时几乎一定会被触发。

5. 运行与验证避坑：五条血泪经验与排查清单

十五种方法在同一套协议下跑起来之后，真正的战斗才刚刚开始——验证结果、排查发散、对比口径，每一步都有坑。下面的五条问题是我在复现和改造这套代码时真实踩过的，按出现频率排序，每条都按「现象 → 原因 → 解决」给你排查路径。

5.1 牛顿法迭代两步就直接飞掉

现象：result.converged 为 false，residual_history 里残差先降一点，第二步突然跳到 1e30 量级，然后全是 NaN。

原因：初值落在牛顿法收敛域之外。牛顿法的局部收敛性质决定了它只在真解附近才保证收敛，Rosenbrock 这类方程在 [-1.2, 1] 附近恰恰有一段「反斜面」，完整的牛顿步直接把你弹到数值溢出区。这是牛顿法最经典的翻车方式，和代码实现无关，是方法本身的数学性质。

解决：先别换方法，把params.method改成newton_descent，线搜索的折半机制会自动把膨胀的步长压回来，大多数初值问题到这里就解决了。如果下山法还在飞，再升级到trust_region或者homotopy。经验数据是：90% 的初值发散靠牛顿下山就能兜住，剩下 10% 交给同伦延拓。

5.2 手写解析 Jacobian 和数值差分结果对不上

现象：你在自己改造时按公式手写了 Jacobian，跑出来结果比纯数值差分的版本差很多，甚至不收敛。

原因：绝大多数情况是手写 Jacobian 的表达式顺序和残差函数定义不一致。比如残差函数里写了x(2) - x(1)^2，但 Jacobian 的第 (1,1) 分量手算时把 -2*x(1) 写成了 -x(1)，这种低级错误在变量多的方程组里很难用眼睛发现。另一个常见原因是残差函数内部做了变量代换，但 Jacobian 没有链式法则。

解决：强制做一致性校验，不要相信眼睛。写一个几分钟的验证脚本，在随机点上计算norm(J_analytic - J_numeric) / norm(J_numeric)，这个比值应该小于 1e-6。如果大于它，逐列比对 Jacobian，把数值差分的结果打印出来和手写结果对齐。从那以后我每次接入新模型都强制走一遍这步校验，再也不信手推的导数。

5.3 残差降到 1e-10 但解还在漂

现象：residual_history 显示残差早就低于 tol，但 result.x 的数值和真解差出两个量级，迭代也停了。

原因：残差范数小不代表解准。Rosenbrock 函数的谷底就像一条平缓的水槽，垂直于谷底的方向残差下降很快，沿着谷底的方向残差几乎不变。这种「香蕉峡谷」地形会让纯残差判据提前触发假收敛，解停在半山腰而不是谷底。

解决：改用双判据，残差范数和步长同时满足阈值才判定收敛。排查时把 x 的每一步变化量norm(x_new - x_old)也记录成一条曲线，和残差曲线画在一起。如果残差平了但步长还在 1e-2 量级，说明还在峡谷里慢慢地挪，给它加迭代次数而不是放宽容差。

5.4 方程有复数解，实数域迭代跑出 NaN

现象：求解过程中 result.x 突然变成 NaN，residual_history 最后几个值直接是 Inf。

原因：方程组本身只有复数解，你给了一个实数初值，数值差分在前向扰动时碰上了需要开负数的中间表达式，直接产出 NaN。更隐蔽的情况是同伦延拓路径穿过了解流形的分叉点，t 推进步长过大时跃迁到复平面又弹不回来。

解决：先用符号计算或画图确定方程组解的性质。如果是复数解，需要改残差函数让 x 允许复数；如果确认实数域里有解，把同伦延拓的 dt 从 0.1 调小到 0.05 再跑一遍，缩小追踪步长能减少跃迁概率。这里的排查信号很明确——看到 NaN 先查问题本身有没有实数解，而不是先怀疑代码。

5.5 方法对比跑出来的结论和论文相反

现象：run_all_demo 跑完，BFGS 的迭代次数比牛顿法多出好几倍，和高引论文里的「BFGS 优于牛顿法」结论矛盾。

原因：对比口径不公平。牛顿法每步要花 O(n^2) 做数值差分求 Jacobian，BFGS 每步只是 O(n^2) 的矩阵更新，但迭代次数确实比牛顿法多。只看迭代次数，牛顿法永远领先；只看单步成本，BFGS 永远领先。更隐蔽的是线搜索参数——BFGS 对线搜索宽容度的敏感度远高于牛顿法，统一参数跑出来的结果不能代表方法的天赋。

解决：把对比指标拆成三个维度：迭代次数、总残差函数调用次数、末端精度。虚线是 BFGS 的迭代次数优势主要体现在工程场景——残差函数求值成本高的时候，少算几次 F 比少算几次迭代更值钱。跑对比前还要检查每个方法有没有被分配合理的内部参数，比如信赖域的初始半径和 LM 的阻尼初值，这些参数对不对直接决定对比有没有意义。

6. 二次开发思路：把方法库接进你自己的工程模型

6.1 换残差函数的三步走

最常做的改造是把自己的工程模型接进来。第一步，写一个和 example_rosenbrock 同样协议的残差函数：

% my_model_residual.m % 工程模型残差模板：x 是待求变量，F 是残差，长度必须一致 function F = my_model_residual(x) F = zeros(size(x)); F(1) = ...; % 你的第一个方程残差 F(2) = ...; % 你的第二个方程残差 end

第二步，在调用脚本里把@example_rosenbrock换成@my_model_residual，其余什么都不用改。第三步，先跑一次homotopy拿到粗糙解，再用粗糙解当初值跑bfgs或者lm精化，这是我推荐的工程接入流程。如果要传外部参数，用@(x) my_model_residual(x, data1, data2)构造匿名函数即可。

6.2 强制走一遍验证三件套

接入新模型后，我的习惯是强制走三个验证动作：第一，把数值 Jacobian 和一个自己推导的解析版本做一致性比对，比值小于 1e-6 才放心；第二，用相差一个量级的两个初值各跑一遍，两次结果一致才说明收敛到了同一个解而非局部陷阱；第三，把 residual_history 用 semilogy 画出来，理想的收敛曲线应该单调下降且尾部陡峭，如果曲线出现平台或反弹，优先怀疑线搜索参数而不是算法本身。

6.3 扩展新方法的最小改动

这套代码的扩展点设计得相当清晰：新增方法时，复制 core 目录下一个结构最接近的文件，比如要加阻尼最小二乘就复制 gauss_newton.m，保留输入输出协议和 result 结构，只改中间的迭代更新公式，然后在 solver_base.m 的 switch 里加一个 case 分支。整个改动就这两个文件，不影响其他十四个方法。从那以后我每次接新工程模型都强制走一遍「残差函数 + 三件套验证 + 对比跑分」的固定流程，省掉的返工时间远比写这几个脚本的时间多。这套源码包在下载页直接拿，跑通 run_single_demo 再改自己的模型，希望帮到你。

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