数值分析矩阵特征值特征向量计算

合集下载

数值分析—矩阵特征值问题的数值计算

( vk +1 )i vk = α1 x1 , lim = λ1 . k k →∞ λ k →∞ ( v ) k i 1
lim
可见，当 k 充分大时， vk 近似于主特征向量（相差一个常数倍）， vk +1 与 vk 的对应非零分量的比值近似于主特征值。在实际计算中，需要对计算结果进行规范化。因为当 λ1 <1 时， vk 趋于零；当 λ1 >1 时， vk 的非零分量趋于无穷，从而计算时会出现下溢或上溢。为此，对向量
λ1 ≥ λ2 ≥ L ≥ λn ,
∑x y
i =1 i
n
i
为向量 x 和 y 的内积。
定理 8.3 设 A 为 n 阶实对称矩阵，其特征值都为实数，排列为对应的特征向量 x1 , x2 ,L , xn 组成正交向量组，则有 1) 对任何非零向量 x ∈ R n ，有 λn ≤ R( x) ≤ λ1 ， 2) λ1 = max R ( x) = R( x1 ) ，
Ζ = ( z1 , z2 ,L, zn )T ∈ R n ，记 max( Ζ) = zi ，其中 zi = Ζ ∞ ，这样，我们有如下乘幂法的实用
的计算公式：任取 v0 = u0 ≠ 0 ，对于 k = 1, 2,L 分别计算 vk = Auk −1 , uk = vk / max(vk ). 求出对应矩阵的主特征向量和特征值的近似值，有下面的定理。定理 8.4
m1 0 M = 0 0 0
称为质量矩阵，而
0 m2 0 0 0
0 0 m3 0 0
0 0 0 m4 0
0 0 0 0 m5 0 0 −k4 k 4 + k5 − k5 0 0 0 − k5 k5

数值分析第3章矩阵特征值与特征向量的计算

第三章矩阵特征值和特征向量计算
工程实践中有许多问题，如桥梁或建筑物的振动，机械
机件、飞机机翼的振动，及一些稳定性分析和相关分析可转化为求矩阵特征值与特征向量的问题。
设A (aij )nn是n阶方阵, 如果数和 n 维非零向量x满足 Ax x，则称为 A 的一个特征值， x称为矩阵A对应于的特征向量。
（2）在正交相似变换下，矩阵元素的平方和不变。设A (aij )nn ,
n
n
U为正交矩阵，记B U T AU (bij )nn ,则
ai2j
bi2j
i, j1
i, j1
Jacobi方法的基本思路：通过一次正交变换,将A中一对非
零的非对角元素化成零,并且使得非对角元素的平方和减少。
反复进行上述过程，使变换后的矩阵的非对角元素的平方和
2a pq
(
/
4)则有a
(1) pq
a(1) qp
1
0 1 0.93
1 0 0 0.93 1
0
0
1
0
0
0.93
1
0 1/ 0.93 1 0
0 0.93 1/ 0.93
按算法迭代3次， 3.0000954,与准确值3的误差小于10－4，u
(1, 0.9992431, 0.9991478)T 与准确值(1, -1,1)T比较，残差 r 0.001.
空间中的二维坐标旋转矩阵。
坐标旋转矩阵U pq ( )是正交矩阵.
设A为实对称矩阵，且apq aqp 0,若记
A(1)
U
T pq
AU
pq
(ai(j1) )
aaq((p1q1p))
a pp a pp
cos2 aqq sin2 sin2 aqq cos2

数值分析--第四章--特征值特征向量计算（乘幂法）

数值分析--第四章--特征值特征向量计算（乘幂法）
摘要：n阶⽅阵A满⾜AX=λx，λ为矩阵A的特征值，x为特征值对应的特征向量。

⼀.乘幂法（求模最⼤特征值及对应特征向量）
设矩阵A有n个相性⽆关的特征向量x1,x2,x3,.....xn，相应的特征值λ1,λ2,λ3,.....λn（由⼤到⼩排列）。

迭代法引⼊：上⼀章学了迭代法求解线性⽅程组Ax=b的解，给定任⼀的初始值v0，不断迭代可以得到Ax=b的解。

同理，给定任⼀⾮零的n维向量v0，不断迭代可以得到矩阵A的特征向量，
对于初始向量v0可以由A的n个线性⽆关的特征向量表⽰：
带⼊迭代⽅程中：
当迭代次数k趋近于⽆穷⼤时，可得到最⼤特征值λ1对应的特征向量a1x1(与x1线性相关)
同理，当迭代次数趋近于⽆穷⼤时，可得到绝对值最⼤的特征值，λ1
其中，m表⽰向量中的绝对值最⼤的那个元素值
如何利⽤迭代法求解按模最⼤特征值和特征向量
说明：
1.初始值：我们给定初始值x0=[1,1,1]^T,取特征值1
2.第⼀次迭代：
对应的近似特征值取：
3.第⼆次迭代：
⼆.改进乘幂法
这个规范化处理的⽬的：防⽌数据溢出或是数据消失
从上⾯可以看出，改进乘幂法即是每次迭代出的特征向量都进⾏⼀次规范化处理。

数值分析QR方法求矩阵特征值和特征向量

四.实验代码: function [H,B]=Hessenberg(A) n=length(A);B=eye(n);for k=1:n-2X=zeros(n-k,1);H=eye(n);for i=1:n-kX(i)=A(i+k,k);enda=max(abs(X));if a==0.0breakendX=X/a;c=X(1);b1=sqrt(sum(X.^2));if X(1)>=0b1=-b1;endX(1)=X(1)-b1;b=b1^2-b1*c;H0=eye(n-k)-X*X'/b;for i=1:n-kfor j=1:n-kH(i+k,j+k)=H0(i,j);end endA=H*A*H;B=B*H;endH=A;一.实验题目:QR方法求矩阵的特征和特征向量二.设计目的:学会利用镜面变换进行矩阵的QR分解及利用将幂法求特征值和特征向量,熟悉Matlab编程环境。

三.设计原理:利用镜像变换将A相似变换为Hessenberg B矩阵。

记录变换矩阵。

运用Householder矩阵进行QR分解,QR方法为:B1=BB1=Q1R1B2=R1Q1....Bm=QmRmBm+1=RmQmBm+1与Bm相似,从而特征值相等。

再利用原点位移的反幂法求B(或A)的特征向量。

反幂法用来计算矩阵按模最小的特征值及其特征向量,也可用来计算对应与一个给定近似特征值的特征向量。

设A∈R n×n为非奇异矩阵,A 的特征值依次记为|λ1|≥|λ2|≥|λ3|≥…≥|λn |,相应的特征向量为x1 ,x2,…,x n,则A-1的特征值为|1/λn|≥|1/λn-1|≥…≥|1/λ1 | ,相应的特征向量为x n ,x.所以计算A的按模最小的特征值λn的问题就是计算n-1,…,x1A-1的按模最大的特征值问题。

对于A-1应用幂法迭代(称为反幂法),可求得矩阵A-1的主特征值1/λn,从而求得A的按模最小的特征值λn。

西安科技大学研究生数值分析课件7章矩阵特征值与特征向量计算

