南开大学计算物理讲义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位有效数字为精确解舍入到位浮点数进行计算用13用高斯顺序消去法求解.)4000. 0 ,09980. 0,4000.

4、0(是一个很坏的解是一个很坏的解解得解得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而交换列时而交换列时解的解的次序发生改变次序发生改变解的次序不变。解的次序不变。交换行时交换行时, ,20；，对；中选绝对值最大者中选绝对值最大者从从选主元选主元、第一步消元、第一步消元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定理定理 设A为非奇异矩阵，方程组的增广矩阵为，如果对C应用

9、高斯-若当方法化为，则2.2高斯-若当方法59=563452231AA1-例例2 用高斯若当消元法求的逆矩阵. =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矩阵特征值和特征向量的数值解法矩阵特征值和特征向量的数值解法设矩阵RA，如果存在数及非零向量 x满足方程 xAx =，则称为矩阵A的一个特征值， 特征值的计算，直接从特征方程0)det(=-AI 当 n 稍大一些，会遇到很大困难。行列式展开本身就很不容易，随后是高次代数。 因此，矩阵特征值的求解，主要是数值解法。 x 为矩阵A的相应于特征值的特征向量。67矩阵特征值和

12、特征向量的数值解法矩阵特征值和特征向量的数值解法幂法幂法逆幂法逆幂法古典古典Jacobi法法Jacobi法的改进法的改进QR算法算法68QR算法是目前求一般矩阵全部特征值和特征算法是目前求一般矩阵全部特征值和特征向量行之有效的一种方法，它适合于对称矩向量行之有效的一种方法，它适合于对称矩阵，也适合于非对称矩阵。阵，也适合于非对称矩阵。QR算法最早在算法最早在1961年由年由J.G.Francis提出的，后来经一系列提出的，后来经一系列数学家进行深入讨论数学家进行深入讨论,并作出了做有成效的改并作出了做有成效的改进与发展。进与发展。QR算法算法69矩阵的矩阵的QR分解分解 QRA =定理定理:

13、设矩阵设矩阵A非奇异，则一定存在正交矩阵非奇异，则一定存在正交矩阵Q，上三角矩阵，上三角矩阵R，使，使且当要求且当要求R的主对角元素均为正数时，则上述分的主对角元素均为正数时，则上述分解式是唯一的。解式是唯一的。QR算法是对算法是对A进行一系列的正交相似变换，达进行一系列的正交相似变换，达到求出矩阵到求出矩阵A的全部特征值和相应的特征向量的全部特征值和相应的特征向量。1T(Q-=）Q分解：kkkRQA = 构造：,.)3 , 2 , 1(1=+kQRAkkk 这里 Qk为正交矩阵，Rk为上三角矩阵，且当 Rk主对角元均为正数时，则上述正交三角分解唯一。 1TQ(Q-=）kkAR）TQ(=,.)

14、3 ,2, 1(1=+kQRQAQAkkkkTkk正交矩阵(1)用用用左乘用左乘(1)式式TQ(2)(3)将将(3)式代入式代入(2)71727374从10A可 以 看 出 ， 已 近 似 接 近 对 角 矩 阵 ， 即 有 特 征 值,2680. 1,0035. 3,7282. 4321与矩阵 A 的三个精确解 2679. 133, 3,7321. 433321-=+= 相比，已有良好精确度。随着迭代次数增加，nA将收敛到矩阵A 的三个精确特征值。 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三对角方程组的追赶法81例例 用追赶法解三对角线方程组=+-=-+-=-+-=-120202124343232121xxxxxxxxxx82835.6

17、 Other approaches to solving linear systemsIteration84)(1)(1)(11122112323121222213132121111-=-=-=nnnnnnnnnnnnnxaxaxaraxxaxaxaraxxaxaxarax1. The Jacobi Iteration系数矩阵系数矩阵A可逆且主对角元素可逆且主对角元素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设设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前三次叠代的结果，要求取前三次叠代的结果，要求取)(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迭代迭代法求解方程组法求解方程组 前三次叠代的结果，要求取前三次叠代的结果，要求取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）迭代法）迭代法求解方求解方程组程组前三次叠代的结果，要求取前三次叠代的结果，要求取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从此例看出从此例看出,高斯高斯塞德尔迭代法比雅可比迭代塞德尔迭代法比雅可比迭代法收敛快法收敛快(达到同样的精度所需迭代次数少达到同样的精度所需迭代次

27、数少),但但这个结论这个结论,在一定条件下才是对的在一定条件下才是对的,甚至有这样的甚至有这样的方程组方程组,雅可比方法收敛，而高斯雅可比方法收敛，而高斯塞德尔迭代塞德尔迭代法却是发散的法却是发散的.104定义定义1：如果矩阵的每一行中，不在主对角线：如果矩阵的每一行中，不在主对角线上的所有元素绝对值之和小于主对角线上元素上的所有元素绝对值之和小于主对角线上元素的绝对值，即的绝对值，即niaaiinijjij, 2, 11=则称矩阵则称矩阵A按行严格对角占优，类似地，也有按行严格对角占优，类似地，也有按列严格对角占优。按列严格对角占优。迭代收敛的充分条件迭代收敛的充分条件105定理：若线性方程

28、组定理：若线性方程组AX = b的系数矩的系数矩阵阵A按行严格对角占优，则雅克比迭代法和按行严格对角占优，则雅克比迭代法和高斯高斯赛得尔迭代法对任意给定初值均收赛得尔迭代法对任意给定初值均收敛。敛。如果矩阵如果矩阵A严格对角占优，那么高斯严格对角占优，那么高斯赛得赛得尔迭代法的收敛速度快于雅克比迭代法的收敛尔迭代法的收敛速度快于雅克比迭代法的收敛速度。速度。106=-36203312362381114321xxxTx)0,0,0(0=用用 Gauss-Seidel迭代迭代法求解时方程组法求解时方程组 前三次叠代的结果，要求取前三次叠代的结果，要求取1073. Successive Over Relaxation (SOR)超松驰法超松驰法使用迭代法的困难是计算量难以估计，有使用迭代法的困难是计算量难以估计，有些方程组的迭代格式虽然收敛，但收敛速些方程组的迭代格式虽然收敛，但收敛速度慢而使计算量变得很大。度慢而使计算量变得很大。超松驰迭代法是一种线性加速方法超松驰迭代法是一种线性加速方法。超松驰迭代法是高斯超松驰迭代法是高