Matlab实现 LU分解法解线性方程组（全选主元&&列选主元）

选主元LU分解

实验内容：列选主元LU分解和全选主元LU分解求解线性方程组
计算方法：

1. 全选主元消元法
   1.1 初始化
   根据参数A、b,记录下矩阵、右端项的尺寸n；
   以得到的尺寸n初始化解向量x；
   同时，以尺寸(n-1)初始化一个记录列交换顺序的向量tcolum;
   1.2 进入循环
   for k=1:n-1
   选主元:|A(p1,p2)|=max{|A(i,j)|:i=k:n,j=k:n}
   记录列交换顺序:tcolum(k)=p2
   换行:A(k,k:n)与A(p1,k:n)互换,同时b(k)与b(p1)互换
   换列:A(k:n,k)与A(k:n,p2)互换
   消元:
       ~~~    for i=k+1:n
          ~~~~~~       A(i,k)=A(i,k)/A(k,k);
          ~~~~~~       for j=k+1:n
             ~~~~~~~~~          A(i,j)=A(i,j)-A(i,k)A(k,j)
          ~~~~~~       end
       ~~~    end
   end
   注：记录列交换顺序，举例：tcolum(2)=5表示第二列和第五列互换。
   寻找第k个主元时，第k个主元原始位置的列数总要大于k。
   最多需要交换n-1次（矩阵A维数nn）,所以向量tcolum是n-1维的。
   1.3 回代求解
   x(n)=b(n)/A(n,n)
   for k=n-1: -1:1
        ~~~~     s=sum(A(k,k+1:n));
        ~~~~     x(k)=(b(k)-s)/A(k,k);
   end
   1.4 理顺解向量
   for i=n-1: -1:1
        ~~~~     交换x(i)与x(tcolum(i))
   end
2. 列选主元消元法
   2.1 初始化
   根据参数A、b,记录下矩阵、右端项的尺寸n；
   以得到的尺寸n初始化解向量x；
   2.2 进入循环
   for k=1:n-1
        ~~~~     选列主元:|A(p1,k)|=max{|A(i,j)|:i=k:n}
        ~~~~     换行:A(k,k:n)与A(p1,k:n)互换,同时b(k)与b(p1)互换
        ~~~~     消元:
        ~~~~     for i=k+1:n
            ~~~~~~~~         A(i,k)=A(i,k)/A(k,k);
            ~~~~~~~~         for j=k+1:n
                ~~~~~~~~~~~~             A(i,j)=A(i,j)-A(i,k)*A(k,j)
            ~~~~~~~~         end
        ~~~~     end
   end
   2.3 回代求解
   x(n)=b(n)/A(n,n)
   for k=n-1: -1:1
        ~~~~     s=sum(A(k,k+1:n));
        ~~~~     x(k)=(b(k)-s)/A(k,k);
   end
   计算结果：
   全选主元计算结果的残量范数为8.4428e-011
   列选主元计算结果的残量范数为1.8175e-010
   源程序：
   1.定义全选主元函数

function x=Lu_Decomplete(A,b)
%定义全选主元--解Ax=b
[n,n2]=size(A);
n3=length(b);
if n~=n2|n~=n3 %判断输入的合法性
    error('wrong input!');
end
x=ones(n,1);%记录解
tcolum=zeros(n-1,1);%记录列变换
for k=1:n-1 %从（1，1）位置到（n-1,n-1）位置
    %进行全选主元
    p1=k;p2=k;p_max=abs(A(k,k));%p1、p2分别记录/将选主元/的行、列
    for i=k:n
        for j=k:n
            if p_max<abs(A(i,j))
                p_max=abs(A(i,j));
                p1=i;p2=j;
            end
        end
    end
    %对A、b换行
    if p1>k 
        for i=k:n
            t=A(k,i);
            A(k,i)=A(p1,i);
            A(p1,i)=t;
        end
        t=b(k);
        b(k)=b(p1);
        b(p1)=t;
    end
    tcolum(k)=p2;
%对A换列
    if p2>k
        for i=1:n
            t=A(i,k);
            A(i,k)=A(i,p2);
            A(i,p2)=t;
        end
    end 
 %消元
   if A(k,k)<eps
     error('Fail solving');
 end
    for i=k+1:n
        A(i,k)=A(i,k)/A(k,k);
        for j=k+1:n
            A(i,j)=A(i,j)-A(i,k)*A(k,j);
        end
        b(i)=b(i)-A(i,k)*b(k);
    end
end
%回代求解
x(n)=b(n)/A(n,n);
for i=n-1:-1:1 
    z=0;
    for j=i+1:n
        z=z+x(j)*A(i,j);
    end 
    x(i)=(b(i)-z)/A(i,i);
end
%理顺解向量!!!重点在此！！
for i=n-1:-1:1
    t=x(i);
    x(i)=x(tcolum(i));
    x(tcolum(i))=t;
end

2.定义列选主元函数

function x=Lu_Decolum(A,b)
%定义列选主元LU分解
[n,n2]=size(A);
n3=length(b);
if n~=n2|n~=n3 %判断输入的合法性
    error('wrong input!');
end
x=ones(n,1);%初始化解向量
for k=1:n-1
    %选出列主元
    p=k;p_max=abs(A(k,k));
    for i=k+1:n
        if p_max<abs(A(i,k))
            p_max=abs(A(i,k));
            p=i;
        end
    end
    %据选出的列主元进行换行
    if p>k
        for i=k:n
           t=A(k,i);
            A(k,i)=A(p,i);
            A(p,i)=t;
        end
        t1=b(k);
        b(k)=b(p);
        b(p)=t1;
    end
    %对矩阵进行LU分解
    if A(k,k)==0
        break;
    end
    for j=k+1:n
        m=A(j,k)/A(k,k);
        for i=k+1:n
            A(j,i)=A(j,i)-m*A(k,i);
        end
        b(j)=b(j)-m*b(k);
    end
end
%回代法解方程组
x(n)=b(n)/A(n,n);
for k=n-1:-1:1
    z=0;
    for i=k+1:n
        z=z+x(i)*A(k,i);
    end
    x(k)=(b(k)-z)/A(k,k);
end

3.调用函数计算结果
在命令行中，
输入A=randn(1000,1000);b=ones(1000,1);
调用全选主元函数：
x1=Lu_Decomplete(A,b);
计算此时结果的残量范数：
norm(Ax1-b,2).
同理，调用列选主元函数：
x2=Lu_Decolum(A,b);
计算此时结果的残量范数：
norm(Ax2-b,2).