7 矩阵特征值与特征向量地计算设A 为n 阶方阵,所谓A 地特征值问题是求数λ和非零向量x ,使x Ax λ=成立.数λ称作A 地一个特征值,非零向量x 称作与特征值λ对应地特征向量.求给定方阵地特征值与特征向量是先求解特征方程()||0E A ϕλλ=-=然后对应于每一个特征值i λ,再求解退化地齐次线性方程组()0i E A x λ-=从而得到A 地特征值i λ及对应地特征向量x .但是这种方法计算机很大,计算过程复杂,因此有必要研究相对简单地数值解法.本章主要介绍三类计算特征值地方法：计算大型（稀疏）矩阵主特征地幂法与反幂法,计算中小型（实对称）矩阵全部特征值地Jacobi 法,计算中小型矩阵全部特征值地QR 法.7.1 特征值估计在矩阵特征值计算中,有时需要对特征值所在范围给出一个估计.这里介绍一种从矩阵地元素出发,运用较简便地运算估计特征值地方法.定义7-1 设()n m ij A a C ⨯=∈,称由不等式||ii i z a R -≤在复平面上确定地区域为矩阵A 地第i 个盖尔圆（Gerschgorin 圆）,并用i G 表示.其中1||ni ij j j i R a =≠=∑称为盖尔圆i G 地半径(1,2,,)i n =.定理7-1 矩阵()n m ij A a C ⨯=∈地一切特征值均落在它地n 个盖尔圆地并集中,即1(1,2,,)ni jj G i n λ=∈=.证明设λ是A 地任一特征值,12(,,,)T n x x x x =是λ对应地特征向量.令01||max ||i i i nx x ≤≤=,则00i x ≠.由Ax x λ=,可得001()ni j j i j a x x λ==∑.即∑≠==-ni j j j j i i i i x a x a 000001)(λ于是有 000000011i i jni j j ji ni j j i jji i i R x x ax x aa ≤≤=-∑∑≠=≠=λ这表明任一特征值0i G λ∈,从而也在A 地第n 个盖尔圆地并集中.例7-1 估计矩阵10.10.20.30.530.10.210.310.50.20.30.14A ⎡⎤⎢⎥⎢⎥=⎢⎥-⎢⎥---⎣⎦地特征值范围. 解 A 地4个盖尔圆为：1:|1|0.6G z -≤ 2:|3|0.8G z -≤ 3:|1| 1.8G z +≤ 4:|4|0.6G z +≤画在复平面上其区域如图7-1所示.图7-1 例7-1盖尔圆分布图于是A 地全部特征值就在这4个盖尔圆地并集中.为了更确切地知道某个特征值落在哪个或哪几个盖尔圆地并集中,给出如下第二盖尔圆盘定理.定理7-2 若A 地n 个盖尔圆中,有m 个盖尔圆构成地一个连通域（所谓连通域,是指其中地任意两点都可以用位于该区域内地一条折线连接起来）,且该连通域与其余n m -个盖尔圆严格分离,则在该连通域中恰好有A 地m 个特征值（重特征值按重数重复计算）.特别地,每个孤立地盖尔圆恰有A 地一个特征值（证明从略）.由定理2可知,在例1中2G 与4G 中各有A 地一个特征值,而1G 与3G 构成地连通部分中有两个特征值,但不能确定这两个特征值具体落在哪个盖尔圆中.例7-2 估计矩阵10.80.50A -⎡⎤=⎢⎥⎣⎦地特征值范围. 解 A 地两个盖尔圆为：1:|1|0.8G z -≤,2:|0|0.5G z -≤在复平面上地区域如图7-2所示.图7-2 例7-2盖尔圆分布图此时只能判断A 地两个特征值落在1G 与2G 地并集中,至于是每个盖尔圆中各有一个特征值还是两个特征值都落在其中一个盖尔圆上则无法确定.实际上,由于1,21(12λ=±,1,2||0.5λ=>,所以两个特征值都不会在盖尔圆2G 中,而是落在盖尔圆1G 中.对于某些矩阵,可利用相似变换矩阵具有相同特征值地性质得到更确切地特征值范围.设()ij n m A a ⨯=,取正数12,,,n d d d 构成对角阵12diag(,,,)n D d d d =,对A 作相似变换,令1()iij n n jd B DAD a d -⨯==,由于B 相似于A ,所以B 与A 地特征值完全相同,又由于B 与A 地主对角线元素对应相等,所以B 与A 地盖尔圆圆心相同.这表明,若适当选取正数12,,,n d d d ,可以改变盖尔圆地半径,从而有可能将相交地盖尔圆分离得到仅含一个特征值地孤立盖尔圆.选取12,,,n d d d 地一般方法是：欲使A 地第i 个盖尔圆i G 地半径大而其余盖尔圆变小,就取1i d >,其余1()j d j i =≠.例7-3 求矩阵2050.841011210A j ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦地特征值范围. 解 A 地3个盖尔圆为：1:|20| 5.8G z -≤,2:|10|5G z -≤,3:|10|3G z j -≤其中1G 与2G 相交,而3G 孤立.记3G 中所含地一个特征值为3λ,如图7-3所示.为分离2G 与1G ,可以让A 地第3行元素绝对值变大,第3列元素绝对值变小.现取diag(1,1,2)D =,则12050.44100.52410B DAD j -⎡⎤⎢⎥==⎢⎥⎢⎥⎣⎦图7-3 例3盖尔圆分布图图7-4 例7-3分离后盖尔圆分布图其3个盖尔圆分别是：1:|20| 5.4G z '-≤,2:|10| 4.5G z '-≤,3:|10|6G z j '-≤ 显然,B 地盖尔圆是3个孤立地盖尔圆,如图7-4,注意,此情况下,3G '地半径变大了.例7-4 设矩阵()ij n n A a ⨯=按行严格对角占优,则A 可逆.证明由线性代数知,A 可逆地充分条件是||0A ≠,而1||nj j A λ==∏（其中j λ是A 地特征值）,所以只要证明0j λ≠即可(1,2,,)j n =. 设λ是A 地任一特征值,则必存在某个盖尔圆i G 使∑≠=≤-ij ij i ii a R a λ.若0j λ=,则有∑≠≤ij ij ii a a ,而这与A 按行严格对角占优矛盾,故应有0λ≠,由λ地任意性,得||0A ≠.7.2 幂法与反幂法在线性代数中,设A 是n 阶方阵,若A 存在n 个线性无关地特征向量,则称这n 个特征向量构成A 地一个完全地特征向量组.例如,对矩阵320230005A -⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦,110430102B -⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦通过求解特征方程,不难求出A 地三个特征值为1231,5λλλ===,B地三个特征值为1232,1λλλ===.方阵A 可以找到三个线性无关地特征向量,而方阵B 找不到三个线性无关地特征向量.我们称方阵A 可对角化,而B 不可对角化. 7.2.1 幂法幂法地基本思想是构造一个向量序列使之逼近主特征值（矩阵地按模最大地特征值）对应地特征向量,然后求出主特征值.该方法简单易行,但收敛速度较慢.现设()ij n n A a ⨯=有一个完全地特征向量组12,,,n x x x ,其对应地特征值是12,,,n λλλ.已知A 地主特征值是单根1λ,即特征值满足条件12||||||n λλλ>≥≥任取一个非零初始向量0u ,由矩阵A 构造向量序列102210110k k k u Au u Au A u u Au A u++=⎧⎪==⎪⎪⎨⎪==⎪⎪⎩由于A 地完全特征向量组可以作为向量空间n R 地一组基,因此0u 可由12,,,n x x x 线性表示,即有01122n n u a x a x a x =+++ （设10a ≠）于是011122211111121()()k k k k k n n nn kk k i i i k i u A u a x a x a x a x a x a x λλλλλλελ===+++⎡⎤=+=+⎢⎥⎣⎦∑ 其中21()nk i k i i i a x λελ==∑.注意到),,2(11n i i=<λλ,故当k →∞时,0k ε→,因此有111k k u a x λ≈由于1x 是主特征值1λ对应地特征向量,其乘上常数因子11k a λ仍为1λ地特征向量,故当k 充分大时,迭代向量k u 是1λ地特征向量地近似向量.为了利用迭代向量求出主特征值1λ地近似值,设()k i u 表示k u 地第i 个分量,则1111111()()()[]()()()k i i k ik i i k iu a x u a x ελε+++=+ 于是 11()lim()k ik k iu u λ+→∞= 这表明两相邻迭代向量对应分量地比值收敛于主特征值,亦即当k 充分大时,可用两相邻迭代向量地分量比作为主特征值地近似值,即11()()k ik iu u λ+≈若主特征值是A 地r 重实特征值,即12(1)r r n λλλ===≤≤,对应地r 个线性无关特征向量为12,,,n x x x .则有01111()r nkk k i k i i i i i i r u A u a x a x λλλ==+⎡⎤==+⎢⎥⎣⎦∑∑当k 充分大时,11rkk i i i u a x λ=≈∑即k u 仍为主特征值对应地特征向量地近似向量,相邻两迭代向量地分量比仍为主特征值地近似值.综上所述,有定理7-3 设A 是n 阶实矩阵,具有完全地特征向量组,主特征值是r 重根,即112||||||||(1)r r n r n λλλλ++>≥≥≥≤≤则对任意非零初始向量0u ,迭代向量0k k u A u = 满足 111lim(0)rki ikk i u a x a λ→∞==≠∑ ,11()lim ()k ik k iu u λ+→∞= 或 11rk k i i i u a x λ=≈∑,11()()k ik iu u λ+≈ 这样用非零初始向量0u 及矩阵A 构造向量序列{}k u 以计算A 地主特征值1λ及相应地特征向量地方法称为幂法.不过从上面地讨论中可以看到,如果1||1λ>或11<λ,迭代向量k u 当k →∞时,其不为零地分量就会趋于无穷大或趋于零.为克服这个缺点,可以在每步迭代中加上对向量规范化地步骤,使迭代向量地数量级保持在一个稳定地量级上,归纳起来,幂法地计算步骤为：步骤 1 给定非零初始向量0u ,精度12,εε,令00v u =；令(0)10max()v λ=,1=k ；步骤 2 迭代1-=k k Av u ,()1max()k k u λ=,其中)max(k u 表示k u 绝对值最大地分量；步骤3 规范化max()kk k u v u =；步骤 4 若11k k v v ε--<且()(1)112||k k λλε--<,则k v 即为1λ地近似特征向量,()1k λ即为1λ地近似值；否则,1+=k k ,转步骤2继续迭代.例7-5 用幂法计算1.0 1.00.51.0 1.00.250.50.252.0A ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦地主特征值和相应地特征向量,结果见表7-1.表7-1而此题地准确值为1 2.5365258λ= 1(0.748221,0.649661,1.000000)T x =7.2.2 幂法地加速幂法地收敛速度由比值21r λλ=来确定,r 越小收敛越快,而当1r ≈时收敛可能很慢.为了克服这一缺点,常采用原点平移法对幂法进行加速.设B A pE =-,其中p 是待定参数.显然,若A 地特征值为12,,,n λλλ,则B 地相应特征值(1,2,,)i k i n =为12,,,n p p p λλλ---,且A .B 地特征向量相同.这是因为对A 有特征方程||0i A E λ-=,而对B 有特征方程|||()|0i i B k E A p k E -=-+=,所以,i i i i p k k p λλ=+=-另一方面,若i x 是A 地对应i λ地特征向量,即i i i Ax x λ=则 ()()i i i i i i Bx A pE x Ax px p x λ=-=-=-原点平移法地思想是引入矩阵B ,恰当地选择参数p ,使11k p λ=-是B 地主特征值,且其速比2211maxB A p r r p λλλλ-=<=-,这样用幂法求B 地主特征值1k 地收敛速度就快于用幂法求A 地主特征值1λ,而一旦1k 求出,由11k p λ+=可得A 地主特征值,达到了加速地目地.但是为了选取恰当地选择参数p ,需要对A 地特征值地分布地大致了解. 例7-6 设4阶方阵A 有特征值15(1,2,3,4)j jj λ=-=其速比210.9A r λλ=≈.作变换 (12)B A pEp =-=则B 地特征值为12k =,21k =,30k =,41k =-,其速比2112B A k r r k ==<. 设A 地实特征值满足121n n λλλλ->≥≥>若2,n λλ地值可大致估计出,若要求1λ,考察待定参数p 地选取. 在原点平移法通过变换pE A B -=后,不论p 如何选取,矩阵地B 主特征值也只能是在n p λ-或 1p λ-.若希望求1λ,就应选择p ,使1p λ-称为B 地主特征值,即1||||n p p λλ->-这时B 地收敛速比B r 是比值21||/||p p λλ--和1||/||n p p λλ--中地较大者,即211||||max ,||||n B p p r p p λλλλ⎧⎫--=⎨⎬--⎩⎭显然B r 依赖于p 地选取,记做()B r p .为了使应用幂法求B 地主特征值地收敛速度尽可能快,我们希望选择最佳参数*p ,使*()min ()B B r p r p =由B r 地表示式（求二者之间地较大值）和)(*p r B 对)(p r B 地最小化要求,只有当2||||n p p λλ-=-时,()B r p 达到最小.由于2n p p λλ-=-会有得到矛盾地结果（2n λλ=）,所以只能是2()n p p λλ-=--即 *22np λλ+=类似地,若用反幂法求最小特征值n λ,若1n λ-,1λ可大致估计,取最佳平移参数*112n p λλ-+=例7-7 取0.75p =,用原点平移法,计算例7-7中矩阵A 地主特征值.解作变换B A pE =-,则0.2510.510.250.250.50.25 1.25B ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦对B 应用幂法,计算结果见表7-2.即1 1.7865914k ≈,则A 地主特征值1λ为110.75 2.5365914k λ=+=与例7-5比较,上述结果比例7-5迭代15次还好.表7-27.2.3 反幂法设方阵A 按模最小地特征值是n λ,且0n λ≠,则A 可逆.于是,由n n n Ax x λ=,可得11n n nA x x λ-=,这表明1nλ是1A -地主特征值.反幂法就是将幂法应用于1A -,通过求出1A -地主特征值得到A 地按模最小地特征值及其对应地特征向量.定理7-4 设A 是n 阶实矩阵,具有完全地特征向量组,其特征值满足12||||||0n λλλ≥≥≥>则对任意非零初始向量00u v =,按下述方式构造地迭代向量11k k u A v --= ,max()kk k u v u =满足lim max()n k k n x v x →∞=, 1lim max()k k nu λ→∞= /max()k n n v x x ≈,1max()k nu λ≈在实际计算中,可先对A 进行LU 分解,通过求解1k k Ly v -= ,k k Uu y =来求解方程组1k k Au v -=.反幂法地计算步骤为：步骤1 预先取定非零向量00u v =；给定精度12,εε；取(0)0m a x ()nu μ=；步骤2 对矩阵A 作LU 分解,A LU =；令1=k ；步骤3 求解方程组1k k Ly v -= ,k k Uu y = 得到迭代向量k u ；步骤4 规范化max()kk k u v u =步骤5 若11k k v v ε--<且()(1)2||k k n n μμε--<,则k v 即为A 地对应于n λ地近似特征向量,()1k nμ即为n λ地近似值；否则,令1+=k k ,转步骤3继续迭代.7.3 矩阵地两种正交变换本节先介绍镜面(初等)反射变换和平面旋转变换,它们是QR 算法和Jacobi 算法地基础.7.3.1 豪斯荷尔德（House holder ）变换定义7-2 设有方阵B ,若当1i j >+时(,1,2,,)i j n =,0ij b =,则称B 是上Hessenberg 矩阵,即1112121222,1n n n n nn b b b b b b B b b -⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦定义7-3 设向量ω满足21ω=,矩阵2T H E ωω=- （ω是列向量）称为初等反射矩阵,又称House holder 矩阵,记为()H ω,即211212212221212222122()2212n n n n n H ωωωωωωωωωωωωωωωω⎡⎤---⎢⎥---⎢⎥=⎢⎥⎢⎥---⎢⎥⎣⎦其中(1,2,,)i i n ω=是ω地分量.可以证明初等反射阵是对称阵()T H H =.正交阵()T H H E =. 例7-8 设向量0α≠,试证矩阵222TH E ααα=- 是一个初等反射阵. 证明令2αωα=,则 222221||||||||1αωααα=== 由定义7-3,2222TTH E E ααωωα=-=-是初等反射阵.定理7-5 设,x y 是两个不相等地n 维列向量,且22||||||||x y =,则存在一个初等反射阵H,使得Hx y =证明令2||||x yx y ω-=-,由例7-8可知22()()22||||T T Tx y x y H E E x y ωω--=-=-- 是一个初等反射阵.由于22||||()()T T T T Tx y x y x y x x y x x y y y -=--=--+ 注意到22||||||||x y =,即T T x x y y =,又()T T T T x y x y y x == ,故22||||2()T Tx y x x y x -=-从而22()()2||||T T x y x x y x Hx x x y --=--y y x x =--=)(. 例7-9 设1(1,2,2),(1,0,0)T T x e ==,用Householder 变换将x 化为与1e 同方向地向量.解因为2||||3x =,可设13y e =,则22||||||||x y = 取21,1,1)||||T x y w x y -==--,构造Householder 矩阵[]11122212111,1,12123311221T H E ww -⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥=-=--=-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦则13Hx e =推论设向量12(,,,)0T n x x x x =≠,12()||||r sign x x =,且1x r ε≠-,则存在初等反射阵1222||||T T uu H E E uu u ρ-=-=- 使1Hx r ε=- .其中,1(1,0,,0)T ε=,1u x r ε=+,22||||/2u ρ=.设12(,,,)T n u u u u =,则12(,,,)T n u x r x x =+22222122222112111||||[()]221(2)2()n n u x r x x r rx x x x r r x ρ==++++=+++++=+引入初等反射阵地目地,是设法用一系列初等反射阵将原始矩阵约化成上Hessenberg 阵.由于约化过程是逐列进行地,我们先给出计算Hx 地算法步骤,该算法算出H 及r ,使Hx r ε=-,u 地分量冲掉x 地分量.（1）1max ||i i nx η≤≤=；（2）(1,2,,)ii i x x u i n η←==,此步规范化是为避免计算r 时产生溢出；（3） 12211()()nii r sign x x ==∑；（4）11u u r ←+；（5） 1ru ρ=；（6） r r η←；于是初等反射阵1T H E uu ρ-=-,1Hx r ε=-.如果要将H 作用于矩阵A ,设i a 是A 地第i 列向量,则12(,,,)n A a a a =,12(,,,)n HA Ha Ha Ha = 其中,11()()(1,2,,)T T i i i i Ha E uu a a u a ui n ρρ--=-=-=.下面讨论用初等反射阵约化原始矩阵A 称为上Hessenberg 阵地步骤.11121(1)(1)2122211121(1)(1)212212n n n n nn a a a a a a a A A A a a a a a ⎡⎤⎢⎥⎡⎤⎢⎥===⎢⎥⎢⎥⎣⎦⎢⎥⎣⎦步骤1 不妨设(1)210a ≠（否则这一步不需约化）,选择初等反射阵1R ,使(1)12111R a r ε=-,其中： 1(1)(1)2212112(1)1211112111121211111()(())(1)1()2ni i T r sign a a u a r n u r r a R E u u εερρ=-⎧=⎪⎪⎪=+-⎨⎪==+⎪⎪=-⎩∑是维单位坐标列向量令11100U R ⎡⎤=⎢⎥⎣⎦则(2)(2)(2)(1)111213111212111(2)(2)(1)(1)222312112210A a A a A R A U AU a A R a R A R ⎡⎤⎡⎤===⎢⎥⎢⎥⎣⎦⎣⎦其中,(2)11A 是21⨯阵,(2)22a 是2n -维列向量,(2)23A 是2n -阶方阵.步骤k 设对A 已进行了1k -步约化,即111(2)()()()()11121,111,11(2)()()()1222,12,2()()()1,1,()()()1,1,11,()()(),1(2,3,,1)k k k k k k k k k k k n k k k k kn k k k k kk k k k n k k k k k k k k nk k k nkn k nnA U A U k n a a a a a a r a a a a r a a a a a a a a a ----+--++++++==-⎡⎢-⎢=-⎣()()()111213()()22230k k k k k A a A a A ⎤⎥⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎦⎡⎤=⎢⎥⎣⎦其中,()11k A 是(1)k k ⨯-阵,()22k a 是n k -维列向量,()23k A 是n k -阶方阵.设()220k a ≠,选初等反射阵()k R n k -阶,使()221k k k R a r ε=-,其中1ε是n k -维单位坐标向量,可得1()()221,1()221()1,1()(())()nk k k k k ik i k k kk k k k kk nT k k k k r sign a a u a r r r a R E u u ερρ+=++-⎧=⎪⎪⎪=+⎨⎪=+⎪⎪=-⎩∑ 令 00k k E U R ⎡⎤=⎢⎥⎣⎦则 ()()()1112131()()2223()()()111213()12300k k k k k k k k k k k k k k k k k k k k k A a A R A U A U R a R A R A a A R r R A R ε+⎡⎤==⎢⎥⎣⎦⎡⎤=⎢⎥-⎣⎦ 可见1k A +地左上角1k +阶子阵为上Hessenberg 阵,从而约化又进了一步.重复此过程,直到122112211(2)122(3)233(1)1n n n n n nn A U U U AU U U a r a r a r a -----=⨯⨯⨯⎡⎤⎢⎥-⨯⨯⎢⎥⎢⎥=-⨯⎢⎥⎢⎥⎢⎥-⎣⎦使原始矩阵A 在一系列初等反射阵地作用下,约化为上Hessenberg 阵.综上所述,有定理7-6.定理7-6 如果A 是n 阶实矩阵,则存在初等反射阵122,,,n U U U -,使221122n n U U U AU U U C --=（上Hessenberg 阵）例7-10 试证矩阵A 与其约化成为地上Hessenberg 阵C 有相同地特征值.证明记221n P U U U -=,由于初等反射阵是正交对称阵,故122T n P U U U -=,且P 是正交阵,故T PAP C =.于是||||||||||||T T C E PAP E P A E P A E λλλλ-=-=-=-其中T PP E =,||||1T P P =.这表明A 与C 具有相同地特征多项式,即两者有相同地特征值.进一步,设y 是C 地对应于特征值λ地特征向量,即Cy y λ=,则有T PAP y y λ= ()()T T A P y P y λ=这表明T P y 为A 地对应于λ地特征向量,于是求原始矩阵A 地特征值与特征向量可转化为求上Hessenberg 阵C 地特征值和特征向量.定理7-7 若A 是实对称矩阵,则存在初等反射阵122,,,n U U U -使2211221112211()n n n n n U U U AU U U c b b c b C b b c ----⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎣⎦对称三对角阵证明由定理7-6,存在初等反射阵可使A 约化为上Hessenberg 阵C ,当A 是对称矩阵时,C 亦为对称阵,即T C C =,且T C 亦为上Hessenberg 阵,故C 是对称三对角阵.例7-11 用豪斯荷尔德方法将下述矩阵化为上Hessenberg 阵.1437232427A A ---⎡⎤⎢⎥==⎢⎥⎢⎥⎣⎦解（1）对1k =,确定变换阵111000U R ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦,(1)2124a ⎡⎤=⎢⎥⎣⎦ 其中1R 为初等反射阵,使(1)121110R a r ⎡⎤=-⎢⎥⎣⎦(1)12124.472136r a ==≈(1)12111 6.472136244u a r ε⎡⎡⎤+=+=≈⎢⎢⎥⎣⎦⎣⎦11121()2)28.94427r r a ρ=+≈[]1111110 6.4721361 6.472136401428.944270.4472070.8944230.8944230.447216TR E u u ρ-=-⎡⎤⎡⎤=-⎢⎥⎢⎥⎣⎦⎣⎦--⎡⎤=⎢⎥-⎣⎦（2）计算(1)122R A .记(1)221232(,)27A a a ⎡⎤==⎢⎥⎣⎦,于是 (1)1221112 3.1304967.155419(,) 1.788855 1.341640R A R a R a --⎡⎤==⎢⎥-⎣⎦其中,111111111()()(1,2)T T i i i i R a E u u a a u a u i ρρ--=-=-=（3）计算(1)121A R 及(1)1221()R A R ,即求 1(1)121211(1)1223373.1304967.1554191.788855 1.341640T T T b A R b R R R A b ⎡⎤--⎡⎤⎡⎤⎢⎥⎢⎥==--⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥-⎣⎦⎣⎦7.6026340.4472127.800030.3999990.399999 2.200000-⎡⎤⎢⎥=-⎢⎥⎢⎥-⎣⎦其中,11111()(1,2,3)T T T Ti i i b R b b u u i ρ-=-=（4）计算2111A U AU =.(1)12121(1)1221447.6026340.4472124.4721367.8000030.39999900.399999 2.2000000A R A r R A R ⎡⎤--⎡⎤⎢⎥⎢⎥⎢⎥==--⎢⎥⎢⎥-⎢⎥-⎢⎥⎣⎦⎢⎥⎣⎦为上Hessenberg 阵.7.3.2 平面旋转变换定义7-4 称矩阵111(,)111i j csi P i j scj ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第列第列第行第行为平面旋转矩阵,又称Givens 矩阵,其中cos c θ=,sin s θ=.平面旋转阵(,)P i j 是一个正交阵,与单位阵只有在(,),(,),(,i i i j j j j i四个位置上地元素不一样,用其左乘矩阵A 只改变A 地第i 行和第j 行元素.设12(,,,)T n x x x x =则平面旋转变换Px y =地结果为⎪⎩⎪⎨⎧≠=+-=+=ji k x y cx sx y sx cx y k kj i j j i i ,若令/i c x =,j s x =, 则平面旋转变换向量y 地第i个分两为22j i x x +,第j 个分量为0,其余分量即为x 对应地分量.和初等反射变换一样,用平面旋转变换也可以将一个方阵化为上Hessenberg 矩阵,也可以将将一个方阵化为上三角矩阵.7.4 QR 算法7.4.1 矩阵地QR 分解定理7-8 设A 是可逆矩阵,则存在正交矩阵121,,,n P P P -使121()n P P P A R -=上三角矩阵且R 地主对角元素0(1,2,,1)ii r i n >=-.证明若10(2,3,,)j a j n ==,则A 地第一列不需约化.若有某个 10(2)j a j n ≠≤≤,则可选择1(1,)j P j P =使A 地第一列中第j 个元素变为零.一般地,可设平面旋转矩阵12131,,,n P P P ,使(2)(2)11121(2)(2)222113122(2)(2)200nn nn nn r a a a a P P P A A a a ⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎣⎦记111312nP P P P =,则12P A A =.同理,若(2)20(3,4,,)j a j n ≠=,可选取23242,,,n P P P 使(2)(2)(2)1112131(3)(3)22232(3)(3)2212323333(3)(3)3nn n n n n nn r a a a r a a P P P A A a a a a -⎡⎤⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎣⎦记2223nP P P =,则213P P A A =.重复上述过程,可得一系列正交阵121,,,n P P P -使11121222121n n n nn r r r r r P P P A R r -⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎣⎦ 定理7-9 （矩阵地QR 分解）如果n 阶实矩阵A 可逆,则A 可分解为一正交阵Q 和上三角阵R 地乘积,即A QR =,且当R 地对角元素都为正数时分解唯一.证明由定理8知存在正交阵11,,n P P -使121n P P P A R -=为上三角阵,记121T n Q P P P -=,于是T Q A R =由于(1,2,,1)i P i n =-是正交阵,则T Q 亦为正交阵,故A QR =. 若A 有两种QR 分解,记为1122A Q R Q R ==其中12,R R 为上三角阵且主对角元素都为正数,12,Q Q 为正交阵,于是12121T Q Q R R -=注意121R R -是上三角阵地乘积,结果仍为上三角阵,而12,TQ Q 是正交阵,所以121R R -也应是正交阵.若记121D R R -=,由其上三角性T D 应是下三角阵,1D -应是上三角阵；由其正交性由1T D D -=,故D 只能是对角阵,且有2T D D D E ==.又因12,R R 地主对角元素都为正数,即有222212diag[,,,]diag[1,1,,1]n D d d d E ===故1(1,2,,)i d i n ==,则D E =,于是12R R =,12Q Q =.例7-12 求矩阵⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=212240130A 地QR 分解. 解方法1：利用初等反射阵进行QR 分解令(0)1(0,0,2)T a =,取(0)112||||2d a ==,则)2,0,2(81211)0(111)0(11-=--=e d ae d a u1110012010100TH E u u ⎡⎤⎢⎥=-=⎢⎥⎢⎥⎣⎦,⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=1302402121A H 再令(0)2(4,3)T a =,取(0)222||||5d a ==,则(1)2212(2)22121,3)||||T a d e u a d e -==--2224312345TH E u u ⎡⎤=-=⎢⎥-⎣⎦令2210014305534055H H ⎡⎤⎢⎥⎢⎥⎡⎤⎢⎥==⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥-⎣⎦于是21212051002H H A R ⎡⎤⎢⎥=-=⎢⎥⎢⎥-⎣⎦故123405521243005155002100T TA H H R ⎡⎤-⎢⎥⎡⎤⎢⎥⎢⎥⎢⎥==-⎢⎥⎢⎥⎢⎥-⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦方法2：利用平面旋转阵进行QR 分解. 取1202,0100221221=+==+=s c ,则130********T ⎡⎤⎢⎥=⎢⎥⎢⎥-⎣⎦,132********T A ⎡⎤⎢⎥=-⎢⎥⎢⎥--⎣⎦再取53)3(43,54)3(44222222-=-+-==-+=s c ,则231004305534055T ⎡⎤⎢⎥⎢⎥⎢⎥=-⎢⎥⎢⎥⎢⎥⎣⎦,2313212051002T T A R ⎡⎤⎢⎥=-=⎢⎥⎢⎥-⎣⎦ 故13233405521243005155002100T T A T T R ⎡⎤-⎢⎥⎡⎤⎢⎥⎢⎥⎢⎥==-⎢⎥⎢⎥⎢⎥-⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦例7-13 求矩阵⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=110133044A 地QR 分解,使得R 地对角线元素为正数.解 A A =1地第一列T x ]0,3,4[1=,521=x .用1x 构造镜面反射阵1H ,使得T y x H ]0,0,5[111==,令T y x u ]0,3,1[111-=-=,有⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-=-=10005453053542221111u u u E H T ,⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-==11054005355112A H A 2A 地第2列对角线以下为T x ]1,0[2=,122=x .用2x 构造镜面反射阵2~H ,使得T y x H ]0,1[~222==,令T y x u ]1,1[222-=-=,易得 ⎥⎦⎤⎢⎣⎡=-=01102~222222u u u E H T,⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎦⎤⎢⎣⎡=010100001~122H H 于是有R A H A =⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡-==54001105355333,⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-==010540535305421H H Q容易验证,QR A =.请读者用平面旋转变换对本例地矩阵A 进行QR 分解.7.5.3 QR 算法QR 算法就是利用QR 分解构造一个矩阵序列{}k A ,当k 充分大时,k A 是近似地上三角矩阵,而该上三角阵地对角元素便是原始矩阵A 地全部特征值.设1()n n ij n n A A a R ⨯⨯==∈,对A 做QR 分解,即A QR =其中R 为上三角阵,Q 为正交阵.利用这个分解可得新矩阵（对QR 交换乘积）2T A RQ Q AQ == 由于2A 是1A 经过正交相似变换得到地,因此2A 与1A 有相同地特征值.再对2A 做QR 分解,按上述方式又可得新矩阵3A ,且3A 与2A 也具有相同地特征值.具体地说,其步骤为：设1A A =,做QR 分解111A Q R =求矩阵211111T A R Q Q A Q ==求得k A 后对k A 作QR 分解k k k A Q R =求矩阵1Tk k k k k k A R Q Q A Q +==只要A 可逆,由定理9可知,按上述方法可唯一确定矩阵序列{}k A ,且序列中任意k A 与原始矩阵有相同特征值.因此只要恰当选择正交相似变换阵12,,,k Q Q Q ,使1111111T T TT TT T k k k k k k k k k k k k k k A Q A Q Q Q A Q Q Q Q Q A Q Q Q +----====当k →∞时,逼近一个上三角阵,便可求出A 地全部特征值（为所逼近上三角阵地主对角元素）.可见,QR 算法地关键在于选择正交变换阵(1,2,)k Q k =.从定理7-8地证明看到,正交变换阵k Q 是一系列平面转换矩阵地乘积,这些平面旋转矩阵是用来将k A 地主对角线以下元素约化为零地.如果将QR 算法直接应用于原始矩阵,计算量很大,所以在实际计算中,总是先将原始矩阵用豪斯赫尔德方法约化为上Hessenberg 阵,而后再对上Hessenberg 阵应用QR 算法.可以证明,由上Hessenberg 阵用QR 算法生成地矩阵序列中地每个矩阵仍为上Hessenberg 阵.7.5 雅可比方法雅可比方法是用来计算实对称矩阵地全部特征值及特征向量地一种有效方法.它地基本思想是,通过一组正交相似变换对称矩阵A 化为对角矩阵,得其全部特征值.定理7-10 设A 为n 阶对称矩阵,T C PAP =,其中P 为正交矩阵,则22||||||||F F C A = 证明一方面2222111||||()()nnnFiji i j i A a tr A A λ======∑∑∑另一方面2221||||()()()nTFi i C tr C C tr C C λ====∑由假设()()i i A C λλ=,故22||||||||F F C A =.设n n A R ⨯∈为对称矩阵,(,)P i j 为一平面旋转矩阵,则T C PAP =（其中()ij n n C c ⨯=）地元素计算公式为：（1）22cos sin 2sin cos ii ii jj ij c a a a θθθθ=++22sin cos 2sin cos jj ii jj ij c a a a θθθθ=+-（2）1()sin 2cos 22ij ji jj ii ij c c a a a θθ==-+ （3）第i 行元素和第j 列元素cos sin (,)ik ki ik jk c c a a k i j θθ==+≠ （4）第j 行元素和第i 列元素 cos sin (,)jk kj jk ik c c a a k i j θθ==+≠（5）(,,)lk lkc a l k i j =≠这说明,当A 经过一初等正交相似变换化为C 时,只需按上述公式计算C 地第i 列.第j 列元素,由对称性可得第i 行和第j 行元素,C 地其余元素与A 地对应元素相同.设A 地非对角元素0ij a ≠,我们可选择平面旋转阵(,)P i j ,使T C PAP =地非对角元素0ij ji c c ==.由定理11可选择(,)P i j ,使sin 2cos 202jj iiij ji ij a a c c a θθ-==+=即选择θ,使22(||)4ij ii jja tg a a πθθ=≤-其中定理7-11 设n n A R ⨯∈为对称阵,0ij a ≠为A 地一个非对角元素,则可选择一平面旋转阵(,)P i j ,使T C PAP =地非对角元素0ij ji c c ==且T C PAP =与A 地元素满足下述关系（1）2222(,)ik jk ik jkc c a a k i j +=+≠（2）222222ii jj ii jj ij c c a a a +=++ （3）22(,,)iklk c a l k i j =≠证明由上面地计算ij c 公式直接计算可知（1）成立.由（1）及定理7-10可证（2）.如果用()S A 表示A 地非对角线元素地平方和,()D A 表示A 地对角线元素平方和,则2()()2ijD C D A a =+ ,2()()2ij S C S A a =- 这说明C 地对角线元素平方和比A 地对角线元素平方和增加了22ij a ,C 地非对角线元素平方和比A 地非对角线元素平方和减少了22ij a .下面介绍雅可比方法.首先在A 地非对角元素中选择绝对值最大地元素（称为主元素）,如11||max ||i j lk l ka a ≠=可设110i j a ≠,否则A 已经对角化了.由定理12,选择一平面旋转矩阵111(,)P i j ,使111TAP AP =地非对角元素11110i j j i c c ==. 再选(1)1()lkn n A a ⨯=地非对角元素中地主元素,如 22(1)(1)||max ||0i j lk l ka a ≠=≠由定理12,又可选择一平面旋转矩阵222(,)P i j ,使2212T A P A P =地非对角元素2222(2)(2)i j j i a a ==（注意上次消除了地主元素这次又可能变为不是零）. 继续这个过程,连续对A 实行一系列平面旋转变换,消去非对角线绝对值最大地元素,直到将A 地非对角元素全化为充分小为止,从而求得A 地全部（近似）特征值.定理7-12 （雅可比方法地收敛性）设()ij n n A a ⨯=为实对称矩阵,对A 施行上述一系列平面旋转变换1(1,2,)Tm m m mA P A P m -==则 lim ()m m A D→∞=对角矩阵证明记()()m m lk n n A a ⨯=,()2()m m lk l kS a ≠=∑由定理7-11地（2）可得()212()m m m ij S S a +=-其中 ()()||max ||m m ijlk l ka a ≠= 又由于()2()2()(1)()m m m lk ij l kS a n n a ≠=≤-∑即()2()(1)m m ij S a n n ≤- 由以上得12(1)(1)m m S S n n +≤-- 反复应用上式,即得1102(1)(2)(1)m m S S n n n ++≤->-故 lim 0m m S →∞= 可以证明()lim m ll m a →∞存在(1,2,,)l n =. 下面介绍特征向量地计算.由雅可比收敛定理知,当m 充分大时2112T TTmm P P P AP P P D ≈记12T T T T m m R P P P =,则T m R 地列向量就是A 地近似特征向量.计算Tm R 可采用累积地办法,用一数组R 保存Tm R ,开始时R E ←,以后对A 每进行一次平面旋转变换,就进行计算Tm R RP ←用初等正交阵T m P 右乘R 只需计算R 地两列元素,若记(,)m m P P i j =,则Tm RP 地计算公式为()()cos ()sin (1,,)()()sin ()cos li li lj li li lj l n θθθθ←+⎧⎪=⎨←+⎪⎩R R R R R R关于sin θ和cos θ地计算如下.由定理7-11知,当0ij a ≠时,可选θ满足2tg2ij ii jja a a θ=-方ii jj a a ≠时,由22tg 1tg21tg dθθθ=≡- 得到tg θ地二次方程2tg 2tg 10d θθ+-=解得tg θ=选取tg 0d d θ>=<由此得 |tg |1θ≤可由集合{},,ii jj ij a a a 来计算sin ,cos θθ,设0,||max ||ij ij lk l ka a a ≠≠=,则210tg ,()10cos sin cos ii jj ija a d a d t s d d c t ct sθθθθ-⎧=⎪⎪⎪≥⎧⎪=≡=⎨-<⎨⎩⎪⎪=≡⎪⎪=⋅=≡⎩如果jj ii ij a a a -<<,则12ij ii jja t d a a ≈=-,将c,s 代入定理7-9地（1）中可得ii ii ij jj jj ij ij ji c a ta c a ta c c ⎧=+⎪=-⎨⎪==⎩ 每迭代一次地主要工作是选m A 地非对角线元素中地主元素与计算T 111m m m +++=A P AP .首先计算sin ,cos ,θθ,只需计算1m +A 地第i 列,第j 列元素,再算对称元素,不用做3个矩阵地乘法.计算机计算时,需要两组工作单元,以便存储A （或m A ）和R .可用()2()m m lk l ka ε≠=<∑S 控制迭代终止,其中ε是要求地精度.例7-14 用雅可比方法计算对称阵210121012⎡⎤-⎢⎥--⎢⎥⎢⎥-⎣⎦A = 地特征值.解第1步0=A A ,选非对角线元素中地主元素121(1,2)a i j =-==0,1,1/0.7071068,1/0.7071068d t c s ======T 111100.7071068030.70710680.70710680.70710682⎡⎤-⎢⎥==-⎢⎥⎢⎥--⎣⎦A P AP第2步在1A 中选非对角元素地主元素(1)130.7071068(1,3)a i j =-==0.7071068,0.5176381,0.8880738,0.4597008d t c s ====T 22120.63397460.325057600.325057630.627963000.62798302.366025-⎡⎤⎢⎥=--⎢⎥⎢⎥-⎣⎦A P A P 第3步在2A 中选非对角元素地主元素(2)230.627930(2,3)a i j =-==0.5047869,0.6153960,0.8516540,0.5241045d t c s =-=-==-T 33230.63397460.27683660.17036420.27683663.38644600.170364201.979579⎡⎤--⎢⎥=-⎢⎥⎢⎥-⎣⎦A P A P 第4步在3A 中选非对角元素地主元素(3)120.2768366(1,2)a i j =-==4.971292,0.09958013,0.9950785,0.09909004d t c s ====T 44340.606407200.169525803.4140130.016881400.16952580.016881401.979579⎡⎤-⎢⎥=⎢⎥⎢⎥-⎣⎦A P A P 第5步在4A 中选非对角元素地主元素(4)130.1695258(1,3)a i j =-==4.050038,0.1216293,0.9926842,0.1207395d t c s ==== 2T 255450.58578790.20382521000.203825210 3.4140130.0167579000.016757902.000198--⎡⎤⨯⎢⎥=⨯⎢⎥⎢⎥⎣⎦A P A P 于是A 地特征值为1233.414013, 2.000198,0.5857879λλλ===A 地精确特征值为12(1 3.414214λ=≈,22λ=,32(10.585786λ=-≈ 且可逐步求出412345T T T T T T R P P P P P =地列向量,即得A 地近似特征向量.雅可比方法是一个求对称矩阵A 地全部特征值及特征向量地迭代方法,精确度较高,但计算量较大,对稀松带状矩阵经过平面旋转变换后其稀松带状将被破坏,所以很少使用.习题71．设911203111(2102113810A j j B ⎡⎤⎡⎤⎢⎥⎢⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦试估计它们地特征值所在地范围.2．编写幂法程序,并求矩阵732341213A -⎡⎤⎢⎥=⎢⎥⎢⎥--⎣⎦地主特征值及对应地特征向量（准确到小数点后3位）.3．若p 是A 地特征值j λ地一个近似值,且||||()j i p p i j λλ-<-≠则1j pλ-是1()A pE --地主特征值.试用反幂法求矩阵134231111A ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦地最接近于6地特征值及对应地特征向量.4．设有向量(2,1,2)Tx=,试构造初等反射阵H,使(3,0,0)THx=.5．设(2,3,0,5)Tx=,(1,0,0,5)Te=,用Householder变换化x为与e同方向向量.6．设031042212A⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦,求其QR分解.7．设221022212A⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦,求其QR分解.8．利用初等反射阵将134312421A⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦正交相似约化为对称三对角阵.9．试用平面旋转变换阵对矩阵A作QR分解,其中111021245A⎡⎤⎢⎥=-⎢⎥⎢⎥-⎣⎦.10．按下列要求编写程序框图.（1）将一般矩阵用豪斯赫尔德方法约化称上Hessenberg阵.（2）对矩阵作QR分解.（3）对上Hessenberg阵应用QR算法求全部特征值及相应地特征向量.11．用QR算法求矩阵120211013A⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦地全部特征值.12．设A是对称矩阵,λ和(1)x x=是A地一个特征值及相应地特征向量.又设p是一个正交阵,使1(1,0,0,,0)Tpx e==证明T=是第一行和第一列除了λ外,其余元素均为零.B PAP。

数值分析第四章矩阵特征值与特征向量的计算

192.9996. 973
12
➢ 幂法的加速—原点移位法
应用幂法计算矩阵A的主特征值的收敛速度主要
由比值 r=|2/1|来决定, 但当r接近于1时, 收敛可能
很慢. 这时可以采用加速收敛的方法.
引进矩阵
B=A－0I
其中0为代选择参数. 设A的特征值为1, 2, …, n, 则B的特征值为1－0, 2－0, …, n－0, 而且A, B
10
2 1 0 例用幂法求矩阵 A 0 2 1
0 1 2
的按模最大的特征值和相应的特征向量.
取 x(0)=(0, 0, 1)T, 要求误差不超过103.
解 y 0 x 0 0 ,0 ,1 T ,
x 1 A 0 0 y , 1 , 2 T , 1 m x ( 1 ) ) a 2 , x
y(1)
x(1)
1
(0,0.5,1)T
x ( 2 ) A ( 1 ) 0 . 5 y , 2 , 2 . 5 T ,2 m x ( 2 ) ) 12 1a . 5 ,
y(2)
x(2) 2
(0.2,0.8,1)T
x ( 3 ) A ( 2 ) 1 . 2 y , 2 . 6 , 2 . 8 T ,3 m x ( 3 ) ) 2 a . 8 ,
x
(
k
1
)
Ax
(k )
A k1 x (0)
在一定条件下, 当k充分大时:
1
x ( k 1) i
x
( i
k
)
相应的特征向量为: x(k1) 4
➢ 幂法的理论依据
n
对任意向量x(0), 有 x(0) tiui ,
i1
x(k1) Ax(k) Ak1x(0)

特征值和特征向量计算的数值方法

特征值和特征向量计算的数值方法在数学和计算机科学领域中，特征值和特征向量是非常重要的概念。

特征值和特征向量的计算有许多不同的数值方法，本文将介绍其中一些常见的数值方法，并分析它们的优劣和适用范围。

一、特征值和特征向量的定义在矩阵理论中，给定一个n×n的矩阵A，如果存在一个非零向量v和一个标量λ，使得Av=λv，那么称v为矩阵A的特征向量，λ为矩阵A的特征值。

特征值和特征向量的计算可以帮助我们理解矩阵的性质以及解决一些实际问题。

二、幂法幂法是计算特征值和特征向量的常用数值方法之一。

幂法的基本思想是通过多次迭代，逐渐逼近矩阵的特征值和特征向量。

具体操作如下：1. 初始化一个非零向量b0；2. 进行迭代计算：bi+1 = A * bi / ||A * bi||；3. 取出近似特征向量的最后一列：v = bn；4. 进行迭代计算特征值：λ = (Av)T * v / (vT * v)。

幂法的主要优点是简单易懂，且只需要进行矩阵向量乘法和内积计算。

然而，幂法仅能求取具有最大特征值的特征向量，而且对于存在多个特征值相等的情况并不适用。

三、反幂法反幂法是幂法的一种改进方法，用于求取矩阵A的最小特征值和对应的特征向量。

反幂法的基本步骤如下：1. 初始化一个非零向量b0；2. 进行迭代计算：bi+1 = (A - μI)^-1 * bi / ||(A - μI)^-1 * bi||；3. 取出近似特征向量的最后一列：v = bn；4. 进行迭代计算特征值：λ = (Av)T * v / (vT * v)。

反幂法的改进之处在于引入了矩阵的逆运算，通过使用矩阵A减去一个合适的常数μ乘以单位矩阵来实现。

反幂法适用于矩阵A的特征值接近于μ的情况。

四、QR方法QR方法也是一种常用的特征值计算方法，它适用于求解所有特征值以及对应的特征向量。

QR方法的基本思想是将一个矩阵分解为正交矩阵Q和上三角矩阵R的乘积，然后迭代地将矩阵A转化为更接近上三角形的形式。

矩阵特征值的数值解法

矩阵特征值的数值解法矩阵的特征值是在矩阵与其特征向量之间的关系中的数值解。

特征值在各个领域中都有广泛应用，包括物理、工程、金融等。

在解决实际问题时，我们经常需要计算矩阵的特征值，因此研究如何求解矩阵特征值的数值方法是非常重要的。

1. 幂迭代法（Power Iteration）幂迭代法是求解矩阵特征值的一种简单而常用的数值方法。

它的基本思想是通过不断迭代矩阵与向量的乘积，使得向量趋近于该矩阵的一个特征向量。

具体步骤如下：（1）初始化一个非零的初始向量x。

（2）进行迭代计算，即$x^{(k+1)}=Ax^{(k)}/，Ax^{(k)}，$。

（3）当向量x的相对误差小于一些预设的精度要求时，停止迭代，此时的x即为矩阵A的一个特征向量。

（4）将x带入特征值的定义式$\frac{Ax}{x}$，计算出特征值。

幂迭代法的优点是简单易实现，计算速度较快，缺点是只能求解特征值模最大的特征向量，而且对于存在特征值模相近的情况，容易收敛到错误的特征值上。

2. QR迭代法（QR Iteration）QR迭代法是一种较为稳定的求解矩阵特征值的数值方法。

它的基本思想是通过不断进行QR分解，使得矩阵的特征值逐渐收敛。

具体步骤如下：（1）将矩阵A进行QR分解，得到正交矩阵Q和上三角矩阵R，令$A_1=RQ$。

（2）将$A_1$再次进行QR分解，得到新的矩阵$A_2=R_1Q_1$。

（3）重复步骤（2），直到得到收敛的矩阵$A_k$，此时$A_k$的对角线上的元素即为矩阵A的特征值。

QR迭代法的优点是对于特征值模相近的情况仍然能够收敛到正确的特征值上。

缺点是每次QR分解都需要消耗大量的计算量，迭代次数较多时计算速度较慢。

3. Jacobi迭代法（Jacobi's Method）Jacobi迭代法是一种通过对称矩阵的对角线元素进行迭代操作，逐步将非对角元素变为零的求解特征值的方法。

具体步骤如下：（1）初始化一个对称矩阵A。

矩阵特征值和特征向量的数值解法

(k )
的常用方法是迭代每一步对向量 u
规范化。引入函数 max（ u
(k )
），它表示取
向量 u (k ) 中按模最大的分量，例如， u (k ) =(2,-5,4)T,则 max( u (k ) )=-5,这样
u (k ) 的最大分量为 1，即完成了规范化。 (k ) max (u )
7.1 幂法
7.1.1 幂法原理及实用幂法幂法主要用于求矩阵按模最大的特征值和相应的特征向量。设矩阵,2,..., n) 满足：
| λ1 |>| λ 2 |≥| λ3 |≥ ... ≥| λ n | （7.1.1）
相应的 n 个特征向量 xi (i = 1,2,..., n) 线性无关。上述假设表明， λ1 为非零单实根， x1 为实特征向量。
k →∞
k →∞
lim v ( k ) =
x1 max( x1 )
事实上，由式（7.1.5）知
v
(k )
=
Ak u ( 0 )
∏m
i =0
k
i
算法 7.1.1 实用幂法 (1) 输入： aij (i, j = 1,2, L n), ui (i = 1,2, L), ε ; (2) k = 1; m0 = max(ui );
7.1 幂法
幂法基本原理是：任取非零实向量 u
(0)
，做迭代
u ( k ) = Au ( k −1) = Ak u ( 0 ) (k = 1,2,...)
则
( 7 .1 . 2 )
λ1 = lim
这里 u j 表示向量 u
(k ) (k )
u (jk +1) u (jk )
k →∞

