数值分析(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);