Hilbert代数方程组的数值解法

1、Hilbert代数方程组的数值解法1 Hilbert矩阵的条件数和矩阵的阶数的关系编制Matlab程序：clear,clcformat long%format shortfor n=1:14;A=hilb(n);%A=rand(n);condA(n尸cond(A,inf);endplot(log10(condA)得出图形条件数的对数-阶次18164 2 0 8 6111)Atfnovo1ngr02468101214n图1: Hilbert矩阵和阶次的关系由图可见a.Hilbert矩阵的2-条件数的对数与阶次近似正比b.阶次超过12后，计算困难，舍入误差会导致计算不准结论：Hilbert矩阵的2

2、-条件数随阶次指数增长，关系大概是：Log10(Cond(A)尸 1.33791720780254(n-1)2 .各种求解方法的对比2.1 直接法：Gauss消去方法（程序见附录，下同）在双精度型变量下，即 eps=2.220446049250313e-016,对于直接法 Gauss消去方法，求解Hilbert矩阵方程组，在阶次 n=13时，解得：X=10.99997 1.00130.97881 1.1906-0.0354414.6171-7.3967 14.089-12.542 9.9162-2.3817 1.5624出现明显错误，可见 Gauss消去方法对病态问题比较敏感。求解阶次不高。2

3、.2 Jacobi迭代法很遗憾，jakobi迭代法的迭代矩阵的谱半径随阶次线性上升（程序见附录2）,普遍大于1,计算结果发送，无法迭代出满意的结果。需放弃Jacobi迭代收敛性分析2520151050径半谱0510152025阶次图2： Jacobi迭代矩阵的谱半径2.3 Gauss-Seidel迭代与SOR迭代求解先分析SOR迭代矩阵的谱半径和松弛因子w的关系径半谱SOR迭代5.5高斯-赛德尔迭代误差曲线4 -3.53 -2.521.50102030405060708090100迭代次数图4: G-S迭代下降曲线取w=1.5时，迭代次数需58798次前100步下降曲线SOR迭代误差曲线876

4、5差4误34打卅 阳计血】：3匕”1hr”“trt:r100102030405060708090100迭代次数w=0.5迭代次数16290图5: SOR (w=1.5)下降曲线图6: SOR (w=0.5)下降曲线2.4 SOR按照寻优算法的迭代求解寻优得出w=0.1时，谱半径相对较小，迭代次数只需3689次，下降曲线也非常快，如图 71.4SOR迭代误差曲线1.2 1 -0.8差误 0.6 0.4 0.2 -n _ I： ,.二”：二：；二II0102030405060708090100迭代次数图7: SOR (w=0.1 )下降曲线2.5 共轲梯度法迭代求解奇迹发生了，迭代只需要 9次就完

5、成了，下降也非常快，见图最速下降迭代误差曲线1234迭代次数2.6 各算法综合对比综上所述，在双精度的数值精度下和e=0.00001下，在笔记本可接受的计算时间里，做各种算法的对比如下：GaussJacobiG-SSORCG能计算到n 阶122120020005500迭代次数/最多中等较少计算时间/慢较快超快稳定性较差最差较好很好3.讨论病态方程组的求解3.1 对于Hilbert病态方程组的求解Gauss直接法和Jacobi迭代法都比较差，G-S法和SOR虽能求解，但由于其条件数都接近于1,收敛太慢，计算时间太长，因此实际可计算的阶次不高。3.2 采用高精度计算我们分别采用单精度和双精度型变量

6、，以Gauss消去方法为求解办法，做数值实验，x2=gauss(A, b)x3=gauss(single(A), single(b)当阶次n=6时，得x2 =0.999999999999081.000000000026300.999999999821721.000000000464350.999999999486911.00000000020235x3 =1.00031300.99069881.06451270.82937911.19058060.9242350可以发现，x3已经明显偏离了正解，得出结论，采用高精度计算所得解有更多有效位。3.3 采用预处理方法4000阶的问题，速度加快很多，迭

7、代预处理共轲梯度法见附录，采用后，可以发现，计算 次数也减少了。附录1 Gauss?肖去方法function x=Gauss(a,b)% Gauss消去法n,m=size(a);nb=length(b);det=1;%存储行列式值x=zeros(n,1);for k=1:n-1for i=k+1:nif a(k,k)=0returnendm=a(i,k)/a(k,k);for j=k+1:na(i,j)=a(i,j)-m*a(k,j);endb(i)=b(i)-m*b(k);enddet=det*a(k,k);enddet=det*a(n,n);for k=n:-1:1% 回代for j=k+

8、1:nb(k)=b(k)-a(k,j)*x(j);endx(k)=b(k)/a(k,k);end附录2 Jacobi迭代的收敛性分析clear,clcformat long%format shortfor n=1:25;A=hilb(n);D = diag(diag(A);B = inv(D) * (D - A);x=ones(n,1);e=0.000001;e=eps;b=A*x;f = inv(D) * b;% %谱半径大于等于1 就不收敛p(n) = max(abs(eig(B);endplot(p)xlabel( 阶次 )ylabel(谱半径)title(Jacobi 迭代收敛性分析)

9、附录3 G-S迭代function x, k, Ve = gauss_seidel(A, b, e,M)NAr, NA = size(A);if NAr = NA, error(A is not a square matrix); end%检查 A ， b 的大小是否匹配Nbr, Nb = size(b);if (Nbr = NA) & (1 = Nb), error(A and b have non-compatible dimensions); end%e 不能小于计算机的最小精度if e = 1sprintf(G-S iteration method not convergent) %r

10、eturnendwhile (r = e) & (k = M)x0 = x;x = B*x0 + f;k = k + 1;r = max(abs(x - x0);%if k = size(Ve, 1)%Ve(k,1) = r;%endend附录4 SORt代function x,n=SOR(A,b,w,e,M)x0 = zeros(size(b);%初始解设置为与b同型的零向量% if(w=2)%error;%return;% endD=diag(diag(A);%求 A 的对角矩阵L=-tril(A,-1);%求 A 的下三角阵U=-triu(A,1);%求 A 的上三角阵B=inv(D-L

11、*w)*(1-w)*D+w*U);f=w*inv(D-L*w)*b;n=0;%迭代次数r=1+e;while r=ex =B*x0+f;n=n+1;r = max(abs(x - x0);% Ve(n)=r;x0=x;if(n=M)disp(Warning: 迭代次数太多，可能不收敛！ );return;endend% figure% stem(Ve(1:100)% title(SOR 迭代误差曲线 );% xlabel(迭代次数)；% ylabel(误差)；.附录5 CG方法function x,n= conjgrad (A,b,x0)r1 = b-A*x0;p = r1;n = 0;for i=1:rank(A)%以下过程可参考算法流程if(dot(p,A*p) 1.0e-50)%循环结束条件break;endalpha = dot(r1,r1)/dot(p,A*p);x = x0+ alpha*p;r2 = r1- alpha*A*p;if(r2 =epsalpha = dot(r1,z1)/dot(p,A*p);x = x0 + alpha*p;r2 = r1 - alpha*A*p;z2 = iM*r2;belta = dot(r2,z2)/dot(r1,z1);p = z2+belta*p;n = n + 1;tol = norm(x-x0);x0 =