数值分析第三章小结

第三章矩阵特征值与特征向量的计算--------学习小结一、本章学习体会本章我们学习了矩阵特征值与特征向量的计算方法即幂法、反幂法、Jacobi方法和QR方法。

下边介绍一下四种方法各自的特点和适用范围。

幂法：主要用于计算矩阵按模最大的特征值及其相应的特征向量；反幂法：主要用于计算矩阵按模最小的特征值及其相应的特征向量；Jacobi法：用于求实对称矩阵的全部特征值和特征向量的方法；QR法：则适用于计算一般实矩阵的全部特征值，尤其适用于计算中小型实矩阵的全部特征值。

归结起来，这四种方法有一个共同的特点，即都是用了迭代的方法来求矩阵的特征值和特征向量。

还有利用用MATLAB自带的解法求解特征值和特征向量，其自带函数Eig即得到结果是虚数也可以算出，并且结果自动正交化。

二、本章知识梳理在工程技术中，计算矩阵的特征值和特征向量主要使用数值解法。

本章将阐述幂法、反幂法、Jacobi 方法、和QR 方法，并且只限于讨论实矩阵的情况。

3.1 幂法和反幂法（1）幂法幂法主要用于计算矩阵的按模为最大的特征值和相应的特征向量，其思想是迭代。

