南開大學(xué)計(jì)算物理講義4

1、15 Linear equations 5.1 Gauss elimination 5.2 Pivoting5.3 LU factorisation 5.4 The determinant and Inverse of a matrix5.5 Banded matrices and Tridiagonal matrices5.6 Other approaches to solving linear systems2In this chapter we deal with basic matrix operations, such as the solution of linear equati

2、ons, calculate the inverse of a matrix, its determinant etc.35.1 Gauss elimination 4例1 用高斯消元法求解方程組=+=-+=+53213614282321321321xxxxxxxxx=+=-+=+53213614282321321321xxxxxxxxx=+-=-+=+-91012414282321321321xxxxxxxxx704=-=-+=+12603114282321321321xxxxxxxxx001522814321=-=xxx111332=+=xx26123-=-=x6matrix operat

3、ions78910 x1 = x2 = x3 = 1.115.2 Pivoting主元交換法主元交換法列主元交換法列主元交換法行主元交換法行主元交換法 Partial pivoting部分主元交換法部分主元交換法12例如例如. .求解方程組求解方程組Tx)3675. 0 ,05104. 0,4904. 0(*-=xxx000. 3000. 2000. 1643. 5072. 1000. 2623. 4712. 3000. 1000. 3000. 2001. 0321=-:4,4位有效數(shù)字為精確解舍入到位浮點(diǎn)數(shù)進(jìn)行計(jì)算用13用高斯順序消去法求解.)4000. 0 ,09980. 0,4000.

4、0(是一個(gè)很壞的解是一個(gè)很壞的解解得解得Tx- - -= =00. 2002. 1000. 100. 500005. 32.0040000. 3000. 2001. 0003. 2002. 1000. 1006. 6001. 40005. 32.0040000. 3000. 2001. 0000. 3000. 2000. 1643. 5072. 1000. 2623. 4712. 3000. 1000. 3000. 2001. 0)|A(.b - - -= =Tx)3675. 0 ,05104. 0,4904. 0(*-=14b6866.0500.0000.38676.1008015.13.1

5、7605.6431.00722.000-0015.1500.0000.30023.30005.208015.13.17605.6431.00722.000-000.1000.2000.3000.3000.2001.0623.43.7121.000-643.5072.1000.2000.3000.2000.1643.5072.1000.2623.4712.3000.1000.3000.2001.0)|A(.3i - - - - -= = =Tx)3676.0,51080.0,88544.0(- - -= =得得用主元交換法求解用主元交換法求解Tx)3675. 0 ,05104. 0,4904.

