c3_3.m

这是我从一本基于MATLAB实现的数值分析算法的书中自带的光盘程序

%*******************************************************
%程序C3_3.M
%选主元后的LU分解

% A: 待分解的系数矩阵；
% P: 选主元时的行变换矩阵；
% PA: A选主元后的矩阵；
%********************************************************

function  C3_3
clear
%A=input('A = ');    
A=[2 1 -3;-1 3 2;3 1 -3]
%首先选主元
m=size(A)
n=m(1)
P=eye(n);
for k=1:n-1 
    for p=k:n
        ekp=zeros(n);
        I=eye(n);
        ekp=I-(I(:,k)-I(:,p))*((I(:,k)-I(:,p))');
       if abs(A(p,k))==max(abs(A(k:n,k)))
          %A=ekp*A;
          P=ekp*P;
       end
    end
end
P
PA=P*A

%下面进行LU分解，PA是待进行LU分解的矩阵：
%计算L第1列和和U的第2行，并且，存放于原系数矩阵PA的相应位置：
for i=2:n
    PA(i,1)=PA(i,1)/PA(1,1);
    PA(2,i)=PA(2,i)-PA(2,1)*PA(1,i);
end
%计算L的第2列，并且，存放于原系数矩阵PA的相应位置：
for i=3:n
    PA(i,2)=(PA(i,2)-PA(i,1)*PA(1,2))/PA(2,2);
end
%计算L和U的剩余行和列，并且，存放于原系数矩阵A的相应位置：
for i=3:n
    for g=i:n
        sum1=0;
        for k=1:i-1 
            sum1=sum1+PA(i,k)*PA(k,g);
        end
        PA(i,g)=PA(i,g)-sum1;
    end
    for p=i+1:n
        sum2=0;
        for g=1:i-1
            sum2=sum2+PA(p,g)*PA(g,i);
        end
        PA(p,i)=(PA(p,i)-sum2)/PA(i,i);
    end
end
%输出矩阵L和U：
M=zeros(n,n);
for i=2:n
    for k=1:i-1
        M(i,k)=PA(i,k);
    end
end   
L=I+M
U=PA-M