Matlab實現——求矩陣的逆(LU分解)

 

Program  ( ：Solve By Factorization with Pivoting )

思路及原理：

就得到：

程式：

function  X=Ni(A)%Input  - A is an N x N matrix%Output - I is an N x N inverse matrix of A %and I(j,:)containing the solution to AX(:,j) =E(:,j).%Initialize X, Y,the temporary storage matrix C, and the row % permutation information matrix R[N,N]=size(A);B=eye(N);   %B is an N x N identity matrixX=zeros(N,N);Y=zeros(N,N);C=zeros(1,N);R=1:N;%the next steps is to find the factorization(factorize for only once)for p=1:N-1%Find the pivot row for column p   [max1, j]=max(abs(A(p:N,p)));%Interchange row p and j      C=A(p,:);      A(p,:)=A(j+p-1,:);      A(j+p-1,:)=C;      d=R(p);      R(p)=R(j+p-1);      R(j+p-1)=d;      if A(p,p)==0      'A is singular.  No unique solution'      break      end   %Calculate multiplier and place in subdiagonal portion of A      for k=p+1:N         mult=A(k,p)/A(p,p);     A(k,p) = mult;         A(k,p+1:N)=A(k,p+1:N)-mult*A(p,p+1:N);      endendfor j=1:N        %when j is fixed then the method is similar to the Program 3.3    %Solve for Y(:,j)    Y(1,j) = B(R(1),j);    for k=2:N        Y(k,j)= B(R(k),j)-A(k,1:k-1)*Y(1:k-1,j);    end    %Solve for X(:,j)    X(N,j)=Y(N,j)/A(N,N);    for k=N-1:-1:1        X(k,j)=(Y(k,j)-A(k,k+1:N)*X(k+1:N,j))/A(k,k);    endend

 

如果運行程式可以看到：Ni(A)和inv(A)運算得到的逆矩陣是相同的

而且 A*Ni(A)=E 所以結果是令人滿意的

此方法中LU非直接三角分解只用了一次，通過增加一個j的迴圈，實現方程組的逐個求解，將得到的N個解向量C(:,j)合到X中得到最終結果。

個人感覺此方法的穩定性不錯，暫時不需要改進了吧。

設計這個程式的時候並沒遇到什麼困難，思路已經想好了：

求解N個方程AXj=Ej;只要通過兩個矩陣的對應列向量來儲存Xj和Ej以及增加一個j迴圈，然後利用原來的Program就能達到預期的目的。