Hilbert代数方程组的数值解法.docx

1、Hilbert代数方程组的数值解法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)得出图形图1：Hilbert矩阵和阶次的关系由图可见a.Hilbert矩阵的2-条件数的对数与阶次近似正比; b.阶次超过12后，计算困难，舍入误差会导致计算不准结论：Hilbert矩阵的2-条件数随阶次指数增长，关系大概是：Log10(Cond(A

2、)= 1.33791720780254(n-1)2.各种求解方法的对比2.1 直接法：Gauss消去方法（程序见附录，下同）在双精度型变量下，即eps=2.220446049250313e-016，对于直接法Gauss消去方法，求解Hilbert矩阵方程组，在阶次n=13时，解得:X=1 0.99997 1.0013 0.97881 1.1906 -0.035441 4.6171 -7.3967 14.089 -12.542 9.9162 -2.3817 1.5624出现明显错误，可见Gauss消去方法对病态问题比较敏感。求解阶次不高。2.2 Jacobi迭代法很遗憾，jakobi迭代法的迭代

3、矩阵的谱半径随阶次线性上升（程序见附录2），普遍大于1，计算结果发送，无法迭代出满意的结果。需放弃图2：Jacobi迭代矩阵的谱半径2.3 Gauss-Seidel迭代与SOR迭代求解先分析SOR迭代矩阵的谱半径和松弛因子w的关系图3：SOR迭代矩阵的谱半径和松弛因子w的关系可以发现SOR迭代和G-S迭代的谱半径都普遍接近于1，因此松弛因子的选择对计算结果的影响不大。我们选择w=1，即G-S迭代来做数值实验：选定e=0.000001；迭代次数：15610次下降曲线：图4：G-S迭代下降曲线取w=1.5时，迭代次数需58798次前100步下降曲线图5：SOR（w=1.5）下降曲线w=0.5 迭代

4、次数16290图6：SOR（w=0.5）下降曲线2.4 SOR按照寻优算法的迭代求解寻优得出w=0.1时，谱半径相对较小，迭代次数只需3689次，下降曲线也非常快，如图7图7：SOR（w=0.1）下降曲线2.5 共轭梯度法迭代求解奇迹发生了，迭代只需要9次就完成了，下降也非常快，见图图8：共轭梯度法下降曲线2.6 各算法综合对比综上所述，在双精度的数值精度下和e=0.00001下，在笔记本可接受的计算时间里，做各种算法的对比如下：GaussJacobiG-SSORCG能计算到n阶122120020005500迭代次数/最多中等较少计算时间/慢较快超快稳定性较差最差一般较好很好3. 讨论病态方程

5、组的求解3.1 对于Hilbert病态方程组的求解Gauss直接法和Jacobi迭代法都比较差，G-S法和SOR虽能求解，但由于其条件数都接近于1，收敛太慢，计算时间太长，因此实际可计算的阶次不高。3.2 采用高精度计算我们分别采用单精度和双精度型变量，以Gauss消去方法为求解办法，做数值实验，x2=gauss(A, b)x3=gauss(single(A), single(b)当阶次n=6时，得x2 = 0.99999999999908 1.00000000002630 0.99999999982172 1.00000000046435 0.99999999948691 1.0000000

6、0020235x3 = 1.0003130 0.9906988 1.0645127 0.8293791 1.1905806 0.9242350可以发现，x3已经明显偏离了正解，得出结论，采用高精度计算所得解有更多有效位。3.3 采用预处理方法预处理共轭梯度法见附录，采用后，可以发现，计算4000阶的问题，速度加快很多，迭代次数也减少了。 附录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-1 for i=k+1:n if a(k,k)=0

7、 return end m=a(i,k)/a(k,k); for j=k+1:n a(i,j)=a(i,j)-m*a(k,j); end b(i)=b(i)-m*b(k); end det=det*a(k,k);enddet=det*a(n,n); for k=n:-1:1 %回代 for j=k+1:n b(k)=b(k)-a(k,j)*x(j); end x(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)

8、 * (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迭代收敛性分析)附录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);

9、if (Nbr = NA) & (1 = Nb), error(A and b have non-compatible dimensions); end%e不能小于计算机的最小精度if e = 1 sprintf(G-S iteration method not convergent) %returnendwhile (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 SOR迭代function x,n=SOR

10、(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*w)*(1-w)*D+w*U);f=w*inv(D-L*w)*b;n=0; %迭代次数r=1+e;while r=e x =B*x0+f; n=n+1; r = max(abs(x - x0);% Ve(n)=r; x0=x; if(n=M) disp(Warning: 迭代次数太多，可能

11、不收敛！); 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; end alpha = dot(r1,r1)/dot(p,A*p); x = x0+ alpha*p; r2 = r1- alpha*A*p; if(r2 =eps alpha = 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 = x; %更新迭代值 r1 = r2; z1 = z2;end