特征值问题数值计算上机实验.docx

特征值问题数值计算上机实验练习6.16：取初始向量为v0=1,1,1,1T,用乘幂法计算主特征值及其相应的特征向量。请绘制主特征值误差下降曲线，以及特征子空间距离的下降曲线。然后采用Atiken加速技巧和Rayleigh 商技术分别对算法进行加速，并完成类似的工作。解：采用教科书上8.1算法编写程序。在n=100时，由于1,1,1,1T在主特征值对应的特征向量上投影为0，如果取v0=1,1,1,1T，若无舍入误差的影响，就会导致乘幂法求解出的特征值为第二主特征值。故这里取初始向量v0为随机向量。为了便于比较，在每一步迭代中，将Atiken加速值和Rayleigh商加速值分别计算出来，并作出如下图像：（下图显示的为从10001步后开始）从上面两图可以看出，初始阶段乘幂法会剧烈震荡，在迭代一段时间后会突然特征值误差变大。在迭代6000次后，震荡现象不再出现。从第二张图中可以看出，采用Atiken和Rayleigh商加速得到的误差下降较快。由于Atiken加速和Rayleigh商加速均是对特征值进行的，对特征向量无影响，故特征子空间的距离的下降情形一致。下图为特征子空间距离的下降曲线。n=101时情况类似。练习6.17用反幂法求解最靠近2的特征值，及其对应的特征向量。观察是否有所谓的“一次收敛”的特性。解：当n=100时，T100的靠近2的特征值有1.9688963761593和2.0311036238407，可以看出特征值并不非常靠近2，预期反幂法的“一次收敛”特性并不明显。进行20次反幂法得到的特征值序列如下： 9.048615164438272.031103055857341.326560550532182.031103815513901.955013485188222.031103560391512.030789846401602.031103645088842.031256812288802.031103616774162.031060088198912.031103626200652.031119631851712.031103623054542.031098582417292.031103624102992.031105362183862.031103623753272.031105362183862.03110362386986从上面序列可以看出，迭代一次后特征值并不靠近真实特征值，但利用反幂法，迭代20步后，精度已达到10-10级别。从特征向量的距离的角度来看：（其中间断处是由于值为0）从上图可以看出特征向量无“一次收敛”特性，但经过6次迭代后，误差已达到10-8数量级。当n=101时，T101的靠近2的特征值为1.9384098828和2以及 2.06159011711234，即恰好2为T101的特征值，取q=2，由于此时矩阵T101-qI奇异，故可对平移量q进行一个非常小的扰动，取=10-15.对T101-(q+)I进行反幂法，预计此时反幂法的“一次收敛”特性非常明显。进行20次反幂法得到计算序列结果为：1.99999999999998222222222222222222从特征向量的角度来看，可得迭代过程中特征向量与用eig函数求出的特征向量的子空间之间的距离如下表：2.1073424255447e-080000000000000000000从上表可以看出，进行一次迭代后计算得到的特征向量与eig函数求出的特征向量子空间的距离已为10-8数量级，故可以认为两者相同，即反幂法确有“一次收敛”到特征向量的特性。不过并未观察到后续迭代会起到反作用的现象（即做多次反幂法后，特征子空间的距离反而可能会变大），这可能是由于MATLAB计算的精度太高导致。练习6.18 用幂法求解第二主特征值及其特征向量；用同时迭代法求解Tn的前两个主特征值。比较两者的计算效果。解：对于用乘幂法求解第二主特征值，先采用降阶法，这里取S1为Householder变换矩阵，先利用乘幂法求出Tn的主特征值1及其对应的特征向量x1， 然后求出降阶阵A2，此处利用MATLAB的eig函数，发现A2与Tn除1外的其余特征值的差距不大。再对A2用乘幂法，可以求出A2的主特征值，即Tn的第二特征值。见下表：而若用同时迭代法则可直接求出Tn的前两个主特征值。采用教科书上算法，可得计算结果如下：当n=100真实值降阶法+乘幂法误差同时迭代法误差3.999032564583983.999032564462331.21668009001041e-103.999032564583983.10862446895044e-153.996131194267193.996131196032091.75433800819746e-093.99613119426726.21724893790088e-15当n=101真实值降阶法+乘幂法误差同时迭代法误差3.999051439426733.999051439426734.44089209850063e-163.999051439426731.77635683940025e-15 3.99620665747409 3.99620665578206 1.69099045876919e-09 3.99620665747409 7.54951656745106e-15从上面两表可以看出，n=100时，利用乘幂法和降阶法求出的特征值误差为10-10和10-9量级，而利用同时迭代法求出的特征值误差为10-15和10-15量级。n=101时，利用乘幂法和降阶法求出的特征值误差为10-16和10-9量级，而利用同时迭代法求出的特征值误差为10-15和10-15量级。由此可以看出降阶法求出的第二个主特征值的精度远低于同时迭代法求出的第二个主特征值精度。理论分析如下，因为使用降阶法，第二个主特征信息的计算精度会受到第一个主特征信息计算精度的影响，由于舍入误差的积累，特征信息的逐次求解过程，势必造成计算精度越来越差。而子空间同时迭代法可以同时求出矩阵的前若干个主特征值信息，避免了舍入误差的积累，故对较靠后的主特征值也有很好的精度。练习6.19 分别用古典Jacobi方法、循环Jacobi方法和阈值Jacobi方法求解Tn的全部特征值；绘制相应的收敛过程。解： （1）用三种Jacobi方法求出的全部特征值对应的误差的下降来刻画收敛过程。在每做完一次Givens变换后，将更新后的矩阵的对角元与Tn的真实特征值进行比较，取误差的绝对值的最大值记录下来作为每一步的误差。由于用三种不同的Jacobi方法得到的对角矩阵，其对角元的排列顺序不确定，故须将对角元按从小到大排列后再进行比较。可以得到如下图像：从上面两图可以看出，使用古典Jacobi方法，矩阵的对角元最快地趋于真实特征值。阈值Jacobi比循环Jacobi稍微快一些。（2）用进行相同次数的Givens变换后，矩阵非对角元的F-范数值的下降来刻画收敛过程。为了比较三种不同的方法的效率，下面将经过相同次数的Givens变换后三种不同的Jacobi方法分别对应的非对角元的F-范数值的大小显示在一张图上，由此来进行比较。由于阈值Jacobi方法的效率也部分依赖于消元容限的选取，这里需要指出消元容限。采用的消元容限为：1=1nE0F,n=1nE0F 从上面两图可以看出，仅从进行相同次数的Givens变化的角度来看，应用古典Jacobi方法效率最高，阈值和循环Jacobi差距不大。从理论上来说也是合理的，这是由于古典Jacobi方法每次都是将模最大的元变为0，所以效率最高。但是，对于较大矩阵的计算机来说，计算是容易的，搜索模最大的元素的代价是巨大的。故经典Jacobi用计算机实际执行起来效率不高。而循环Jacobi算法扫描整个矩阵，按预定的次序将元素变为0，省去了搜索最大元的代价。阈值Jacobi则是每一次消元时只对大于容限的值进行处理。练习6.20 用二分法加原点位移反幂法求解在(1,2)之间的特征值。实验步骤：1.strum序列pii=0n的符号相同数计算采用讲义给出的如下算法：q1=p1()qi=i-i2qi-1,i=2,3,.符号相同数恰好为序列q1,q2,中非负数的个数。2.利用循环依次计算出(0,1)区间上的全部特征值的区间估计。3.对上面的区间估计在利用反幂法得到特征值的近似。实验数据：n=100和101时，直接利用MATLAB eig函数，得到T100、T101位于(1,2)区间上的特征值如下表：（左边为n=100时的情形，右边为n=101时的情形）1.018011838053361.05381288632798 1.072672936029351.108523288446921.128231163049261.164079310226431.184632770116621.220428253414641.241823192319241.27751666762569 1.299747101617221.335290401040681.358348460369271.393694651773911.417570575545501.452674019855841.477356153574281.512172559783251.537647355770061.572133833587011.598385854288561.632500964366861.659512888555201.69321669024263 1.720969322112151.754223418670571.782695699829051.81546328107341.844632305421991.876878187732111.906719219225161.938409882887661.96889637615930设置区间误差为10-5,利用二分法，得到的特征值区间估计如下：（左边两列值为n=100时的情形,右边两列值为n=101时的情形）1.0180051.0180131.053810121.053817749 1.072671.0726781.1085205081.1085281371.1282271.1282351.1640777591.1640853881.1846311.1846391.2204208371.2204284671.2418211.2418291.2775115971.2775192261.2997441.2997511.3352890011.3352966311.3583451.3583531.3936920171.3936996461.4175641.4175721.4526672361.4526748661.4773561.4773641.5121688841.5121765141.5376431.5376511.5721282961.5721359251.5983811.5983891.6324996951.6325073241.6595081.6595151.6932144171.6932220461.7209631.720971.7542190551.7542266851.7826921.78271.8154602051.8154678341.8446271.8446351.8768768311.876884461.9067151.9067231.9384078981.9384155271.9688951.968903对上面得到的区间估计值，利用原点位移反幂法得到的特征值近似为：1.0180118381.0538128861.0726729361.1085232881.1282311631.164079311.184632771.2204282531.2418231921.2775166681.2997471021.3352904011.358348461.3936946521.4175705761.452674021.4773561541.512172561.5376473561.5721338341.5983858541.6325009641.6595128891.693216691.7209693221.7542234191.78269571.8154632811.8446323051.8768781881.9067192191.9384098831.968896376二分法可以很快地缩小特征值的存在区间，即得到特征值的近似值。在利用原点位移反幂法，即可得到很精确的特征值。练习6.21：用对称隐式QR方法求解全部特征值：解：利用Wilkinson位移量和“驱逐出境”策略，得到隐式对称QR方法。隐式对称QR算法:（1）输入T（n*n的对称三对角矩阵）（2）收敛性判定：（i）把所有满足条件Ti+1,1=Ti,1+1Ti,i+Ti+1,i+1u的Ti+1,1和Ti,i+1置零。其中u为机器精度。（ii）确定最大的非负整数m和最小的非负整数l,使得T=T11000T22000T33其中T33为（m*m）的对角阵，T22为(n-l-m*n-l-m)不可约的三对角阵。（3）如果m=n,则输出相关信息，结束；否则，进行下一步。（4）QR迭代，对T22用带Wilkinson