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

-8
主特征值误差的下降曲线 乘幂法 Atiken加 速 Rayleigh加 速
3 2.5 2 1.5 1 0.5 0
1.1
1.2
1.3
1.4
1.5 1.6 迭代步数
1.7
1.8
1.9 x 10
2
4
从上面两图可以看出， 初始阶段乘幂法会剧烈震荡， 在迭代一段时间后会突然特征值误差变 大。在迭代 6000 次后，震荡现象不再出现。从第二张图中可以看出，采用 Atiken 和 Rayleigh 商加速得到的误差下降较快。 由于 Atiken 加速和 Rayleigh 商加速均是对特征值进行的， 对特征向量无影响， 故特征子空间 的距离的下降情形一致。下图为特征子空间距离的下降曲线。
练习 6.20 用二分法加原点位移反幂法求解在(1,2)之间的特征值。 实验步骤： 1.strum 序列 pi μ ������ i=0 的符号相同数计算 采用讲义给出的如下算法： q1 = ������1 (������) ������������2 q i ������ = ������������ − ������ − , ������ = 2,3, …. ������������−1 ������ 符号相同数恰好为序列q1 ������ , ������2 ������ , …中非负数的个数。
从上面序列可以看出，迭代一次后特征值并不靠近真实特征值，但利用反幂法，迭代 20 步 后，精度已达到10−10 级别。 从特征向量的距离的角度来看：
10 10 10
特征子空间距离大小
0
n=100, 特 征 子 空 间 距 离 的 下 降 曲 线
-1
-2
10 10 10 10 10 10
-3
-4
-5
10
-250
阈值 古典 循环 0 0.5 1 1.5 2
Givens 变 换 次 数 2.5 3 3.5 4 4.5 x 10 5
4
10
-300
10
50
n=101, 阈 值 、 循 环 及 古 典 Jacobi方 法 中 非 对 角 线 元 F范 数 值 的 下 降
10
0
10
非 对 角 线 元 F范 数 的 大 小
3.9962066557820 6
1.69099045876919e-0 9
3.9962066574740 9
7.54951656745106e-1 5
从上面两表可以看出，n=100 时，利用乘幂法和降阶法求出的特征值误差为10−10 和10−9量 级，而利用同时迭代法求出的特征值误差为10−15 和10−15 量级。n=101 时，利用乘幂法和降 阶法求出的特征值误差为10−16 和10−9 量级，而利用同时迭代法求出的特征值误差为10−15 和10−15 量级。 由此可以看出降阶法求出的第二个主特征值的精度远低于同时迭代法求出的第二个主特征 值精度。 理论分析如下， 因为使用降阶法， 第二个主特征信息的计算精度会受到第一个主特征信息计 算精度的影响，由于舍入误差的积累，特征信息的逐次求解过程，势必造成计算精度越来越 差。 而子空间同时迭代法可以同时求出矩阵的前若干个主特征值信息， 避免了舍入误差的积 累，故对较靠后的主特征值也有很好的精度。
3.9990325645839 8 3.9961311942671 9
降阶法+乘幂 法
3.9990325644623 3 3.9961311960320 9 0
误差
1.21668009001041e-1
同时迭代法
3.9990325645839 8 3.9961311942672
误差
3.10862446895044e-1 5 6.21724893790088e-1 5
特征值误差
阈值 古典 循环
1.2 1 0.8 0.6 0.4 0.2 0
0
1
2
3 Givens 变 换 次 数
4
5 x 10
6
4
从上面两图可以看出，使用古典 Jacobi 方法，矩阵的对角元最快地趋于真实特征值。阈值 Jacobi 比循环 Jacobi 稍微快一些。 （2）用进行相同次数的 Givens 变换后，矩阵非对角元的 F-范数值的下降来刻画收敛过程。
解：采用教科书上 8.1 算法编写程序。在 n=100 时，由于 1,1,1, … ,1 T 在主特征值对应的特 征向量上投影为 0，如果取v0 = 1,1,1, … ,1 T ，若无舍入误差的影响，就会导致乘幂法求解 出的特征值为第二主特征值。故这里取初始向量v0 为随机向量。为了便于比较，在每一步迭 代中，将 Atiken 加速值和 Rayleigh 商加速值分别计算出来，并作出如下图像：
从特征向量的角度来看，可得迭代过程中特征向量与用 eig 函数求出的特征向量的子空间之 间的距离如下表： 2.1073424255447e-08 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
从上表可以看出，进行一次迭代后计算得到的特征向量与 eig 函数求出的特征向量子空间的 距离已为10−8 数量级，故可以认为两者相同，即反幂法确有“一次收敛”到特征向量的特 性。 不过并未观察到后续迭代会起到反作用的现象 （即做多次反幂法后， 特征子空间的距离反而 可能会变大） ，这可能是由于 MATLAB 计算的精度太高导致。
特征值问题数值计算上机实验
练习 6.16： 取初始向量为v0 = 1,1,1, … ,1 ������ ,用乘幂法计算主特征值及其相应的特征向量。 请绘制主特征值误差下降曲线，以及特征子空间距离的下降曲线。然后采用 Atiken 加速 技巧和 Rayleigh 商技术分别对算法进行加速，并完成类似的工作。
为了比较三种不同的方法的效率，下面将经过相同次数的 Givens 变换后三种不同的 Jacobi 方法分别对应的非对角元的 F-范数值的大小显示在一张图上，由此来进行比较。由于阈值 Jacobi 方法的效率也部分依赖于消元容限的选取，这里需要指出消元容限。采用的消元容限 为：σ1 = n ������0
练习 6.18 用幂法求解第二主特征值及其特征向量；用同时迭代法求解Tn 的前两个主特
征值。比较两者的计算效果。
解：对于用乘幂法求解第二主特征值，先采用降阶法，这里取S1 为 Householder 变换矩阵， 先利用乘幂法求出Tn 的主特征值λ1 及其对应的特征向量x1 ，然后求出降阶阵 A2，此处利用 MATLAB 的 eig 函数，发现 A2 与Tn 除λ1 外的其余特征值的差距不大。 再对A2 用乘幂法，可以求出 A2 的主特征值，即Tn 的第二特征值。见下表： 而若用同时迭代法则可直接求出Tn 的前两个主特征值。 采用教科书上算法， 可得计算结果如 下： 当 n=100 真实值
n=100, 阈 值 、 循 环 及 古 典 Jacobi方 法 中 特 征 值 误 差 的 下 降 2 1.8 1.6 1.4
特征值误差
阈值 古典 循环
1.2 1 0.8 0.6 0.4 0.2 0
0
0.5
1
1.5
2 2.5 3 Givens 变 换 次 数
3.5
4
4.5 x 10
5
4
n=101, 阈 值 、 循 环 及 古 典 Jacobi方 法 中 特 征 值 误 差 的 下 降 2 1.8 1.6 1.4
2.利用循环依次计算出(0,1)区间上的全部特征值的区间估计。 3.对上面的区间估计在利用反幂法得到特征值的近似。 实验数据： n=100 和 101 时，直接利用 MATLABeig 函数，得到T100 、T101 位于(1,2)区间上的特征值如下 表： （左边为 n=100 时的情形，右边为 n=101 时的情形）
1.3265605505321 8 1.95501348518822
2.03078984640160 2.03125681228880 2.03106008819891 2.03111963185171 2.03109858241729 2.03110536218386 2.03110536218386
-6
-7
-8
0
2
4
6
8
10 12 迭代步数
14
16
18
20
（其中间断处是由于值为 0） 从上图可以看出特征向量无“一次收敛”特性，但经过 6 次迭代后，误差已达到10−8 数量 级。 当 n=101 时，T101 的靠近 2 的特征值为 1.9384098828 和 2 以及 2.06159011711234，即恰好 2 为T101 的特征值，取 q=2，由于此时矩阵T101 − ������������ 奇异，故可对平移量 q 进行一个非常小 的扰动，取ε = 10−15 .对T101 − (������ + ������)������进行反幂法，预计此时反幂法的“一次收敛”特性非 常明显。进行 20 次反幂法得到计算序列结果为： 1.99999999999998 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
-50
10
-100
10
-150
10
-200
10
-250
10
-300
阈值 古典 循环 0 1 2
Givens 变 换 次 数 3 4 5 x 10 6
4
从上面两图可以看出，仅从进行相同次数的 Givens 变化的角度来看，应用古典 Jacobi 方法 效率最高，阈值和循环 Jacobi 差距不大。从理论上来说也是合理的，这是由于古典 Jacobi 方法每次都是将模最大的元变为 0，所以效率最高。但是，对于较大矩阵的计算机来说，计 算是容易的， 搜索模最大的元素的代价是巨大的。 故经典 Jacobi 用计算机实际执行起来效率 不高。而循环 Jacobi 算法扫描整个矩阵，按预定的次序将元素变为 0，省去了搜索最大元的 代价。阈值 Jacobi 则是每一次消元时只对大于容限的值进行处理。