Why doesn't the L in my partial pivot LU decomposition switch rows?

6 次查看（过去 30 天）
Jessica Arroyo
Jessica Arroyo 2023-2-19
回答： Piyush Patil 2023-3-1
Ran in:
Everything else works well besides my L, please help.
A= [1 2 3 1; 4 5 6 2; 7 8 0 4; 0 1 3 1];
GE = GE_partial(A);
L = 4×4
1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1
L = 4×4
1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1
L = 4×4
1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1
%%
function[L,U,P] = GE_partial(A)
[n,n] = size(A); %nxn matrix A
L = eye(n,n);
P = eye(n,n);
U = A;
for j =1:n-1 % looping over the columns of U
%Finding the row width with the largest magnitude
[~, r] = max(abs(U(j:n,j)));
% Calling the index of this row k
k = r + j - 1;
% Swap rows of j and k for U and P
U([j,k],:) = U([k,j],:); %Swaping the first and the third row.
P([j, k],:) = P([k, j],:); % swaping the first ad third row
% Swap rows j and k in L for coulmns 1 to j-1
for c = 1:j-1
L([j,k],c) = L([k,j],c)
end
end
end
  3 个评论
显示 1更早的评论
Torsten
Torsten 2023-2-20
编辑：Torsten 2023-2-20
j=1, k=3 means: your loop to modify L is not entered. Do you understand ?

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Piyush Patil
Piyush Patil 2023-3-1
Hello Jessica,
As @Torsten correctly mentioned, you should check the values of variables j and k in the outer “for” loop.
You will notice that –
For j = 1, value of k = 3 but since value of j is 1, so it won’t enter the for loop that modifies the matrix L
For j = 2, value of k = 2. Therefore, both j and k are pointing to same row and hence swapping of rows won’t happen.
Similarly, for j = 3, value of k = 3 and hence swapping won’t take place.
Now, to add on to this discussion, since you are trying to perform LU Decomposition with partial pivoting, you need to modify matrix L and U by computing multipliers and eliminating the elements below pivot element.
You can do this by adding the following logic to the code –
% Compute multipliers and eliminate elements below pivot
for i=j+1:n
L(i,j) = U(i,j)/U(j,j);
U(i,:) = U(i,:) - L(i,j)*U(j,:);
end
Also, you can replace the for loop that switches rows of matrix L by a single line –
L([j, k], 1:j-1) = L([k, j], 1:j-1);
Both will give the same result.
So, finally your entire function should look like –
function[L,U,P] = GE_partial(A)
[n,n] = size(A); %nxn matrix A
L = eye(n,n);
P = eye(n,n);
U = A;
for j =1:n-1 % looping over the columns of U
%Finding the row width with the largest magnitude
[~, r] = max(abs(U(j:n,j)));
% Calling the index of this row k
k = r + j - 1;
% Swap rows of j and k for U, P, and L
U([j, k],:) = U([k, j],:);
P([j, k],:) = P([k, j],:);
L([j, k], 1:j-1) = L([k, j], 1:j-1);
% Compute multipliers and eliminate elements below pivot
for i=j+1:n
L(i,j) = U(i,j)/U(j,j);
U(i,:) = U(i,:) - L(i,j)*U(j,:);
end
end
end

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Polynomials 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by