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

实验一病态线性代数方程组的求解
1.估计Hilbert矩阵2-条件数与阶数的关系
运行tiaojianshu.m 输入m=10 可以得到如下表的结果
2.选择不同维数，分别用Guass消去（LU分解），Jacobi迭代，GS 迭代，SOR迭代求解，比较结果。

说明：Hx=b，H矩阵可以由matlab直接给出，为了设定参考解，我们先设x为分量全1的向量，求出b,然后将H和b作为已知量，求x，与设定的参考解对比。

对于Jacobi迭代，GS迭代，SOR迭代，取迭代初值x0为0向量，迭代精度eps=1.0e-6，迭代次数<100000, SOR迭代中w=1.2和0.8分别计算。

a. n=5
b. n=8
c. n=10
d. n=15
取不同的n值，得到如下结果:
对于Guass法，可以看出来，随着n的增大，求解结果误差变大，这是因为随着n增大，系数矩阵的条件数变大，微小的扰动就容易造成很大的误差。

最后得不到精确解。

对于Jacobi迭代，计算结果为Inf，说明是发散的。

对于GS迭代和SOR迭代，结果是收敛的，但是可以看出迭代次数比较多，并且对于不同维数GS和SOR收敛速度不一样，有时候GS快，有时SOR快。

对SOR取不同的w迭代速度也不一样，存在一个最优的松弛因子w。

并且可以知道，迭代次数多少跟初值x0也有关系。

3.讨论病态问题求解的算法。

通过上面的实验分析，可以看出，求解病态矩阵的时候要小心，否则可能得不到所要求的精确度。

可以采用高精度运算，用双倍多倍字长，使得由于误差放大而损失若干有效数字位之后，还能保留一些有效位。

另外可以通过对原方程作某些预处理，降低系数矩阵的条件数，因为cond(aA)=cond(A),所以不能通过将每一个方程乘上相同的常数来达到这个目标，可考虑将矩阵的每一行和每一列分别乘上不同的常数，亦即找到可逆的对角阵D1和D2将方程组化为
D1AD2y=D1b，x=D2y
这称为矩阵的平衡问题，但是这样计算量比原问题本身要多。

或者通过变分原理将求解线性方程组的问题转化为等价的求解无约束函数最优化问题的极小值等等，可以参考
[1]郑洲顺，黄光辉，杨晓辉求解病态线性方程组的混合算法
实验一所编程序如下：
1.求条件数
tiaojianshu.m
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
2.Guass法
①Guass.m
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
②Doolittle.m
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);%调用上三角系数矩阵求解函数
③solveDownTriangle.m
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
④solveUpTriangle.m
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=0;
end
x(i,1)=(b(i)-s)/A(i,i);
end
3.Jacobi法
function [x,n]=jacobi(a,x0)
a=input('请输入系数矩阵的维数n=');
x0=input('请输入迭代初始向量x0');
eps=1.0e-6;%解的精度控制
m=10000;%迭代步数控制
H=hilb(a);%生成h矩阵
x_exact=ones(a,1);%求出x精确值
b=H*x_exact;
D=diag(diag(H));
L=-tril(H,-1);
U=-triu(H,1);
B=D\(L+U);
f=D\b;
x=x0;
n=0;
tol=1;
while tol>=eps
x=B*x0+f;
n=n+1;
tol=norm(x-x0);
x0=x;
if (n>=m)
disp('迭代次数太多，可能不收敛');
return;
end
end
disp(x)
disp(n)
4.GS法
function [x,n]=gauseidel(a,x0)
a=input('请输入系数矩阵的维数n=');
x0=input('请输入迭代初始向量x0');
eps=1.0e-6;%解的精度控制
m=100000;%迭代步数控制
H=hilb(a);%生成h矩阵
x_exact=ones(a,1);%求出x精确值
b=H*x_exact;
D=diag(diag(H));
L=-tril(H,-1);
U=-triu(H,1);
G=(D-L)\U;
f=(D-L)\b;
x=x0;
n=0;
tol=1;
while tol>=eps
x=G*x0+f;
n=n+1;
tol=norm(x-x0);
x0=x;
if (n>=m)
disp('迭代次数太多，可能不收敛');
return;
end
end
disp(x)
disp(n)
5.SOR法
function [x,n]=SOR(a,x0,w)
a=input('请输入系数矩阵的维数n=');
x0=input('请输入迭代初始向量x0=');
w=input('请输入w=');
eps=1.0e-6;%解的精度控制
m=100000;%迭代步数控制
H=hilb(a);%生成h矩阵
x_exact=ones(a,1);%求出x精确值
b=H*x_exact;
if (w<=0||w>=2)
error;
return;
end
D=diag(diag(H));
L=-tril(H,-1);
U=-triu(H,1);
B=inv(D-L*w)*((1-w)*D+w*U);
f=w*inv((D-L*w))*b;
x=x0;
n=0;
tol=1;
while tol>=eps
x=B*x0+f;
n=n+1;
tol=norm(x-x0);
x0=x;
if (n>=m)
disp('迭代次数太多，可能不收敛');
return;
end
end
disp(x)
disp(n)
欢迎您的下载，资料仅供参考！。