设n ⨯n 实矩阵A 具有n 个线性无关的特征向量,,...,,321n x x x x 其相应的特征值n λλλ...21，，满足如下不等式 n λλλλ≥≥≥> (321)其中i i i x Ax λ= ）。

（n i ,...2,1=现在要求出1λ和相应的特征向量。

任取一n 维非零向量0u ，从0u 出发，按照如下的递推公式 1-=k k Au u ），，（...21=k 因n 维向量组n x x x ,...,21线性无关，故对于向量0u ，必存在唯一的不全为零的数组n ααα,...,21，使得n n x x x u ααα...22110++=n k n k k k k k k x A x A x A u A u A Au u ααα+++=====--......22110221=⎥⎥⎦⎤⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛++⎪⎪⎭⎫ ⎝⎛+=+++n kn n k kn k n n k k x x x x x x 12122111222111......λλαλλααλλαλαλα 设01≠α。

矩阵的特征值和特征向量的计算

矩阵的特征值和特征向量的计算矩阵的特征值和特征向量是线性代数中比较重要的概念。

在机器学习、信号处理、图像处理等领域都有着广泛的应用。

本文将会介绍矩阵的特征值和特征向量的概念、意义以及计算方法。

一、特征值和特征向量的定义对于一个n阶方阵A，如果存在一个n维向量v和一个常数λ，使得下面的等式成立：Av=λv那么称λ为矩阵A的特征值，v为矩阵A的特征向量。