6、0(*-=155.2.1 Partial pivoting部分主元交換法部分主元交換法16行主元交換法行主元交換法175.2.2 Full pivoting全主元交換法全主元交換法1819而交換列時(shí)而交換列時(shí)解的解的次序發(fā)生改變次序發(fā)生改變解的次序不變。解的次序不變。交換行時(shí)交換行時(shí), ,20；，對(duì)；中選絕對(duì)值最大者中選絕對(duì)值最大者從從選主元選主元、第一步消元、第一步消元jtilathenaifdonjnitlanjiaijijij= = =j. The diagonal is 1 and the upper matrix elements are zero. We solve thisequ

7、ation column by column (increasing order of j). 55In a similar way we can define an equationwhich gives us the inverse of the matrix C, labelled E in the equation below. with the following general equation5611111111/ 11cece=2212111222121211/0cceecece-=+22222222/ 11cece=332312211113331323121311/ )(0c

8、ceceececece-=+57A calculation of the inverse of a matrix could then be implemented in the following way:Set up the matrix to be inverted.Call the LU decomposition function.Check whether the determinant is zero or not.Then solve column by column eqs.58IAXn=IACn=BIn.BA=-1定理定理 設(shè)A為非奇異矩陣，方程組的增廣矩陣為，如果對(duì)C應(yīng)用

9、高斯-若當(dāng)方法化為，則2.2高斯-若當(dāng)方法59=563452231AA1-例例2 用高斯若當(dāng)消元法求的逆矩陣. =1005630104520012313I |AC60 9103130012010001231-10133100012010035201- 11133100012010231001- -=-1330122311A61Finally, an important property of hermitian and symmetric matrices is that they have realeigenvalues.62Real orthogonal matrices is1TA(A-

10、=）IA(A(AA(1TTT=-）IA(A(A(AT1TT=-）ijkkkaa=jiijkkkiaa=j63Eigenvalue problemAx= x 64Eigenvalue problemAx= x(A- I)x=0Solutions for the system (1) exist if, and only if, the determinant of the matrix A- I is zero65Eigen-vectors Ax= x(A- I)x=0Once we have determined the eigenvalues,we may solve the system

11、of linear equations to find a set of solutions for each value of . These solutions are called the eigen-vectors. 66矩陣特征值和特征向量的數(shù)值解法矩陣特征值和特征向量的數(shù)值解法設(shè)矩陣RA，如果存在數(shù)及非零向量 x滿足方程 xAx =，則稱為矩陣A的一個(gè)特征值， 特征值的計(jì)算，直接從特征方程0)det(=-AI 當(dāng) n 稍大一些，會(huì)遇到很大困難。行列式展開本身就很不容易，隨后是高次代數(shù)。 因此，矩陣特征值的求解，主要是數(shù)值解法。 x 為矩陣A的相應(yīng)于特征值的特征向量。67矩陣特征值和

12、特征向量的數(shù)值解法矩陣特征值和特征向量的數(shù)值解法冪法冪法逆冪法逆冪法古典古典Jacobi法法Jacobi法的改進(jìn)法的改進(jìn)QR算法算法68QR算法是目前求一般矩陣全部特征值和特征算法是目前求一般矩陣全部特征值和特征向量行之有效的一種方法，它適合于對(duì)稱矩向量行之有效的一種方法，它適合于對(duì)稱矩陣，也適合于非對(duì)稱矩陣。陣，也適合于非對(duì)稱矩陣。QR算法最早在算法最早在1961年由年由J.G.Francis提出的，后來經(jīng)一系列提出的，后來經(jīng)一系列數(shù)學(xué)家進(jìn)行深入討論數(shù)學(xué)家進(jìn)行深入討論,并作出了做有成效的改并作出了做有成效的改進(jìn)與發(fā)展。進(jìn)與發(fā)展。QR算法算法69矩陣的矩陣的QR分解分解 QRA =定理定理:

13、設(shè)矩陣設(shè)矩陣A非奇異，則一定存在正交矩陣非奇異，則一定存在正交矩陣Q，上三角矩陣，上三角矩陣R，使，使且當(dāng)要求且當(dāng)要求R的主對(duì)角元素均為正數(shù)時(shí)，則上述分的主對(duì)角元素均為正數(shù)時(shí)，則上述分解式是唯一的。解式是唯一的。QR算法是對(duì)算法是對(duì)A進(jìn)行一系列的正交相似變換，達(dá)進(jìn)行一系列的正交相似變換，達(dá)到求出矩陣到求出矩陣A的全部特征值和相應(yīng)的特征向量的全部特征值和相應(yīng)的特征向量。1T(Q-=）Q分解：kkkRQA = 構(gòu)造：,.)3 , 2 , 1(1=+kQRAkkk 這里 Qk為正交矩陣，Rk為上三角矩陣，且當(dāng) Rk主對(duì)角元均為正數(shù)時(shí)，則上述正交三角分解唯一。 1TQ(Q-=）kkAR）TQ(=,.)

14、3 ,2, 1(1=+kQRQAQAkkkkTkk正交矩陣(1)用用用左乘用左乘(1)式式TQ(2)(3)將將(3)式代入式代入(2)71727374從10A可 以 看 出 ， 已 近 似 接 近 對(duì) 角 矩 陣 ， 即 有 特 征 值,2680. 1,0035. 3,7282. 4321與矩陣 A 的三個(gè)精確解 2679. 133, 3,7321. 433321-=+= 相比，已有良好精確度。隨著迭代次數(shù)增加，nA將收斂到矩陣A 的三個(gè)精確特征值。 755.5 Banded matrices and Tridiagonal matricesIn Banded matrices, the on

15、ly non-zero elements fall within some distance of the leading diagonal. LU Factorisation may readily be modified to account for banded structure. For example, if elements outside the range ai,ib to ai,i+b are all zero, then the summations in the LU Factorisation algorithm need be performed only from

16、 k=i or k=i+1 to k=i+b. Moreover, the factorisation loop FOR q=i+1 TO n can terminate at i+b instead of n.76=-nnnnnnnnnnaaaaaaaaaaA111121232221121177Tridiagonal matrices78=-nnnnnbacbacbacbA1112221179-nnnnnbacbacbacb1112221180三對(duì)角方程組的追趕法81例例 用追趕法解三對(duì)角線方程組=+-=-+-=-+-=-120202124343232121xxxxxxxxxx82835.6

17、 Other approaches to solving linear systemsIteration84)(1)(1)(11122112323121222213132121111-=-=-=nnnnnnnnnnnnnxaxaxaraxxaxaxaraxxaxaxarax1. The Jacobi Iteration系數(shù)矩陣系數(shù)矩陣A可逆且主對(duì)角元素可逆且主對(duì)角元素nna,.,a ,a2211均不為零均不為零.85)(1)(1)(1)3)232)131333)13)2)323)121222)12)1)313)212111)11knnkkkknnkkkknnkkkxaxaxaraxxaxaxa

