第3章 线性代数计算方法

第3章 线性代数计算方法,1 高斯消去法 2 高斯约当消去法3 解实三对角线性方程组的追赶法4 矩阵的三角分解5 行列式和逆矩阵的计算6 迭代法 7 迭代法的收敛性8 矩阵的特征值与特征向量的计算,1 高斯消去法,1.1 顺序消去法 如果线性方程组(31)的系数矩阵A具备某特殊形式,例如其为上三角矩阵,且aii0,i=1,2,n,这时方程组(31)实际为,(34),由方程组(34)的最后一个方程直接可得,将其代入倒数第二个方程可求得,如此再解出xn-2,x2,x1,一般有,(35),其中规定,以上讨论告诉我们,对具有上三角形系数矩阵的方程组(34)求解极为方便。当然,若方程组(31)的系数矩阵为下三角形,则求解也很方便。于是对于一般形式的方程组(31),我们总设法把它化为系数矩阵呈上(或下)三角形的方程组来求解。为了达到目的,可利用消去法进行。现举例如下: 解方程组,(36),作-消去中的x1,作-4消去中的x1,则方程组(36)化为,(36),对方程组(36)作- ,得到三角形方程组,(36),从方程组(36“)的方程解出x3,将所得的结果代入方程求出x2,再把x3、x2同时代入方程解出x1。这样可求出方程组的解为 上述求解方程组的方法就是高斯(Gauss)消去法。从式(36)到 (36)的过程称为消元过程而由(36)求出x3、x2、x1的过程称为回代过程。因此用高斯消去法求解性方程组要经过消元和回代两个过程。,在计算机上实现时,我们常把方程组右端的常数项排于系数矩阵的第n+1列,这样顺序高斯消去法的计算步骤为: 1.消元过程 对于k=1,2,n-1,若按顺序有某一ark0,rk,则交换k与r行,然后计算 ,(311),2. 回代过程 对于k=n,n-1,2,1,计算,(312),1.2 主元素消去法 前述顺序消去法是按序通过用a11,a(1)22,a(n-2)n-1 (a(k-1)kk0)作为除数来达到消元目的的。在实际计算时,由于舍入误差的影响,计算结果会改变很大,甚至于完全失真。例如用高斯消去法求解下列方程组(用四位有效数字计算):,-105得 (1.000-1.000105)x2=1.000-6.000104化简可得 x2=0.6000回代求得 x1=105(0.6-0.6000)=0而方程组的解应为 x1=04000 x2=06000,显然用上述方法求出的解x1与方程组的实际解相差很大。若改变两个方程的顺序,即 x1+x2=1 10-5 x1+x2=0.6 -10-5得 (1000-100010-5)x2=0.6-1.00010-5 化简得 x2=0.6000 回代求得 x1=(1-0.6000)=0.4000,高斯主元素消去法是顺序消去法的一种改进。它的基本思想是在逐次消元时总是选绝对值最大的元素(称之为主元)做除数,按顺序消去法的步骤消元。 这里主要介绍求解线性方程组最常用的列主元素消去法和全主元素消去法。 1.列主元素消去法 所谓列主元素消去法就是在每一步消元过程中取系数子矩阵的第一列元素中绝对值最大者作主元。对线性方程组(31)进行n-1次消元后,可得到上三角形方程组,必须指出的是方程组(313)中的系数aij(ij)和右端的bi已经改变了,并非与原来相同。这样就可对方程组(313)回代求解。 例1 用列主元素消去法解方程组,(313),取四位有效数字计算。 解 中-18为主元,交换和得, , ,+12/18,+1/18得, ,第二列消元时,主元为1167,交换方程和得, ,+1/1167得, ,回代求得 x1=1000,x2=2000,x3=3001方程组的实际解 x1=1,x2=2,x3=3,在计算机上求解时,前面已经讲过,常把右端常数项 b 作为系数矩阵A的第n+1列,从而得到增广矩阵(A, b), 仍记为A,于是, 列主元素消去法的计算过程为: (1)消元过程。对于k=1,2,n-1进行下述运算: 选主元,确定r,使得 若ark=0,说明系数矩阵为奇异,则停止计算否则进行下一步。 交换A的r、k两行 对i=k+1,k+2,n,j=k+1,k+2,n+1计算,(314),aij-aikakj/akkaij (315) (2) 回代过程。对于k=n,n-1,1计算,(316),2. 全主元素消去法 所谓全主元素消去法,就是每步消元时选取系数子矩阵中绝对值最大的元素作主元。经过n-1次消元后,方程组(31)可化为上三角形方程组,(317),这里方程组(317)中的系数aij(ij)及bi一般改变了。特别是未知数的排列顺序,由于全主元素的消元过程中,其系数矩阵可能进行了列对调,那么未知数也相应地作了对调。例如,若矩阵第i列与第j列对调,则未知数xi与xj也相应地对调了,xi的结果实质上为xj的结果。,图 3.1,图 3.1,例2 用全主元素消去法求解方程组, ,解这里主元为-18,交换方程与方程得, ,再全选主元,主元为2333,交换x2和x3所在的两列,同时改变两未知数的排列号得, ,-0944/2333得, ,已经化为三角方程组,回代求解 x1=1000,x2=3000,x3=2000 这里未知数x2与x3已对调,所以应恢复解的顺序,方程组的实际精确解为 x1=1000,x2=2000,x3=3000 在计算机上求解时,我们仍把增广矩阵(A, b) 记为A,具体计算过程为 (1)消元过程。对k=1,2,n-1进行下列运算:,图 3.2,图 3.2, 选主元,确定r,t使得 若art=0,则系数矩阵为奇异的,停止计算否则进行下一步。 交换A中的r、t两行及t、k两列,并记下交换的足码t、k。 对i=k+1,k+2,n,j=k+1,k+2,n+1计算 (2) 回代过程。对k=n,n-1,1计算,(318),(319),(320),2 高斯约当消去法,前面所述的消去法均要进行两个过程,即消元过程和回代过程。但对消元过程稍加改变可以把方程组(31)化为对角形,(321),此时求解就不要回代了。这种无回代过程的主元素消去法称为 高斯约当(Jordan)消去法。 特别是方程组(321)还可化为,(322),显然等号右端即为方程组的解。 对于n阶线性方程组(31),其增广矩阵为,首先把主元素(按列选主元或全选主元)调换到主对角线上,并化为1,再将主元素所在列的其它元素消为0,则第一次消元后增广矩阵化为,显见经n步消元后即得方程组的解。 其具体计算步骤为: 对k=1,2,n 选主元(列选或全选主元) 对j=k,k+1,n+1计算,(323), 对i=1,2,n(ik),j=k+1,k+2,n+1计算 aij-aikakj aij (324) 计算 xk=a kn+1 还得指出,若用到全选主元,最后应注意恢复解的顺序。,3 解实三对角线性方程组的追赶法,实三对角线性方程组是一类特殊的方程组。我们将利用所谓“追赶法”解决。 设所给的实三对角线性方程组为,其中 a1c10 akbk+ck, bkck0,k=2,3,n-1 anbn0 则三对角矩阵A的各阶主子式都不等于0,也即A为非奇异。用数学归纳法证明: 令A的k阶主子式为k, 2=a1a2-b2c10 其中,因为,所以方程式(326)右端的k-1阶行列式满足假设的条件,由此可得k0,即矩阵A的各阶主子式都不等于0。 由于实三对角线性方程组的系数矩阵A为非奇异,则有唯一的一组解。下面来求解。 从方程组(325)的第一个方程解得,令,(327),(328),将式(328)代入(325)式的第二个方程解出x2,(329),令,则,再将(329)式代入(325)式的第三个方程解出x3,一般有 xk=uk-vkx k+1,k=1,2,n-1 (330) 其中,k=2,3,n-1 (331),当k=n时,因为 bnx n-1+anxn=rn又 x n-1=u n-1-v n-1 xn于是,(332),例3 解线性方程组,解 依公式(327)、(331)求得,v3=0(一般v3不必算)再由式(332)、(330)求得方程组的解,图 3.3,图 3.3,4 矩阵的三角分解,对于线性方程组(31),若其系数矩阵A能够分解成两个三角形矩阵相乘,譬如, A=LU L为下三角矩阵,U为上三角矩阵,则求解方程组(31)就化为求解 LUx=b,令,Ux=y,(334),(335),则 Ly=b (335)、(336)都为三角形方程组,先求出(336)中的 y,然后代入(335)就可以求出原方程组(31)的解x。,(336),4.1 消去法与矩阵的初等变换 从矩阵的初等变换来看,前面所述的主元素消去法(假定主元素经过行列调换以后均在对角线上),就是通过一系列的初等变换将方程组(31)的系数矩阵A化为三角矩阵。其每一步的消元过程相当于对增广矩阵(A,b)作一次初等变换。 令 在消去第一列中的ai1(i=2,3,n)时,左乘的矩阵为,这样就把A化为一个上三角矩阵U，即,在需要进行行交换时，还得左乘某些排列矩阵Pij,A=LU (337) 称式(337)为矩阵A的LU分解或三角分解。 由上述讨论,只要系数矩阵A非奇异,经过适当的行交换后总可以将A分解为一个下三角矩阵L和一个上三角矩阵U之积。下面首先讨论矩阵三角分解的唯一性。,4.2 矩阵三角分解的唯一性 设非奇异矩阵A存在一种三角分解LU,一般这种分解未必唯一,这是因为任意一个同阶的非奇异对角矩阵D总有 (LD)(D-1 U)=L1U1 它也是矩阵A的一种三角分解。为了确保三角分解的唯一性,需对分解加以规格化。 为使LU分解规格化,可以把矩阵U换成DU,其中L为单位下三角矩阵,U为单位上三角矩阵,D为对角矩阵,即,我们称 A=LDU(338) 为矩阵A的一种LDU分解。关于矩阵A的LDU分解有如下定理: 定理1对于nn阶矩阵A存在唯一的LDU分解的必要充分条件是A的各阶主子矩阵A1,A2,An均非奇异。 证先证必要性。设A有唯一的LDU分解,即 A=LDU 由于 Ak=LkDkUk,其中Ak、Lk、Dk、Uk分别表示A、L、D、U的k阶主子矩阵。因为L、U是单位三角矩阵,D为非奇异对角矩阵,所以Lk、Uk、Dk均非奇异,从而知Ak必为非奇异矩阵。 再证充分性。用数学归纳法证明: 当k=1时,A1显然有这种三角分解且是唯一的,即 A1=L1D1U1=l11d1u11 其中 l11=1,d1=a11,u11=1,设对k=n-1时结论成立,将A按下面方式分块:,这里 s=(a1n,a2n,a n-1n)T rT=(an1,an2,a nn-1) a=ann再将A、D、U仿照A分块:,这里,设,若通过A的元素能唯一地定出L中的IT,U中的u及D中的d,就可确认A唯一地有LDU分解了。因为,其中r、s、a均为已知,于是有 L n-1 D n-1 u=s (D n-1 U n-1 )TI=r ITD n-1 u+d=a,由于A的各阶主子矩阵非奇异,因而Ln-1、Dn-1、Un-1 均为非奇异,于是可以唯一地求出 所以,A存在唯一的LDU分解。 证毕。由定理1不难得到定理2(证明留给读者)。定理2 nn阶矩阵A非奇异的必要充分条件为A(或A经行、列变换后)存在LDU分解。,推论1 奇异矩阵不能进行LDU分解。 推论2 若方阵A有奇异主子矩阵,则A不能直接进行LDU分解。 证 用反证法。若A能直接进行LDU分解,即A=LDU。设A的k(kn)阶主子矩阵Ak奇异,则按分块矩阵的乘法有 Ak=LkDkUk 而Lk为单位下三角矩阵,Dk为非奇异对角矩阵,Uk为单位上三角矩阵,显然与推论1矛盾,所以结论成立。,可验证A为非奇异矩阵,但其左上角二阶主子矩阵奇异,故由推论2不能直接进行LDU分解。因A为非奇异,由定理2,经适当行变换后存在LDU分解。若经行变换,例 设,则LDU分解为,A1=LDU,若经列变换,则A2的LDU分解为,4.3 LU分解方法 我们已经知道,若矩阵A的各阶主子矩阵非奇异,则存在唯一的LDU分解,即 A=LDU 常常称 (LD)U=L1U,L(DU)=LU1 (339) 分别为矩阵A的克劳特 (Crout)分解和 杜利特尔(Doolittle)分解。仍记为 A=LU,只不过若是克劳特分解,则L为下三角矩阵,U为单位上三角矩阵若是杜利特尔分解,此时L为单位下三角矩阵,U为上三角矩阵。显然这两种分解均是唯一的。现分别介绍求解线性方程组的两种LU分解方法。 1.克劳特分解方法 设A为nn阶非奇异矩阵,且各阶主子矩阵为非奇异,则矩阵A的克劳特分解为 A=LU (340) 其中,L、U中元素的确定,可不用消去法,而直接由式(340)来求得。按矩阵乘法规则有 因为u11=1,所以 a i1=li1u11=li1, i=1,2,n 于是可求得L的第一列元素为 li1=ai1,i=1,2,n (342) 又因为 a1j=l11u1j,j=2,3,n 所以 u1j=a1j/l11,j=2,3,n (343),(341),于是U的第一行元素也就由式(343)计算出来了。假定L的前k-1列以及U的前k-1行元素都已算出,由于ukk=1,于是从式(341)有,(344),(345),这样,L、U中的元素都已求出。计算L的各列与U的各行的次序如图3。4所示 。,图 3.4,对方程组(31)的系数矩阵A作出LU分解后,方程组便化为 LUx=b 则求解上列方程组就化为依次解方程组 Ly=b (346) Ux=y (347) 由于L为下三角矩阵,U为单位上三角矩阵,故方程组(346)、(347)的求解极为方便。先由(346)式解出y,然后代入(347)式求得x,亦即求得了方程组(31)的解。而求解(346)、(347)式的计算公式分别为,用克劳特分解求解线性方程组(31)的计算过程为: LU分解过程:对于k=1,2,n依次计算,(348),(349),从上述计算过程还可看出,求解yk的公式完全可以插在计算L、U元素的两个公式之间。其计算框图见图3.5。,图 3.5,图 3.5,例4 用克劳特分解方法求解下列方程组,解 令,利用矩阵乘法可得到,这样原方程组就化为依次求下列两个三角形方程组,代入第二个方程组可求得原方程组的解为,2. 杜利特尔分解方法 杜利特尔分解方法是令 A=LU 这里,关于L、U的元素可按以下公式计算: 对于k=1,2,n进行,(350),(351),在L与U的元素算出后就可以按Ly=b(352) Ux=y(353)解方程组LUx=b,从而求得方程组(31)的解。,例5 用杜利特尔分解方法求解下列方程组的解:,解 利用式(350)、(351)可求得,再由式(352)、(353)求得原方程组的解为,从上面的讨论我们可以看到,用LU分解方法解线性方程组(31),就是将其化为依次求解下面两个三角形方程组: Ly=b Ux=y,显然求解Ly=b的过程即为消元过程,它只不过是对 LUx=b 两端左乘L-1而得,即 Ux=L-1 b=y 而求解Ux=y的过程即为回代过程。因此三角分解方法实质上还是高斯消去法。而克劳特分解方法中的元素lii(i=1,2,n)和杜利特尔分解方法中的元素uii(i=1,2,n)实际上就是高斯消去法中的各次主元素。,4.4 乔累斯基(Cholesky)分解方法 设A为nn阶正定矩阵,记为A0。根据前节的讨论,其有唯一的分解式 A=LDU 其中L为单位下三角矩阵,D为非奇异对角矩阵,U为单位上三角矩阵。由于A是对称的,于是有 U=LT 从而A的LDU分解式可写成 A=LDLT (354) 又因A为正定矩阵,故D的主对角元素均为正数。记,则A可以表示成 为简便起见,我们去掉L1的下标,仍把A的这种分解记作 A=LLT (355) 其中L是一个下三角矩阵。分解式(355)便是通常所说的正定矩阵A的乔累斯基分解。,下面介绍两种具体的乔累斯基分解方法。 1. LLT分解方法(平方根法) 对于分解式,由矩阵乘法规则,从而得,(356),(357),当i=j时,也即,于是,这样由式(356)、(357)确定L的元素后,就可按 Ly=b (358) LTx=y (359)解方程组 LLTx=b (360)其解x即为方程组 Ax=b的解。,由式(358)、(359)依次计算y和x可按下列递推公式进行:,(361),(362),例6 用LLT分解方法解线性方程组,取四位小数计算。 解用公式(357)、(356)依次计算,由式(361)求得 y11.3416,y2-0.4474,y31.7320再由式(362)求得原方程组的解为 x10.9993,x20.9989,x30.9999实际上,原方程组的准确解为 x1=1,x2=1,x3=1LLT分解方法的计算框图见图36。,2. LDLT分解方法 LLT分解方法的计算中含有开方运算,而LDLT分解方法就避免了这一点。 设A0,则其有唯一的分解式 A=LDLT 其中,图 3.6,图 3.6,类似于前面的讨论,由矩阵乘法规则可得,(363),也即,因为ljj=1,从而有,故有,当i=j时,(364),根据式(363)、(364)求出矩阵L、D的元素后,就可按照 Ly=b (365) LTx=D-1 y (366),解方程组 LDLTx=b 其解x就是线性方程组(31)的解。 具体计算x和y的递推公式为,(367),(368),图 3.7,图 3.7,5 行列式和逆矩阵的计算,5.1 行列式的计算 行列式的数值计算,只要将行列式化为三角形式(上三角或下三角),求解就很方便了。因此,前面所述的解线性方程组的方法都可以用来化行列式detA为三角形式。 1. 高斯消去法 设 detA为n阶行列式,即,按矩阵A经若干次初等变换以后,其主元素分别为r11,r22,rnn,于是,2. LU分解法克劳特分解法:detA=detLdetU=detL=l11l22lnn杜利特尔分解法:detA=detLdetU=detU=u11u22unn乔累斯基分解法:这里A为正定矩阵。detA=detLdetLT=l211l222l2nn或者detA=detLdetDdetLT=detD=d1d2dn,5. 逆矩阵的计算设A为非奇异矩阵,则A可逆,即有AA-1=I令A-1=X=(x1,x2,xn)于是有AA-1=A(x1,x2,xn)=(e1,e2,en),从而求逆矩阵A-1的问题就化为求解下列几个有相同系数矩阵A的线性方程组问题 Axi=ei, i=1,2,n (369) 因此,前面介绍的求解线性方程组的方法均可用来计算逆矩阵。,例7 设,求A-1。解将方程组写成同一个增广矩阵,即,先对第一列消元,将a11化为1,a21、a31化为0,有,再将第二列a32化为0,得,将a33化为1,并将a23、a13化为0,得,最后将a13化为0,得,所以,6 迭代法,6.1 简单迭代法 设所给的线性方程组为 Mx=g (370) 其中,且系数矩阵M为非奇异,g0。将方程组(370)改写成等价形式 x=Ax+f (371) 这种改写的方法很多,例如将M分解为两个矩阵之差 M=B-C 其中矩阵B可逆,于是方程组(370)成为 Bx=Cx+g x=B-1 Cx+B-1g 令 A=B-1 C,f=B-1 g,即得(371)式。当然选取的B应该便于求逆,如B为单位矩阵或对角矩阵等。 对任意给定的初始向量x(0),根据(371)构造迭代公式 x(k+1)=Ax(k)+f, k=0,1,2, (372) 这里A称为迭代矩阵。用分量形式可写成,(373),算出迭代公式(372)的解的各次近似值x(1),x(2), ,x(k),。这种迭代求解的方法称为简单迭代法。式(372)称为简单迭代公式。特别当M的对角元素均不为零且按绝对值来说较大时,常取,则 A=D-1 C, f=D-1 g,例8 用简单迭代法解下列方程组,解将方程组写成等价形式,取初始值x(0)= 0,按迭代公式,表 31,图 3.8,6.2 赛德尔(Seidel)迭代法 从简单迭代法看到,已知x(k)计算x(k+1)时,需要保留x(k)和x(k+1)两个分量。逐次用前面算出的新分量来计算下一个分量,这就是赛德尔迭代法的基本思想。 设方程组(370)的等价形式为 x=Ax+f 将矩阵A分解为 A=B+C 其中,于是 x=Bx+Cx+f (374) x(k+1)=Bx(k+1)+Cx(k)+f, k=0,1,2, (375),例9 用赛德尔迭代法解方程组,解 将原方程组写成等价形式并按(375)构造赛德尔迭代公式,表 32,由式(375)可得 (I-B)x(k+1)=Cx(k)+f 因为矩阵I-B是非奇异矩阵,则I-B可逆,所以迭代公式(375)可写成 x(k+1)=(I-B)-1 Cx(k)+(I-B)-1 f,k=0,1,2, (377) 这里迭代矩阵 A=(I-B)-1 C 由此看出,对方程组(370)用公式(375)作赛德尔迭代相当于对方程组(370)用(377)式作简单迭代。因此,赛德尔迭代法实际上是某种简单迭代法。,6.3 松驰法 赛德尔迭代公式(375)为 x(k+1)=Bx(k+1)+Cx(k)+f 现令 x=x(k+1)-x(k)=Bx(k+1)+Cx(k)+f-x(k)于是 x(k+1)=x(k)+x (378),x(k+1)可以看作在向量x(k)上加修正项x而得到。若在修正项的前面加上一个参数,便得到松驰法的迭代公式 x(k+1)=x(k)+x=(1-)x(k)+(Bx(k+1)+C x(k)+f) k=0,1,2, (379) 或者用分量形式写成,(380),其中叫做松驰因子,当1时叫超松驰,1时叫低松驰,=1时就是赛德尔迭代法。松驰因子的选取直接影响到松弛法的收敛性。 因为(I-B)-1存在,所以还可以把(379)改写成 x(k+1)=(I-B)-1(1-)I+C)x(k)+(I-B)-1f (381) 这是简单迭代公式,迭代矩阵为 A=(I-B)-1 (1-)I+C) 实际上,赛德尔迭代法或松弛法都是某种简单迭代法,它们仅是迭代矩阵A不同。,7 迭代法的收敛性,前面介绍的例7、例8分别用简单迭代法和赛德尔迭代法都得到了较好的迭代效果。若将原方程组改写成如下等价形式:,表 33,可见当k越大时,其迭代值与准确解 x1=1.1, x2=1.2, x3=1.3 相差越大,即向量序列x(k)不收敛。此例如用赛德尔迭代法迭代求解,得到的结果也令人失望。为此,很有必要讨论迭代法的收敛性。也就是研究等价方程 x=Ax+f 的收敛条件。为了达到此目的,我们首先介绍向量范数、矩阵范数等基本知识,然后再讨论迭代法的收敛性。,7.1 向量范数 在三维空间中,常用三种方法来衡量一个向量r=(x,y,z)T的“长度”,即,这种衡量的方法可推广到n维向量 x=(x1,x2,xn)T也即我们下面介绍的向量范数。,定义1设n维向量xRn,记对应向量x的一个实数为x,若x满足下面三个性质: (1)非负性,即x0,当且仅当x=0时,x=0 (2)齐次性,x=x,为任意实数 (3)三角不等式性,x+yx+y,yRn。 则称该实数x为向量x的范数。 可验证,对于不小于1的正数p,具有向量范数的三条基本性质,称之为向量x中的p范数,记为xp,即,在Rn中,常用的几种向量范数有:,7.2 矩阵范数 在数值计算中,为了进行某种估计,常常要比较不同向量的范数,向量有时以Ax的形式出现,其中A为nn阶矩阵 A=(aij)nn 这就需要寻求x和Ax之间的某种关系。我们希望有 AxCx 其中C为某一非负实数,它与x无关,只与A有关。为此只要取,即可。因为,定义2 设A为nn阶矩阵,定义,为矩阵A的范数。可以证明,矩阵范数满足下列的几个性质:(1)非负性,A0,当且仅当A= 0时,A=0(2)齐次性,A=A,为任意实数(3)三角不等式性,A+BA+B,B与A为同阶矩阵(4)AxAx(5)ABAB,7.3 迭代法的收敛性 定义3若,则称向量序列x(k)收敛于向量x,记为,定义4若,则称矩阵序列A(k)收敛于矩阵A,记为,谱半径设nn阶矩阵A的特征值为i(i=1,2,n),则称,(382),为矩阵A的谱半径。 定理3设(A)是矩阵A的谱半径,则对于A的任何一种范数A,都有 (A)A (383) 设矩阵A的任一特征值i与其对应的特征向量x有关系式 Axi=ixi,两端取范数,再利用矩阵范数的性质 ixiAxi 因为xi 0,则xi0,所以 iA 故 (A)A,定理4设A是nn阶矩阵,由A的各次幂所组成的矩阵序列Ak收敛于零,即,的必要充分条件是 (A)1,定理5 对于任意初始向量x(0)和常数项f,由迭代公式 x(k+1)=Ax(k)+f, k=0,1,2, 收敛的充要条件是 (A)1,定理6 若迭代矩阵A的范数A1,则迭代公式x(k+1)=Ax(k)+f收敛,且有误差估计式,(387),(388),8 矩阵的特征值与特征向量的计算,求矩阵的特征值与特征向量是自然科学许多分支中一个重要的问题。nn阶矩阵A的特征值是其特征方程 det(A-I)=0 (396) 的根,而属于特征值的特征向量是线性方程组 (A-I)x= 0 (397) 的非零解。,8.1 乘幂法 在实际问题中,往往不需要求矩阵A的全部特征值,而只要计算其按模最大或最小的特征值,也即只需求矩阵A的部分特征值。乘幂法便是一种行之有效的方法。 设n阶实矩阵A有n个线性无关的特征向量,其特征值按模排列如下: 12n(398) 其相应的特征向量分别记为 x1,x2,xn 现求按模最大的特征值1及其相应的特征向量x1。,取非零初始向量v0,作迭代构成一向量序列,(399),因为特征向量x1,x2,xn线性无关,故v0可以表示成这n个线性无关的向量的线性组合,即 v0=1x1+2x2+nxn (3100),于是 vk=Ak(1x1+2x2+nxn) =1k1x1+2k2x2+nknxn (3101) 设1为实数,将式(3101)改写成,(3102),若12,由式(398)知,因此,当k充分大时,只要10,就有 vkk11x1 (3103) (A-1I)vk1k1(A-1I)x1=0 (3104) 这说明vk是相应于特征值1的特征向量的近似向量。 由式(3103),当k充分大时, v k+1k+111x1=k11x11vk (3105),例10 求矩阵,的按模最大的特征值与相应的特征向量。,表 34,(A6v0)1/(A5v0)13.4 (A6v0)2/(A5v0)23.4 (A6v0)3/(A5v0)33.4 (A7v0)1/(A6v0)13.4 (A7v0)2/(A6v0)23.4 (A7v0)3/(A6v0)33.4,于是有 13.4相应的特征向量为 x1(792,-1120,792),若把x1的第一个分量化为1,则有 x1(1,-1.4,1) 事实上,矩阵A的特征值与相应的特征向量为,在计算时,为避免特征值的模过大或过小而造成的溢出,现对算法改进如下:,(3106),其中mk为vk中按模最大的分量。式(3106)实际上是每次对向量vk标准化成uk,然后再迭代。容易验证,这说明uk收敛到1所对应的标准化了的特征向量x1。此外,由于 maxv k+1=maxAuk,于是当k时,(3107),因此有当k时 m k+1=maxv k+11 (3108) 这样当k充分大时 1 m k+1 (3109) 1的符号可用uk的第一个非零分量在相邻两次计算中的符号来决定。若相邻两次的符号相同,则 1 m k+1 否则 1- m k+1 此时uk即为相应按模最大的特征值1的特征向量的近似向量。,图 3.9,图 3.9,设r为矩阵A按模最小的特征值,x(r)是相应的崐特征向量,若detA0,则r0,由于 Ax(r)=rx(r) 于是 (3110) 可见1/r是A-1的按模最大的特征值,且x(r)是A-1的相应于1/r的特征向量。 任取一非零向量v0,按乘幂法作迭代 v k+1=A-1 vk, k=0,1,2, (3111),从理论上讲可求得矩阵A的按模最小的特征值r及相应的特征向量x(r)。但A-1往往是未知的,所以实际求解时