默认情况下,MATLAB的lu始终执行旋转.如果您在尝试执行传统的LU分解算法时,例如对角线系数等于0,则在执行高斯消元法创建上三角矩阵U时需要对角系数,因此无法工作除以零误差.需要旋转以确保分解稳定.
但是,如果您可以保证矩阵的对角线系数不为零,那么它非常简单,但您必须自己编写.您所要做的就是在矩阵上执行高斯消元法,并将矩阵缩减为梯形缩减形式.结果减少的梯形形式矩阵是U,而去除高斯消除中L的下三角形部分所需的系数将被放置在下三角形半部以制造U.
假设你的矩阵存储在A中,这样的东西可以工作.记住我在这里假设一个方阵.非旋转LU分解算法的实现放在名为lu_nopivot的MATLAB函数文件中:
function [L, U] = lu_nopivot(A)
n = size(A, 1); % Obtain number of rows (should equal number of columns)
L = eye(n); % Start L off as identity and populate the lower triangular half slowly
for k = 1 : n
% For each row k, access columns from k+1 to the end and divide by
% the diagonal coefficient at A(k ,k)
L(k + 1 : n, k) = A(k + 1 : n, k) / A(k, k);
% For each row k+1 to the end, perform Gaussian elimination
% In the end, A will contain U
for l = k + 1 : n
A(l, :) = A(l, :) - L(l, k) * A(k, :);
end
end
U = A;
end
作为一个运行的例子,假设我们有以下3 x 3矩阵:
>> rng(123)
>> A = randi(10, 3, 3)
A =
7 6 10
3 8 7
3 5 5
运行算法给我们:
>> [L,U] = lu_nopivot(A)
L =
1.0000 0 0
0.4286 1.0000 0
0.4286 0.4474 1.0000
U =
7.0000 6.0000 10.0000
0 5.4286 2.7143
0 0 -0.5000
将L和U相乘得出:
>> L*U
ans =
7 6 10
3 8 7
3 5 5
…这是原始矩阵A.