18、raxxaxaxarax（-=-=-=+-+)1()1()(kikikixxx),()1()1(2)1(1+knkkxxx其中其中 Tnx,.x ,xx002010=為初始向量為初始向量.86下面介紹迭代格式的矩陣表示：下面介紹迭代格式的矩陣表示：)(1)(1)(1)3)232)131333)13)2)323)121222)12)1)313)212111)11knnkkkknnkkkknnkkkxaxaxaraxxaxaxaraxxaxaxarax（-=-=-=+其中 Tnx,.x ,xx002010=為初始向量.), 2, 1(1)(1)1(nixabaxkjnijjijiiiki=-=+8

19、7設(shè)設(shè)D = diag (a11, a22, , ann)，將，將AX = b改寫為：改寫為：FBXXkk+=+)()1(稱為雅克比迭代矩陣。稱為雅克比迭代矩陣。), 2, 1(1)(1)1(nixabaxkjnijjijiiiki=-=+AX = (D (D - A) X = bDX = (D - A) x + bX = (I D-1A) x + D-1b記記 B = I D-1A F = D-1 b則迭代格式的向量表示為則迭代格式的向量表示為88bAx =2 .02 .1)0(x前三次疊代的結(jié)果，要求取前三次疊代的結(jié)果，要求取)(1)(1)2)323)121222)12)1)313)212

20、111)11knnkkkknnkkkxaxaxaraxxaxaxarax（-=-=+89)1 (21)22(4131)3(31)1)1)12)2)2)11kkkkkkxxxxxx（-=-=-=-=+1 .0)2 .11 (21)1 (21933.032 .0131)01)12)02)11-=-=-=-=-=（xxxxbAx =2 .02 .1)0(x901 .0)2 .11 (21)1 (211514933.032 .0131)01)12)02)11-=-=-=-=-=（xxxx301033.0)933.01 (21)1 (213031033.131 .0131)11)22)12)21=-=-

21、=-=-=（xxxx601017.0)033.11 (21)1 (219089989.03033.0131)21)32)22)31-=-=-=-=-=-=（xxxx91-=017.0989.0)3(x-=100.0933.0)1(x=000.0000.1)(x=033.0033.1)2(x=180/1180/181)4(x92=1342ApbAx =32bp6 .02 .1*33336 .02 .0*2121)01)12)02)11-=-=-=-=-=（xxxx=2 .02 .1)0(x936 .02 .1*33336 .02 .0*2121)01)12)02)11-=-=-=-=-=（xxx

22、x2 .16 .0*33332 .26 .0*2121)11)22)12)21=-=-=+=-=（xxxx6 .32 .2*33334 .12 .1*2121)21)32)22)31=-=-=-=-=-=（xxxx94=-36332012361114238321xxxTx)0,0,0(0=用用JacobiJacobi迭代迭代法求解方程組法求解方程組 前三次疊代的結(jié)果，要求取前三次疊代的結(jié)果，要求取952. The Gauss-Seidel Iteration)(1)(1)(1)3)232)131333)13)2)323)121222)12)1)313)212111)11knnkkkknnkkk

