《数值分析希尔伯特病态线性方程组.docx》由会员分享，可在线阅读，更多相关《数值分析希尔伯特病态线性方程组.docx（16页珍藏版）》请在三一文库上搜索。
1、病态线性代数方程组的求解理论的分析表明，求解病态的线性代数方程组是困难的。考虑方程组Hx =1. . . _b 的求解，其中 H 为 Hilbert 矩阵，H=(hj)" , hj =, i, j =1,2,.,ni j -11. 估计Hilbert矩阵2-条件数与阶数的关系；2. 选择问题的不同维数，分别用 Gaus或肖去法，Jacob迭代，GS迭代和SOR迭代 求解，比较结果；3. 讨论病态I可题求解的算法。解：1、取Hilbert矩阵阶数最高分另U为n=20和n=100。采用Hilbert矩阵的2-条件数作为实验的比较对象，画出的曲线如下图所示：lg(cond (Hn)nlg(
2、cond(Hn)n 关系图2020从图中可以看出，在n< 13之刖，图像接近直线，在n> 13之后，图像趋于 平缓，在一定的范围内上下波动。为了比较图像的线性部分，作出一条直线与已 知曲线进行比较。比较直线的关系式为： lg(cond(Hn)=1.519n-1.833,结果下 图所示。lg(cond(Hn)n 关系图30 .,.IIIIIIIBI25 一_20 ._-g10 -5 -0 -5 1111111102468101214161820n从图2中可以看出，当n较小时，lg(cond(Hn) n之间近似满足线性关系。当n继续增大到100时，lg(cond(Hn) n关系图下图所
3、示：lg(cond(Hn)n 关系图25 ,!11(11!120 -. 寸一 .顶厂/ A d'rn 15.-Hog 10 -5-r0 111111110102030405060708090100n从图中可以看出，图像的走势符合在n=20时的猜想，在n大丁一定的值之后，图像趋丁平缓，且在一定范围内震荡，同时乂有一定上升趋势，但上升速度 很慢。2、选择不同的阶数n,设方程组的精确解为xz=(1,1,-；1)T进行计算，用四种方法解 x_Guass1 x_Jacobi1 x_GS1 x_SOR侗比表如下表所小。四种方法解 x_Guass> x_Jacobi、x_GS x_SOR对比表
4、nx_Guassx_Jacobix_GSx_SOR31-8.22E+3070.9999911-1.552e+30810.999991-Inf0.999971411.0172e+30810.999991Inf0.999651.00011Inf1.00080.999741Inf0.999481.000251-6.07E+3070.998791.00011-1.2685e+3081.02230.999041-1.6592e+3080.905111.0041-Inf1.14210.994111-Inf0.930951.002861-6.64E+3070.9994911-1.4322e+3081.005
5、90.999541-Inf0.995181.00111-Inf0.951491.00091-Inf1.10210.996211-Inf0.945361.002371-9.20E+3071.000711-Inf0.981921.00021-Inf1.10220.99791-Inf0.808241.00621-Inf1.0620.994641-Inf1.14480.998981-Inf0.899561.0021811.4378e+3081.00190.999971Inf0.958751.0011Inf1.19620.994081Inf0.730241.01141Inf0.954490.994641
6、Inf1.18240.998061Inf1.14660.996681Inf0.828041.0042915.09E+3071.00310.9999211.1719e+3080.939461.001911.6249e+3081.24890.990291Inf0.754461.01540.99999Inf0.838240.995491Inf1.09850.998760.99998Inf1.22470.994051Inf1.11480.999691Inf0.775611.0046101-7.82E+3071.00360.999941-Inf0.93761.00111-Inf1.21420.99555
7、1-Inf0.889821.00450.99992-Inf0.779991.00111.0002-Inf0.958251.00030.99964-Inf1.16090.997481.0003-Inf1.22070.997760.99982-Inf1.08170.999271-Inf0.750191.0031117.97E+3071.00290.999951Inf0.958571.00070.99998Inf1.09130.998091.0002Inf1.09090.99980.99894Inf0.802911.0031.0037Inf0.828491.0010.99187Inf1.02360.
8、999361.0111Inf1.18860.997760.99081Inf1.21280.997991.0042Inf1.06010.999610.99916Inf0.736721.00281218.11E+3071.00120.999971Inf0.995921.00030.99993Inf0.919041.00021.0009Inf1.29740.996210.99316Inf0.887531.0041.0298Inf0.743091.00140.91765Inf0.874511.0011.1481Inf1.08370.998340.82734Inf1.22150.997821.1258I
9、nf1.21480.998030.9479Inf1.04310.999931.0094Inf0.715151.00291314.90E+3070.998980.999990.999971.1863e+3081.03880.99991.00121.6981e+3080.74141.00220.98043Inf1.46940.993251.1771Inf1.0031.00450.033234Inf0.706161.00154.3909Inf0.746491.0023-6.8988Inf0.953940.9991313.35Inf1.16010.99821-11.81Inf1.26350.99739
10、9.4534Inf1.21940.99831-2.2127Inf1.01981.00021.5352Inf0.675951.00321415.94E+3070.9967710.999991.4525e+3081.07990.999511.0004Inf0.58371.00410.99479Inf1.59150.990681.0373Inf1.12841.00460.85719Inf0.711571.00151.2465Inf0.650911.00351.176Inf0.824241-0.85802Inf1.06130.99895.3287Inf1.24210.99733-4.4541Inf1.
11、3020.997485.0263Inf1.21660.99841-0.64129Inf0.985591.00061.2863Inf0.621961.0035211Inf0.9894910.99978Inf1.16240.999291.0084Inf0.559161.00170.86586Inf0.938321.00142.1174Inf1.53570.9961-4.2693Inf1.44540.9985715.159Inf1.05070.99966-18.28Inf0.698661.00176.2524Inf0.526541.002113.831Inf0.538331.002110.008In
12、f0.681851.0014-49.838Inf0.892471.000445.871Inf1.11260.99942-27.1Inf1.29840.9986154.555Inf1.41970.99808-34.976Inf1.45850.99793-31.788Inf1.40570.9981824.491Inf1.25930.9988541.386Inf1.02150.99993-43.821Inf0.697761.001413.526Inf0.295371.0032Gauss消去法求解：选择问题的阶数为38时，用Gauss消去法求得的解与精确解一致，当阶数 为914时，解开始出现偏差，而且
13、n越大，偏差越大。用迭代法求解：取初始向量为1.2（1,1,；1）T.无论n为多少阶，用Jacobi迭代方法迭代出现发散的不稳定现象，无法求解；用GS迭代方法迭代不发散，能求得解，但收敛非常缓慢，当迭代次数取得 相当大（ 20000次）时解仍在精确解附近波动；取w=1.5，用SO双代方法迭代不发散，能求得解，收敛速度较GS迭代快一些，但仍非常缓慢。选择问题的阶数为21时：用Gauss消去法求得的解与精确解相差很大，相差 10的量级。用迭代法求解时，取初始向量为1.2（1,1,；1）t用Jacobi迭代方法迭代发散，无法求解；用GS迭代方法迭代不发散，能求得解，但收敛非常缓慢，迭代20000次后
14、， 算得的值与精确值1相差很大。用SOR4代方法迭代不发散，能求得解，收敛速度还行。从上面的结果可以看出当病态问题的阶数升高时作为直接法的Gauss消去法乂能求解变成不能求解。用Jacobi迭代方法迭代发散，无法求解；而GS和SOR 迭代法在阶数升高时仍能求解，选择合适的w,可使SOR迭代法收敛速度加快。 但在阶数较低时直接法能求得精确解而迭代法却总存在一定的误差。可见直接法与迭代法各有各的优势与不足。3、关丁病态问题的求解，主要的方法是对原方程作某些预处理，以降低系数矩阵的条件数。例如选择非奇异矩阵 P,Qr,使方程组Ax = b化为等价方程组(PAQ)y=Pb ,原方程的解x=Qy.原则上
15、应使矩阵PAQ的条件数比A有所改善。 一般P和Q可选择为三角形矩阵或对角矩阵。 理论上最好选择对角阵P和Q，满 足：cond(PAQ) =min cond(PAQ)。以Hilbert矩阵为例。lg(cond(hn)/cond(Hn)n 关系图2XJnwnuanocrnb001I5012503003504005045011上图中的hn=Dn HnDn ( Dn为Hn的对角线兀素开万构成的对角矩阵)，n 为希尔伯特矩阵的阶数。我们观察到，随希尔伯特矩阵阶数的增大，函数 lgcond(hn)/cond(Hn)在-3,1.5区间波动，主要集中在-1.5,0.5区间。我们知道 当 cond(hn)壬co
16、nd(Hn)时,有 lgcond(hn)/cond(Hn)壬0 ,在上图中,我们可以 很容易的观察到，对丁大部分n,函数值都是小丁或者等丁零的，这说明 Hn经过 预处理后的hn的条件数较小，在一定程度上改善了原希尔伯特矩阵的性质。由丁希尔伯特矩阵病态的性质，对丁对原希尔伯特矩阵进行预处理后的新希 尔伯特矩阵的条件数在一定的区间呈波动的变化规律。从整体上观察，对丁大多数的n,进行预处理后的希尔伯特矩阵的条件数有明显的降低，这就说明预处理 改善了大多数阶数的希尔伯特矩阵的性质。第1小题Matlab编程如下:aln6li= (9r(DNsluoJi-roQ6)os (寸CXIO;Z0UO二®
17、;<u(uh)puoo)6_¥_m OCXISNSCO二(UH)puoo)6_-)oqB_AOCXIO;Z0UO二 u-)oqB_x (=二£)>1£MO_d coco8vF6Tgur:PUCD (Uwo80ll(u)>i (H)PU8"U二 >1(u)qznHwuu OVA-ss qgel/M 照 w zw(9r(DNsluoJi-roQ6)os (寸CXIO;Z0UO二®<u (UH)PU8)6_¥_M OCXISNSCO二(UH)puoo)6_-)oqB_AOCXI由Z0UO二 u-)oqB_x (
18、cxlmzmodPUCD(u 二 x)08oll(u)x(H)pu?(u)m(u)qznHCXIrlrCXI一 。&clearfor m=3:25n=m % Hilbert矩阵的阶数n1=n+1H=hilb(n)K=cond(H)u=1.2x0=u*ones(n,1)xz=ones(n,1)b=H*xza=HD=diag(diag(a)U=-triu(a,1)L=-tril(a,-1)w=1.5d=1.0e-6max=1000c=a,b% Gaus前肖去法for k=1:nbmax=0; % 选主元，并存放在变量bmax中for i=k:nif bmax-abs(c(i,k)<0
19、bmax=abs(c(i,k);l=i; %记下绝对值最大者所在位置endendif bmax<1.0e-20'stop' pause; %如果主元绝对值小于1.0e-20,则A奇异，计算终止end%进行行交换if l=k %如果l=k,则不需要交换for j=k:n1t=c(l,j); c(l,j)=c(k,j); c(k,j)=t;endend%进行消元t=1/c(k,k);c(k,j)=c(k,j)*t; c(k+1:n,j)=c(k+1:n,j)-c(k+1:n,k)*c(k,j);for j=k+1:n1end endC的第n+1列中%回代法求解上三角形方程组，
20、解存放在矩阵for k=2:ni=n1-k;for j=i+1:nc(i,n1)=c(i,n1)-c(i,j)*c(j,n1);endendxGuass=c(1:n,n1)% Jacobi迭代法BJ=D(L+U)fJ=Dbx2=x0xJacobi=BJ*x2+fJkJacobi=1while (norm(xJacobi-x2)>=d)&(kJacobi<max)x2=xJacobixJacobi=BJ*x2+fJkJacobi=kJacobi+1end% GS迭代法BG=(D-L)UfG=(D-L)bx3=x0xGS=BJ*x3+fGkGS=1while (norm(xGS
21、-x3)>=d)&(kGS<max)x3=xGSxGS=BG*x3+fGkGS=kGS+1end% SOR迭代法lw=(D-w*L)(1-w)*D+w*U)fS=(D-w*L)b*wx4=x0xSOR=lw*x4+fSkSOR=1while (norm(xSOR-x4)>=d)&(kSOR<max)x4=xSORxSOR=lw*x4+fSkSOR=kSOR+1end for p=1:mxGuass1(p,m)=xGuass(p,1)xJacobi1(p,m)=xJacobi(p,1)xGS1(p,m)=xGS(p,1)xSOR1(p,m)=xSOR(p,
22、1)x(p,4*m-3)=xGuass1(p,m)x(p,4*m-2)=xJacobi1(p,m)x(p,4*m-1)=xGS1(p,m)x(p,4*m)=xSOR1(p,m)endend第3小题Matlab编程如下：cleark=500i=1:kfor n=1:kH=hilb(n)K(n)=cond(H)D1=diag(sqrt(diag(H)D=inv(D1)Hy=D*H*DKy1(n)=cond(Hy)Ky(n)=log10(Ky1(n)/K(n)endplot(i,Ky)xlabel('n','fontsize',20)ylabel('lg(cond(hn)','fontsize',20)title('lg(cond(hn)/cond(Hn)n 关系图,'fontsize',24)set(gca,'fontsize',16)