数值分析第五章 矩阵特征值与特征向量的计算
特征值与特征向量的计算方法
特征值与特征向量的计算方法特征值与特征向量是矩阵理论中的重要概念，用于解决矩阵特征与变换特性的相关问题。
在本文中，将介绍特征值与特征向量的定义和计算方法，以及它们在实际问题中的应用。
一、特征值与特征向量的定义在矩阵理论中，对于一个n阶方阵A，如果存在一个非零向量x，使得Ax=kx（k为标量），那么k称为矩阵A的特征值，x称为对应于特征值k的特征向量。
特征向量可以理解为在矩阵变换下保持方向不变的向量，而特征值则表示特征向量在变换中的伸缩比例。
二、要计算特征值和特征向量，可以使用以下步骤：1. 首先，由于特征值和特征向量的定义基于方阵，所以需要确保矩阵A是方阵，即行数等于列数。
2. 接下来，根据特征值和特征向量的定义方程Ax=kx，将其改写为(A-kI)x=0（I为单位矩阵）。
3. 为了求解此方程组的非零解，需要求出(A-kI)的零空间（核）。
4. 将(A-kI)的零空间表示为Ax=0的齐次线性方程组，采用高斯消元法或其它线性方程组求解方法，求得方程的基础解系，即特征向量。
5. 特征向量已找到，接下来通过将每个特征向量代入原方程式Ax=kx中，计算出对应的特征值。
值得注意的是，特征值是一个多重属性，即一个特征值可能对应多个线性无关的特征向量。
此外，方阵A的特征值计算方法存在多种，如幂迭代法、QR迭代法等。
三、特征值与特征向量的应用特征值与特征向量在物理、工程、经济等领域具有广泛的应用。
1. 物理学中，特征值与特征向量可用于解析力学、量子力学等领域中的问题，如研究振动系统的固有频率、粒子的角动量等。
2. 工程学中，特征值与特征向量可用于电力系统的稳定性分析、机械系统的振动模态分析等。
3. 经济学中，特征值与特征向量可用于描述经济模型中的平衡点、稳定性等重要特征。
此外，特征值与特征向量在图像识别、数据降维、网络分析等领域也有重要的应用。
总结：特征值和特征向量在矩阵理论中有着重要的地位和应用价值。
通过计算特征值和特征向量，可以揭示矩阵在变换中的性质和特点，并应用于各个学科领域，为问题求解提供了有效的工具和方法。
第五章矩阵的特征值与特征向量
λ1λ2 ⋯λn = A；
及 题 有 以 依 意 f (λ) = (λ − 4)(λ −1)(λ + 2).
α = (1 , 1) 为A的属于特征值5的特征向量.
T
李正兴2011-11-10 3
说明： 说明：
(1)如果α1，α 2都是A的属于特征值λ0的特征向量，则 k1α1 + k2α 2 (k1α1 + k2α 2 ≠ 0)也是A的属于特征值λ0的特征向量.
特别地
(2)如果α 是A的属于特征值λ0的特征向量， 则kα (k ≠ 0)也是A的属于特征值λ0的特征向量. (3)如果α1，α 2都是A的属于特征值λ0的特征向量，则
所以
2E − A = (2 − λ1 )(2 − λ2 )⋯ (2 − λn ).
例7 设A，B均为n阶矩阵，且B = P −1 AP P为n阶可逆矩阵）， （
证明: tr ( B ) = tr ( P −1 AP ) = tr ( APP −1 ) = tr ( A ) .
李正兴2011-11-10 19
李正兴2011-11-10 18
4. tr ( AB ) = tr ( BA )
例6 已知n阶矩阵 的n个特征值是 λ1 , λ2 ,⋯ , λn , 已知 阶矩阵A的 个特征值是 阶矩阵 求 2E − A . 依题意， 解 依题意，2E − A 的特征值是
2 − λ1 , 2 − λ2 ,⋯ , 2 − λn ,
第五章 矩阵的特征值与特征向量 §5.1 矩阵的特征值与特征向量 §5.2 相似矩阵与矩阵的对角化 §5.3 实对称矩阵的对角化
矩阵特征值与特征向量计算
矩阵特征值与特征向量计算在数学中，矩阵是一种非常基础而且重要的概念，它可以被看做是一种线性变换的表示。
在矩阵中，特征值和特征向量是两个非常重要的概念，它们在运用矩阵进行计算、测量和定量分析时扮演着至关重要的角色。
一、矩阵特征值的计算方法特征值是一个矩阵的固有属性，它表示在进行线性变换时，各个方向上对应的比例因子，具有很重要的几何意义。
计算一个矩阵的特征值需要使用到线性代数的基础知识和运算。
对于一个n阶方阵A，如果存在一个非零向量x和一个标量λ，使得Ax=λx，则λ是矩阵A的一个特征值，而x是对应的特征向量。
在实际计算中，我们首先需要求解方程det(A-λI)=0，其中I是指n阶单位矩阵。
这个方程的解即为矩阵A的特征值，它们可以是实数或复数。
当然，在计算特征值时，使用一些优化的方法可以更快地得出结果，例如使用特征值分析法或雅可比方法。
二、矩阵特征向量的计算方法在获得了矩阵的特征值之后，我们可以通过简单的代数运算来计算它们对应的特征向量。
设λ为矩阵A的一个特征值，x为一个对应的特征向量，我们有以下等式：(A-λI)x=0这可以被看做是一个齐次线性方程组，将它转化成矩阵形式，我们得到以下方程：(A-λI)X=0其中X=[x1,x2,...,xn]为特征向量的矩阵形式。
对于特征向量矩阵X，我们需要求解出它的非零解。
这需要使用到线性代数的基本技巧，例如高斯消元法或LU分解等。
三、矩阵特征值和特征向量的应用矩阵特征值和特征向量的应用非常广泛，从计算机科学到物理学、化学、经济学、金融学等各个领域都有它们的应用。
以下是几个主要的应用领域：1. 机器学习和人工智能在机器学习和人工智能中，特征值和特征向量经常用于降维和数据分析。
通过分析一个数据矩阵的特征值和特征向量，我们可以找到它们对应的主要特征，从而对大型数据进行有效的分析和处理。
2. 物理学和化学在物理学和化学中，特征值和特征向量可以用于计算量子力学、分析分子结构、电子轨道等问题。
矩阵特征值与特征向量的计算方法
矩阵特征值与特征向量的计算方法矩阵是一个广泛应用于线性代数、微积分和物理学等领域的数学对象。
在许多问题中，矩阵和线性变换起着重要作用，并且特征值与特征向量是矩阵理论中的两个核心概念。
本文将介绍矩阵特征值与特征向量的定义、性质以及计算方法。
一、特征值与特征向量的定义给定一个n阶矩阵A，如果存在一个非零向量x，使得A与x的线性组合仍然是x的倍数，即有Ax = λx其中λ为常数，称λ为A的特征值，x为对应于λ的特征向量。
从几何意义上理解，特征向量是不被矩阵变换影响方向，只被影响长度的向量。
特征值则是描述了矩阵变换对于特定方向上的伸缩倍数。
二、特征值与特征向量的性质1. 特征向量构成的向量空间没有零向量。
证明：设x为A的特征向量，有Ax=λx，则A(cx) =cAx=cλx=λ(cx),即A的任意常数倍(cx)仍是x的倍数，因此cx也是A的特征向量。
特别地，对于λ≠0时，x/λ也是A的特征向量。
2. A的特征值的个数不超过n个。
证明：考虑特征值λ1, λ2,…,λt，对应于各自的特征向量x1,x2,…,xt。
利用向量线性无关性可得，至少存在一个向量y不属于x1,x2,…,xt的张成空间内，此时Ay不能被表示成λ1x1,λ2x2,…,λtxt的线性组合，因此Ay与y方向没有重合部分，由此可得λ1, λ2,…,λt最多就是n个。
3. 如果特征向量x1,x2,…,xt彼此不共线，则它们就可以作为Rn空间的一组基。
证明：设x1,x2,…,xt是不共线的特征向量，考虑它们张成的向量空间V，在此空间中，A的作用就是对向量做伸缩变换，且Λ(xj) = λj。
对于每个向量y ∈ V，y可以表示成如下形式：y = c1x1 + c2x2 + ··· + ctxt由于x1,x2,…,xt构成V的基，因此c1,c2,…,ct唯一确定了向量y。
因此，对于任意的向量y，可以得到:Ay = A(c1x1 + c2x2 + ··· + ctxt)= c1Ax1 + c2Ax2 + ··· + ctAxt= λ1c1x1 + λ2c2x2 + ··· + λtctxt由于{x1,x2,…,xt}是V的一组基，c1,c2,…,ct是唯一确定的，因此Ay也被唯一确定了。
矩阵的特征值与特征向量的简易求法
矩阵的特征值与特征向量的简易求法特征值与特征向量对于矩阵的性质和变换有着重要的意义。
矩阵的特征值可以帮助我们判断矩阵的相似性、可逆性以及矩阵的对角化等；而特征向量可以帮助我们理解矩阵的线性变换、寻找矩阵的基矢量等。
求解矩阵的特征值与特征向量可以采用多种方法。
下面介绍两种常见的简易求法：特征多项式法和幂迭代法。
特征多项式法是求解矩阵特征值与特征向量的一种常见方法。
其步骤如下：步骤1：对于n阶方阵A，求解其特征多项式，即特征方程det(A-λI)=0。
其中，I为单位矩阵，λ为未知数。
步骤2：将特征多项式化简，得到一个关于λ的方程，如λ^n+c1λ^(n-1)+c2λ^(n-2)+...+cn=0。
步骤3：解这个n次方程，得到n个特征值λ1,λ2,...,λn。
步骤4：将每个特征值λi带入原方程(A-λI)X=0，求解对应的特征向量。
特征多项式法适用于任意阶数的方阵，但是对于高阶矩阵，其计算过程可能比较复杂，需要借助数值计算工具。
幂迭代法是一种迭代求解特征值与特征向量的方法，适用于对于方阵的特征值为实数且相近的情况。
其步骤如下：步骤1：选取一个初始向量X(0)，通常是一个n维非零向量。
步骤2：迭代计算：X(k+1)=A*X(k)，其中k为迭代次数，A为待求特征值与特征向量的方阵。
步骤3：计算迭代步骤2中得到的向量序列X(k)的模长，即，X(k)。
步骤4：判断，X(k)-X(k-1)，是否满足预定的精度要求，如果满足，则作为矩阵A的近似特征向量；否则，返回步骤2继续进行迭代。
步骤5：将步骤4得到的近似特征向量作为初始向量继续迭代，直至满足精度要求。
幂迭代法的优点是求解简单、易于操作，但由于其迭代过程，只能得到一个特征值与特征向量的近似解，且只适用于特征值为实数的情况。
在实际应用中，根据具体问题的要求，可以选择适合的方法来求解矩阵的特征值与特征向量。
除了特征多项式法和幂迭代法，还有QR分解法、雅可比迭代法等其他方法。
矩阵特征值与特征向量的求法
矩阵特征值与特征向量的求法一、矩阵特征值与特征向量的定义矩阵特征值（eigenvalue）是指一个矩阵在某个非零向量上的线性变换结果等于该向量的常数倍，这个常数就是该矩阵的特征值。
而对应于每个特征值，都有一个非零向量与之对应，这个向量就是该矩阵的特征向量（eigenvector）。
二、求解矩阵特征值与特征向量的方法1. 特征多项式法通过求解矩阵A减去λI（其中λ为待求解的特征值，I为单位矩阵）的行列式det(A-λI)=0来求解其特征值。
然后将每个特征值代入到(A-λI)x=0中，即可求得对应的特征向量x。
2. 幂法幂法是一种迭代方法，通过不断地将A作用于一个初始向量x上，并将结果归一化，最终得到收敛到最大（或最小）特征值所对应的特征向量。
具体步骤如下：(1) 选取任意一个非零初始向量x；(2) 将Ax除以x中最大元素得到新的向量y=A*x/max(x)；(3) 将y归一化得到新的向量x=y/||y||；(4) 重复步骤2-3，直到收敛。
3. QR分解法QR分解是将矩阵A分解为Q和R两个矩阵的乘积，其中Q是正交矩阵（即Q^T*Q=I），R是上三角矩阵。
通过不断地对A进行QR分解，并将得到的Q和R相乘，最终得到一个上三角矩阵T。
T的对角线元素就是A的特征值，而对应于每个特征值，都可以通过反推出来QR分解中的Q所对应的特征向量。
4. Jacobi方法Jacobi方法也是一种迭代方法，通过不断地施加相似变换将A转化为对角矩阵D。
具体步骤如下：(1) 选取任意一个非零初始矩阵B=A；(2) 找到B中绝对值最大的非对角元素b(i,j)，记其位置为(i,j)；(3) 构造Givens旋转矩阵G(i,j,k)，使其作用于B上可以消去b(i,j)，即B=G^T*B*G；(4) 重复步骤2-3，直到所有非对角元素均趋近于0。
三、总结以上介绍了求解矩阵特征值与特征向量的四种方法：特征多项式法、幂法、QR分解法和Jacobi方法。
第五章矩阵的特征值与特征向量的计算doc
在自然科学和工程设计中的许多问题，如电磁振荡、桥梁振动、机械振动等，常归结为求矩阵的特征值和特征向量.求矩阵的特征值和特征向量的问题是代数计算中的重要课题.本章着重介绍直接计算矩阵的特征值和特征向量的MATLAB 程序、间接计算矩阵的特征值和特征向量的幂法、反幂法、雅可比方法、豪斯霍尔德方法和QR 方法及其它们的MATLAB 计算程序.最后我们还讨论广义特征值问题.直接计算特征值和特征向量的MATLAB 程序计算特征值和特征向量的MATLAB 程序从以上的讨论可以看到，有许多问题归结为求矩阵的特征值和特征向量，而用手工计算高阶矩阵的特征值与特征向量的难度较大，但是，计算机软件MATLAB 提供了直接计算特征值与特征向量的MATLAB 函数 （见表5–1），下面介绍这些函数的使用方法.命 令功 能 b = eig(A)输入方阵A ，运行后输出b 为由方阵A 的全部特征值构成的列向量 [V,D] = eig (A) 输入对称矩阵A ，运行后输出D 为由A 的全部特征值构成的对角矩阵，V 的各列为对应于特征值的特征向量构成的矩阵，使得AV = DV[V,D] = eig (A,'nobalance')输入方阵A ，运行后输出D 为由A 的全部特征值构成的对角矩阵，V 的各列为对应于特征值的特征向量构成的矩阵，使得AV = DV ；如果A 是对称矩阵，则输出的结果与程序 [V,D] = eig (A)的运行结果相同幂法及其MATLAB 程序幂法是求实矩阵A 的主特征值（即实矩阵A 按模最大的特征值）及其对应的特征向量的一种迭代方法.幂法的MATLAB 程序设n 阶实矩阵A 的n 个特征值为n λλλ,,,21 ，且满足021>≥≥≥n λλλ ，A 的主特征值1λ对应的特征向量为1X ，则我们可以用下面的MATLAB 程序计算1λ和1X 的近似用幂法计算矩阵A 的主特征值和对应的特征向量的MATLAB 主程序输入的量：n 阶实矩阵A 、n 维初始实向量V 0、计算要求的精度jd 、迭代的最大次数max1；输出的量：迭代的次数k 、A 的主特征值1λ的近似值lambda 、1λ对应的特征向量1X 的近似向量V k 、相邻两次迭代的误差W c .如果迭代次数已经达到最大的迭代次数max1，则给出提示的相关信息.第五章 矩阵的特征值与特征向量的计算根据迭代公式（），现提供用幂法计算矩阵A 的主特征值和对应的特征向量的MATLAB 主程序如下：function [k,lambda,Vk,Wc]=mifa(A,V0,jd,max1)lambda=0;k=1;Wc =1; ,jd=jd*;state=1; V=V0;while ((k<=max1)&(state==1))Vk=A*V; [m j]=max(abs(Vk)); mk=m;tzw=abs(lambda-mk); Vk=(1/mk)*Vk;Txw=norm(V-Vk); Wc=max(Txw,tzw); V=Vk;lambda=mk;state=0;if (Wc>jd)state=1;endk=k+1;Wc=Wc;endif (Wc<=jd)disp('请注意：迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下：')elsedisp('请注意：迭代次数k 已经达到最大迭代次数max1,主特征值的迭代值lambda,主特征向量的迭代向量Vk,相邻两次迭代的误差Wc 如下：')endVk=V;k=k-1;Wc;例 用幂法计算下列矩阵的主特征值和对应的特征向量的近似向量，精度510-=ε.并把（1）和（2）输出的结果与例中的结果进行比较.(1)⎪⎪⎭⎫ ⎝⎛-=4211A ； (2)⎪⎪⎪⎭⎫ ⎝⎛=633312321B ；(3)⎪⎪⎪⎭⎫ ⎝⎛--=1124111221C ；（4）⎪⎪⎪⎭⎫ ⎝⎛---=20101350144D . 解 （1）输入MATLAB 程序>> A=[1 -1;2 4];V0=[1,1]';[k,lambda,Vk,Wc]=mifa(A,V0,,100),[V,D] = eig (A), Dzd=max(diag(D)), wuD= abs(Dzd- lambda), wuV=V(:,2)./Vk,运行后屏幕显示结果请注意：迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下：k = lambda = Wc =33Vk = V = wuV =Dzd = wuD =3 由输出结果可看出，迭代33次，相邻两次迭代的误差W c ≈ 19e-007，矩阵A 的主特征值的近似值≈ 00和对应的特征向量的近似向量V k ≈（ 00， 00T), lambda 与例中A 的最大特征值32=λ近似相等，绝对误差约为 37e-006，V k 与特征向量X =T22k T )1,21(- )0(2≠k 的第1个分量的绝对误差约等于0，第2个分量的绝对值相同.由wuV 可以看出，2λ的特征向量V (:,2) 与Vk 的对应分量的比值近似相等.因此，用程序计算的结果达到预先给定的精度510-=ε.(2) 输入MATLAB 程序>>B=[1 2 3;2 1 3;3 3 6]; V0=[1,1,1]';[k,lambda,Vk,Wc]=mifa(B,V0,,100), [V,D] = eig (B),Dzd=max(diag(D)), wuD= abs(Dzd- lambda), wuV=V(:,3)./Vk,运行后屏幕显示结果请注意：迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下：k = lambda = Wc = Dzd = wuD =3 9 0 9 0Vk = wuV =V =由输出结果可看出，迭代3次，相邻两次迭代的误差W c =0，实对称矩阵B 的主特征值的近似值lambda =9和对应的特征向量的近似向量V k =（ 00， 00， 00T ),lambda 与例中B 的最大特征值93=λ相同，V k 与特征向量X =T33k T )2,1,1( )0(3≠k 的对应分量成比例.从wuV 的每个分量的值也可以看出，3λ的特征向量V (:,3) 与Vk 的对应分量的比值相等.因此，用程序计算的结果达到预先给定的精度510-.此例说明，幂法对实对称矩阵B 的迭代速度快且计算结果精度高，(3) 输入MATLAB 程序>> C=[1 2 2;1 -1 1;4 -12 1];V0=[1,1,1]';[k,lambda,Vk,Wc]=mifa(C,V0,,100), [V,D] = eig (C),Dzd=max(diag(D)), wuD= abs(Dzd- lambda), Vzd=V(:,1),wuV=V(:,1)./Vk,运行后屏幕显示请注意：迭代次数k 已经达到最大迭代次数max1,主特征值的迭代值lambda,主特征向量的迭代向量Vk,相邻两次迭代的误差Wc 如下：k = lambda = Wc =100Dzd = wuD =Vk= Vzd = wuV =由输出结果可见，迭代次数k 已经达到最大迭代次数max 1=100，并且lambda 的相邻两次迭代的误差≈ 58>2,由wuV 可以看出，lambda 的特征向量Vk 与真值Dzd 的特征向量Vzd 对应分量的比值相差较大，所以迭代序列发散.实际上，实数矩阵C 的特征值的近似值为i ,i ,010*********.000321=-==λλλ ，并且对应的特征向量的近似向量分别为X T1=1k （，，）T， X =T 22k （, ）T ， X =T33k ( , + ,+ T 0,0(21≠≠k k , 03≠k 是常数）.此例说明，当n 阶实矩阵有复数特征值时，不宜用幂法计算它的主特征值1λ对应的特征向量1X .（4）输入MATLAB 程序>> D=[-4 14 0;-5 13 0;-1 0 2]; V0=[1,1,1]';[k,lambda,Vk,Wc]=mifa(D,V0,,100), [V,Dt] = eig (D),Dtzd=max(diag(Dt)), wuDt= abs(Dtzd- lambda), Vzd=V(:,2),wuV=V(:,2)./Vk,运行后屏幕显示结果请注意：迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下：k = lambda = Wc =19Dtzd = wuDt =Vk = Vzd = wuV =由输出结果可见，迭代19次，相邻两次迭代的误差≈ 52e-006，矩阵D 的主特征值的近似值≈ 01和对应的特征向量的近似向量为V k ≈（ 40， 29， 00T ).用eig(A)计算，lambda 与对应的特征值的真值Dtzd 的绝对误差wuDt 510-<，二者的特征向量V (:,2) 与Vk 的对应分量的比值近似相等（请对比wuV 的每个分量的值）.因此，用程序计算的结果达到预先给定的精度510-=ε.由例的计算结果可见，用幂法计算实对称矩阵B 的主特征值对应的特征向量时，得到的迭代序列的收敛速度最快且计算结果精度也最高；非实对称矩阵对应的迭代序列的敛散性不定，有时发散（例如矩阵C ）,有时收敛（例如矩阵A 和D ），且收敛速度较慢，比较（1）和（4）的计算结果可知，用幂法得到的迭代序列的收敛速度与矩阵的阶数无关；当实矩阵有复数特征值时，不宜用幂法计算它的主特征值对应的特征向量.反幂法和位移反幂法及其MATLAB 程序反幂法是求非奇异矩阵A 的按模最小特征值及其对应的特征向量的一种迭代方法.原点位移反幂法的MATLAB 程序设n 阶实矩阵A 的n 个特征值n λλλ,,,21 都不相同，且满足n j n n λλλλ~~-<- )(n j <，且A 的按模最小特征值n λ对应的特征向量为n X ，则我们可以根据迭代公式（）和（）编写两种MATLAB 程序分别计算n λ和n X 的近似值和近似向量.（一） 原点位移反幂法的MATLAB 主程序1根据原点位移反幂法的迭代公式（），现提供用原点位移反幂法计算矩阵A 的按模最用原点位移反幂法计算矩阵A 的特征值和对应的特征向量的MATLAB 主程序1输入的量：n 阶实矩阵A 、n 维初始实向量V 0、特征值的近似值jlamb 、计算的精度jd 、迭代的最大次数max1；输出的量：迭代的次数k 、A 的特征值n λ的近似值lambdan 、与n λ对应的特征向量n X 的近似向量V k 、相邻两次迭代的误差Wc .如果迭代次数已经达到最大的迭代次数max1，则给出提示的相关信息.[n,n]=size(A); A1=A-jlamb*eye(n); jd= jd*;RA1=det(A1);if RA1==0disp('请注意：因为A-aE 的n 阶行列式hl 等于零，所以A-aE 不能进行Lu 分解.')returnendlambda=0;if RA1~=0for p=1:nh(p)=det(A1(1:p, 1:p));endhl=h(1:n);for i=1:nif h(1,i)==0disp('请注意：因为A-aE 的r 阶主子式等于零，所以A-aE 不能进行Lu 分解.')returnendendif h(1,i)~=0disp('请注意：因为A-aE 的各阶主子式都不等于零，所以A-aE 能进行Lu 分解.')k=1;Wc =1;state=1; Vk=V0;while ((k<=max1)&(state==1))[L U]=lu(A1); Yk=L\Vk;Vk=U\Yk; [m j]=max(abs(Vk));mk=m;Vk1=Vk/mk; Yk1=L\Vk1;Vk1=U\Yk1;[m j]=max(abs(Vk1));mk1=m;Vk2=(1/mk1)*Vk1;tzw1=abs((mk-mk1)/mk1);tzw2=abs(mk1-mk);Txw1=norm(Vk)-norm(Vk1);Txw2=(norm(Vk)-norm(Vk1))/norm(Vk1);Txw=min(Txw1,Txw2); tzw=min(tzw1,tzw2); Vk=Vk2;mk=mk1; Wc=max(Txw,tzw); Vk=Vk2;mk=mk1;state=0;if (Wc>jd)state=1;endk=k+1;%Vk=Vk2,mk=mk1,endif (Wc<=jd)disp('A-aE 的秩R(A-aE)和各阶顺序主子式值hl 、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下：')elsedisp('A-aE 的秩R(A-aE)和各阶顺序主子式值hl 、迭代次数k 已经达到最大迭代次数max1,按模最小特征值的迭代值lambda,特征向量的迭代向量Vk,相邻两次迭代的误差Wc 如下：')endhl,RA1endend[V,D]=eig(A,'nobalance'),Vk;k=k-1;Wc;lambdan=jlamb+1/mk1;例 用原点位移反幂法的迭代公式（），根据给定的下列矩阵的特征值n λ的初始值n λ~，计算与n λ对应的特征向量n X 的近似向量,精确到 1.（1）⎪⎪⎪⎭⎫ ⎝⎛----210242011,2.0~2=λ；（2）⎪⎪⎭⎫ ⎝⎛-4211，001.2~2=λ；（3）⎪⎪⎪⎭⎫ ⎝⎛--3315358215211，8.26~3=λ. 解 （1）输入MATLAB 程序>> A=[1 -1 0;-2 4 -2;0 -1 2];V0=[1,1,1]';[k,lambda,Vk,Wc]=ydwyfmf(A,V0,,,10000)运行后屏幕显示结果请注意：因为A-aE 的各阶主子式都不等于零，所以A-aE 能进行Lu 分解.A-aE 的秩R(A-aE)和各阶顺序主子式值hl 、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下：k = lambda = Wc = hl =3Vk = V = D =0 00 00 0（2）输入MATLAB 程序>> A=[1 -1;2 4];V0=[20,1]';[k,lambda,Vk,Wc]=ydwyfmf(A,V0,,,100)运行后屏幕显示结果请注意：因为A-aE 的各阶主子式都不等于零，所以A-aE 能进行Lu 分解.A-aE 的秩R(A-aE)和各阶顺序主子式值hl 、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下：k = lambda = Wc = hl =2Vk = V = D =2 00 3（3）输入MATLAB 程序>> A=[-11 2 15;2 58 3;15 3 -3];V0=[1,1,-1]';[k,lambdan,Vk,Wc]=ydwyfmf(A,V0,, ,100)运行后屏幕显示结果请注意：因为A-aE 的各阶主子式都不等于零，所以A-aE 能进行Lu 分解.A-aE 的秩R(A-aE)和各阶顺序主子式值hl 、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下：k = lambdan= Wc = hl =2Vk = V = D =0 00 00 0例 用原点位移反幂法的迭代公式（），计算⎪⎪⎪⎭⎫ ⎝⎛-----=1026471725110A 的分别对应于特征值 1.001~11=≈λλ，.001 2~22=≈λλ， 001.4~33=≈λλ的特征向量1X ，2X ，3X 的近似向量，相邻迭代误差为.将计算结果与精确特征向量比较，其中=1λ 999 999 999 97…，=2λ 999 999 999 99…，=3λ 000 000 000 02…，分别对应特征向量为 =1X ( 000 000 000 00， 000 000 000 00， 000 000 000 00T )，=2X ( 999 999 999 99， 000 000 000 00， 000 000 000 00T )，=3X ( 000 000 000 00， 000 000 000 00， 000 000 000 00T ).解 （1）计算特征值 1.001~11=≈λλ对应的特征向量1X 的近似向量.输入MATLAB程序>> A=[0 11 -5;-2 17 -7;-4 26 -10];V0=[1,1,1]';[k,lambda,Vk,Wc]= ydwyfmf(A,V0,, ,100),[V,D]=eig(A);Dzd=min(diag(D)), wuD= abs(Dzd- lambda),VD=V(:,1),wuV=V(:,1)./Vk,运行后屏幕显示结果请注意：因为A-aE 的各阶主子式都不等于零，所以A-aE 能进行Lu 分解.A-aE 的秩R(A-aE)和各阶顺序主子式值hl 、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下：hl =k = lambda = RA1 =5Vk = VD = wuV =Wc = Dzd = wuD =从输出的结果可见，迭代5次，特征向量1X 的近似向量1~X 的相邻两次迭代的误差≈e-009，由wuV 可以看出，1~X = Vk 与VD 的对应分量的比值相等.特征值1λ的近似值lambda ≈与初始值=1~λ的绝对误差为，而与 1λ的绝对误差为，其中 =1X T )000000000001.000 , 000000000000.500- , 000000000000.500( -, =1~X T )000000000001.000 , 000000000000.500- , 000000000000.500(-. （2）计算特征值.001 2~22=≈λλ对应特征向量2X 的近似向量.输入MATLAB 程序>> A=[0 11 -5;-2 17 -7;-4 26 -10];V0=[1,1,1]';[k,lambda,Vk,Wc]=ydwyfmf(A,V0,, ,100) ,[V,D]=eig(A); WD=lambda-D(2,2),VD=V(:,2),wuV=V(:,2)./Vk,运行后屏幕显示结果请注意：因为A-aE 的各阶主子式都不等于零，所以A-aE 能进行Lu 分解.A-aE 的秩R(A-aE)和各阶顺序主子式值hl 、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下：hl =k = Wc = lambda = WD =2Vk = VD = wuV =从输出的结果可见，迭代2次，特征向量2X 的近似向量2~X 的相邻两次迭代的误差≈，2~X 与2X 的对应分量的比值近似相等.特征值2λ的近似值≈与初始值=2~λ的绝对误差约为，而lambda 与2λ的绝对误差约为，其中 =2~X T )00000000000000.1,99999999999499.0,99999999999249.0(---, =2X T ) 000000000001.000- ,000000000000.500- ,99999999999-0.249( . （3）计算特征值 001.4~33=≈λλ对应特征向量3X 的近似向量.输入MATLAB 程序>> A=[0 11 -5;-2 17 -7;-4 26 -10];V0=[1,1,1]';[k,lambda,Vk,Wc]=ydwyfmf(A,V0,, ,100)[V,D]=eig(A); WD=lambda-max(diag(D)),VD=V(:,3),wuV=V(:,3)./Vk,运行后屏幕显示结果请注意：因为A-aE 的各阶主子式都不等于零，所以A-aE 能进行Lu 分解.A-aE 的秩R(A-aE)和各阶顺序主子式值hl 、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下：hl =k = lambda = Wc = WD =2Vk = VD = wuV =从输出的结果可见，迭代2次，特征向量3X 的近似向量3~X 的相邻两次迭代的误差≈，3~X 与3X 的对应分量的比值近似相等.特征值3λ的近似值 4.001~4.0022=≈λ与初始值lambda 的绝对误差近似为001.0，而lambda 与3λ的绝对误差约为，其中 =3X ( 000 000 000 00， 000 000 000 00， 000 000 000 00T ),=3~X T )000000000001.000 ,100000000000.600 ,10000000000.400(.综上所述，用原点位移反幂法的迭代公式（）求矩阵的全部特征值分别对应特征向量的收敛速度快，且精确度高.但是求矩阵的全部特征值的结果不理想.（二）、原点位移反幂法的MATLAB 主程序2根据迭代公式（），我编写用原点位移反幂法计算矩阵A 的按模最小特征值和对应的用原点位移反幂法计算矩阵A 的特征值和对应的特征向量的MATLAB 主程序2输入的量：n 阶实矩阵A 、n 维初始实向量V 0、特征值的近似值jlambn 、计算的精度jd 、迭代的最大次数max1.输出的量：迭代的次数k 、A 的特征值n λ的近似值lambdan 、n λ对应的特征向量n X 的近似向量V k 、相邻两次迭代的误差Wc .如果迭代次数已经达到最大的迭代次数max1，则给出提示的相关信息.[n,n]=size(A); jd= jd*;A1=A-jlamb*eye(n);nA1=inv(A1);lambda1=0;k=1;Wc =1;state=1; U=V0;while ((k<=max1)&(state==1))Vk=A1\U; [m j]=max(abs(Vk)); mk=m; Vk=(1/mk)*Vk; Vk1=A1\Vk;[m1 j]=max(abs(Vk1)); mk1=m1,Vk1=(1/mk1)*Vk1;U=Vk1,Txw=(norm(Vk1)-norm(Vk))/norm(Vk1);tzw=abs((lambda1-mk1)/mk1);Wc=max(Txw,tzw); lambda1=mk1;state=0;if (Wc>jd)state=1;endk=k+1;endif (Wc<=jd)disp('请注意迭代次数k,特征值的近似值lambda,对应的特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下：')elsedisp('请注意迭代次数k 已经达到最大迭代次数max1, 特征值的近似值lambda,对应的特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下：')end[V,D] =eig(A,'nobalance'), Vk=U;k=k-1;Wc;lambdan=jlamb+1/mk;例 用原点位移反幂法的迭代公式（），计算例题，并且将这两个例题的计算结果进行比较.再用两种原点位移反幂法的MATLAB 主程序，求979999999990.999~1=λ对应的特征向量.解 （1）计算特征值 1.001~11=≈λλ对应特征向量1X 的近似向量.输入MATLAB 程序>> A=[0 11 -5;-2 17 -7;-4 26 -10];V0=[1,1,1]';[k,lambda,Vk,Wc]=wfmifa1(A,V0,,,100)运行后屏幕显示结果请注意迭代次数k,特征值的近似值lambda,对应的特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下：k = lambda = Wc =5Vk’ =同理可得，另外与两个特征值对应的特征向量，将计算结果列表5–3.表5–3 例题和例题计算结果比较由（）式计算结果 由（）式计算结果 迭代次数k 1X 和其近似向量V k1~X 的相邻迭代误差W c 1λ和其近似值lambda k = 51X =, , T V k =, , T W c = 344 154 436 924e-006 lambda = 000 000 001 38 979999999990.9991 =λk = 51X = , , T V k =, , TW c = 794 763 695 562e-009lambda = 000 000 000 00 979999999990.9991 =λ 迭代次数k 2X 和其近似向量V k 2~X 的相邻迭代误差W c 2λ和其近似值lambda k = 2 2X =, , T V k =, , TW c = 343 455 491 521e-004 lambda = 999 999 687 03 999999999991.9992=λk =22X =, , TV k =, , T W c = 258 787 820 493e-007lambda = 000 000 000 16999999999991.9992=λ 迭代次数k 3X 和其近似向量V k 3~X 的相邻迭代误差W c 3λ和其近似值lambda k = 2 3X =, , T V k =, , TW c = 600 240 971 801e-004 lambda = 999 999 800 30 020*********.0003=λ k = 23X =, , TV k =, , T W c = 005 395 108 913e-007lambda = 999 999 999 90020*********.0003=λ（）的迭代序列的收敛速度比迭代公式（）稍慢些.（2）再用两种原点位移反幂法的MATLAB 主程序，求979999999990.999~1=λ对应的特征向量.输入MATLAB 程序>> A=[0 11 -5;-2 17 -7;-4 26 -10];V0=[1,1,1]';[k,lambda,Vk,Wc]=ydwyfmf(A,V0,,,100)运行后屏幕显示结果请注意：因为A-aE 的各阶主子式都不等于零，所以A-aE 能进行Lu 分解.A-aE 的秩R(A-aE)和各阶顺序主子式值hl 、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下：hl =RA1 = k = 2 lambda =输入MATLAB 程序>> A=[0 11 -5;-2 17 -7;-4 26 -10];V0=[1,1,1]';[k,lambda,Vk,Wc]=wfmifa1(A,V0, ,,100)运行后屏幕显示结果请注意迭代次数k,特征值的近似值lambda,对应的特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下：k = 3 lambda = Wc =Vk = Wc = Vk =雅可比(Jacobi)方法及其MATLAB 程序雅可比方法是用来计算实对称矩阵的全部特征值和对应的特征向量的一种迭代的方法，最早由雅可比给出.自从计算机出现以后，古典的雅可比方法已有了不少的改进和推广.雅可比方法的MATLAB 程序设n 阶实对称矩阵A 的n 个特征值为n λλλ,,,21 ，特征值i λ对应的特征向量为i X ，则我们可以用下面的MA TLAB 程序计算A 的全部特征值和对应的特征向量的近似用雅可比方法计算对称矩阵A 的特征值和对应的特征向量的MATLAB 主程序输入的量：n 阶实对称矩阵A 、计算的精度jd 、迭代的最大次数max1；输出的量：计算全过程中每次迭代的序号k , 选取的旋转主元mk 及其所在的行数i 和列数j ，c ，t ，pii 和pij 的值，旋转矩阵P k ，正交矩阵V k =P 1P 2…P k ，对称矩阵Bk ,判断是否满足控制迭代终止的条件Wc ，以特征向量为列向量的矩阵V ,特征值为对角元的对角矩阵D .如果迭代次数已经达到最大的迭代次数max1，则给出提示的相关信息.根据雅可比迭代求阶实对称矩阵A 的特征值i 和对应的特征向量i 的近似向量，精度为ε的一般步骤，现提供MATLAB 主程序如下：function [k,Bk,V,D,Wc]=renyujjacobite(A,jd,max1)[n,n]=size(A);Vk=eye(n);Bk=A;state=1;k=0;P0=eye(n);Aij=abs(Bk-diag(diag(Bk)));[m1 i]=max(Aij);[m2 j]=max(m1);i=i(j);while ((k<=max1)&(state==1))k=k+1,aij=abs(Bk-diag(diag(Bk)));[m1 i]=max(abs(aij));[m2 j]=max(m1);i=i(j),j,Aij=(Bk-diag(diag(Bk)));mk=m2*sign(Aij(i,j)),Wc=m2,Dk=diag(diag(Bk));Pk=P0;c=(Bk(j,j)-Bk(i,i))/(2*Bk(i,j)),t=sign(c)/(abs(c)+sqrt(1+c^2)),pii=1/( sqrt(1+t^2)), pij=t/( sqrt(1+t^2)),Pk(i,i)=pii;Pk(i,j)=pij;Pk(j,j)=pii; Pk(j,i)=-pij;Pk,B1=Pk'*Bk;B2=B1*Pk; Vk=Vk*Pk,Bk=B2,if (Wc>jd)state=1;elsereturnendPk;Vk;Bk=B2;Wc;endif (k>max1)disp('请注意迭代次数k 已经达到最大迭代次数max1,迭代次数k,对称矩阵Bk,以特征向量为列向量的矩阵V,特征值为对角元的对角矩阵D 如下：')elsedisp('请注意迭代次数k,对称矩阵Bk,以特征向量为列向量的矩阵V,特征值为对角元的对角矩阵D 如下：')endWc;k=k; V=Vk;Bk=B2;D=diag(diag(Bk));[V1,D1] =eig(A,'nobalance')例 用雅可比方法的MATLAB 程序计算例中矩阵A 的特征值i λ和对应的特征向量i X （4,3,2,1=i ）.解 （1）保存名为为M 文件；（2）输入MATLAB 程序>> A=[12 -56 3 -1;-56 7 2 0;3 2 5 1;-1 0 1 12];[k,B,V,D,Wc]=renyujjacobite(A,,100)（3）运行后屏幕显示计算的全过程中每次迭代的序号k, 选取的旋转主元mk及其所在的行数i和列数j，c，t，pii和pij的值，旋转矩阵P k，正交矩阵V k=P1P2…P k，对称矩阵B k,判断是否满足控制迭代终止的条件W c，以特征向量为列向量的矩阵V,特征值为对角元的对角矩阵D和相关信息如下：k =i = j = mk = Wc =1 2 1 -56 56c = t =pii = pij =Pk =0 00 00 0 00 0 0Vk =0 00 00 0 00 0 0Bk =k =i = j = mk = Wc =2 3 2c = t =pii = pij =Pk =0 0 00 00 00 0 0Vk =0 00 0 0Bk =k = i = j = mk = Wc =3 4 3c = t =pii = pij =Pk =0 0 00 0 00 00 0Vk =0 0Bk =k =i = j = mk = Wc =4 1 3c = t =pii = pij =Pk =0 00 0 00 00 0 0Vk =Bk =k = i = j = mk = Wc =5 4 2c = t =pii = pij =Pk =0 0 00 00 0 00 0Vk =Bk =k =i = j = mk = Wc =6 4 1c = t =pii = pij =Pk =0 00 0 00 0 00 0Vk =Bk =k =i = j = mk = Wc =7 3 2c = t =+002pii = pij =Pk =0 0 0 0 0 0 00 0 0 …………………………………………………………………………请注意迭代次数k,对称矩阵Bk,以特征向量为列向量的矩阵V,特征值为对角元的对角矩阵D如下：V1 =D1 =0 0 0 0 0 0 0 0 0 0 0 0 k = 10 B =V =D =0 0 0 0 0 0 0 0 0 0 0 0 Wc =豪斯霍尔德(Householder)方法及其MATLAB 程序在雅可比方法中，每一次迭代产生两个值为零的非对角元，但接下来的迭代使它们变成非零.因此需要许多次迭代才能使得非对角元足够接近零.豪斯霍尔德(Householder)方法可以在一次迭代中产生多个值为零的非对角元，而且在接下来的迭代中使它们保持为零.豪斯霍尔德方法及其MATLAB 程序设向量O X n ≠⨯1，确定一个n 阶初等反射矩阵n P ，使得1⨯n n X P 的后分量1-n 个求初等反射矩阵P ，使得PX 的第一个分量以外的其余的分量都为零的MATLAB 主程序输入的量：n 维向量O X ≠；输出量：迭代次数k ,对称矩阵B ,以特征向量为列向量的矩阵V ,特征值为对角元的对角矩阵D .如果迭代次数已经达到最大的迭代次数max1，则给出提示的相关信息.function [xigema,rou,miou,P,PX]=renHouseholder(X) n=size(X);nX=norm(X,2); xigema=nX*sign(X(1));rou=xigema*(xigema+X(1)); miou=[xigema,zeros(1,n-1)]'+X,E=eye(n,n); C=2*miou*(miou)'; P=E-C/(norm(miou,2)^2); PX=P*X;例 设向量=X ()T1,2,2，确定一个初等反射矩阵P ，使得PX 的后两个分量为零. 解 输入MATLAB 程序>> X=[2 2 1]'; [xigema,rou,miou,P ,PX]=renHouseholder(X)运行后屏幕显示结果P = PX =矩阵约化为上豪斯霍尔德（Householder ）矩阵及其MATLAB 程序我们讨论的问题是用正交相似变换（即豪斯霍尔德(Householder) 变换）约化一般的矩阵A 为上豪斯霍尔德(Householder) 矩阵.用2-n 步正交相似（即豪斯霍尔德(Householder)）变换将n 阶矩阵A 规约成上豪斯用豪斯霍尔德变换将n 阶矩阵A 规约成上豪斯霍尔德矩阵的MATLAB 主程序 输入的量：n 阶实对称矩阵A ；输出的量：计算全过程中计算次数k ，Sk =)()(,12)(22k k k k a sign A +，uk =1)(22e A k k σ+，ck =2221k μ，Pk =k Tk k c E μμ-，Uk =⎪⎪⎭⎫⎝⎛k k P O O E ,上豪斯霍尔德矩阵Ak . n=size(A); Ak=A; for k=1:n-2k,Sk=norm(Ak(k+1:n,k))*sign(Ak(k+1,k)), uk= Ak(k+1:n,k)+ Sk*eye(n-k,1), ck=(norm(uk,2)^2)/2,Pk= eye(n-k,n-k)-uk*uk'/ck,Uk=[eye(k,k),zeros(k,n-k);zeros(n-k, k),Pk], A1=Uk*Ak;Ak=A1, end例 用初等反射矩阵正交相似约化实矩阵A 为上Householder 矩阵.其中⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=34 19- 37 78- 41- 31 11 72- 98 10.2- 78- 32-94- 21 12 1 0 1- 63- 72 1 5 2 3 17- 32 0 2 7 56- 51- 1712- 34 52- 12A . 解 输入MATLAB 程序>> A=[12 -52 34 -12 17 -51;-56 7 2 0 32 -17;3 2 5 1 72 -63;-1 0 1 12 21 -94;-32-78 98 -72 11;31 -41 -78 37 -19 34]; [k,Sk,uk,ck,Pk,Uk,Ak]=Householdrer1(A) 运行后屏幕显示结果k = Sk = ck = 1 +003 uk = Pk =Uk =0 0 0 0 0Ak =k = Sk = ck =2 +003uk = Pk =Uk =0 0 0 0 00 0 0 0 00 00 00 00 0Ak = …………………………………………………………………………k = Sk = ck =4uk = Pk =Uk =0 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 00 0 0 0Ak =实对称矩阵的三对角化及其MATLAB程序用2-n 步正交相似变换将n 阶实对称矩阵A 规约成三对角形式，可以用下面的MATLAB 程序完成.将n 阶实对称矩阵A 规约成三对角形式的MATLAB 主程序 输入的量：n 阶实对称矩阵A ； 输出的量：三对角矩阵T .[n,n]=size(A); for k=1:n-2s=norm(A(k+1:n,k),2); if (A(k+1,k)<0) s=-s; endr=sqrt(2*s*(A(k+1,k)+s)); U(1:k)=zeros(1,k); U(k+1)=(A(k+1,k)+s)/r; U(k+2:n)=A(k+2:n,k)'/r; V(1:k)=zeros(1,k);V(k+1:n)=A(k+1:n,k+1:n)*U(k+1:n)'; C=U(k+1:n)*V(k+1:n)'; P(1:k)=zeros(1,k);P(k+1:n)=V(k+1:n)-C*U(k+1:n); A(k+2:n,k)=zeros(n-k-1,1); A(k,k+2:n)=zeros(1,n-k-1); A(k+1,k)=-s; A(k,k+1)=-s;A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-2*U(k+1:n)'*P(k+1:n)-2*P(k+1:n)'*U(k+1:n); end T=A;例 用初等反射矩阵正交相似约化实对称矩阵A 为三对角矩阵.其中⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛------------------=5261215121416134237299021237312611451721253233219612371564901435612A 解 输入MATLAB 程序>> A=[12 -56 3 -14 -90 -4;-56 71 23 61 -9 -21;3 23 53 12 -72 51;-14 61 12 7323 21;-90 -9 -72 23 -34 -61;-41 -21 51 21 -61 -52];T=house(A)运行后屏幕显示结果T =0 0 0 0 0 0 00 0 0 0 0 0 0 0 00 0 0 0QR 方法及其MATLAB 程序QR 方法是一种变换方法，它是计算中小型矩阵的全部特征值问题的最有效的方法之一.这种方法主要计算上Householder 矩阵和对称三对角矩阵的全部特征值问题，并且收敛速度快，算法稳定.对于一般的矩阵A ，首先用正交相似变换（即豪斯霍尔德(Householder)方法）约化为上豪斯霍尔德(Householder) 矩阵或对称三对角矩阵，然后用QR方法计算该上Householder矩阵或对称三对角矩阵的全部特征值问题.QR方法的MATLANB程序QR方法也称作正交三角分解（Orthogonal-triangular decomposition）.在MATLAB函数库中有编好的程序，我们只要直接调用这些程序即可.调用方法见表5–4.命令功能[Q,R] = qr (A)（1）如果输入的n阶矩阵A是非奇异的，运行后输出n阶非奇异的上三角形矩阵R和n阶正交矩阵Q，使得A = QR；（2）如果输入的n阶矩阵A是奇异矩阵，运行后输出n阶奇异的上三角形矩阵R和n阶正交矩阵Q，使得A = QR；（3）如果输入的矩阵A是nm⨯阶，nm>且，则运行后输出nm⨯阶上三角形矩阵R （其中R的nm-行和这行以后的元全为零）和m阶正交矩阵Q，使得A = QR，其中实际上仅仅是Q的前n列与R计算.请参考例题例和例= qr (A)（1）如果输入的n阶矩阵A是非奇异的，运行后输出n阶置换矩阵E，n阶非奇异的上三角形矩阵R和n阶正交矩阵Q，使得AE = QR.（2）如果输入的n阶矩阵A是奇异矩阵，运行后输出n阶置换矩阵E，n阶奇异的上三角形矩阵R和n阶正交矩阵Q，使得AE = QR；（3）如果输入的矩阵A是nm⨯阶，nm>且，则运行后输出nm⨯阶上三角形矩阵R （其中R的nm-行和这行以后的元全为零）和m阶正交矩阵Q，使得AE = QR，其中实际上仅仅是Q的前n列与R计算.说明：输出的n阶列置换矩阵E使得上三角形矩阵R的主对角线上的元的绝对值按列从大到小排列.请参考例题例和例= qr (A,0)[Q,R] = QR(A,0)产生"economy size"分解.（1）如果输入的n阶方阵A，则运行后输出的结果与[Q,R] = qr(A)相同；（2）如果输入的矩阵A是nm⨯阶，nm>且，则运行后输出n阶上三角形矩阵R和nm⨯阶矩阵Q （其中Q的列向量两两正交），使得Q的前n列与R的积等于A.请参考例题例= qr (A,0)[Q,R,E] = QR(A,0)产生"economy size" 分解，使得QR = A (:,E).其中E是置换向量，E中的元表示元1所在的行数.例E=[1 3 2]表示置换矩阵为⎪⎪⎪⎭⎫⎝⎛=1111E其它功能与上类似.请参考例题例如果n 阶实矩阵A 的特征值为),,2,1(n i i =λ，则),,2,1(n i s k i =-λ是-=A B k s E 的特征值.将上述思想融入到下面的用带原点位移的QR 方法计算实矩阵的特征值.带原点位移的QR 方法特别适用于计算实对称矩阵的特征值.在这里为了介绍方便，我把用最末元位移法选择位移量k s 的带原点位移的QR 方法求A 的所有近似特征值的方法称作最末元位移QR 方法.如果n 阶矩阵A 是实对称矩阵，则用带原点位移的QR 方法中的最末元位移QR 方法计算A 的所有近似特征值可以用下面的MA TLAB 程序完成，并且给出全部计算过程用最末元位移QR 方法求实对称矩阵A 全部特征值的MATLAB 主程序 输入的量：A 是n 阶实对称矩阵，t 是精度t-=10ε的指数，max1是最大迭代次数. 输出的量：i 表示求第i 个特征值，k 是迭代次数，s k 是原点位移量，Q k 和R k 是B k的QR 分解，A t 迭代矩阵，tzgk 是矩阵A 的特征值的近似值，矩阵B 是m 阶矩阵A t 的（m -1）阶主子矩阵及其相关信息.当迭代矩阵A k 降一阶时，给出相关的提示.QR 方法求实对称矩阵A 全部特征值的MATLAB 主程序如下,取名为.function tzg=qrren4(A,t,max1)[n,n]=size(A); k=0;Ak=A;tzg=zeros(n); state=1; for i=1:n;while ((k<=max1)&(state==1)&(n>1))b1=abs(Ak(n,n-1)); b2=abs(Ak(n,n)); b3=abs(Ak(n-1,n-1)); b4=min(b2, b3); jd=10^(-t); jd1=jd*b4; if (b1>=jd1)sk=Ak(n,n); Bk=Ak-sk*eye(n); [Qk,Rk]=qr(Bk);At=Rk*Qk+sk*eye(n); k=k+1;tzgk=Ak(n,n);disp('请注意：下面的i 表示求第i 个特征值，k 是迭代次数，sk 是原点位移量，')disp(' Bk=Ak-sk*eye(n),Qk 和Rk 是Bk 的QR 分解，At=Rk*Qk+sk*eye(n)迭代矩阵：')i,state=1;k,sk,Bk,Qk,Rk,At,Ak=At; elsedisp('请注意：i 表示求第i 个特征值，tzgk 是矩阵A 的特征值的近似值，k是迭代次数，')disp(' 下面的矩阵B 是m 阶矩阵At 的（m-1）阶主子矩阵，即At 降一阶.')i,tzgk=Ak(n,n),tzg(n,1)=tzgk;k=k,sk,Ak;B=Ak(1:n-1,1:n-1), Ak=B;n=n-1;state==1; i=i+1; end end endtzg(1,1)=Ak;tzg=sort(tzg(:,1));tzgk=Akdisp('请注意：n 阶实对称矩阵A 的全部真特征值lamoda 和至少含t 个有效数字的近似特征值tzg 如下：')lamoda=sort(eig(A))例 用最末元位移QR 方法求下列实对称矩阵的全部近似特征值，并将计算结果与A 全部真特征值比较.其中,2 1 1 1 13 1 21 1 4- 212 2 5)1(⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=A 精度为=ε510-；,52612151214161342372990212373126114517212532332196123715641901435612)2(⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛------------------=A 精度为=ε410-.解 （1）首先保存用最末元位移QR 方法求实对称矩阵A 全部特征值的MATLAB 主程序为M 文件，取名为.在MATLAB 工作窗口输入程序>> A=[5 2 2 1;2 -4 1 1;2 1 3 1;1 1 1 2]; tzg=qrren4(A,5,100)运行后屏幕显示结果请注意：下面的i 表示求第i 个特征值，k 是迭代次数，sk 是原点位移量，Bk=Ak-sk*eye(n),Qk 和Rk 是Bk 的QR 分解，At=Rk*Qk+sk*eye(n)迭代矩阵： i = 1 k = 1 sk = 2 Bk =3 2 2 1 2 -6 1 1 2 1 1 1 1 1 1 0 Qk =0 Rk =0 00 0 0 At =请注意：i 表示求第i 个特征值，tzgk 是矩阵A 的特征值的近似值，k 是迭代次数， 下面的矩阵B 是m 阶矩阵At 的（m-1）阶主子矩阵，即At 降一阶. i = 1 tzgk = k = 1 sk = 2 B =请注意：下面的i 表示求第i 个特征值，k 是迭代次数，sk 是原点位移量，Bk=Ak-sk*eye(n),Qk 和Rk 是Bk 的QR 分解，At=Rk*Qk+sk*eye(n)迭代矩阵： i = 2 k = 2 sk = Bk =。
矩阵特征值计算矩阵的特征值和特征向量
矩阵特征值计算矩阵的特征值和特征向量矩阵是线性代数中的重要概念之一，它在众多学科领域中都有广泛的应用。
而矩阵的特征值和特征向量则是矩阵分析与应用中的核心内容之一。
本文将详细介绍矩阵特征值的计算方法，以及如何求解矩阵的特征向量。
1. 特征值和特征向量的定义首先，我们来了解一下什么是矩阵的特征值和特征向量。
给定一个n阶方阵A，如果存在一个数λ以及一个非零n维列向量X，使得满足下述条件：AX = λX那么，λ就是矩阵A的一个特征值，而X则是对应于特征值λ的特征向量。
特征值和特征向量的求解在很多应用中都具有重要的意义。
2. 特征值的计算方法接下来，我们介绍几种常见的特征值计算方法。
2.1 特征多项式法特征多项式法是求解特征值的一种常用方法。
它利用方阵A减去λ乘以单位矩阵I之后的行列式为零的性质，构造出特征多项式，并求解多项式的根即可得到特征值。
举个例子，对于二阶方阵A = [a, b; c, d]，其特征多项式为：| A - λI | = | a-λ, b; c, d-λ | = (a-λ)(d-λ) - bc = 0解这个方程可以得到A的特征值。
2.2 幂迭代法幂迭代法也是一种常见的特征值计算方法。
它利用特征向量的性质，通过迭代计算来逼近矩阵的特征值。
其基本思想是，给定一个初始向量X0，不断迭代计算：Xk+1 = AXk然后对得到的向量序列进行归一化处理，直到收敛为止。
最后得到的向量X就是对应的特征向量，而特征值可以通过如下公式计算：λ = X^TAX / X^TX2.3 QR方法QR方法是一种数值稳定性较好的特征值计算方法。
它利用矩阵的QR分解的性质来逐步逼近矩阵的特征值。
首先，对矩阵A进行QR分解，得到一个正交矩阵Q和一个上三角矩阵R。
然后，将分解后的矩阵R与矩阵Q逆序相乘，得到一个新的矩阵A'。
重复进行QR分解和相乘的操作，直到收敛为止。
最后，得到的矩阵A'的对角线上的元素即为矩阵A的特征值。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
会达到加速收敛的目的.
如把Aitken加速方法用于例3,则有
k
0 1 2 3 … 10 11 12
k
10 7.2 6.5 … 6.003352 6.001675 6.000837
u(k)
(1,1,1)T (1,0.8,0.1)T (1,0.75,-0.111)T (1,0.730769,-0.188034)T ………………….. (1,0.714405,-0.249579)T (1,0.714346,-0.249790)T (1,0.714316,-0.249895)T
适当选取非奇异对角矩阵D=diag(d1,d2,…,dn),则矩
阵D-1AD与矩阵A有相同的特征值,且对角元素相同. 而矩阵
D-1AD对应的n个圆盘为
n a d ij j Ri | a ii j 1 di j i , i 1,2, , n
v Au max( v (k ) ) k u ( k ) v ( k ) k
(k ) ( k 1)
, k 1,2,3,
n
i 1
可得
u
(k )
A v i 2 n k (0) max( A v ) max( a1 x 1 a i ( ) k x i )
i
k
(0)
a1 x 1 a i ( ) k x i
i2
1
u
所以
(k )
A v i 2 n k (0) max( A v ) max( a1 x 1 a i ( ) k x i )
1 i
k
(0)
a1 x 1 a i ( ) k x i
i
n
i2
( k 1) i
( k 1)
a1 x 1 1 v
k 1 1
(k )
1 v
或取
/v
(k ) i
i 1,2, , n
( ( v nk 1) 1 v1( k 1) v 2k 1) 1 ( ( k ) ( k ) ( k ) ) n v1 v 21 vn
0 1
2 3 4
1 0.25
0.10250 0.042292 0.017451
0 0.2
0.083333 0.034389 0.014190 0.41 0.41260 0.41263 0.41665 0.41267 0.41263
可取0.41263 ,x1(0.017451,0.014190)T .
k 0 1 2 3 … 10 11 12
k 10 7.2 6.5 … 6.003352 6.001675 6.000837
u(k) (1,1,1)T (1,0.8,0.1)T (1,0.75,-0.111)T (1,0.730769,-0.188034)T ………………….. (1,0.714405,-0.249579)T (1,0.714346,-0.249790)T (1,0.714316,-0.249895)T
1
lim u
k
(k )
x1 max( x1 )
max[ a1 x 1 a i ( ) k x i ]
i
其收敛速度由比值|2/1|来确定. 又由于
n
k max( v ) 1
(k )
max[ a1 x 1 a i ( ) k 1 x i ]
i
R(x) (Ax, x) (x, x)
设n阶实对称矩阵A的特征值为12…n ,则有
1 R(x)n
1 max R(x)
x0
, n min R(x)
x0
( A R(x)E)x 2 min ( A E)x 2
可见,如果x是特征向量,则R(x)是对应的特征值. 如果x仅 是一近似特征向量,则R(x)是对应的特征值的一个近似值.
r 1 1 n 1
若a1,a2,…,ar不全为零, 则对充分大的k有
k v ( k ) 1 [a1 x 1 a 2 x 2 a r x r ]
由于a1x1+a2x2+…+arxr 是对应1的特征向量, 若仍记为x1 , 则有: v(k) 1kx1 ,故前面的结论仍然成立. 3. 设1=-2,且 |1=|2|>3 n ,这时,(5.1) 式可写成
§1 乘幂法和反幂法
§1.1 乘幂法
乘幂法是用来求矩阵A按模最大的特征值和相应的特 征向量的方法. 设A是单构矩阵, 即A有n个线性无关的特征向量.
A的n个特征值为
|1 2 n
对应的特征向量为 x1,x2,…xn 线性无关. 我们要求1 和 x1 . 乘幂法的基本思想是取初始向量v(0)Rn,作迭代 v(k+1) =Av(k) =Ak+1v(0) , k=0,1,2,… 产生迭代序列v(k). 由于x1,x2,…xn 线性无关, 从而有 v(0) =a1x1+a2x2+…+anxn 故有 v(k) = Akv(0) =a11kx1+a22kx2+…+annkxn (5.1)
对非零向量v,用max(v)表示v的按绝对值最大的分量,
称向量u=v/max(v)为向量v的规范化向量. 例如, 设向量
v=(2,1,-5,-1)T,则max(v)=-5,u=(-0.4,-0.2,1,0.2)T.可
见规范化向量u总满足‖u‖=1.
乘幂法的规范化计算公式为: 任取初始向量u(0)=v(0) 0,计算
可取 1 6.000837, x1 (1,0.714316,-0.249895)T.
实际上,A的3个特征值分别为1=6,2=3,3=2.
2. 设1=2==r,且 |1>r+1n ,这时,(5.1) 式可写成
k v ( k ) 1 [a1 x1 a 2 x 2 a r x r a r 1 ( ) k x r 1 a n ( ) k x n ]
Ri=aiiri,i=1,2,…,n
则 (1)A的任一特征值至少位于其中一个圆盘内;
(2)在m个圆盘相互连通(而与其余n-m个圆盘互不连通)
的区域内,恰有A的m个特征值(重特征值按重数记).
例1 设矩阵
4 1 0 A 1 0 1 1 1 4
而特征向量 x1 v(k). 乘幂法的收敛速度取决于|2/1|的大小.
例2
求矩阵A的按模最大的特征值
1 4 A 1 5
解
k v1(k)
1 6
1 5
取v(0)=(1,0)T ,计算v(k)=Av(k-1), 结果如下
v2(k) v1(k)/v1(k-1) v2(k)/v2(k-1)
12 13
3.454288 4.631924
(0.539495, 1, 1)T (0.465893, 1, 1)T
2 4
1 12 13 3.454288 4.631924 4,
x1=13u(13)+1u(12)=(4.315961, 8.631924, 8.631924)T,
试讨论A的特征值的分布. 解 由A确定的3个圆盘分别为
R1=-41, R2=2, R3=+42 y 所以
315 -2<22 -63<-2
-6 -4 -2 0 2 3 4 5
x
实际上， 1=4.20308 , 2=-0.442931 , 3=-3.76010
i 2 n
1
所以
i2
1
lim k 1
k
因此,当k充分大时可取:
1 k , x1 u(k) .
如用规范化乘幂法解例2,仍取u(0)=v(0)=(1,0)T,则有 k k u1(k) u2(k) 0 1 0.25 2 0.41 3 4 0.412602 0.412627
~ k
6.266667 ……… 6.000017 6.000003 6.000000
k 0 1 2 3 4 5 6 7 8 9 10 11
k 11 3.553628 4.679204 3.461124 4.635465 3.452655 4.632116 3.454315 4.631929 3.454291 4.631920
u(k) (1,1,2)T (0.454545, 0.909091, (0.537222, 0.972116, (0.465201, 0.994041, (0.539392, 0.998269, (0.465721, 0.999627, (0.539487, 0.999892, (0.465890, 0.999975, (0.539495, 0.999993, (0.465893, 0.999999, (0.539495, 1, 1)T (0.465893, 1, 1)T 1)T 1)T 1)T 1)T 1)T 1)T 1)T 1)T 1)T
1
0
1
0.8
1
36 0.813138
故可取 1 0.412627, x1 (1,0.813138)T.
例3 设
4 14 0 A 5 13 0 1 0 2
用乘幂法求A的按模最大的特征值和相应特征向量. 解 取初值u(0)=v(0)=(1,1,1)T,计算得
1. 设|1>2n , 这时,(5.1)式可写成
k v ( k ) 1 [a1 x1 a 2 ( ) k x 2 a n ( ) k x n ]
2 n 1 1
若a10, 则对充分大的k有
v
因而有
(k )
a1 x 1 , v
k 1
1 vi( k 2) / vi( k )