23、knnkkkxaxaxaraxxaxaxaraxxaxaxarax（-=-=-=+In the Jacobi Iteration96)(1)(1)(1)3)1232)1131333)13)2)323)1121222)12)1)313)212111)11knnkkkknnkkkknnkkkxaxaxaraxxaxaxaraxxaxaxarax（-=-=-=+In the Gauss-Seidel Iteration-+)1()1()(kikikixxx),()1()1(2)1(1+knkkxxx Tnx,.x ,xx002010=97高斯高斯賽得爾（賽得爾（Gauss-Seidel）迭代法）迭代

24、法-=+=+-=+)(1)1(11)1(1kjnijijkjijijiiikixaxabax （i = 1, 2, n） FUxLxxkkk+=+)()1()1( FUxxLIkk+=-+)()1()(FLIUxLIxkk1)(1)1()()(-+-+-=98FLIUxLIxkk1)(1)1()()(-+-+-=)()1(FXBXkk+=+FBXXkk+=+)() 1(JacobiJacobi迭代法迭代法 B = I D-1A F = D-1 bULIB1)(-=FLIF1)(-=Gauss-Seidel迭代法迭代法99bAx =2 .02 .1)0(x用用高斯高斯賽得爾（賽得爾（Gauss-

25、Seidel）迭代法）迭代法求解方求解方程組程組前三次疊代的結(jié)果，要求取前三次疊代的結(jié)果，要求取100)1 (21)22(4131)3(31)11)11)12)2)2)11+-=-=-=-=kkkkkkxxxxxx（033.0)933.01 (21)1 (21933.032 .0131)11)12)02)11=-=-=-=-=（xxxxbAx =2 .02 .1)0(x101033.0)933.01 (21)1 (21933.032 .0131)11)12)02)11=-=-=-=-=（xxxx0055.0)989.01 (21)1 (21989.03033.0131)21)22)12)21=

26、-=-=-=-=（xxxx001.0)998.01 (21)1 (21998.030055.0131)31)32)22)31=-=-=-=-=（xxxx102-=017.0989.0)3(x-=100.0933.0)1(x=000.0000.1)(x=033.0033.1)2(x=001.0998.0)3(x=033.0933.0)1(x=006.0989.0)2(xGauss-Seidel IterationJacobi Iteration103從此例看出從此例看出,高斯高斯塞德爾迭代法比雅可比迭代塞德爾迭代法比雅可比迭代法收斂快法收斂快(達(dá)到同樣的精度所需迭代次數(shù)少達(dá)到同樣的精度所需迭代次

27、數(shù)少),但但這個(gè)結(jié)論這個(gè)結(jié)論,在一定條件下才是對(duì)的在一定條件下才是對(duì)的,甚至有這樣的甚至有這樣的方程組方程組,雅可比方法收斂，而高斯雅可比方法收斂，而高斯塞德爾迭代塞德爾迭代法卻是發(fā)散的法卻是發(fā)散的.104定義定義1：如果矩陣的每一行中，不在主對(duì)角線：如果矩陣的每一行中，不在主對(duì)角線上的所有元素絕對(duì)值之和小于主對(duì)角線上元素上的所有元素絕對(duì)值之和小于主對(duì)角線上元素的絕對(duì)值，即的絕對(duì)值，即niaaiinijjij, 2, 11=則稱矩陣則稱矩陣A按行嚴(yán)格對(duì)角占優(yōu)，類似地，也有按行嚴(yán)格對(duì)角占優(yōu)，類似地，也有按列嚴(yán)格對(duì)角占優(yōu)。按列嚴(yán)格對(duì)角占優(yōu)。迭代收斂的充分條件迭代收斂的充分條件105定理：若線性方程

28、組定理：若線性方程組AX = b的系數(shù)矩的系數(shù)矩陣陣A按行嚴(yán)格對(duì)角占優(yōu)，則雅克比迭代法和按行嚴(yán)格對(duì)角占優(yōu)，則雅克比迭代法和高斯高斯賽得爾迭代法對(duì)任意給定初值均收賽得爾迭代法對(duì)任意給定初值均收斂。斂。如果矩陣如果矩陣A嚴(yán)格對(duì)角占優(yōu)，那么高斯嚴(yán)格對(duì)角占優(yōu)，那么高斯賽得賽得爾迭代法的收斂速度快于雅克比迭代法的收斂爾迭代法的收斂速度快于雅克比迭代法的收斂速度。速度。106=-36203312362381114321xxxTx)0,0,0(0=用用 Gauss-Seidel迭代迭代法求解時(shí)方程組法求解時(shí)方程組 前三次疊代的結(jié)果，要求取前三次疊代的結(jié)果，要求取1073. Successive Over Relaxation (SOR)超松馳法超松馳法使用迭代法的困難是計(jì)算量難以估計(jì)，有使用迭代法的困難是計(jì)算量難以估計(jì)，有些方程組的迭代格式雖然收斂，但收斂速些方程組的迭代格式雖然收斂，但收斂速度慢而使計(jì)算量變得很大。度慢而使計(jì)算量變得很大。超松馳迭代法是一種線性加速方法超松馳迭代法是一種線性加速方法。超松馳迭代法是高斯超松馳迭代法是高