Hilbert代数方程组数值解法.docxVIP

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))= 1.33791720780254(n-1)2.各种求解方法的对比2.1 直接法：Gauss消去方法（程序见附录，下同）在双精度型变量下，即eps=2.220446049250313e-016，对于直接法Gauss消去方法，求解Hilbert矩阵方程组，在阶次n=13时，解得:X=[10.999971.00130.978811.1906-0.0354414.6171-7.396714.089-12.5429.9162-2.38171.5624]出现明显错误，可见Gauss消去方法对病态问题比较敏感。求解阶次不高。2.2 Jacobi迭代法很遗憾，jakobi迭代法的迭代矩阵的谱半径随阶次线性上升（程序见附录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 迭代次数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. 讨论病态方程组的求解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.00000000020235x3 = 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 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