特征向量是非零向量，因为如果v为0向量，等式就无法成立。

另外，特征向量不唯一，如果v是A的特征向量，k是任意一个非零常数，那么kv也是A的特征向量。

但特征值是唯一的。

二、特征值和特征向量的意义矩阵的特征值和特征向量有着重要的物理和数学含义。

对于一个矩阵A，它的特征向量v和特征值λ描述的是矩阵A对向量v的作用和量变化。

当一个向量v与矩阵A相乘时，向量v的方向可能会发生变化，而特征向量v就是那些方向不变的向量，仅仅发生了缩放，这个缩放的倍数就是特征值λ。

也就是说，特征向量v在被矩阵A作用后仍保持了原来的方向，并且只发生了缩放。

从物理角度理解，矩阵的特征值和特征向量可以描述线性系统的固有特性。

在某些情况下，如机械振动、电路等自然界现象中，系统本身就带有某种特有的振动频率或固有响应。

而这些系统在一些特殊的情况下可以通过线性代数描述，正是因为它们具有特征值和特征向量。

三、特征值和特征向量的计算矩阵的特征值和特征向量可以通过求解特征方程来计算。

特征方程的形式为det(A-λI)=0，其中det(A-λI)表示A-λI的行列式，I是单位矩阵。

求解特征方程可以得到矩阵A的n个特征值λ1,λ2,…,λn。

