MATLAB LU分解法实现：从线性方程组求解到可复用工程代码

简介：这份资源聚焦数值计算中的经典算法，提供LU分解法求解线性方程组的MATLAB实现，面向正在学习数值分析、计算方法课程的高校学生，以及需要将算法移植到C语言环境的开发者。压缩包共2个文件，包含1个m脚本与1个txt说明文档，整体约1KB，体量轻巧，便于快速阅读与二次修改。脚本以最基础的步骤演示LU分解的完整流程，未引入动态算法，逻辑直白，方便读者逐行理解矩阵分解与回代求解的核心思路；说明文档则对程序用途与运行方式作出补充，降低上手门槛。对于C程序员而言，其中涉及的函数调用与流程组织方式具有参考价值，稍加改写即可迁移到C实现。目前已有2351人学习下载，适合作为课堂练习、算法入门或代码移植的起步材料。

1. LU 分解法的 MATLAB 实现：从线性方程组到可复用代码的落地路径

很多做数值计算的同学第一次接触 LU 分解，都是在求解线性方程组Ax = b的时候。直接求逆矩阵inv(A)*b在课堂上能跑通，但一旦矩阵维度上到几百上千，数值稳定性和计算量都会让你后悔。LU 分解法的 MATLAB 实现，本质上就是把系数矩阵A拆成一个下三角矩阵L和一个上三角矩阵U，让原本一次昂贵的求解变成「一次分解、多次回代」，特别适合那些系数矩阵不变、右端项b反复变化的场景，比如结构力学里的多工况加载、电路网络的多次激励分析、热传导方程的隐式时间推进。这篇文章面向的是已经会用 MATLAB 做基本矩阵运算、但想把 LU 分解真正写进工程代码的从业者和学生，我会从原理、选型、手写实现、内置函数对比、避坑到进阶技巧，把这条路走一遍。

2. LU 分解的数学骨架与 MATLAB 里的三种落地形态

在动手写代码之前，得先把「分解到底在分解什么」这件事说清楚，否则后面调参数、看残差都会变成玄学。

2.1 从高斯消元到 Doolittle 分解：为什么 L 的对角线默认是 1

高斯消元法求解Ax = b的过程，其实就是在对A做一系列初等行变换，把它变成一个上三角矩阵。把这些行变换的「逆操作」按顺序收集起来，就得到了一个单位下三角矩阵L，而对角线之所以是 1，是因为每一步消元只做「减去某一行的倍数」，没有对主元行做缩放。这就是 Doolittle 分解：A = L*U，其中L是单位下三角，U是上三角。

对应的数学表达是：

A(i,j) 消元后：U(i,j) = A(i,j) - sum(L(i,k)*U(k,j)), k = 1..i-1 L(i,j) = (A(i,j) - sum(L(i,k)*U(k,j))) / U(j,j), k = 1..j-1

MATLAB 里实现这个双重循环，最直观的写法是三层嵌套，但那样在解释型语言里会慢得让人怀疑人生。所以实际落地时，我一般会走「向量化外积更新」的路线，把第三层循环压掉。

2.2 三种实现形态：手写循环、向量化、内置 lu

在 MATLAB 里做 LU 分解，常见有三条路：

内置lu底层调的是 LAPACK 的dgetrf，带部分主元选取，数值稳定性和速度都不是手写能比的。但手写版本的价值在于：你能在分解过程中插入自己的逻辑，比如只对某个稀疏结构做特殊处理，或者在做分解的同时统计主元增长情况。选型的原则很简单——如果只是解方程，直接用lu；如果要研究算法行为或做定制化预处理，才自己写。

2.3 部分主元选取：不选主元为什么会在对角占优矩阵上翻车

考虑矩阵：

A = [1e-20, 1; 1, 1];

如果不选主元，第一步消元时用1e-20做除数，乘数会变成1e20，浮点误差被放大到不可接受。部分主元选取就是在每一列消元前，把当前列绝对值最大的元素所在行换到主元位置，用行交换矩阵P记录这些交换，最终得到P*A = L*U。

MATLAB 内置lu默认就是带部分主元的，返回三个输出时[L,U,P] = lu(A)满足L*U = P*A。注意这个顺序，很多人会记反，导致回代时b没有同步置换，结果全错。

3. 手写 LU 分解：向量化版本与完整求解流程

这一章给出一份可以直接抄作业的代码，包含分解和回代两部分，并逐段说明参数和逻辑。

3.1 向量化 Doolittle 分解函数

