简介这是一份关于Doolittle分解法MATLAB实现的文档资源面向数值分析课程学生、计算数学爱好者及相关算法开发者帮助读者快速理解并上手矩阵三角分解的编程过程。资源只包含1个Word文档压缩包大小20KB文档内提供了完整的MATLAB函数源码函数接收方阵A和向量b输出下三角矩阵L、上三角矩阵U以及解向量X并包含方阵性质判断、顺序主子式是否全非零的检查以及按行列逐步计算L、U并回代求解的完整流程。代码带有中文注释逻辑清晰读者可直接复制使用也能在此基础上改造成自己的求解函数。文档已有269人学习内容简短但覆盖了从输入校验到结果输出的全部关键环节十分适合线性代数实验、MATLAB课程设计或数值分析复习时查阅。1. Doolittle 分解法matlab程序的实现不选主元的 LU 分解能用到哪里第一次跑 Doolittle 分解法 MATLAB 程序时很多人会认为它不过是高斯消元的另一种写法。真正动手把 A 拆成 L 和 U 之后才会意识到它的价值在于左端项分解一次可以反复用于多个右端项。本资源给出的是一个完整的函数实现输入方阵 A 和右端项 b输出单位下三角矩阵 L、上三角矩阵 U、中间向量 Y 以及最终解 X并内置了顺序主子式检查不符合分解条件时直接在命令行给出中文提示。适合数值分析课程设计、需要封装线性方程组求解的仿真脚本以及想对照 Doolittle 分解公式逐行理解程序的从业者。对于已经用惯A\b的人也可以借这个程序看清不选主元做法在实际计算中的边界。2. Doolittle 分解的算法原理与三层循环结构2.1 为什么 L 对角线上全是 1Doolittle 分解是 LU 分解的一种不带主元实现的约定。它把方阵 A 分解为 ALU其中 L 是单位下三角矩阵U 是普通上三角矩阵。单位下三角意味着 L 的对角线元素全部为 1所以第一列不需要额外求解直接用 a21/a11 就能得到 l21消元过程中的乘数都被记录在 L 的下三角位置。这样 A 的 n² 个有效元素正好对应 L 和 U 中除对角线 1 以外的所有位置存储上和原矩阵等量。程序里Leye(n)就是先把对角线全部初始化为 1之后只填下三角部分Uzeros(n)先占位再填上三角部分。如果换用 Crout 分解则 U 对角线为 1公式和代码的除法位置都会不同所以在引用网上下载的 matlab 文档资料前先确认代码采用的是哪一种三角分解约定。对后面回代过程也有影响。求解 Lyb 时因为 L 对角线为 1理论上 Y(1)b(1)不需要做除法。原程序为了保持循环统一在回代第一段里仍然执行/L(i,i)由于 L(i,i)1结果不受影响。这个细节不影响正确性但阅读时容易让人误以为 L 不是单位下三角。2.2 分解公式与三层循环的对应关系Doolittle 分解的核心公式可以写成两行。第一步初始化 U 的第一行和 L 的第一列u1j a1j, j1,…,nli1 ai1 / u11, i2,…,n之后对 k2,…,n逐行计算 U 的第 k 行再逐列计算 L 的第 k 列u(kj) a(kj) − ∑(r1 to k−1) l(kr)u(rj), jk,…,nl(mk) (a(mk) − ∑(r1 to k−1) l(mr)u(rk)) / u(kk), mk1,…,n程序里的三层循环就是按照这两套公式逐项累加。看下面这段与资源中实际处理方式一致的核心代码for k 2:n % 先计算 U 的第 k 行 for j k:n s 0; for r 1:k-1 s s L(k,r) * U(r,j); end U(k,j) A1(k,j) - s; end % 再计算 L 的第 k 列 for m k1:n s 0; for r 1:k-1 s s L(m,r) * U(r,k); end L(m,k) (A1(m,k) - s) / U(k,k); end end这段代码中s 是累加器分别保存公式中的求和项。第一个内层循环的 j 从 k 开始因为 U 的第 k 行前 k−1 列已经在之前的计算中确定为 0不需要重写。第二个内层循环的 m 从 k1 开始因为 L 的对角线已经固定为 1且第 k 列对角线位置由单位下三角约定决定。注意U(k,k)是在前一个内层循环中已经算出的所以计算L(m,k)时可以安全使用它作为除数。如果U(k,k)恰好为 0 或小到接近 0这里就会出现除零或数值放大这是不选主元分解的天然弱点。2.3 det 顺序主子式检查的数学意义Doolittle 分解能够唯一进行的充分必要条件是 A 的所有顺序主子式都不为 0。程序在分解前用det(A1(1:i,1:i))逐个检查左上角子矩阵只要存在一个为 0就把 YY 置为 1最终提示矩阵不能进行 Doolittle 分解并返回。这个判断在数学定义上是准确的但在数值计算中有两个问题。第一det 的结果是标量它只能判断严格为 0 的情况无法反映接近 0 带来的数值风险第二行列式的计算代价偏高对 n 阶矩阵做 n 次 det 调用累计复杂度约为 O(n⁴)比后续分解本身的 O(n³) 更慢。n 比较小的时候影响不大但如果把这个函数直接套到几百阶的稀疏矩阵上det 检查反而可能成为整个脚本最耗时的部分。更实用的做法是在分解过程中直接监控 U(k,k) 是否接近 0这一点放到最后一部分专门说明。3. 函数落地重命名 doolittle.m、路径配置与参数调用3.1 文件命名与函数名同步原始代码第一行写的是function [L,U,Y,X]qw2014210705(A,b)而注释要求将文件重命名为doolittle.m。MATLAB 规定当文件名与函数名不一致时调用以文件名为准函数名本身会被忽略。为了让命令行调用更直观第一行最好同步改成function [L,U,Y,X] doolittle(A,b)然后把文件保存为doolittle.m。修改后原函数名前缀qw2014210705就不再使用。常见问题有两个一是把文件名存成doolitle.m少写一个字母 o二是保存时从 Word 或网盘中复制出来后文件名带.txt后缀MATLAB 无法识别为函数文件。保存完成后在命令行输入which doolittle如果能看到完整路径说明函数已经被正确识别。3.2 放入 bin 还是加入搜索路径“存入 bin 文件”是压缩包作者自己环境里的习惯并不是 MATLAB 的强制要求。MATLAB 只会从当前工作目录和搜索路径中查找函数。如果你把解压后的文件放在D:\numerical\doolittle\doolittle.m推荐直接执行addpath(D:/numerical/doolittle); savepath;第一行将目录加入当前会话的搜索路径第二行把路径写入pathdef.m这样下次启动 MATLAB 时仍然生效。如果不想永久修改环境也可以直接cd(D:/numerical/doolittle)切到函数所在目录。使用addpath的好处是工作目录可以保留在其他位置同时能调用该目录下的其他数值实验脚本。路径配置后重新执行which doolittle验证返回空结果说明路径配置有问题。3.3 输入输出参数与调用签名函数的输入输出参数设计如下参数维度含义备注An×n输入系数矩阵必须为方阵且顺序主子式非零bn×1右端项应为列向量行数与 A 一致Ln×n单位下三角矩阵对角线元素固定为 1Un×n上三角矩阵与 L 满足 AL*UYn×1Lyb 的中间解前代过程的结果Xn×1原方程组的解向量回代过程的最终结果函数签名要求同时传入 A 和 b。即使只关心分解结果也需要写成[L,U] doolittle(A,b);不能省略 b。如果传入的 b 是行向量后续代码中的b(i,1)会越界因为 MATLAB 矩阵索引是按列优先的行向量的第二行并不存在。因此在实际调用时建议用b(:)把向量统一成列向量。3.4 完整调用示例与验证下面用 Pascal 矩阵做一次完整调用。Pascal 矩阵的顺序主子式全部非零能直接通过程序开头的 det 检查A pascal(4); b [1; 2; 3; 4]; [L, U, Y, X] doolittle(A, b); disp(L ); disp(L); disp(U ); disp(U); disp(Y ); disp(Y); disp(X ); disp(X); fprintf(||A-L*U|| %e\n, norm(A - L*U)); fprintf(||A*X-b|| %e\n, norm(A*X - b));代码中的norm(A-L*U)计算分解误差衡量 L 和 U 相乘后能否还原原始矩阵 Anorm(A*X-b)计算残差衡量计算出的解代入原方程后是否成立。理论上这两个值都应该接近 0。查看 Y 也有意义因为 L 是单位下三角矩阵Y 的第一项应当严格等于 b(1)这是快速判断前代过程是否正确的一种人工检查方法。此时不要拿[L2,U2]lu(A);直接和 L、U 对比因为内置lu默认带列主元返回的矩阵需要经过置换才能对应到原矩阵 A。4. 数值实验从四阶方程组到 Hilbert 矩阵的稳定性对比4.1 使用已知解构造测试用例用未知解测试很难判断程序是对是错。更可靠的做法是先生成 A再指定一个 x_true令 bA*x_true这样最终解就是已知的。测试脚本如下A [4 -1 0 0; -1 4 -1 0; 0 -1 4 -1; 0 0 -1 4]; x_true [1; 2; 3; 4]; b A * x_true; [L, U, Y, X] doolittle(A, b); disp(计算解 X:); disp(X); disp(真实解 x_true:); disp(x_true); fprintf(误差范数: %e\n, norm(X - x_true));这个三对角矩阵来自一维热传导问题的离散化顺序主子式全部非零并且不需要选主元。四阶规模下程序输出的 X 应该与 x_true 一致到浮点精度。如果某个分量差得很远可以对照 X 和 x_true 的数值快速定位是分解过程出错还是回代方向出错。这种构造方式比随机给一个 b 更能体现程序内部的计算逻辑。4.2 残差、分解误差与前代误差的隔离数值验证时不同残差的指向并不相同。下面的表格列了常用检查项检查项命令说明求解残差norm(A*X-b)解代入原方程后的误差理想为 0分解误差norm(A-L*U)反映 L 和 U 是否能还原 A前代误差norm(L*Y-b)隔离 L 和 Y 的配合是否正确回代误差norm(U*X-Y)隔离 U 和 X 的配合是否正确如果norm(A-L*U)很大说明分解过程有误如果它很小但norm(A*X-b)很大问题通常出在回代阶段。norm(L*Y-b)和norm(U*X-Y)分别是前代和回代两个阶段的直接检查。原程序在顺序主子式检查之后直接进行分解中间没有设置任何阶段性输出因此这四条命令可以帮你快速判断错误到底在哪一层。4.3 Hilbert 矩阵压力测试Doolittle 分解的适用边界可以用 Hilbert 矩阵来观察。Hilbert 矩阵的元素是 1/(ij−1)各阶顺序主子式都非零能通过程序开头的 det 检查但条件数随矩阵阶数指数增长。测试脚本如下for n [5, 8, 10] A hilb(n); b ones(n, 1); [L, U, Y, X] doolittle(A, b); fprintf(n%d, cond(A)%.3e, 残差%.3e\n, ... n, cond(A), norm(A*X - b)); endcond(A)是 MATLAB 内置条件数函数。逻辑上条件数越大方程组越敏感舍入误差越容易被放大。n5 时残差通常还能接受n 增大到 10 时U 的主元会因为浮点舍入变得很小程序虽然不报错但解向量的有效数字可能严重丢失。这个实验不是说明函数写错了而是说明不选主元的分解法在数学条件满足时依然可能因数值条件不佳而失效。4.4 复杂度的观察方法如果关心大规模表现可以在调用前后加tic和tocn 100; A pascal(n); b ones(n, 1); tic; [L, U, Y, X] doolittle(A, b); toc;但更值得关注的不是单次运行时间而是 det 检查在总时间里的占比。n100 时循环内进行 n 次行列式计算累计计算量明显高于后续分解本身。若测试时发现 n100 还很快、n300 突然变慢优先怀疑 det 检查。想改进的话可以直接把最前面的 det 循环删掉改成在分解过程中判断 U(k,k) 是否为 0具体思路放到最后一部分。5. 验证 Doolittle 分解程序的三种自检方法与主元扩展思路5.1 对照内置 lu 做整体验证最可靠的参照是 MATLAB 内置lu但要注意它默认带列主元返回的 L、U 和 P 满足 PALU。比较脚本如下A pascal(6); b ones(6, 1); [L, U, Y, X] doolittle(A, b); [L2, U2, P] lu(A); X_builtin U2 \ (L2 \ (P * b)); fprintf(与内置 lu 的解误差: %e\n, norm(X - X_builtin));内置lu返回的 L2 是经过行交换后的下三角矩阵不能直接把L2*U2和 A 相减。需要先用 P 对 b 做行置换得到P*b再依次左除 L2 和 U2。如果这个误差在 1e-10 级别说明自定义的 Doolittle 分解在可分解矩阵上的求解行为和内置函数一致。5.2 用条件数判断何时不该用这个资源程序开头的 det 检查只能判断严格奇异的情况无法识别接近奇异。这时候可以在调用前先看条件数A hilb(10); fprintf(条件数 %e\n, cond(A));当cond(A)大于 1e12 时即使程序能跑完解的有效数字也可能只有一两位。这属于算法选型问题不是修改几个小数位能解决的。实际工程中如果矩阵来自数值仿真建议先运行condest(A)这类估计函数做快速判断条件数过大时直接改用A\b或带主元的lu(A)不要在本资源上继续调整。5.3 把 det 检查替换为 U(k,k) 监控如果保留 Doolittle 分解但想让它更实用一个低成本的改造是把 det 循环删掉在每一轮算出 U(k,k) 后立即判断它是否接近 0。改造后的关键片段如下for k 2:n for j k:n s 0; for r 1:k-1 s s L(k,r) * U(r,j); end U(k,j) A1(k,j) - s; end if abs(U(k,k)) 1e-12 warning(第 %d 个主元接近零建议改用列主元 lu, k); return; end for m k1:n s 0; for r 1:k-1 s s L(m,r) * U(r,k); end L(m,k) (A1(m,k) - s) / U(k,k); end end这段代码中U(k,k)在前一个 j 循环结束时已经算出判断它是否接近 0比一开始计算全部顺序主子式更快。用warning而不是disp可以携带具体的 k 值方便定位是哪一层消元出现问题。真正的列主元扩展需要在主元过小时把当前行与后续行交换同时调整 L 中已计算的前 k−1 列和 b 的分量改造量比监控U(k,k)更大。先把上述验证方法跑通再考虑主元扩展会更稳妥。本文还有配套的精品资源点击获取 Doolittle分解法MATLAB实现：不选主元LU分解的算法与稳定性分析
简介这是一份关于Doolittle分解法MATLAB实现的文档资源面向数值分析课程学生、计算数学爱好者及相关算法开发者帮助读者快速理解并上手矩阵三角分解的编程过程。资源只包含1个Word文档压缩包大小20KB文档内提供了完整的MATLAB函数源码函数接收方阵A和向量b输出下三角矩阵L、上三角矩阵U以及解向量X并包含方阵性质判断、顺序主子式是否全非零的检查以及按行列逐步计算L、U并回代求解的完整流程。代码带有中文注释逻辑清晰读者可直接复制使用也能在此基础上改造成自己的求解函数。文档已有269人学习内容简短但覆盖了从输入校验到结果输出的全部关键环节十分适合线性代数实验、MATLAB课程设计或数值分析复习时查阅。1. Doolittle 分解法matlab程序的实现不选主元的 LU 分解能用到哪里第一次跑 Doolittle 分解法 MATLAB 程序时很多人会认为它不过是高斯消元的另一种写法。真正动手把 A 拆成 L 和 U 之后才会意识到它的价值在于左端项分解一次可以反复用于多个右端项。本资源给出的是一个完整的函数实现输入方阵 A 和右端项 b输出单位下三角矩阵 L、上三角矩阵 U、中间向量 Y 以及最终解 X并内置了顺序主子式检查不符合分解条件时直接在命令行给出中文提示。适合数值分析课程设计、需要封装线性方程组求解的仿真脚本以及想对照 Doolittle 分解公式逐行理解程序的从业者。对于已经用惯A\b的人也可以借这个程序看清不选主元做法在实际计算中的边界。2. Doolittle 分解的算法原理与三层循环结构2.1 为什么 L 对角线上全是 1Doolittle 分解是 LU 分解的一种不带主元实现的约定。它把方阵 A 分解为 ALU其中 L 是单位下三角矩阵U 是普通上三角矩阵。单位下三角意味着 L 的对角线元素全部为 1所以第一列不需要额外求解直接用 a21/a11 就能得到 l21消元过程中的乘数都被记录在 L 的下三角位置。这样 A 的 n² 个有效元素正好对应 L 和 U 中除对角线 1 以外的所有位置存储上和原矩阵等量。程序里Leye(n)就是先把对角线全部初始化为 1之后只填下三角部分Uzeros(n)先占位再填上三角部分。如果换用 Crout 分解则 U 对角线为 1公式和代码的除法位置都会不同所以在引用网上下载的 matlab 文档资料前先确认代码采用的是哪一种三角分解约定。对后面回代过程也有影响。求解 Lyb 时因为 L 对角线为 1理论上 Y(1)b(1)不需要做除法。原程序为了保持循环统一在回代第一段里仍然执行/L(i,i)由于 L(i,i)1结果不受影响。这个细节不影响正确性但阅读时容易让人误以为 L 不是单位下三角。2.2 分解公式与三层循环的对应关系Doolittle 分解的核心公式可以写成两行。第一步初始化 U 的第一行和 L 的第一列u1j a1j, j1,…,nli1 ai1 / u11, i2,…,n之后对 k2,…,n逐行计算 U 的第 k 行再逐列计算 L 的第 k 列u(kj) a(kj) − ∑(r1 to k−1) l(kr)u(rj), jk,…,nl(mk) (a(mk) − ∑(r1 to k−1) l(mr)u(rk)) / u(kk), mk1,…,n程序里的三层循环就是按照这两套公式逐项累加。看下面这段与资源中实际处理方式一致的核心代码for k 2:n % 先计算 U 的第 k 行 for j k:n s 0; for r 1:k-1 s s L(k,r) * U(r,j); end U(k,j) A1(k,j) - s; end % 再计算 L 的第 k 列 for m k1:n s 0; for r 1:k-1 s s L(m,r) * U(r,k); end L(m,k) (A1(m,k) - s) / U(k,k); end end这段代码中s 是累加器分别保存公式中的求和项。第一个内层循环的 j 从 k 开始因为 U 的第 k 行前 k−1 列已经在之前的计算中确定为 0不需要重写。第二个内层循环的 m 从 k1 开始因为 L 的对角线已经固定为 1且第 k 列对角线位置由单位下三角约定决定。注意U(k,k)是在前一个内层循环中已经算出的所以计算L(m,k)时可以安全使用它作为除数。如果U(k,k)恰好为 0 或小到接近 0这里就会出现除零或数值放大这是不选主元分解的天然弱点。2.3 det 顺序主子式检查的数学意义Doolittle 分解能够唯一进行的充分必要条件是 A 的所有顺序主子式都不为 0。程序在分解前用det(A1(1:i,1:i))逐个检查左上角子矩阵只要存在一个为 0就把 YY 置为 1最终提示矩阵不能进行 Doolittle 分解并返回。这个判断在数学定义上是准确的但在数值计算中有两个问题。第一det 的结果是标量它只能判断严格为 0 的情况无法反映接近 0 带来的数值风险第二行列式的计算代价偏高对 n 阶矩阵做 n 次 det 调用累计复杂度约为 O(n⁴)比后续分解本身的 O(n³) 更慢。n 比较小的时候影响不大但如果把这个函数直接套到几百阶的稀疏矩阵上det 检查反而可能成为整个脚本最耗时的部分。更实用的做法是在分解过程中直接监控 U(k,k) 是否接近 0这一点放到最后一部分专门说明。3. 函数落地重命名 doolittle.m、路径配置与参数调用3.1 文件命名与函数名同步原始代码第一行写的是function [L,U,Y,X]qw2014210705(A,b)而注释要求将文件重命名为doolittle.m。MATLAB 规定当文件名与函数名不一致时调用以文件名为准函数名本身会被忽略。为了让命令行调用更直观第一行最好同步改成function [L,U,Y,X] doolittle(A,b)然后把文件保存为doolittle.m。修改后原函数名前缀qw2014210705就不再使用。常见问题有两个一是把文件名存成doolitle.m少写一个字母 o二是保存时从 Word 或网盘中复制出来后文件名带.txt后缀MATLAB 无法识别为函数文件。保存完成后在命令行输入which doolittle如果能看到完整路径说明函数已经被正确识别。3.2 放入 bin 还是加入搜索路径“存入 bin 文件”是压缩包作者自己环境里的习惯并不是 MATLAB 的强制要求。MATLAB 只会从当前工作目录和搜索路径中查找函数。如果你把解压后的文件放在D:\numerical\doolittle\doolittle.m推荐直接执行addpath(D:/numerical/doolittle); savepath;第一行将目录加入当前会话的搜索路径第二行把路径写入pathdef.m这样下次启动 MATLAB 时仍然生效。如果不想永久修改环境也可以直接cd(D:/numerical/doolittle)切到函数所在目录。使用addpath的好处是工作目录可以保留在其他位置同时能调用该目录下的其他数值实验脚本。路径配置后重新执行which doolittle验证返回空结果说明路径配置有问题。3.3 输入输出参数与调用签名函数的输入输出参数设计如下参数维度含义备注An×n输入系数矩阵必须为方阵且顺序主子式非零bn×1右端项应为列向量行数与 A 一致Ln×n单位下三角矩阵对角线元素固定为 1Un×n上三角矩阵与 L 满足 AL*UYn×1Lyb 的中间解前代过程的结果Xn×1原方程组的解向量回代过程的最终结果函数签名要求同时传入 A 和 b。即使只关心分解结果也需要写成[L,U] doolittle(A,b);不能省略 b。如果传入的 b 是行向量后续代码中的b(i,1)会越界因为 MATLAB 矩阵索引是按列优先的行向量的第二行并不存在。因此在实际调用时建议用b(:)把向量统一成列向量。3.4 完整调用示例与验证下面用 Pascal 矩阵做一次完整调用。Pascal 矩阵的顺序主子式全部非零能直接通过程序开头的 det 检查A pascal(4); b [1; 2; 3; 4]; [L, U, Y, X] doolittle(A, b); disp(L ); disp(L); disp(U ); disp(U); disp(Y ); disp(Y); disp(X ); disp(X); fprintf(||A-L*U|| %e\n, norm(A - L*U)); fprintf(||A*X-b|| %e\n, norm(A*X - b));代码中的norm(A-L*U)计算分解误差衡量 L 和 U 相乘后能否还原原始矩阵 Anorm(A*X-b)计算残差衡量计算出的解代入原方程后是否成立。理论上这两个值都应该接近 0。查看 Y 也有意义因为 L 是单位下三角矩阵Y 的第一项应当严格等于 b(1)这是快速判断前代过程是否正确的一种人工检查方法。此时不要拿[L2,U2]lu(A);直接和 L、U 对比因为内置lu默认带列主元返回的矩阵需要经过置换才能对应到原矩阵 A。4. 数值实验从四阶方程组到 Hilbert 矩阵的稳定性对比4.1 使用已知解构造测试用例用未知解测试很难判断程序是对是错。更可靠的做法是先生成 A再指定一个 x_true令 bA*x_true这样最终解就是已知的。测试脚本如下A [4 -1 0 0; -1 4 -1 0; 0 -1 4 -1; 0 0 -1 4]; x_true [1; 2; 3; 4]; b A * x_true; [L, U, Y, X] doolittle(A, b); disp(计算解 X:); disp(X); disp(真实解 x_true:); disp(x_true); fprintf(误差范数: %e\n, norm(X - x_true));这个三对角矩阵来自一维热传导问题的离散化顺序主子式全部非零并且不需要选主元。四阶规模下程序输出的 X 应该与 x_true 一致到浮点精度。如果某个分量差得很远可以对照 X 和 x_true 的数值快速定位是分解过程出错还是回代方向出错。这种构造方式比随机给一个 b 更能体现程序内部的计算逻辑。4.2 残差、分解误差与前代误差的隔离数值验证时不同残差的指向并不相同。下面的表格列了常用检查项检查项命令说明求解残差norm(A*X-b)解代入原方程后的误差理想为 0分解误差norm(A-L*U)反映 L 和 U 是否能还原 A前代误差norm(L*Y-b)隔离 L 和 Y 的配合是否正确回代误差norm(U*X-Y)隔离 U 和 X 的配合是否正确如果norm(A-L*U)很大说明分解过程有误如果它很小但norm(A*X-b)很大问题通常出在回代阶段。norm(L*Y-b)和norm(U*X-Y)分别是前代和回代两个阶段的直接检查。原程序在顺序主子式检查之后直接进行分解中间没有设置任何阶段性输出因此这四条命令可以帮你快速判断错误到底在哪一层。4.3 Hilbert 矩阵压力测试Doolittle 分解的适用边界可以用 Hilbert 矩阵来观察。Hilbert 矩阵的元素是 1/(ij−1)各阶顺序主子式都非零能通过程序开头的 det 检查但条件数随矩阵阶数指数增长。测试脚本如下for n [5, 8, 10] A hilb(n); b ones(n, 1); [L, U, Y, X] doolittle(A, b); fprintf(n%d, cond(A)%.3e, 残差%.3e\n, ... n, cond(A), norm(A*X - b)); endcond(A)是 MATLAB 内置条件数函数。逻辑上条件数越大方程组越敏感舍入误差越容易被放大。n5 时残差通常还能接受n 增大到 10 时U 的主元会因为浮点舍入变得很小程序虽然不报错但解向量的有效数字可能严重丢失。这个实验不是说明函数写错了而是说明不选主元的分解法在数学条件满足时依然可能因数值条件不佳而失效。4.4 复杂度的观察方法如果关心大规模表现可以在调用前后加tic和tocn 100; A pascal(n); b ones(n, 1); tic; [L, U, Y, X] doolittle(A, b); toc;但更值得关注的不是单次运行时间而是 det 检查在总时间里的占比。n100 时循环内进行 n 次行列式计算累计计算量明显高于后续分解本身。若测试时发现 n100 还很快、n300 突然变慢优先怀疑 det 检查。想改进的话可以直接把最前面的 det 循环删掉改成在分解过程中判断 U(k,k) 是否为 0具体思路放到最后一部分。5. 验证 Doolittle 分解程序的三种自检方法与主元扩展思路5.1 对照内置 lu 做整体验证最可靠的参照是 MATLAB 内置lu但要注意它默认带列主元返回的 L、U 和 P 满足 PALU。比较脚本如下A pascal(6); b ones(6, 1); [L, U, Y, X] doolittle(A, b); [L2, U2, P] lu(A); X_builtin U2 \ (L2 \ (P * b)); fprintf(与内置 lu 的解误差: %e\n, norm(X - X_builtin));内置lu返回的 L2 是经过行交换后的下三角矩阵不能直接把L2*U2和 A 相减。需要先用 P 对 b 做行置换得到P*b再依次左除 L2 和 U2。如果这个误差在 1e-10 级别说明自定义的 Doolittle 分解在可分解矩阵上的求解行为和内置函数一致。5.2 用条件数判断何时不该用这个资源程序开头的 det 检查只能判断严格奇异的情况无法识别接近奇异。这时候可以在调用前先看条件数A hilb(10); fprintf(条件数 %e\n, cond(A));当cond(A)大于 1e12 时即使程序能跑完解的有效数字也可能只有一两位。这属于算法选型问题不是修改几个小数位能解决的。实际工程中如果矩阵来自数值仿真建议先运行condest(A)这类估计函数做快速判断条件数过大时直接改用A\b或带主元的lu(A)不要在本资源上继续调整。5.3 把 det 检查替换为 U(k,k) 监控如果保留 Doolittle 分解但想让它更实用一个低成本的改造是把 det 循环删掉在每一轮算出 U(k,k) 后立即判断它是否接近 0。改造后的关键片段如下for k 2:n for j k:n s 0; for r 1:k-1 s s L(k,r) * U(r,j); end U(k,j) A1(k,j) - s; end if abs(U(k,k)) 1e-12 warning(第 %d 个主元接近零建议改用列主元 lu, k); return; end for m k1:n s 0; for r 1:k-1 s s L(m,r) * U(r,k); end L(m,k) (A1(m,k) - s) / U(k,k); end end这段代码中U(k,k)在前一个 j 循环结束时已经算出判断它是否接近 0比一开始计算全部顺序主子式更快。用warning而不是disp可以携带具体的 k 值方便定位是哪一层消元出现问题。真正的列主元扩展需要在主元过小时把当前行与后续行交换同时调整 L 中已计算的前 k−1 列和 b 的分量改造量比监控U(k,k)更大。先把上述验证方法跑通再考虑主元扩展会更稳妥。本文还有配套的精品资源点击获取