数值分析(Hilbert矩阵)病态线性方程组的求解Matlab程序

（Hilbert 矩阵）病态线性方程组的求解
理论分析表明，数值求解病态线性方程组很困难。

考虑求解如下的线性方程组的求解Hx = b ，期中H 是Hilbert 矩阵，()ij n n H h ⨯=，11
ij h i j =
+-，i ，j = 1,2，…，n 1. 估计矩阵的2条件数和阶数的关系 2. 对不同的n ，取(1,1,,1)n
x =∈K ?，分别用Gauss 消去，Jacobi 迭代，Gauss-seidel 迭
代，SOR 迭代和共轭梯度法求解，比较结果。

3. 结合计算结果，试讨论病态线性方程组的求解。

第1小题：
condition.m %第1小题程序
t1=20;%阶数n=20
x1=1:t1;
y1=1:t1;
for i=1:t1
H=hilb(i);
y1(i)=log(cond(H));
end
plot(x1,y1);
xlabel('阶数n');
ylabel('2-条件数的对数（log(cond(H)）');
title('2-条件数的对数（log(cond(H)）与阶数n 的关系图');
t2=200;%阶数n=200
x2=1:t2;
y2=1:t2;
for i=1:t2
H=hilb(i);
y2(i)=log(cond(H));
end
plot(x2,y2);
xlabel('阶数n');
ylabel('2-条件数的对数（log(cond(H)）');
title('2-条件数的对数（log(cond(H)）与阶数n 的关系图');
画出Hilbert 矩阵2-条件数的对数和阶数的关系
n=200时
n=20时
从图中可以看出，
1）在n小于等于13之前，图像近似直线
log（cond（H））~1.519n-1.833
2）在n大于13之后，图像趋于平缓，并在一定范围内上下波动，同时随着n的增加稍有上升的趋势
第2小题：
solve.m%m第2小题主程序
N=4000;
xGauss=zeros(N,1);
xJacobi=zeros(N,1);
xnJ=zeros(N,1);
xGS=zeros(N,1);
xnGS=zeros(N,1);
xSOR=zeros(N,1);
xnSOR=zeros(N,1);
xCG=zeros(N,1);
xnCG=zeros(N,1);
for n=1:N;
x=ones(n,1);
t=1.1;%初始值偏差
x0=t*x;%迭代初始值
e=1.0e-8;%给定的误差
A=hilb(n);
b=A*x;
max=100000000000;%可能最大的迭代次数w=0.5;%SOR迭代的松弛因子
G=Gauss(A,b);
[J,nJ]=Jacobi(A,b,x0,e,max); [GS,nGS]=G_S(A,b,x0,e,max);
[S_R,nS_R]=SOR(A,b,x0,e,max,w); [C_G,nC_G]=CG(A,b,x0,e,max);
normG=norm(G'-x);
xGauss(n)=normG;
normJ=norm(J-x);
nJ;
xJacobi(n)=normJ;
xnJ(n)=nJ;
normGS=norm(GS-x);
nGS;
xGS(n)=normGS;
xnGS(n)=nGS;
normS_R=norm(S_R-x);
nS_R;
xSOR(n)=normS_R;
xnSOR(n)=nS_R;
normC_G=norm(C_G-x);
nC_G;
xCG(n)=normC_G;
xnCG(n)=nC_G;
end
Gauss.m%Gauss消去法
function x=Gauss(A,b)
n=length(b);
l=zeros(n,n);x=zeros(1,n);
%消去过程
for i=1:n-1
for j=i+1:n
l(j,i)=A(j,i)/A(i,i);
for k=i:n
A(j,k)=A(j,k)-l(j,i)*A(i,k);
end
b(j)=b(j)-l(j,i)*b(i);
end
end
%回代过程
x(n)=b(n)/A(n,n);
for i=n-1:-1:1
c=A(i,:).*x;
x(i)=(b(i)-sum(c(i+1:n)))/A(i,i);
end
Jacobi.m%Jacobi迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m 可能最大的迭代次数
function [x,n]=Jacobi(A,b,x0,e,m)
n=length(A);
D=diag(diag(A));
U=-triu(A,1);
L=-tril(A,-1);
B=D\(L+U);
f=D\b;
x=B*x0+f;
n=1;
while norm(x-x0)>e
x0=x;
x=B*x0+f;
n=n+1;
if n>m
disp('Jacobi迭代次数过多，迭代可能不收敛');
break;
end
end
G_S.m%Gauss-Seidel迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数
function [x,n]=G_S(A,b,x0,e,m)
n=length(A);
D=diag(diag(A));
U=-triu(A,1);
L=-tril(A,-1);
B=(D-L)\U;
f=(D-L)\b;
x=B*x0+f;
n=1;
while norm(x-x0)>e
x0=x;
x=B*x0+f;
n=n+1;
if n>m
disp('Gauss-Seidel迭代次数过多，迭代可能不收敛');
break;
end
end
SOR.m%SOR超松弛迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数,w松弛因子
function [x,n]=SOR(A,b,x0,e,m,w)
n=length(A);
D=diag(diag(A));
U=-triu(A,1);
L=-tril(A,-1);
B=(D-w*L)\((1-w)*D+w*U);
f=(D-w*L)\b*w;
x=B*x0+f;
n=1;
while norm(x-x0)>e
x0=x;
x=B*x0+f;
n=n+1;
if n>m
disp('SOR超松弛迭代次数过多，迭代可能不收敛');
break;
end
end
CG.m%CG共轭梯度法,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数
function [x,n]=CG(A,b,x0,e,m)
r=b-A*x0;
p=r;
alpha=(r'*r)/(p'*(A*p));
x=x0+alpha*p;
r1=b-A*x;
n=1;
while norm(r1)>e
belta=(r1'*r1)/(r'*r);
p=r1+belta*p;
r=r1;
x0=x;
alpha=(r'*r)/(p'*(A*p));
x=x0+alpha*p;
r1=b-A*x;
n=n+1;
if n>m
disp('CG共轭梯度法迭代次数过多，迭代可能不收敛');
break;
end
end。