function [L, U] = lu_doolittle(A) % LU_DOOLITTLE 不带主元选取的 Doolittle 分解 % 输入：A —— n x n 数值矩阵，要求各阶顺序主子式非零 % 输出：L —— 单位下三角矩阵，U —— 上三角矩阵，满足 A = L*U n = size(A, 1); if size(A, 2) ~= n error('输入矩阵必须为方阵'); end L = eye(n); % 单位下三角，对角线预置为 1 U = zeros(n); for k = 1:n % 计算 U 的第 k 行：用 A 的第 k 行减去已算出的 L(k,1:k-1)*U(1:k-1,:) U(k, k:n) = A(k, k:n) - L(k, 1:k-1) * U(1:k-1, k:n); if abs(U(k,k)) < 1e-14 error('主元接近零，矩阵可能需要主元选取，第 %d 步失败', k); end % 计算 L 的第 k 列：用 A 的第 k 列减去已算出的 L(1:k-1,1:k-1)*U(1:k-1,k) if k < n L(k+1:n, k) = (A(k+1:n, k) - L(k+1:n, 1:k-1) * U(1:k-1, k)) / U(k,k); end end end

逻辑说明：外层循环k逐列推进，每一轮先算出U的第k行，再算出L的第k列。U(k, k:n)这一段用向量化方式一次性更新整行，避免了内层对j的循环；L(k+1:n, k)同理，一次性更新整列。参数上，1e-14是我常用的主元阈值，低于它基本可以判定矩阵接近奇异，继续算下去结果没有意义。

3.2 前代与回代：把分解结果用起来

function x = solve_lu(L, U, b) % SOLVE_LU 利用已分解的 L、U 求解 A*x = b % 输入：L、U —— lu_doolittle 的输出；b —— n x 1 右端向量 % 输出：x —— 解向量 n = length(b); % 前代：解 L*y = b，L 是单位下三角 y = zeros(n, 1); for i = 1:n y(i) = b(i) - L(i, 1:i-1) * y(1:i-1); end % 回代：解 U*x = y x = zeros(n, 1); for i = n:-1:1 x(i) = (y(i) - U(i, i+1:n) * x(i+1:n)) / U(i,i); end end

逻辑说明：前代利用L的对角线为 1，不需要除法；回代时U的对角线元素参与除法，所以前面分解阶段必须保证U(k,k)不为零。参数上，b可以是列向量，也可以是多列矩阵，此时y和x的初始化要相应改成zeros(n, size(b,2))，这样一次分解就能处理多个右端项，这正是 LU 分解相对直接求逆的核心优势。

3.3 一个可复现的验证脚本

% 验证脚本：随机矩阵上对比手写 LU 与内置 lu rng(0); % 固定随机种子，保证可复现 n = 200; A = rand(n) + n * eye(n); % 对角占优，保证无主元选取也稳定 b = rand(n, 1); % 手写版本 tic; [L1, U1] = lu_doolittle(A); x1 = solve_lu(L1, U1, b); t1 = toc; % 内置版本 tic; [L2, U2, P2] = lu(A); x2 = U2 \ (L2 \ (P2 * b)); t2 = toc; fprintf('手写残差：%e，耗时：%.4f s\n', norm(A*x1 - b), t1); fprintf('内置残差：%e，耗时：%.4f s\n', norm(A*x2 - b), t2); fprintf('解的差异：%e\n', norm(x1 - x2));

逻辑说明：rng(0)固定随机种子，保证每次运行结果一致，这是做数值实验的基本习惯。A = rand(n) + n*eye(n)构造对角占优矩阵，避免手写版本因主元过小报错。残差norm(A*x - b)是检验解质量最直接的指标，正常情况下应该在1e-10量级以下。耗时对比能直观看到手写与 LAPACK 的差距，200 阶时手写大约慢一个数量级，维度越大差距越明显。

4. 避坑与排查：LU 分解在 MATLAB 里最容易翻车的五个地方

这一章记录的是我在实际项目里踩过的坑，每一条都按「现象 → 原因 → 解决」写，方便你对照排查。

4.1 现象：残差正常但解完全不对

原因：使用[L,U,P] = lu(A)后，回代时忘记对b做置换，直接用了原始的b。内置lu返回的P满足L*U = P*A，所以正确的求解顺序是先算P*b，再前代回代。

解决：写成y = L \ (P*b); x = U \ y;，或者用[L,U,P] = lu(A, 'vector')拿到置换向量，用b(P)做索引置换，后者在多次求解时更快。