接下来，针对每个特征值λi，都可以通过求解线性方程组(A-λiI)v=0来得到一个特征向量vi。

需要注意的是，一个矩阵的特征值和特征向量并不一定都能够求出来，只有在某些情况下才可以求出。

例如，对于一个非方阵，就不存在特征值和特征向量。

另外，如果矩阵的特征值出现重复，那么对应于这些特征值的特征向量可能无法确定，可以使用广义特征向量来处理。

矩阵特征值和特征向量的数值解

风险管理
在风险管理模型中，可以使用风险矩阵的特征值和特征向量来分析风险的分布和相关性，从而制定有效的风险管理策略。
感谢您的观看
THANKS
稳定性分析通常通过比较不同数值解法的计算结果，观察其误差随舍入精度的变化情况来进行。
稳定性好的算法能够在不同舍入精度下保持一致的计算结果，而稳定性差的算法则可能导致计算结果的较大偏差。
数值解法的收敛性分析
01
收敛性分析是评估数值解法求解特征值和特征向量问题的有效性的关键步骤。

02
收敛性分析主要关注算法是否能够收敛到正确的解，以及收敛
通过求解矩阵的特征值和特征向量，可以找到线性变换下的不变量，从而更好地理解和分析线性变换的性质和行为。
特征值和特征向量在矩阵的奇异值分解和QR分解等矩阵分解方法中也有着重要的应用，这些分解方法在许多科学计算和工程领域中都有广泛的应用。
在微分方程中的应用
01
矩阵特征值和特征向量在解决微分方程问题中也有着重要的应用。
速度的快慢。
收敛速度的快慢通常用收敛阶数来衡量，收敛阶数越高，收敛
03
速度越快。
数值解法的误差估计
01
误差估计是对数值解法计算结果的精度进行量化的重要手段。
02
误差估计通常通过比较数值解法的计算结果与精确解之间的差异来进行。
03
误差估计可以帮助我们了解算法的精度，从而在实际应用中选择合适的算法和舍入精度。
在研究热传导问题时，热传导矩阵的特征值和特征向量可以用来确定温度场的分布和变化。
在工程问题中的应用实例
结构分析
在结构分析中，结构的质量矩阵和刚度矩阵的特征值和特征向量可以用来确定结构的固有频率和振型，从而评估结构的稳定性和安全性。

