MATLAB矩陣的LU分解及在解線性方程組中的應用

作者：凱魯嘎吉 - 博客園
http://www.cnblogs.com/kailugaji/

三、實驗程序

五、解答（按如下順序提交電子版）

1.(程序)

(1)LU分解源程序：

function [l,u]=lu12(a,n)
for k=1:n-1
    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
l=eye(n);
u=zeros(n,n);
for k=1:n
    for i=k:n
        u(k,i)=a(k,i);
    end
end
for k=1:n
    for j=1:k-1
        l(k,j)=a(k,j);
    end
end      

(2)直接三角分解法源程序：

function [a,l,u,y,x]=direct_triangle(a,b,n)
%a為N*N矩陣，b為n*1列向量
for k=1:n-1
    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
l=eye(n);
u=zeros(n,n);
for k=1:n
    for i=k:n
        u(k,i)=a(k,i);
    end
end
for k=1:n
    for j=1:k-1
        l(k,j)=a(k,j);
    end
end
y=ones(n,1);
x=ones(n,1);
y(1,1)=b(1,1);
for i=2:n
    s=0;
    for k=1:i-1
        s=s+l(i,k)*y(k,1);
     end
     y(i,1)=b(i,1)-s;
 end
 
 x(n,1)=y(n,1)/u(n,n);
 for j=n-1:-1:1
     s1=0;
     for k1=j+1:n
         s1=s1+u(j,k1)*x(k1,1);    
     end   
     x(j,1)=(y(j,1)-s1)/u(j,j);
 end

2.(運算結果)

（1）求一個4階矩陣的LU分解。

>> a=[10,7,8,7;7,5,6,5;8,6,10,9;7,5,9,10];
>> [l,u]=lu12(a,4)

l =

    1.0000         0         0         0
    0.7000    1.0000         0         0
    0.8000    4.0000    1.0000         0
    0.7000    1.0000    1.5000    1.0000


u =

   10.0000    7.0000    8.0000    7.0000
         0    0.1000    0.4000    0.1000
         0         0    2.0000    3.0000
         0         0         0    0.5000

>>  a=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10];b=[32 23 33 31]';
>> [a,l,u,y,x]=direct_triangle(a,b,4)

a =

   10.0000    7.0000    8.0000    7.0000
    0.7000    0.1000    0.4000    0.1000
    0.8000    4.0000    2.0000    3.0000
    0.7000    1.0000    1.5000    0.5000


l =

    1.0000         0         0         0
    0.7000    1.0000         0         0
    0.8000    4.0000    1.0000         0
    0.7000    1.0000    1.5000    1.0000


u =

   10.0000    7.0000    8.0000    7.0000
         0    0.1000    0.4000    0.1000
         0         0    2.0000    3.0000
         0         0         0    0.5000


y =

   32.0000
    0.6000
    5.0000
    0.5000


x =

    1.0000
    1.0000
    1.0000
    1.0000

比如，希爾伯特矩陣就是一個病態矩陣，在方程組問題求解之前，可以先判斷其條件數是否較大。

源程序：hilbert.m:

function [A,cond1]=hilbert(k)
format rat
A=zeros(k,k);
for m=1:k
    for n=1:k
        A(m,n)=1/(m+n-1);
    end
end
cond1=cond(A,inf);

運行結果：

>> [A,cond1]=hilbert(3)

A =

       1              1/2            1/3     
       1/2            1/3            1/4     
       1/3            1/4            1/5     


cond1 =

     748      
>> [A,cond1]=hilbert(4)

A =

       1              1/2            1/3            1/4     
       1/2            1/3            1/4            1/5     
       1/3            1/4            1/5            1/6     
       1/4            1/5            1/6            1/7     


cond1 =

   28375       
>> [A,cond1]=hilbert(5)

A =

       1              1/2            1/3            1/4            1/5     
       1/2            1/3            1/4            1/5            1/6     
       1/3            1/4            1/5            1/6            1/7     
       1/4            1/5            1/6            1/7            1/8     
       1/5            1/6            1/7            1/8            1/9     


cond1 =

  943656

從結果可見希爾伯特矩陣是一個病態矩陣，用一般的直接法和迭代法會有較大的誤差，甚至嚴重失真。