寻求MATLAB/TI-Nspire CX CAS实现上三角对角为1的LU分解方法

你想要的是Crout型LU分解（上三角矩阵U的对角元全为1，下三角矩阵L为非单位下三角），而MATLAB默认的lu函数实现的是Doolittle型分解（下三角矩阵L的对角元全为1，上三角矩阵U为非单位上三角），直接转置矩阵的方法之所以不对，是因为转置后的分解结果无法直接转换回你需要的原矩阵的Crout分解形式。

先验证一下你给出的目标分解是正确的：

L2 = [2 0 0; -0.5 3.5 0; 3 1 -3];
U = [1 1 -3; 0 1 2; 0 0 1];
L2*U

运行后会得到和你的矩阵A完全一致的结果：[2 2 -6; -0.5 3 8.5; 3 4 -10]，说明这个分解是成立的。

下面提供两种可行的解决方法：

方法1：手动实现Crout分解算法 #

直接按照Crout分解的规则编写代码，能精准得到你需要的L和U：

clc
A = [2, 2, -6; -0.5, 3, 8.5; 3, 4, -10];
n = size(A, 1);
L = zeros(n);
U = eye(n); % 初始化U为单位上三角矩阵

% Crout分解核心步骤
for j = 1:n
    % 计算L的第j列元素
    for i = j:n
        L(i,j) = A(i,j) - L(i,1:j-1)*U(1:j-1,j);
    end
    % 计算U的第j行元素（最后一行无需计算，已为[0...0 1]）
    if j < n
        for k = j+1:n
            U(j,k) = (A(j,k) - L(j,1:j-1)*U(1:j-1,k))/L(j,j);
        end
    end
end

disp('L矩阵（对应你的L2）：');
disp(L);
disp('U矩阵（对角元全为1）：');
disp(U);

运行这段代码后，你会得到和目标完全一致的分解结果。

方法2：利用MATLAB默认lu函数转换得到 #

MATLAB默认的lu函数返回[L, U, P]满足P*A = L*U，其中L是单位下三角矩阵，U是非单位上三角矩阵。我们可以把U拆分为对角矩阵D和单位上三角矩阵U1的乘积：U = D * U1（D的对角元就是U的对角元，U1是U除以D的对角元后得到的单位上三角矩阵）。

推导后可得：A = P' * L * D * U1，其中L*D是下三角矩阵（对角元为D的元素），U1是单位上三角矩阵，正好符合你的需求。对应的代码如下：

clc
A = [2, 2, -6; -0.5, 3, 8.5; 3, 4, -10];
[L_default, U_default, P] = lu(A);

% 拆分U_default为对角矩阵D和单位上三角矩阵U1
D = diag(diag(U_default));
U1 = U_default / D; % 等价于inv(D)*U_default，因D是对角矩阵

% 计算目标L2矩阵
L2 = P' * L_default * D;

disp('L2矩阵：');
disp(L2);
disp('U矩阵（对角元全为1）：');
disp(U1);

注意：由于浮点计算精度问题，结果可能会出现极小量级的误差（如1e-16），可以用round函数或调整显示精度来处理。

为什么你之前的转置方法无效？ #

你对A'做LU分解得到A' = L2*U，转置后得到A = U' * L2'，此时U'是下三角矩阵，L2'是上三角矩阵且对角元为1——这其实是A的Doolittle分解，而非你需要的Crout分解，所以结果不符合预期。

内容的提问来源于stack exchange，提问作者user9546999