1

我正在尝试创建一个程序，该程序将方形（n×n）矩阵作为输入，如果它是可逆的，LU 将使用高斯消元法分解矩阵。

这是我的问题：在课堂上，我们了解到最好更改行，以便您的枢轴始终是其列中的最大数字（绝对值）。例如，如果矩阵当时正在A = [1,2;3,4]切换行[3,4;1,2]，那么我们可以继续进行高斯消元。

我的代码适用于不需要行更改的矩阵，但对于需要行更改的矩阵，则不需要。这是我的代码：

function newgauss(A)
    [rows,columns]=size(A);
    P=eye(rows,columns); %P is permutation matrix
    if(det(A)==0) %% determinante is 0 means no single solution
        disp('No solutions or infinite number of solutions')
        return;
    end
    U=A;
    L=eye(rows,columns);
    pivot=1;
    while(pivot<rows)
        max=abs(U(pivot,pivot));
        maxi=0;%%find maximum abs value in column pivot
        for i=pivot+1:rows
            if(abs(U(i,pivot))>max)
                max=abs(U(i,pivot));
                maxi=i;
            end
        end %%if needed then switch
        if(maxi~=0)
            temp=U(pivot,:);
            U(pivot,:)=U(maxi,:);
            U(maxi,:)=temp;
            temp=P(pivot,:);
            P(pivot,:)=P(maxi,:);
            P(maxi,:)=temp;
        end %%Grade the column pivot using gauss elimination
        for i=pivot+1:rows
            num=U(i,pivot)/U(pivot,pivot);
            U(i,:)=U(i,:)-num*U(pivot,:);
            L(i,pivot)=num;
        end
        pivot=pivot+1;
    end
    disp('PA is:');
    disp(P*A);
    disp('LU is:');
    disp(L*U);
end

澄清：由于我们正在切换行，我们正在寻找分解P（置换矩阵）时间，而不是我们作为输入A的原始数据。A

代码说明：

  1. 首先我检查矩阵是否可逆，如果不是，停止。如果是，则枢轴为 (1,1)
  2. 我在第 1 列中找到最大的数字，然后切换行
  3. 使用高斯消元对第 1 列进行评分，将除点 (1,1) 之外的所有部分都归零
  4. 枢轴现在是 (2,2)，在第 2 列中找到最大的数字...冲洗，重复
4

3 回答 3

3

据我所知，您的代码似乎运行良好，至少对于基本示例A=[1,2;3,4]或A=[3,4;1,2]. 将您的函数定义更改为：

function [L,U,P] = newgauss(A)

所以你可以输出你的计算值（比使用好得多disp，但这也显示了正确的结果）。然后你会看到P*A = L*U。也许您期望直接L*U等于A？lu您还可以通过 Matlab 的函数确认您是正确的：

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

置换矩阵是正交矩阵，所以 P -1 = P T。如果你想找回A你的代码，你可以这样做：

P'*L*U

同样，使用lu带有置换矩阵输出的 Matlab，您可以执行以下操作：

[L,U,P] = lu(A);
P'*L*U

（在检查行列式时，您还应该使用errororwarning而不是使用方式disp，但他们可能没有教过。）

于 2013-07-21T20:17:22.187 回答
1

请注意，该det函数是使用 LU 分解本身来实现的，以计算行列式...递归任何人:)

除此之外，页面末尾有一个提醒，建议使用cond而不是det测试矩阵奇点：

不推荐使用测试奇点abs(det(X)) <= tolerance，因为很难选择正确的容差。该函数cond(X)可以检查奇异矩阵和近奇异矩阵。

COND 使用奇异值分解（参见其实现edit cond.m：）

于 2013-07-22T00:54:15.143 回答
0

对于将来发现此问题并需要有效解决方案的任何人：

OP 的代码不包含在L创建 permutation matrix 时切换元素的逻辑P。给出与 Matlab 函数相同的输出的调整后的代码lu(A)是：

function [L,U,P] = newgauss(A)
    [rows,columns]=size(A);
    P=eye(rows,columns); %P is permutation matrix
    tol = 1E-16; % I believe this is what matlab uses as a warning level
    if( rcond(A) <= tol) %% bad condition number
        error('Matrix is nearly singular')
    end
    U=A;
    L=eye(rows,columns);
    pivot=1;
    while(pivot<rows)
        max=abs(U(pivot,pivot));
        maxi=0;%%find maximum abs value in column pivot
        for i=pivot+1:rows
            if(abs(U(i,pivot))>max)
                max=abs(U(i,pivot));
                maxi=i;
            end
        end %%if needed then switch
        if(maxi~=0)
            temp=U(pivot,:);
            U(pivot,:)=U(maxi,:);
            U(maxi,:)=temp;
            temp=P(pivot,:);
            P(pivot,:)=P(maxi,:);
            P(maxi,:)=temp;

            % change elements in L-----
            if pivot >= 2
                temp=L(pivot,1:pivot-1);
                L(pivot,1:pivot-1)=L(maxi,1:pivot-1);
                L(maxi,1:pivot-1)=temp;
            end
        end %%Grade the column pivot using gauss elimination
        for i=pivot+1:rows
            num=U(i,pivot)/U(pivot,pivot);
            U(i,:)=U(i,:)-num*U(pivot,:);
            L(i,pivot)=num;
        end
        pivot=pivot+1;
    end
end

希望这对将来偶然发现此问题的人有所帮助。

于 2016-02-11T14:12:51.010 回答