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(順列行列) timesを分解しようとしています。AA

コードの説明:

  1. 最初に行列が可逆かどうかを確認し、そうでない場合は停止します。そうであれば、ピボットは (1,1) です
  2. 1列目の最大数を見つけて、行を入れ替えます
  3. ガウス消去法を使用して列 1 を等級付けし、スポット (1,1) 以外をすべてゼロにします。
  4. ピボットは (2,2) になり、列 2 で最大の数を見つけます... すすぎ、繰り返します
4

3 に答える 3

0

将来これを見つけて、実用的な解決策が必要な人のために：

LOP のコードには、順列行列を作成するときに要素を切り替えるためのロジックが含まれていませんP。lu(A)Matlab の関数と同じ出力を与える調整されたコードは次のとおりです。

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 に答える