简介：面向数值计算、科学计算与数学建模学习者的一套 MATLAB 非线性方程组求解代码合集，聚焦多元非线性方程组的求根与迭代算法实现。内容涵盖不动点迭代、牛顿法、离散牛顿法、牛顿-SOR 与牛顿下山法、两点割线法、拟牛顿类（对称秩1、D-F-P、B-F-S）、数值延拓、参数微分法、最速下降、高斯牛顿、共轭梯度及阻尼最小二乘等十余种经典与改进算法，每种方法对应独立函数文件，便于按名称定位并横向比较收敛速度与适用条件。压缩包共 23 个文件，以 22 个 .m 函数脚本为主体，另附 1 个说明文本，整体约 11KB，导入即可运行调试。已有 17052 人学习，可作为数值分析课程实践、课程设计或建模竞赛的求解工具，也便于对照不同迭代格式的收敛性与稳定性。

1. 从 200 步不收敛到 11 步收敛：十多种解法到底在解决什么

同一组三元方程，有人用 fsolve 默认设置跑满 200 次迭代，只换来一句Solver stopped prematurely，残差卡在 1e-2；把变量做一次尺度变换、雅可比换成解析式、算法换成带线搜索的阻尼牛顿，11 步残差就掉到 1e-11。MATLAB 非线性方程组的求解质量，基本都取决于方法选型、初值和缩放这三件事，敲代码反而是最简单的部分。

这份资料把十多种方法收在一个包里：牛顿、修正牛顿、弦截、Broyden 拟牛顿、阻尼牛顿、信赖域 dogleg、Levenberg-Marquardt、高斯-牛顿、不动点迭代、二分与试位、同伦延拓、符号求解，每一类都配了能直接运行的 MATLAB 实现和对照算例。它解决的不是“能不能解”，而是“fsolve 报错时手里还剩几张牌”——做课程设计要对比收敛阶的、工程里被小规模方程组卡住的，都能直接拿去改。

2. 建模先于编码：残差、雅可比与收敛判据的三道关口

同一组方程，写成什么形式交给求解器，结果可能差出两个数量级。这一步做错，后面换再多算法也只是在错误的条件下比较。下面用一个三元方程组贯穿：三个未知量分别代表长度、比例和温度系数，量级差得比较远，正好能暴露缩放问题。

2.1 把工程问题写成 F(x)=0：缩放比公式更重要

先把方程整理成残差向量形式，所有分量都移到等号左边，让“解”对应残差全零。这一步看似机械，但两个细节直接决定后面能不能收敛：一是每个残差分量的量纲要尽量接近，二是未知量的量级要拉平。前者影响TolFun的判定，后者影响雅可比矩阵的条件数。

function F = eqsys(x) % 三元非线性方程组的残差向量，输入输出都用列向量 F = zeros(3,1); F(1) = x(1)^2 + x(2)^2 + x(3)^2 - 4; % 球面约束，量级 ~1 F(2) = x(1)*x(2) - x(3) - 0.5; % 交叉项，量级 ~1 F(3) = exp(-x(1)) - x(2) + x(3)^2 - 1; % 指数项，对 x(1) 极敏感 end

三个分量的系数都控制在 1 附近，是为了让norm(F)有可比性。如果某个分量天然带 1e6 的系数（比如应力平衡式），缩放时就要整体除以这个量级，否则 fsolve 的默认TolFun=1e-6会被这个分量单独决定，另外两个分量实际上根本没收敛到位。

变量缩放用一句代换就能完成，把待求量写成x = s.*y，让y的各分量量级接近：

s = [1; 1; 10]; % 第三个变量真实量级约 10，缩到 1 附近 f = @(y) eqsys(s.*y); % 缩放后的残差函数句柄 Jwrap = @(y) Jfun(s.*y) * diag(s); % 链式法则：dF/dy = dF/dx * dx/dy

s的取法有讲究：按解的预估量级取，不知道就先跑一次看output里的迭代轨迹。Jwrap那行是链式法则，缩放后雅可比要右乘diag(s)，漏掉这一步，牛顿方向会系统性偏小，表现为“每步都下降但下降得极慢”。我一般习惯在函数里留一个scale参数，默认全 1，出问题再调，比事后改代码省事。

2.2 雅可比矩阵的三档来源：解析、符号、有限差分

