简介面向数值计算、科学计算与数学建模学习者的一套 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 的默认TolFun1e-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/dys的取法有讲究按解的预估量级取不知道就先跑一次看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 表示步长小于TolX3 表示残差变化小于TolFun4 表示搜索方向幅值过小0 表示迭代超限负值表示失败。问题在于 2、3、4 这三个“看起来成功”的返回值只说明“不再变化”不说明“变对了”。变量量级 1e6 时默认TolX1e-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. 十多种方法怎么选从牛顿的二次收敛到同伦的兜底方法多不代表都要用真正需要记住的是每种方法的“代价”和“失效条件”。下面这张表是按实际调试顺序排的从上往下先试最省事的不行再往下走。方法每步需要 J局部收敛阶单步成本典型适用常见翻车点牛顿法是2二次高初值在解附近、J 可逆J 奇异或病态时方向爆炸修正牛顿固定 J只用一次1线性低J 计算极贵、解附近够近远离解时收敛慢到不可接受弦截法一维否1.618低单变量、导数写不出来初值两点同侧导致外推跑飞Broyden 拟牛顿否秩一更新超线性中变量多、J 难求更新矩阵逐渐失真需重启阻尼牛顿 线搜索是2近解中高初值一般、需要全局收敛性线搜索参数过严导致步长被砍没信赖域 dogleg是2近解高方阵、J 可靠fsolve 默认信赖域半径收缩到 0 提前退出Levenberg-Marquardt是1~2 自适应高欠定/超定、J 奇异阻尼因子初值不当前期几乎不动高斯-牛顿最小二乘是2残差为零时中矛盾方程组求最小残差残差不为零时退化为线性收敛不动点迭代 xG(x)否1线性极低方程天然可整理成压缩映射谱半径 ≥1 时直接发散二分 / 试位一维否1二分极低只需一个根、要求必收敛只适用于变号区间多维无对应同伦延拓是内层取决于内层极高初值极差、多解、强非线性步长过大跟踪路径断裂vpasolve / 符号自动—高需要高精度或全部根高次多变量时返回空或极慢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 在正规方程上加一个阻尼项(JJ λI)d -JFλ 大时退化为梯度下降稳但慢λ 小时退化为高斯-牛顿快但不稳。它是欠定或超定方程组的首选尤其适合“残差根本无法降到零”的矛盾方程组——此时它会老老实实给你最小二乘解而牛顿族会一直在那里抖动。3.4 同伦、二分与符号不讲收敛阶的兜底手段同伦延拓的思路是造一个“容易解的方程”再慢慢变形到目标方程。固定点同伦H(x,t) F(x) - (1-t)F(x0)在t0时解就是x0t1时就是原方程把 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 endJfun(x) \ F用的是 MATLAB 的左除内部走 LU 分解并带奇异检测比inv(J)*F快一倍以上且数值更稳。x x0(:)强制列向量避免有人传行向量进来导致norm和矩阵乘法悄悄广播这类错误在小规模算例里极难发现——残差照样下降只是下降的是别的东西。hist预分配成nan(maxit,1)而不是空数组追加是因为追加会触发反复内存重分配迭代几百次时能差出可观的耗时预分配还有一个好处没跑满的部分留NaN画收敛曲线时自动断开不会把没执行的迭代当成 0 误差画进去。4.2 fsolve 的算法开关与 options 参数fsolve 是基线先把它的选项摸清楚再决定哪些场景要自己写。常用参数和调整时机如下选项常用取值什么时候改Algorithmtrust-region-dogleg默认/trust-region/levenberg-marquardt方阵且 J 可逆用默认J 接近奇异用 LM矛盾方程用 trust-regionJacobianoff/on/Jfun有解析式就传句柄省掉 n 次差分调用JacobPattern稀疏逻辑矩阵大规模稀疏问题只对非零位置做差分FiniteDifferenceStepSizesqrt(eps)≈1.5e-8残差含开方/指数时放大到 1e-6 量级MaxIterations400默认延拓内层压到 30~50避免单步死磕MaxFunctionEvaluations100*n默认差分雅可比时真实消耗是 n 倍要留余量TolFun1e-6默认残差量级本身小于 1e-3 时必须收紧TolX1e-6默认变量量级超过 1e4 时按比例放大ScaleProblemnone/jacobian变量量纲差三个数量级以上时开jacobianDisplayoff/iter/final排查时开iter看残差是否单调一个典型配置长这样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 是有意的单收紧一个另一个会先触发返回exitflag2或 3看起来成功了其实还差得远。MaxFunctionEvaluations给到 2000 而不是默认的100*n300是因为延拓和多重启时函数调用增长很快卡在这个上限上会导致明明能收敛的算例被判失败这是很常见的一类误判。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 条件新残差要比“按当前下降率的线性预测”还小一点c11e-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)); endB0一般传fdjac(f, x0, 1e-7)得到的初始差分雅可比也可以传单位阵收敛慢很多但更稳。(s*s)是标量除之前不用求逆这一点让 Broyden 单步成本远低于牛顿。reset_every10是我自己的经验值纯 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 (ab)/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; % t0 时 x0 是精确解t1 时是原方程 [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这一项保证t0时x0就是g的精确零点所以第一段几乎不可能失败随着 t 增大方程连续变形到目标形式每一步的初值都来自上一步的解天然在收敛域内。steps取 20 是折中太少则每步变形过大、内层失败太多则函数调用次数暴涨。失败时缩半步重试t2那段比直接返回更有用因为断裂往往发生在某一小段斜率突变的地方。track输出了整条解路径这个不只是调试用——如果路径在某步突然折返或者跳变说明方程本身在那附近有分支点或者多解这时候该做的是加密步长而不是怀疑求解器。5. 收敛失败排查五个现场的现象、成因与止损动作以下五个场景基本覆盖了小规模方程组九成以上的“不收敛”报告每个都按现象、成因、解决三步走。5.1 现场一exitflag 是正的解却是错的现象exitflag2或 3output.iterations也不多但把解代回原方程残差有 1e-2 量级甚至解在物理上明显不合理比如负的长度。成因TolX或TolFun相对变量尺度太松。变量量级 1e-3 时默认TolX1e-6相当于只要相对精度 0.1% 就收工更隐蔽的是残差分量量级不一致某一个分量天然是 1e-4 量级TolFun1e-6实际上只约束了它别的分量根本没动。解决先算后向误差norm(F)/norm(J)这个数比exitflag可信得多。然后把TolFun、TolX按变量和残差的真实量级重设或者干脆做 2.1 节的缩放让所有量回到 O(1)默认值就够用了。5.2 现场二NaN 从第一帧就出现现象Displayiter打开第一行残差就是NaN或Infexitflag-2.1之类的负值没有任何有效迭代。成因初值落在定义域外。方程里有log(x)、sqrt(x)、acos(x)、x^(-0.5)这类项时初值给到负数或者零第一次函数求值就产生复数或NaNfsolve 拿到NaN后无法判断方向直接退出。解决给初值加保护或者做变量代换把约束吸收掉——正数变量取x exp(u)落在 [-1,1] 的变量取x sin(u)这样任意u都能映射回合法域求解器在无约束空间里随便跑。代价是变换后的雅可比可能在某些点上退化需要配合前面的缩放一起用。另一个更省事的做法是包一层try在残差函数内部把越界输入改成边界值但这会让求解器看到虚假的平台只适合调试不适合交付。5.3 现场三雅可比条件数 1e12牛顿方向乱飞现象残差在 1e-1 到 1e-2 之间来回震荡迭代十几步后x突然跳到 1e8 量级rcond(J)输出1e-16J\F弹出接近奇异的警告。成因变量量纲差异过大或者方程之间存在近似线性相关。三个变量分别是 1e-3、1、1e6 量级时雅可比的三列量级差出九个数量级条件数自然爆炸J\F算出来的方向在数值上毫无意义。解决先做 2.1 节的尺度变换把变量拉回同一量级这一步通常能把条件数压下去五六个数量级。变换后仍然病态说明方程本身有冗余或者接近分支点此时换 LM 算法它对JJ加阻尼天然抗病态或者干脆减少方程数量、把问题转成最小二乘。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)); enduniquetol的ByRows, true是关键默认按元素比较会把本来相同的解判成不同阈值 1e-6 要比TolX宽松一点否则同一个解因为末位差异会被重复计数。扫描区间不是随便取的第三个变量如果有物理下界比如温度必须为正就按真实范围取盲目对称扫描既浪费时间又可能扫出一个物理解释不了的数学解。cond(J)那一列是第二重校验。||F||小只说明找到零点cond(J)大说明这个零点是“擦边”过去的——方程在那里几乎是平的微小的参数扰动就能让解跑很远。这类解能不能用要看你的实际工况允不允许这么大敏感度不能只看残差。第三步用第 4.5 节的solve_homotopy从扫描到的解反向走到t0检查路径是否连续。如果路径在中间断裂或者跳变说明这个解和初值之间不存在连续过渡它可能是被数值噪声“撞”上的伪解也可能是真解但孤立的——两种情况处理方式完全不同必须分清楚。提示扫描网格的疏密建议按解的预估间距取一般取变量范围的 1/10 到 1/20太密收益递减太疏会漏掉窄吸引域的解。从那以后我每次交付前都强制走一遍这个三步验证宁愿多花两分钟跑扫描也不想再被“多解只找到一个”这种问题返工。这份 MATLAB 非线性方程组方法包里的每个实现都留了统一的[x, hist]输出改一改就能接进上面的流程省掉自己搭框架的时间。希望帮到你。本文还有配套的精品资源点击获取