matlab 为什么透视代码会给LU分解带来失败？

sczxawaw 于 2023-04-06 发布在 Matlab

关注(0)|答案(1)|浏览(214)

在我发言之前，我是和翻译一起写的，而不是一个会说英语的人。因此，即使我的英语很糟糕，我也希望你能理解我。
在学习了如何编写一个非旋转的LU分解代码之后，我决定尝试编写一个使用旋转的LU分解代码作为家庭作业。由于我的非旋转代码进行得很顺利，所以我认为制作一个旋转的LU分解代码很容易。
但事情变糟了真的很糟。。

% non-pivoting code
function x = mylu_np(A)
    [m,~] = size(A);
    L = eye(m); U = A(:,:);
    for i = 1:m-1
        for j = i+1:m
            L(j,i) = U(j,i)/U(i,i);
            U(j,i) = 0;
            U(j,i+1:m) = U(j,i+1:m) - L(j,i)*U(i,i+1:m);
        end
    end
    x = [L U];
end

我的代码如下。请告诉我这段代码失败的地方。我不知道为什么这么严重的失败发生，即使我只把在枢轴部分。

function x = mylu_pp(A)
    [m,~] = size(A);
    P = eye(m); L = eye(m); U = A(:,:);
    for i = 1:m-1
        for j = i+1:m
            if abs(U(i,i)) < abs(U(j,i))
                U([i j],:) = U([j i],:);
                P([i j],:) = P([j i],:);
            end
        end
        for j = i+1:m
            L(j,i) = U(j,i)/U(i,i);
            U(j,i) = 0;
            U(j,i+1:m) = U(j,i+1:m) - L(j,i)*U(i,i+1:m);
        end
    end
    x = [P L U];
end

我用4*4矩阵A测试了一下，像这样

A = [0.0001 -5.0300 5.8090 7.8320;
    2.2660 1.9950 1.2120 8.0080;
    8.8500 5.6810 4.5520 1.3020;
    6.7750 -2.2530 2.9080 3.9700];
W = mylu_pp(A)
P = W(1:4,1:4); L = W(1:4,5:8); U = W(1:4,9:12);
P*A - L*U

答案是这样的。非常糟糕。

ans = [0.0000 0.0000 0.0000 0.0000;
    6.7749 4.3489 3.4847 0.9967;
    -2.2659 -7.0250 -1.6521 2.1754
    -4.5090 2.6761 -1.8326 -3.1721]

告诉我是哪部分的问题，欢迎各种帮助，我会提前谢谢你的。

matlab

来源：https://stackoverflow.com/questions/75925305/why-does-pivoting-code-bring-failure-to-my-lu-decomposition

1条答案

按热度按时间

w3nuxt5m1#

问题与矩阵L有关。
当交换矩阵U中的行i和j时，我们还必须交换矩阵L中的行i和j中的元素，但仅交换先前更新的元素（而不是整个行）：

if i > 1
    L([i j], 1:i-1) = L([j i], 1:i-1);
end

我思考这个问题的方式相当于重新命名方程组中的变量。
我不是数学家，我的解释不会很好...
假设我们有方程组：A*v = b。
使用LU分解，我们将方程求解为：L*U*v = b。
假设：

[x
v =  y
     z
     w]

在交换矩阵U中的两行之后，我们必须保持与原始方程组的等价性。
保持等价性的方法是对v中的元素进行重命名。
当将行2与行4交换时，我们可以交换y和w的名称，并得到新的v：

[x
v0 =  w
      z
      y]

与名称交换等价的是：

交换矩阵U中的两行。
交换矩阵P中的两行。

P*A将A的行重新排序为新的“交换顺序”。

我们还必须交换矩阵L的行中的元素，但只交换先前更新的元素（而不是整行）。

交换L的第2行和第4行中的相关元素与交换y和w的名称是等价的。
代码示例：

A = [0.0001 -5.0300 5.8090 7.8320;
     2.2660  1.9950 1.2120 8.0080;
     8.8500  5.6810 4.5520 1.3020;
     6.7750 -2.2530 2.9080 3.9700];

[P, L, U] = mylu_pp(A);
P*A - L*U

function [P, L, U] = mylu_pp(A)
    m = size(A, 1);
    P = eye(m); L = eye(m); U = A;
    for i = 1:m-1
        for j = i+1:m
            if abs(U(i,i)) < abs(U(j,i))
                U([i j],:) = U([j i],:);
                P([i j],:) = P([j i],:);
                
                if i > 1
                    % Swap the relevant elements "1:i-1" in the row of matrix L (only the elements that were updated until now).
                    L([i j], 1:i-1) = L([j i], 1:i-1);
                end
            end
        end
        for j = i+1:m
            L(j,i) = U(j,i)/U(i,i);
            U(j,i) = 0;
            U(j,i+1:m) = U(j,i+1:m) - L(j,i)*U(i,i+1:m);
        end
    end
    %x = [P L U];
end

输出P*A - L*U：

ans =

     0     0     0     0
     0     0     0     0
     0     0     0     0
     0     0     0     0

赞(0）回复(0）举报 2023-04-06

我来回答

matlab 为什么透视代码会给LU分解带来失败？

1条答案

相关问题

热门标签

最新问答