MATLAB LU 分解部分旋转
MATLAB LU Decomposition Partial pivoting
我正在尝试主要基于 LU decomposition with partial pivoting Matlab
来处理我的 lu 分解
function [L,U,P] = lup(A)
n = length(A);
L = eye(n);
U = zeros(n);
P = eye(n);
for k=1:n-1
% find the entry in the left column with the largest abs value (pivot)
[~,r] = max(abs(A(k:end,k)));
r = n-(n-k+1)+r;
A([k r],:) = A([r k],:);
P([k r],:) = P([r k],:);
L([k r],:) = L([r k],:);
% from the pivot down divide by the pivot
L(k+1:n,k) = A(k+1:n,k) / A(k,k);
U(k,1:n) = A(k,1:n);
A(k+1:n,1:n) = A(k+1:n,1:n) - L(k+1:n,k)*A(k,1:n);
end
U(:,end) = A(:,end);
end
它似乎适用于大多数矩阵(等同于 matlab lu 函数),但是以下矩阵似乎会产生不同的结果:
A = [
3 -7 -2 2
-3 5 1 0
6 -4 0 -5
-9 5 -5 12
];
我只是不知道出了什么问题。它似乎在链接 post
中提到的矩阵上运行良好
你非常接近。我总共改了三行
for k=1:n-1
变成了 for k=1:n
我们不做 -1 是因为我们还想用你的方法得到 L(n,n)=u(n,n)/u(n,n)=1
我们把它排除在外
L(k+1:n,k) = A(k+1:n,k) / A(k,k);
变成了 L(k:n,k) = A(k:n,k) / A(k,k);
因为你遗漏了 L(k,k)=A(k,k)/A(k,k)=1
因为 k+1
的变化我们不需要从 L 的单位矩阵开始,因为我们现在在对角线上复制 1,所以 L=eyes(n);
变成了 L=zeros(n);
和完成的代码
function [L,U,P] = lup(A)
% lup factorization with partial pivoting
% [L,U,P] = lup(A) returns unit lower triangular matrix L, upper
% triangular matrix U, and permutation matrix P so that P*A = L*U.
n = length(A);
L = zeros(n);
U = zeros(n);
P = eye(n);
for k=1:n
% find the entry in the left column with the largest abs value (pivot)
[~,r] = max(abs(A(k:end,k)));
r = n-(n-k+1)+r;
A([k r],:) = A([r k],:);
P([k r],:) = P([r k],:);
L([k r],:) = L([r k],:);
% from the pivot down divide by the pivot
L(k:n,k) = A(k:n,k) / A(k,k);
U(k,1:n) = A(k,1:n);
A(k+1:n,1:n) = A(k+1:n,1:n) - L(k+1:n,k)*A(k,1:n);
end
U(:,end) = A(:,end);
end
我正在尝试主要基于 LU decomposition with partial pivoting Matlab
来处理我的 lu 分解function [L,U,P] = lup(A)
n = length(A);
L = eye(n);
U = zeros(n);
P = eye(n);
for k=1:n-1
% find the entry in the left column with the largest abs value (pivot)
[~,r] = max(abs(A(k:end,k)));
r = n-(n-k+1)+r;
A([k r],:) = A([r k],:);
P([k r],:) = P([r k],:);
L([k r],:) = L([r k],:);
% from the pivot down divide by the pivot
L(k+1:n,k) = A(k+1:n,k) / A(k,k);
U(k,1:n) = A(k,1:n);
A(k+1:n,1:n) = A(k+1:n,1:n) - L(k+1:n,k)*A(k,1:n);
end
U(:,end) = A(:,end);
end
它似乎适用于大多数矩阵(等同于 matlab lu 函数),但是以下矩阵似乎会产生不同的结果:
A = [
3 -7 -2 2
-3 5 1 0
6 -4 0 -5
-9 5 -5 12
];
我只是不知道出了什么问题。它似乎在链接 post
中提到的矩阵上运行良好你非常接近。我总共改了三行
for k=1:n-1
变成了 for k=1:n
我们不做 -1 是因为我们还想用你的方法得到 L(n,n)=u(n,n)/u(n,n)=1
我们把它排除在外
L(k+1:n,k) = A(k+1:n,k) / A(k,k);
变成了 L(k:n,k) = A(k:n,k) / A(k,k);
因为你遗漏了 L(k,k)=A(k,k)/A(k,k)=1
因为 k+1
的变化我们不需要从 L 的单位矩阵开始,因为我们现在在对角线上复制 1,所以 L=eyes(n);
变成了 L=zeros(n);
和完成的代码
function [L,U,P] = lup(A)
% lup factorization with partial pivoting
% [L,U,P] = lup(A) returns unit lower triangular matrix L, upper
% triangular matrix U, and permutation matrix P so that P*A = L*U.
n = length(A);
L = zeros(n);
U = zeros(n);
P = eye(n);
for k=1:n
% find the entry in the left column with the largest abs value (pivot)
[~,r] = max(abs(A(k:end,k)));
r = n-(n-k+1)+r;
A([k r],:) = A([r k],:);
P([k r],:) = P([r k],:);
L([k r],:) = L([r k],:);
% from the pivot down divide by the pivot
L(k:n,k) = A(k:n,k) / A(k,k);
U(k,1:n) = A(k,1:n);
A(k+1:n,1:n) = A(k+1:n,1:n) - L(k+1:n,k)*A(k,1:n);
end
U(:,end) = A(:,end);
end