牛顿族的收敛速度全靠雅可比。来源分三档，成本从低到高、可靠性从高到低，选哪档取决于方程规模和迭代次数预算。

第一档是手写解析式，最可靠。三项方程的雅可比手推十几分钟，但每一项都是精确的，后面无论怎么改算法都不会因为导数不准而翻车。第二档是符号工具箱自动求导，适合方程项数多、手推容易错但表达式还算规整的场景：

syms x1 x2 x3 F = [x1^2 + x2^2 + x3^2 - 4; x1*x2 - x3 - 0.5; exp(-x1) - x2 + x3^2 - 1]; J = jacobian(F, [x1 x2 x3]); % 3x3 符号矩阵 Jfun = matlabFunction(J, 'Vars', {[x1 x2 x3]}); % 输入整向量，输出 3x3 Ffun = matlabFunction(F, 'Vars', {[x1 x2 x3]}); % 顺手把残差也生成了

'Vars', {[x1 x2 x3]}这个写法让生成的函数接收一个向量而不是三个标量，正好对上 fsolve 的调用约定；不加这行，你会得到一个需要三个入参的函数，套不进统一的调用壳。符号求导的代价在生成阶段，一旦生成为普通函数，运行时和手写没有区别，这是性价比最高的一档。

第三档是有限差分，fsolve 内部默认就在用，但自己写一份能看清步长怎么影响精度：

function J = fdjac(f, x, h) % 前向差分雅可比；h 可为标量或与 x 同长的向量 if isscalar(h), h = h * ones(size(x)); end n = numel(x); f0 = f(x); J = zeros(numel(f0), n); for k = 1:n xp = x; xp(k) = xp(k) + h(k); J(:,k) = (f(xp) - f0) / h(k); % 第 k 列 = 对第 k 个变量的偏导 end end

步长h不是越小越好。前向差分的截断误差是 O(h)、舍入误差是 O(eps/h)，两者平衡在h ≈ sqrt(eps) ≈ 1.5e-8，这也正是 fsolve 的默认值。中心差分精度高一阶，最优步长变成eps^(1/3) ≈ 6e-6，代价是函数调用次数翻倍。方程里含exp、sqrt、log这类在解附近变化剧烈的项时，固定 1e-8 的步长会让差商里全是舍入噪声，雅可比看起来“差不多对”，牛顿方向却完全跑偏——这是后面第 5 章要重点说的一类坑。

2.3 收敛判据与假收敛：别只看 exitflag 正不正

fsolve 的exitflag常见取值里，1 表示一阶最优性满足，2 表示步长小于TolX，3 表示残差变化小于TolFun，4 表示搜索方向幅值过小，0 表示迭代超限，负值表示失败。问题在于 2、3、4 这三个“看起来成功”的返回值，只说明“不再变化”，不说明“变对了”。变量量级 1e6 时，默认TolX=1e-6相当于相对精度只有 1e-12，其实很严；反过来变量量级 1e-3 时，同样 1e-6 就是个极松的条件，求解器可能在一堆残差还有 1e-2 的点上宣布收敛。

稳妥的做法是不信任任何一个判据，自己复核一遍残差和后向误差：

opts = optimoptions('fsolve', 'Display', 'off'); [x, ~, exitflag, output] = fsolve(f, x0, opts); r = norm(f(x)); % 前向误差：残差本身多大 Jx = Jwrap(x); % 后向误差：距离真解还有多远 db = r / max(norm(Jx), eps); % 一阶估计：一步能挪多远 fprintf('exitflag=%d, iter=%d, feval=%d, ||F||=%.3e, 后向误差=%.3e\n', ... exitflag, output.iterations, output.funcCount, r, db);

output.funcCount这一项经常被忽略。用有限差分时，一次迭代内部至少调用 n 次函数，如果funcCount已经几万而迭代数只有几十，说明时间全花在差商上，这时候该做的是给JacobPattern或改解析雅可比，而不是把MaxIterations往上加。后向误差db的意义是：如果这个数已经到 1e-12，前向残差还有 1e-4，那多半是雅可比接近奇异，解本身对扰动极其敏感，继续迭代没有意义，得回头改方程或者换变量。

注意：norm(f(x))要针对缩放后的函数算，不要在缩放坐标和原始坐标之间混着比，量级差一个s，结论会完全反过来。