现代数值分析07 矩阵特征值与特征向量计算(1)

3 1 2
3 2 7
的最靠近 p = -13 的特
Jacobi方法
Jacobi 方法是计算实对称矩阵全部特征值和特征向量的一种古典方法，它在理论和使用上至今还是有一定的意义。我们知道，在相似变换下，矩阵的特征值不变。Jacobi 方法的基本思想是选取一系列正交矩阵 J(pk,qk,θk),对矩阵 A 作正交相似变换，即
1 2 n 1 n
对应的特征向量为 x1, x2, ..., xn ；则 A-1 的特征值为： ,
1
1 1
A 可对角化
1
2
, ,
n1
,
1
n
对应的特征向量仍然为 x1, x2, ..., xn 。
计算A的按模最小特征值 n 计算A-1的按模最大特征值
2 1
的大小。
任取非零向量 v0，计算
v k 1 j vk j
是否近似于某个
常数，如果是，则停机；否则转第 1 步。
练习题求矩阵
12 A 6 6
6 16 2
6 2 16
按模最大的特征值。
要求自己编写Matlab 程序,进行计算并给出结果.
k 1
m a x ( v k 1 ) 1
k m a x 1 k m a x 1
1 x1 1 x1
i i i2 1
n
xi
i i xi i2 1
按模最大的特征值和
相应的特征向量。
要求自己编写Matlab 程序,进行计算并给出结果.
乘幂法的加速
乘幂法的收敛速度取决于 r

