作者：凯鲁嘎吉 - 博客园
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
从结果可见希尔伯特矩阵是一个病态矩阵，用一般的直接法和迭代法会有较大的误差，甚至严重失真。