3. 十多种方法怎么选：从牛顿的二次收敛到同伦的兜底

方法多不代表都要用，真正需要记住的是每种方法的“代价”和“失效条件”。下面这张表是按实际调试顺序排的，从上往下：先试最省事的，不行再往下走。

3.1 牛顿族的二次收敛，代价是一整个雅可比

二次收敛的含义是误差按平方缩小：这次 1e-2，下次 1e-4，再下次 1e-8。这个性质在解附近极强，但它是局部的——初值离解足够远时，牛顿法可能第一步就把点甩到定义域外。判断“够近”没有严格门槛，我的经验是看第一步的残差下降比：||F(x1)||/||F(x0)|| < 0.5说明方向大体可信，接近 1 甚至变大，就该上阻尼或者换初值了。

修正牛顿法把雅可比固定在第一帧，省掉每步的求导和分解。它把二次收敛降到线性，但单步成本可能只有原来的十分之一，在 J 计算极贵（比如每列要跑一次仿真）的场景里反而更快。什么时候用：迭代次数不超过 20 次、每次求 J 超过 1 秒，那基本上固定 J 更划算。

3.2 弦截、Broyden 与不动点：省导数的三条路

弦截法是单变量的割线法，用最近两个点的差商代替导数，收敛阶 1.618，比线性快、比二次慢，是“导数写不出来又想要速度”的折中。它最典型的失败是两点函数值同号且接近，差商趋近于零，下一步直接外推到无穷远。工程上的做法是加一个区间保护：一旦新点跑出已知变号区间，就退回区间内做二分。

Broyden 是多变量的弦截推广，用秩一更新维护一个近似雅可比B，每次迭代只做一次矩阵向量乘，不重新求导。更新公式是B ← B + (y - B*s)*s'/(s'*s)，其中s是步长、y是残差差。这个“好 Broyden”更新在解附近超线性收敛，但它有个隐蔽问题：B只会沿着迭代方向变准，垂直方向的信息永远是初值带来的，跑上三五十步后B可能严重失真。我的做法是每 10 步用一次真实差分雅可比重置B，代价很小，稳定性提升明显。

不动点迭代x = G(x)是门槛最低的方法，把方程整理成x = x - ω*F(x)就能跑。它的收敛条件是G的谱半径小于 1，等价于|1 - ω*λ| < 1，λ 是雅可比的特征值。这解释了为什么松弛因子 ω 不能乱取：λ 很大时 ω 必须很小，而 ω 小又会让收敛慢得难以忍受。这个方法真正的价值在调试阶段——它每一步都简单到可以手算，用来验证残差函数本身写得对不对。

3.3 全局化三板斧：阻尼、信赖域与 LM

阻尼牛顿在牛顿方向上乘一个步长 α，用线搜索（Armijo 或 Wolfe 条件）保证每步残差范数确实下降。它解决的是“方向对但步长太大”这一类失败，代价是要多算几次函数值。α 从 1 开始反复减半，通常两三次就能满足下降条件，所以额外开销不大。

信赖域把“先定方向再定步长”反过来：先定一个半径 Δ，在这个球内求模型的最优步，再看实际下降和预测下降的比值决定放大还是缩小 Δ。dogleg 是其中的经典实现——把高斯-牛顿步和梯度下降步拼成一条折线，在折线上找最优点。fsolve 的默认算法trust-region-dogleg就是它，特点是几乎没有参数要调，对中等规模方阵很稳。

