学步园

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 matrix
X=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);
      end
end
for 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);
    end
end

如果运行程序可以看到：Ni(A)和inv(A)运算得到的逆矩阵是相同的

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

此方法中LU非直接三角分解只用了一次，通过增加一个j的循环，实现方程组的逐个求解，将得到的N个解向量C(:,j)合到X中得到最终结果。

个人感觉此方法的稳定性不错，暂时不需要改进了吧。

设计这个程序的时候并没遇到什么困难，思路已经想好了：

求解N个方程AXj=Ej;只要通过两个矩阵的对应列向量来存储Xj和Ej以及增加一个j循环，然后利用原来的Program就能达到预期的目的。