Hilbert矩阵病态线性代数方程组的求解.doc

实验一 病态线性代数方程组的求解
Hilbert矩阵2-条件数与阶数的关系
运行 输入m=10 可以得到如下表的结果
阶数
1
2
3
4
5
条件数
1
+4
+5
阶数
6
和D2将方程组化为
D1AD2y=D1b，x=D2y
这称为矩阵的平衡问题，但是这样计算量比原问题本身要多。
或者通过变分原理将求解线性方程组的问题转化为等价的求解无约束函数最优化问题的极小值等等，可以参考
[1]郑洲顺，黄光辉，杨晓辉 求解病态线性方程组的混合算法
实验一所编程序如下：
m=input ('input m:=') ;
N=[1:m];
for i=1:m
n=N(i);
H=hilb(n);
k=cond(H);
disp('矩阵的阶数')
disp(n)
disp('矩阵')
disp(H)
disp('矩阵的条件数')
disp(k)
end
①
function guass(n)
n=input('请输入系数矩阵的维数n=');
H=hilb(n);
x_exact=ones(n,1);
b=H*x_exact;
x=Doolittle(H,b)
a=input('是否继续,继续请按1,结束请按0:');
if a==1
guass(n);
else end
②
function x=Doolittle(A,b)
% LU分解求Ax=b 的解
N=size(A);
n=N(1);
L=eye(n,n);%L的对角元素为1
U=zeros(n,n); %U为零矩阵
U(1,1:n)=A(1,1:n)%U第一行
L(1:n,1)=A(1:n,1)/U(1,1)%L第一列
for k=2:n
for i=k:n
U(k,i)=A(k,i)-L(k,1:(k-1))*U(1:(k-1),i)
end
for j=(k+1):n
L(j,k)=(A(j,k)-L(j,1:(k-1))*U(1:(k-1),k))/U(k,k)
end
end
y=solveDownTriangle(L,b);%调用下三角系数矩阵求解函数
x=solveUpTriangle(U,y);%调用上三角系数矩阵求解函数
③
function x=solveDownTriangle(A,b)
%求下三角系数矩阵的线性方程组Ax=b
N=size(A);
n=N(1);
for i=1:n
if (i>1)
s=A(i,1:(i-1))*x(1:(i-1),1);
else
s=0;
end
x(i,1)=(b(i)-s)/A(i,i);
end
④
function x=solveUpTriangle(A,b)
% 求上三角系数矩阵的线性方程组Ax=b
N=size(A);
n=N(1);
for i=n:-1:1
if (i<n)
s=A(i,(i+1):n)*x((i+1):n,1);
else
s=