特征向量和特征值计算

特征向量和特征值计算
特征向量和特征值是线性代数中非常重要的概念，它们在多个领域都有广泛的应用。

本文将介绍如何计算特征向量和特征值。

首先，我们需要知道特征向量和特征值的定义。

对于一个n x n 的矩阵A，如果存在非零向量v和标量λ，使得Av=λv成立，则v
称为A的特征向量，λ称为A的特征值。

计算特征向量和特征值的步骤如下：
1. 求出A-λI的行列式det(A-λI)，其中I是n阶单位矩阵，det表示行列式。

2. 解出方程det(A-λI)=0，得到λ的值。

3. 将λ的值代入(A-λI)x=0，求出x的解，x就是对应的特征向量。

需要注意的是，特征向量是非零向量，所以解出的x值不能为0向量。

在实际计算中，通常使用特征值分解来求解特征向量和特征值。

特征值分解是将一个矩阵分解为特征向量和特征值的形式，可以用于矩阵对角化等问题。

具体步骤如下：
1. 计算矩阵A的特征值λ1, λ2, ..., λn。

2. 对于每个特征值λi，求出对应的特征向量xi。

3. 将所有特征向量组成一个矩阵X=[x1, x2, ..., xn]。

4. 将所有特征值组成对角矩阵Λ=diag(λ1, λ2, ..., λn)。

5. 则A可以表示为A=XΛX^-1。

特征向量和特征值的计算在机器学习、信号处理、量子力学等领域都有广泛的应用。

掌握其计算方法可以帮助我们更好地理解和应用相关的理论和算法。

第8章矩阵特征值与特征向量的计算

也可将上式改写成
Ax k = y k −1 m = max( x ) k k y = x / m , k = 1,2,3, ⋯ k k k
式(8.8)称为反幂法. 显然有 1 lim mk =
k →∞
(8.8)
λn
, lim y k = x n / max( x n )
k →∞
每一步求xk需要求解线性方程组, 可采用LU分解法求解.
i =2
n
其收敛速度由比值|λ2/λ1|来确定.
又由于
A k x0 mk = max( y k ) = max( Axk −1 ) = max( ) k −1 max( A x0 )
所以
lim mk = λ1
k →∞
= λ1
λ max[β 1ξ1 + ∑ β i ( λ1i ) k ξ i ] λ max[ β1ξ1 + ∑ β i ( λ1i ) k −1 ξ i ]
可取 λ1 ≈6.000837, ξ1 ≈(1,0.714316,-0.249895)T.
研究生学位课程数值分析
8.1.2 加速技术
由于
mk = max( x k ) = λ1 + o(
λ2 k λ1
)
所以,乘幂法收敛速度取决于比值|λ2/λ1|,当|λ2/λ1|≈1时,收敛是很慢的.
1.Aitken 加速方法 Aitken
可得
y k = Ax k −1 m = max( y ) k k x = y / m , k = 1,2,3, ⋯ k k k
A x0 xk = = max( A k x0 )
k
λ β1ξ1 + ∑ β i ( λ ) k ξ i

数值分析-第四章矩阵特征值与特征向量的计算2017秋

1 9 2.9996952.
12
幂法的加速—Aitken加速法
若 {ak} 线性收敛于a, 即
lim ak1 a C 0 k ak a 当 k 充分大时，有
ak1 a ak2 a ak a ak1 a

a

ak

(ak1 ak ak2 2ak1
涉及知识：乘幂法
22
上机作业
• 第121页第2题.
23
如果x(0)的选取恰恰使得t1=0, 幂法仍能进行. 因为计算过程中会有舍入误差, 迭代若干次后, 必然会产生一个向量x(k), 它在u1方向上的分量不为零, 这样以后的计算就满足所设条件.
因为 x(k) 1kt1u1, 计算过程中可能会出现上溢 (|1|>1)或下溢成为0 (|1|<1). 为避免出现这一情形, 实

Ak1 x(0)

在一定条件下, 当k充分大时:
1
x ( k 1) i xi(k )
相应的特征向量为: x(k1) 4
幂法的理论依据
n
对任意向量x(0), 有 x(0) tiui ,
i 1
x(k1) Ax(k ) Ak1 x(0)
设t1不为零.
n
)2 ak
13
把上式右端记为 aˆk
aˆk

ak

(ak1 ak )2 ak2 2ak1 ak
可以证明
lim k
aˆk

a,
klim
aˆk ak

a a

0.
即 aˆk a 比 ak a 快.
用 aˆk 逼近a, 这就是Aitken加速法.

矩阵的特征值与特征向量算法

矩阵的特征值与特征向量算法矩阵特征值和特征向量定义A为n阶矩阵，若数λ和n维非0列向量x满足Ax=λx，那么数λ称为A的特征值，x称为A的对应于特征值λ的特征向量。

式Ax=λx也可写成( A-λE)x=0，并且，λE-A，叫做A 的特征多项式。

当特征多项式等于0的时候，称为A的特征方程，特征方程是一个齐次线性方程组，求解特征值的过程其实就是求解特征方程的解。

依据普通线性代数中的概念，特征值和特征向量能够用传统的方法求得，可是实际项目中一般都是用数值分析的方法来计算。

这里介绍一下雅可比(Jacobi)迭代法求解特征值和特征向量。

雅可比(Jacobi)迭代法雅克比方法用于求实对称阵的所有特征值、特征向量。

Jacobi算法计算简单、稳定性好、精度高、求得的特征向量正交性好。

但当A为稀疏阵时，Givens旋转变换将破坏其稀疏性，且只能适用于实对称矩阵。

相关知识•矩阵A与相似矩阵 B = P A P-1的特征值相同。

•若矩阵Q满足QT Q = I,则称Q为正交矩阵。

显然Q-1 = QT,且正交阵的乘积仍为正交阵。

•若A为实对称矩阵，则存在正交阵Q，使Q A QT = diag(λ1,λ2,...,λn)，且QT 的列是相应的特征向量。

•实对称矩阵的特征值均为实数，且存在标准正交的特征向量系。

•Givens 旋转矩阵R(p,q,θ)是正交阵，其中Givens 旋转矩阵R原理：•Jacobi 方法用平面旋转对实对称矩阵 A 做一系列旋转相似变换，从而将A 约化为对角阵，进而求出特征值与特征向量。

•R A RT 与A元素之间的关系：为使R A RT 为对角矩阵，可选择θ为：当A为n阶实对称矩阵时，设A有非对角元，apq ≠ 0 ，设Givens 旋转矩阵R(p,q,θ)为：令C = R A RT,则有：若令C的非对角元素cpq = cqp = 0,则：C与A的元素满足下列关系：说明经旋转变换C = R A RT后，C的对角线元素平方和比A的对角线元素平方和增加了2apq2、而C的非对角元素平方和比A的非对角元素平方和减少了2apq2、如果不断的变换下去，则最后非对角元素可趋于0，即可通过一系列旋转变换，使A与一对角阵相似。

数值分析——矩阵特征值问题计算

17
vk 1k a1x1
即为矩阵 A 的对应特征值 1 的近似特征向量。
且
vk 1 Avk 1k 1a1x1 1vk
用 (vk)i 表示 vk 的第 i 个分量,则当k充分大时,有
vk1 i
vk
i
1
即为主特征值的近似值。
18
定理设 A Rnn 有 n 个线性无关的特征向量，
主特征值 1 满足
a11k x1 a2k2 x2 ankn xn
1k
a1
x1
a2
2 1
k
x2
an
n 1k Fra bibliotek xn 1k a1x1 k
16
其中
k
a2
2 1
k
x2
an
n 1
k
xn
由假设条件从而
j 1
1 j 2, ,n, 所以
lim
k
vk
1k
a1x1
lim
k
k
0
所以当k充分大时，有
vk 1k a1x1
9
0.5, 1, 0.8611
3 5.7222, 11.4444, 8.361 4 5.4621, 10.9223, 8.2306 5 5.5075, 11.0142, 8.2576
11.4444 10.9223 11.0142
0.5, 1, 0.7360 0.5, 1, 0.7536 0.5, 1, 0.7494
n
( Ax, x) (x, x)
1
（2）
n
min x0
( Ax, x) (x, x)
（3）
1
max x0
( Ax, x) (x, x)