Levenberg-Marquardt 在正规方程上加一个阻尼项(J'J + λI)d = -J'F：λ 大时退化为梯度下降（稳但慢），λ 小时退化为高斯-牛顿（快但不稳）。它是欠定或超定方程组的首选，尤其适合“残差根本无法降到零”的矛盾方程组——此时它会老老实实给你最小二乘解，而牛顿族会一直在那里抖动。

3.4 同伦、二分与符号：不讲收敛阶的兜底手段

同伦延拓的思路是造一个“容易解的方程”再慢慢变形到目标方程。固定点同伦H(x,t) = F(x) - (1-t)F(x0)在t=0时解就是x0，t=1时就是原方程，把 t 从 0 分 20 步走到 1，每步用上一步的解热启动，等于人工给了一条连续的解路径。它是初值极差时最有效的兜底，代价是函数调用次数翻十几倍。

二分法只对单变量有意义，但它的价值不可替代：只要区间两端变号，就一定收敛，没有“初值不好”这回事。多变量里没有严格对应的“保收敛”方法，能做的是把问题降维——固定其中 n-1 个变量扫一遍，或者沿某条参数化曲线化成一维求根。

符号求解（solve/vpasolve）适合需要高精度或全部根的场景。vpasolve可以指定搜索区间，返回该区间内的一个解；配合'Random', true多次调用可以凑出多个根，但这不是遍历，不能保证找全。

4. 落到 MATLAB 里：统一接口、fsolve 选项与三个手写求解器

十多种方法如果每个都写一套调用方式，对照实验根本做不下去。先定一个统一的函数签名[x, hist] = solver(f, Jfun, x0, opts)，hist是每步的残差范数，后面对比收敛曲线时直接用，不用改代码。

4.1 统一调用壳：让所有方法共用一个签名

function [x, hist] = solve_newton(f, Jfun, x0, opts) % 统一的牛顿类求解器壳：f 残差句柄，Jfun 雅可比句柄，x0 初值 if nargin < 4, opts = struct(); end tol = getopt(opts, 'tol', 1e-10); % 残差范数收敛阈值 maxit = getopt(opts, 'maxit', 100); % 最大迭代次数 x = x0(:); hist = nan(maxit,1); for k = 1:maxit F = f(x); hist(k) = norm(F); if hist(k) < tol, break; end % 先判收敛，再算方向，省一次求导 d = -Jfun(x) \ F; % 左除比 inv(J)*F 稳定且快 x = x + d; end hist = hist(1:min(k, maxit)); end function v = getopt(s, name, def) if isfield(s, name), v = s.(name); else, v = def; end end

Jfun(x) \ F用的是 MATLAB 的左除，内部走 LU 分解并带奇异检测，比inv(J)*F快一倍以上且数值更稳。x = x0(:)强制列向量，避免有人传行向量进来导致norm和矩阵乘法悄悄广播，这类错误在小规模算例里极难发现——残差照样下降，只是下降的是别的东西。

hist预分配成nan(maxit,1)而不是空数组追加，是因为追加会触发反复内存重分配，迭代几百次时能差出可观的耗时；预分配还有一个好处：没跑满的部分留NaN，画收敛曲线时自动断开，不会把没执行的迭代当成 0 误差画进去。

4.2 fsolve 的算法开关与 options 参数

fsolve 是基线，先把它的选项摸清楚，再决定哪些场景要自己写。常用参数和调整时机如下：

一个典型配置长这样：

opts = optimoptions('fsolve', ... 'Algorithm', 'trust-region-dogleg', ... 'Jacobian', Jwrap, ... % 解析雅可比句柄，缩放链式已含 'TolFun', 1e-12, 'TolX', 1e-12, ... % 默认 1e-6 对这类算例太松 'MaxIterations', 200, ... 'MaxFunctionEvaluations', 2000, ... 'Display', 'iter'); [x, fval, exitflag, output] = fsolve(f, x0, opts);

TolFun和TolX同时收紧到 1e-12 是有意的：单收紧一个，另一个会先触发，返回exitflag=2或 3，看起来成功了其实还差得远。MaxFunctionEvaluations给到 2000 而不是默认的100*n=300，是因为延拓和多重启时函数调用增长很快，卡在这个上限上会导致明明能收敛的算例被判失败，这是很常见的一类误判。

4.3 手写阻尼牛顿：Armijo 线搜索那二十行

阻尼牛顿是“默认牛顿不收敛”时的第一张牌，核心就是每步让步长满足充分下降条件：

function [x, hist] = solve_damped_newton(f, Jfun, x0, opts) % 阻尼牛顿 + Armijo 回溯线搜索，步长必须让残差范数充分下降 if nargin < 4, opts = struct(); end tol = getopt(opts,'tol',1e-10); maxit = getopt(opts,'maxit',200); c1 = getopt(opts,'c1',1e-4); amin = getopt(opts,'amin',1e-8); x = x0(:); hist = nan(maxit,1); for k = 1:maxit F = f(x); r = norm(F); hist(k) = r; if r < tol, break; end J = Jfun(x); if rcond(J) < 1e-14 % 条件数保护：J 快奇异就别硬解 warning('第 %d 步雅可比接近奇异，rcond=%.2e', k, rcond(J)); end d = -J \ F; % 牛顿方向 if any(~isfinite(d)), break; end % 方向含 Inf/NaN，立即退出 alpha = 1; % 从整步开始，不满足就减半 while alpha > amin if norm(f(x + alpha*d)) <= (1 - c1*alpha)*r, break; end alpha = alpha * 0.5; end x = x + alpha*d; end hist = hist(1:min(k, maxit)); end

(1 - c1*alpha)*r就是 Armijo 条件：新残差要比“按当前下降率的线性预测”还小一点，c1=1e-4是很宽松的要求，只排除那些几乎没下降的步。alpha从 1 开始减半而不是从小往大加，是因为纯牛顿步在近解处几乎总能满足条件，一开始就砍半反而拖慢二次收敛阶段。

rcond(J) < 1e-14这个检查看着多余，实际上救过很多次：J\F在接近奇异时不会报错，只会返回一个巨大的方向，然后x + alpha*d直接飞到 1e12 量级，残差变成Inf，后面所有判断全部失效，最后只看到一句“没收敛”，完全不知道从哪崩的。加上这行至少能在命令行里看到是第几步、条件数多少。

4.4 Broyden 与弦截：把雅可比换成一个矩阵

Broyden 适合变量数上到几十、解析雅可比写不出来、差分又太贵的场景。核心就一条秩一更新：

function [x, hist] = solve_broyden(f, x0, B0, opts) % Broyden 拟牛顿：用秩一更新维护近似雅可比，每 10 步差分重置一次 if nargin < 4, opts = struct(); end tol = getopt(opts,'tol',1e-10); maxit = getopt(opts,'maxit',200); reset_every = getopt(opts,'reset_every',10); x = x0(:); B = B0; hist = nan(maxit,1); for k = 1:maxit F = f(x); hist(k) = norm(F); if hist(k) < tol, break; end if mod(k, reset_every) == 0 B = fdjac(f, x, 1e-7); % 定期用真实差分校正，防止漂移 end s = -B \ F; % 拟牛顿方向 xn = x + s; y = f(xn) - F; B = B + (y - B*s) * (s') / (s'*s); % 好 Broyden 更新 x = xn; end hist = hist(1:min(k, maxit)); end

B0一般传fdjac(f, x0, 1e-7)得到的初始差分雅可比，也可以传单位阵（收敛慢很多但更稳）。(s'*s)是标量，除之前不用求逆，这一点让 Broyden 单步成本远低于牛顿。reset_every=10是我自己的经验值，纯 Broyden 跑三五十步后方向质量会明显下降，定期重置把这点补回来，代价是每 10 步多 n 次函数调用。

一维情形用弦截法更直接，注意加区间保护：

function [x, hist] = solve_secant(f, a, b, opts) % 一维弦截/试位：a、b 需满足 f(a)*f(b) < 0，失败时退回二分 if nargin < 4, opts = struct(); end tol = getopt(opts,'tol',1e-12); maxit = getopt(opts,'maxit',100); fa = f(a); fb = f(b); if fa*fb > 0, error('区间两端必须异号'); end hist = nan(maxit,1); for k = 1:maxit x = b - fb*(b-a)/(fb-fa); % 试位法插值点 if x <= min(a,b) || x >= max(a,b) % 插值跑出区间就退回二分 x = (a+b)/2; end fx = f(x); hist(k) = abs(fx); if hist(k) < tol, break; end if fa*fx < 0, b = x; fb = fx; else, a = x; fa = fx; end end hist = hist(1:min(k, maxit)); end

那个if x <= min || x >= max的判断是血泪经验：试位法在函数形状偏斜时，插值点可能长期贴在区间一端，收敛速度退化到比二分还慢，甚至因为浮点误差跑出区间导致f(x)定义域错误。加三行保护，最坏情况退化成二分，保证一定收敛。

4.5 同伦延拓：给难算例铺一条台阶

初值差到牛顿、LM 全部失败时，用它。固定点同伦的写法极简，不需要自己实现路径跟踪：

function [x, ok, track] = solve_homotopy(f, x0, steps, opts) % 固定点同伦：H(x,t) = f(x) - (1-t)*f(x0)，t 从 0 走到 1 if nargin < 3 || isempty(steps), steps = 20; end if nargin < 4, opts = optimoptions('fsolve','Display','off'); end x = x0(:); F0 = f(x); ok = true; track = zeros(steps, numel(x)); for k = 1:steps t = k/steps; g = @(z) f(z) - (1-t)*F0; % t=0 时 x0 是精确解，t=1 时是原方程 [xn, ~, flag] = fsolve(g, x, opts); % 上一步的解作热启动 if flag <= 0 % 该步没收敛，缩步重试一次 t2 = (k-0.5)/steps; g2 = @(z) f(z) - (1-t2)*F0; [xn, ~, flag] = fsolve(g2, x, opts); if flag <= 0, ok = false; break; end end x = xn; track(k,:) = x(:); end track = track(1:min(k, steps), :); end

(1-t)*F0这一项保证t=0时x0就是g的精确零点，所以第一段几乎不可能失败；随着 t 增大，方程连续变形到目标形式，每一步的初值都来自上一步的解，天然在收敛域内。steps取 20 是折中：太少则每步变形过大、内层失败，太多则函数调用次数暴涨。失败时缩半步重试（t2那段）比直接返回更有用，因为断裂往往发生在某一小段斜率突变的地方。

track输出了整条解路径，这个不只是调试用——如果路径在某步突然折返或者跳变，说明方程本身在那附近有分支点或者多解，这时候该做的是加密步长，而不是怀疑求解器。

5. 收敛失败排查：五个现场的现象、成因与止损动作

以下五个场景基本覆盖了小规模方程组九成以上的“不收敛”报告，每个都按现象、成因、解决三步走。

5.1 现场一：exitflag 是正的，解却是错的

现象：exitflag=2或 3，output.iterations也不多，但把解代回原方程，残差有 1e-2 量级，甚至解在物理上明显不合理（比如负的长度）。

成因：TolX或TolFun相对变量尺度太松。变量量级 1e-3 时，默认TolX=1e-6相当于只要相对精度 0.1% 就收工；更隐蔽的是残差分量量级不一致，某一个分量天然是 1e-4 量级，TolFun=1e-6实际上只约束了它，别的分量根本没动。

解决：先算后向误差norm(F)/norm(J)，这个数比exitflag可信得多。然后把TolFun、TolX按变量和残差的真实量级重设，或者干脆做 2.1 节的缩放，让所有量回到 O(1)，默认值就够用了。

5.2 现场二：NaN 从第一帧就出现

现象：Display='iter'打开，第一行残差就是NaN或Inf，exitflag=-2.1之类的负值，没有任何有效迭代。

成因：初值落在定义域外。方程里有log(x)、sqrt(x)、acos(x)、x^(-0.5)这类项时，初值给到负数或者零，第一次函数求值就产生复数或NaN，fsolve 拿到NaN后无法判断方向，直接退出。

解决：给初值加保护，或者做变量代换把约束吸收掉——正数变量取x = exp(u)，落在 [-1,1] 的变量取x = sin(u)，这样任意u都能映射回合法域，求解器在无约束空间里随便跑。代价是变换后的雅可比可能在某些点上退化，需要配合前面的缩放一起用。另一个更省事的做法是包一层try，在残差函数内部把越界输入改成边界值，但这会让求解器看到虚假的平台，只适合调试，不适合交付。

5.3 现场三：雅可比条件数 1e12，牛顿方向乱飞

现象：残差在 1e-1 到 1e-2 之间来回震荡，迭代十几步后x突然跳到 1e8 量级，rcond(J)输出1e-16，J\F弹出接近奇异的警告。

成因：变量量纲差异过大，或者方程之间存在近似线性相关。三个变量分别是 1e-3、1、1e6 量级时，雅可比的三列量级差出九个数量级，条件数自然爆炸，J\F算出来的方向在数值上毫无意义。

解决：先做 2.1 节的尺度变换，把变量拉回同一量级，这一步通常能把条件数压下去五六个数量级。变换后仍然病态，说明方程本身有冗余或者接近分支点，此时换 LM 算法（它对J'J加阻尼，天然抗病态），或者干脆减少方程数量、把问题转成最小二乘。

5.4 现场四：只有一组初值，就断定“只有这一个解”

现象：fsolve收敛到一个解，报告写完，评审或者后续计算发现还有另一组解，而且那组才是物理上合理的。

成因：非线性方程组的解本来就不唯一，牛顿族只会收敛到初值所在吸引域内的那个解。一组初值等于只采样了一个点，没有任何理由认为它是唯一解。

解决：做网格多初值扫描，把收敛且残差小于阈值的解收集起来去重。三维以下用三重循环扫几百个点足够了，每次调用几十毫秒到几百毫秒，总时间可控。去重用uniquetol按行比较，阈值取 1e-6 量级，这一步在第 6 章有完整代码。扫描完还要对每个解检查雅可比的行列式，|det(J)|接近零的解是重根或者奇异解，物理意义要单独讨论。

5.5 现场五：符号 solve 返回空，vpasolve 又慢得离谱

现象：solve对一个看起来不复杂的方程组返回Empty sym: 0-by-1，换成vpasolve能出结果，但耗时几十秒甚至几分钟，稍微改个参数又变成空。

成因：多项式次数稍高（比如三次以上）且变量多于两个时，符号解法基本无能为力，会退化成隐式根表示或者直接放弃；vpasolve走的是数值搜索加高精度迭代，每次调用都要重新扫搜索空间，不做初值提示就会慢。

解决：给vpasolve传初值或区间，vpasolve(eqn, x, [0 10])，命中率会高很多；如果只要一个数值解，直接用fsolve加解析雅可比，速度快一两个数量级。真正需要全部根并且变量数多时，符号法不划算，改用多初值数值扫描更实际。另外，多项式单变量情形用roots比任何迭代法都稳——它就是特征值方法，一次算出所有根，包括复数根。

6. 多初值扫描与延拓验证：把“解出来了”变成“解对了”

单一初值跑通只说明存在一条收敛路径，不说明解找全了、更不说明解是对的。我现在的固定流程是三步：网格扫描找全部候选解、对每个候选解做残差和雅可比双重校验、再用延拓确认解的连续性。

扫描的代码很短，但有两个细节决定成败：

opts = optimoptions('fsolve','Algorithm','levenberg-marquardt', ... 'Display','off','TolFun',1e-12,'TolX',1e-12); sols = []; for a = -3:0.5:3 for b = -3:0.5:3 for c = 0.5:0.5:10 % 第三个变量按物理范围取，不用对称区间 [x, ~, ef] = fsolve(f, [a;b;c], opts); if ef > 0 && norm(f(x)) < 1e-9 sols = [sols, x(:)]; % 列向量拼接，后面按行去重 end end end end U = uniquetol(sols', 1e-6, 'ByRows', true); % 去重后的候选解集 for k = 1:size(U,1) Jk = Jwrap(U(k,:)'); fprintf('解 %d: ||F||=%.2e, cond(J)=%.2e\n', k, norm(f(U(k,:)')), cond(Jk)); end

uniquetol的'ByRows', true是关键，默认按元素比较会把本来相同的解判成不同；阈值 1e-6 要比TolX宽松一点，否则同一个解因为末位差异会被重复计数。扫描区间不是随便取的：第三个变量如果有物理下界（比如温度必须为正），就按真实范围取，盲目对称扫描既浪费时间又可能扫出一个物理解释不了的数学解。

cond(J)那一列是第二重校验。||F||小只说明找到零点，cond(J)大说明这个零点是“擦边”过去的——方程在那里几乎是平的，微小的参数扰动就能让解跑很远。这类解能不能用，要看你的实际工况允不允许这么大敏感度，不能只看残差。

第三步用第 4.5 节的solve_homotopy从扫描到的解反向走到t=0，检查路径是否连续。如果路径在中间断裂或者跳变，说明这个解和初值之间不存在连续过渡，它可能是被数值噪声“撞”上的伪解，也可能是真解但孤立的——两种情况处理方式完全不同，必须分清楚。

提示：扫描网格的疏密建议按解的预估间距取，一般取变量范围的 1/10 到 1/20；太密收益递减，太疏会漏掉窄吸引域的解。

从那以后我每次交付前都强制走一遍这个三步验证，宁愿多花两分钟跑扫描，也不想再被“多解只找到一个”这种问题返工。这份 MATLAB 非线性方程组方法包里的每个实现都留了统一的[x, hist]输出，改一改就能接进上面的流程，省掉自己搭框架的时间。希望帮到你。

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