4.2 现象：手写分解在某一步报「主元接近零」

原因：矩阵本身接近奇异，或者虽然非奇异但某阶顺序主子式为零，不选主元的 Doolittle 分解无法继续。

解决：改用带部分主元的分解，或者直接调用内置lu。如果必须手写，可以在每一步消元前扫描当前列，找到绝对值最大的行做交换，同时维护一个置换向量。判断矩阵是否真的奇异，用cond(A)看条件数，超过1e12基本就要警惕了。

4.3 现象：大矩阵上手写版本慢到无法接受

原因：MATLAB 的for循环在旧版本里开销很大，即使做了向量化，外层循环次数多时仍然明显。

解决：优先用内置lu；如果确实需要自定义逻辑，把最内层循环彻底向量化，或者用pagefun、dlarray做批量处理。另一个思路是先用sparse转稀疏存储，MATLAB 对稀疏矩阵的 LU 有专门的luinc和ichol预处理，适合大规模稀疏系统。

4.4 现象：复数矩阵分解结果与预期不符

原因：MATLAB 的lu对复数矩阵默认做的是「部分主元」，但转置关系变成共轭转置，手写代码里如果用了'而不是.'，会引入共轭，导致L*U不再等于P*A。

解决：复数场景下统一用.'做非共轭转置，或者直接用内置lu，它对复数矩阵的处理是经过验证的。手写版本在复数域下要额外注意主元选取时比较的是模长，不是实部。

4.5 现象：分解成功但内存爆掉

原因：对大规模稀疏矩阵直接做稠密 LU，L和U的填充元（fill-in）会让存储量急剧膨胀，远超原始矩阵。

解决：用lu(A, 'vector')配合colamd或symamd做列近似最小度重排序，能显著减少填充。MATLAB 里[L,U,P,Q] = lu(S)会自动做重排序，返回的Q是列置换，满足L*U = P*S*Q。如果内存仍然不够，考虑迭代法或luinc做不完全分解作为预条件。

5. 进阶技巧：把 LU 分解嵌进工程代码的几个习惯

走到这里，基本的分解和求解已经能跑通了。最后这一章聊几个我在实际项目里养成的习惯，能让 LU 分解从「能跑」变成「好用」。

第一个习惯是把分解和求解分离。很多新手会把lu和回代写在一个函数里，每次换b都重新分解一遍，这是最大的浪费。正确的做法是分解一次，把L、U、P存下来，后续所有右端项复用。在结构力学里，这意味着一次分解可以处理几十个荷载工况；在隐式时间推进里，只要时间步长不变，系数矩阵不变，分解就能一直复用。

第二个习惯是用decomposition对象管理。MATLAB 从 R2017b 开始提供了decomposition函数，可以创建一个可复用的分解对象：

dA = decomposition(A, 'lu'); % 创建 LU 分解对象 x1 = dA \ b1; % 第一次求解 x2 = dA \ b2; % 复用分解，直接换右端项

这个对象的好处是，MATLAB 会自动选择最优的分解策略，并且在多次求解时避免重复分解。对于需要反复求解的场景，代码更干净，性能也更稳。

第三个习惯是监控主元增长。在数值稳定性要求高的场景，我会在分解过程中记录max(abs(U(k,k))) / max(abs(A(:)))的比值，如果这个比值超过1e3，说明主元增长过快，解的精度可能已经受损。这时候要么改用全主元选取，要么对矩阵做预处理（比如行平衡），把每行的量级拉到同一水平。

第四个习惯是用残差而不是解的差异来判断正确性。很多人验证时喜欢对比手写解和内置解的差异，但两个都不准的时候差异也可能很小。更可靠的做法是算norm(A*x - b) / norm(b)的相对残差，这个指标直接反映解是否满足方程。我一般要求相对残差在1e-12以下才认为通过。

最后一个习惯是对稀疏矩阵优先考虑迭代法。LU 分解再快也是直接法，填充元的问题在大规模稀疏场景下无法回避。如果矩阵维度上万、稀疏度低于 1%，先试试gmres或bicgstab配合ilu预条件，往往比直接 LU 更省内存。判断标准很简单：如果nnz(L) + nnz(U)超过nnz(A)的 10 倍，就该考虑换方法了。

这些习惯没有一条是教科书上会强调的，都是被残差打脸、被内存爆掉、被复用场景逼出来的。我现在的做法是，任何涉及 LU 分解的代码，第一版就按「分解与求解分离 + 残差验证」的模板写，后面无论怎么改都不会乱。希望帮到你。

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