矩阵的求解(完整版)实用资料.doc

《矩阵的求解(完整版)实用资料.doc》由会员分享，可在线阅读，更多相关《矩阵的求解(完整版)实用资料.doc（19页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、矩阵的求解（完整版）实用资料(可以直接使用，可编辑 完整版实用资料，欢迎下载）矩阵方程的求解问题白秀琴（平顶山工业职业技术学院，基础部，河南 平顶山 467001）摘要：主要考察了矩阵方程的求解问题，给出了一般矩阵方程当系数矩阵满足不同条件时的两种求解方法。关键词：矩阵；矩阵的逆；矩阵方程矩阵是线性代数中的最重要的部分，它贯穿于线性代数的始终，可以说线性代数就是矩阵的代数，矩阵是处理高等数学很多问题的有力工具。矩阵方程是矩阵运算的一部分，这里我们主要讨论如何求解矩阵方程的问题。掌握简单的矩阵方程的求法，对于求解复杂的矩阵方程有很大帮助。 简单的矩阵方程有三种形式：如果这里的、都是可逆矩阵，则求

2、解时需要找出矩阵的逆，注意左乘和右乘的区别。它们的解分别为：例如，求解方程先考察是否可逆，如果可逆时，方程两边同时左乘，得即这里要注意只能左乘不能右乘，因为矩阵的乘法不满足交换律。同样，对于方程只能右乘，得即而对于方程只能是左乘而右乘，得即看下面解矩阵方程例题：例1： 解：先求出，则则 例2： 解：先求出，则则例3： 解：先求出，则则 例4：解矩阵方程其中，是三阶单位方阵。解：移项，将矩阵方程化为标准形式：由于可逆，两边同时左乘，得注：如果按计算，需要先求，再求，最后相乘，计算量大且易出错。因此应先尽量化简矩阵方程，再计算求解。当矩阵方程中的、不是方阵或者是不可逆的方阵时，前面的方法就不能用了

3、。这时，我们需要用待定元素法来求矩阵方程。设未知矩阵的元素为，即，然后由所给的矩阵方程列出所满足的线性方程组，通过解线性方程组求出所有元素，从而得到所求矩阵。例5：解矩阵方程 解：利用元素法，先确定的行数等于左边矩阵的行数的列数等于积矩阵的列数，则是的矩阵。 设，则 即，于是得方程组解得，所以，其中为任意实数。例6：解矩阵方程其中，由于，所以是不可逆矩阵，需要用元素法求解。设则，即比较第一列元素得，解得同样，比较第二、三列元素可得对应方程组，分别解得，所以可得，其中是任意实数。总之，对于矩阵方程，当系数矩阵是方阵时，先判断是否可逆。如果可逆，则可以利用左乘或右乘逆矩阵的方法求未知矩阵，如果方阵

4、不可逆或是系数矩阵不是方阵，则需要用待定元素法通过解方程确定未知矩阵。参考文献：1赵树塬。线性代数M。北京：中国人民大学出版社，19972李君文，线性代数理论与解题方法M。长沙 ：湖南大学出版社，2002有限元刚度矩阵的压缩存贮组集及快速求解姚松田红旗1中南大学轨道交通安全教育部重点实验室,湖南长沙,410075摘要:基于“细胞元”索引存贮方案提出了一种仅组集有限元刚度矩阵中“非零元素”的方法,该方法最突出的特点是计算所需内存空间与有限元网格节点和单元的编号模式无关,适于进行“自适应网格细化”有限元分析。针对刚度矩阵的“一维压缩存贮”格式,对稀疏矩阵直接解法和预处理共轭梯度法进行了探讨,并编制

5、了相应的计算机程序对某地铁车辆有限元模型进行了分析,计算结果与ANSYS5.7的计算结果一致,相差不超过2%,说明提出的存贮方案和求解方法是正确可靠的。关键词:细胞元;组集; 非零元素;稀疏直接求解法;预处理共轭梯度法1 概 述有限元方法自1956年首次应用于飞机结构分析以来,由于其具有几何适应性强、易于处理非线性等优点,已成为科学和工程计算领域应用最为广泛的数值方法之一1。在结构有限元分析中,对求解域进行空间离散以后,经常产生大型或巨型的稀疏线性方程组AX=F,对称正定的系数矩阵A通常称为结构总体刚度矩阵,其大小为N DN D,N D为结构的总自由度数,刚度矩阵A所需的存贮空间随结构规模增大

6、而增加极快,以至于限制了求解规模。目前,刚度矩阵的存贮方案2-4多采用二维等带宽存贮或Profile/Envelope结构(即一维变带宽。采用二维等带宽存贮时将有可能由于局部带宽过大而使整体刚度系数矩阵的存贮量大大增加;一维变带宽存贮虽然会比等带宽存贮节省内存空间,但是对于Skyline下的零元素还是要加以存贮,且存贮的零元素数量取决于单元和节点的编号模式。毫无疑问,仅对整体刚度矩阵A中的非零元素进行存贮是最节省存贮空间的方法5,由于仅对其中的非零元素进行操作,从而提高了计算效率。另外该存贮格式具有如下特点:无论是对存贮还是方程的求解来说,所需要的内存空间与有限元网格划分时节点和单元的编号模式

7、无关。因此该方法高度适合“h-自适应网格细化”有限元分析,由于没有必要对新产生的节点和单元进行重新编号从而降低了计算要求。本文提出了基于“一维压缩存贮”的整体刚度矩阵组集方法,并针对刚度矩阵的压缩存贮格式,讨论了大型有限元方程组的快速求解方法。2 结构稀疏刚度矩阵的存贮对于某地铁车辆的车体有限元模型,结构总自由度N D为266862。由于刚度矩阵的对称正定性,仅需存贮整体刚度矩阵的下三角(或上三角矩阵。整体刚度矩阵的最大半带宽为263982,表1列出了采用不同存贮格式时所需要存贮的元素总数以及零元素的个数:1教育部博士点基金(20 533007项目资助表1 采用不同存贮格式时所需存贮元素的数目

8、存贮格式零元素个数存贮元素总数一维变带宽635136972 640022307本文方法0 4885335 由上可知,无论是二维等带宽还是一维变带宽格式,都存贮了大量多余的零元素,而本文提出的一维压缩存贮格式将所有非零元素压缩成一维向量存贮,最大可能地节省了存贮空间。同时采用辅助向量来描述非零元素在结构矩阵A中的位置等信息6。在本文中采取三个向量来存贮整体刚度矩阵。分别是:Values:长度为N Z的双精度向量,存贮稀疏矩阵下三角矩阵中的N Z个非零元素,元素按“所在行”的顺序依次存贮,每一行以对角元素结束。Column,长度为N Z的整数向量,存贮对应Values向量中每个元素所在的列号。Ro

9、wIndex,长度为(N D+1的整数向量,存贮每一行起始元素的指针,根据RowIndex(I+1和RowIndex(I的差值可以计算出第I行非零元素的个数。3 整体刚度矩阵的组集记节点自由度为D OF,通常D OF1;考虑到结构有限元分析中节点具有多个自由度的特征(注意:每个节点的自由度是连续编号的,我们将整体刚度矩阵按照节点(而不是节点自由度形成分块子矩阵7-10,在本文中称该子矩阵为“细胞元”。根据有限元法的基本原理,如果编号为I,J的两节点不共存于同一个单元中,那么在整体刚度矩阵中对应于(I,J位置的“细胞元”(大小为D OFD OF的子矩阵必然为零矩阵。它们占据整体刚度矩阵中位置如下

10、:D OF(I-1+M, D OF(J-1+N,其中1M, ND OF。当且仅当I,J两节点共存于同一单元中时,整体刚度矩阵中对应于(I,J位置的“细胞元”才不为零。由上述原理,根据节点-节点之间的拓扑连接关系就可大致推断出整体刚度矩阵中非零元素的分布情况;由于仅仅需要存贮对称整体刚度矩阵的下三角矩阵(或上三角矩阵就可以了。那么只要对编号为I的节点,找出所有与I节点共处于同一单元且编号比I小的相关节点即可。图1中列出了地铁有限元模型中相关节点数最大的节点39022; 图1 39022节点的相关节点图可以看到在39022的14个相关节点中,有12个节点的编号小于39022,说明在下三角矩阵中,与

11、39022节点对应需要存贮12个“细胞元”,同时下三角矩阵与39022节点对应的对角线位置的“细胞元”必然不为零,在这个位置需要存贮的“细胞元”大小为(D OFD OF/2。在本文中,由单元刚度矩阵组集成整体刚度矩阵按照“细胞元”索引存贮方案进行,通过确定“细胞元”在一维数组Values中起始位置,就可将“细胞元”中的各个元素依次存放到一维数组中。由于非对角位置的“细胞元”所占据的存贮空间与对角位置“细胞元”所占据空间不一样,在确定“细胞元”的起始位置时,考虑先按节点编号顺序存放所有非对角位置的“细胞元”,最后再按节点编号顺序存放对角位置的“细胞元”。程序实现如下:ELEM (N E, N P

12、E:逆时针方向存放的单元I E的节点编号,N E为单元总数,N PE为每单元节点数;TEMP(N PE:临时数组,记录单元I E的节点编号信息;CONNECT(N N,M AX:节点I N的相关节点编号,N N为节点总数,M AX为指定变量;NUM(N N:节点I N的相关节点总数;NONZERO(N N:记录与节点I N对应的第一个非对角位置的“细胞元”在一维数组Values中的起始位置;DIAG(N N:记录节点I N的对角“细胞元”在一维数组Values中的起始位置;(1单元号I E从1N E进行循环,TEMP(J=ELEM(I E,J,(J=1N PE。对TEMP(N PE 按从小到大

13、进行排序,排序后第J(J=2N PE个节点TEMP(J在I E单元中其相关节点数目为(J-1,将TEMP(J的相关节点编号置入数组CONNECT(TEMP(J,M AX中;从图1中可以看出,当节点TEMP(J与其相关节点是多个单元公共边上的节点时, CONNECT(TEMP(J,M AX数组中相关节点编号就会重复出现;(2节点号I N从1N N进行循环,对CONNECT(I N,M AX从小到大进行排序,如果相邻的节点编号相同,剔除其中的一个编号;得到没有重复节点编号的CONNECT数组和没有重复计数的NUM数组。(3计算与节点I N对应的第一个非对角位置的“细胞元”在一维数组Values中的

14、起始位置NONZERO(I N和节点I N的对角“细胞元”在一维数组Values中的起始位置DIAG(N N。(4单元号I E从1N E进行循环,求解单元刚度矩阵,其大小为(N PEN PE个“细胞元”,判断在单元I E中的(N PED OF个自由度中是否有被约束的自由度,若有,则采用所在行、所在列“置0”,对角元素“置1”的方法对单元刚度矩阵进行修正。“细胞元”的组集方法如下:DO J1=1, N PEDO J2=1, N PE若 TEMP(J1TEMP(J2确定节点TEMP(J2在节点TEMP(J1的相关节点序列中的相对位置,根据相对位置和NONZERO(N N将(J1,J2位置“细胞元”

15、组集到一维数组中若 TEMP( J1=TEMP(J2根据DIAG(N N 将(J1,J1位置“细胞元”组集到一维数组中ENDDOENDDO(5对整体刚度矩阵组集完成以后,由于施加了约束等原因,“细胞元”中的不少元素为零元素,通过一次循环将其中的零元素剔除,并将所有非零元素按照“行压缩存贮”在三个一维数组中,Values(N Z,Column(N Z,Rowindex(N D+1。4 有限元方程组的快速求解大型有限元方程组的求解是结构分析中耗时最多的环节11,对结构整体刚度矩阵的压缩存贮最大程度节约了存贮空间,使得在普通个人计算机上对更大更复杂的结构进行有限元分析成为可能,同时必须发展与该格式相

16、适应的快速求解方法,在本文中介绍了稀疏直接求解法和预条件共轭梯度迭代求解法,并利用FORTRAN95编制了相应的程序,结合某地铁车体有限元模型,对这两种解法进行了比较;4.1 稀疏矩阵直接求解法(Sparse Direct Solver是对规模不太大结构的较为有效的一种有限元直接求解方法,它充分利用了刚度矩阵的稀疏性。运用以图论为基础的稀疏矩阵重排序技术,改变其稀疏结构(相当于换行换列。本文中利用Cuthill-McKee排序算法11-12,图2列出了经过排序以后的稀疏刚度矩阵非零元素的分布图,阴影位置表示非零元素。这样的结构使得三角分解A=LL T 的过程中“填充元”尽可能少;所需存贮量及运

17、算量大大减少,较好地实现了求解速度和节约存贮空间的结合。 图2 经过Cuthill-McKee排序的刚度矩阵其求解过程分为四步:(1 对节点号进行编排,使得总体刚度矩阵A三角分解以后,下三角矩阵中的非零元素尽可能少。(2 对总体刚度矩阵A进行符号分解(Symbol Factorization以确定L的稀疏结构。(3 对总体刚度矩阵A进行数值分解(Numeric Factorization以求得L。(4 回代求解。4.2 预处理共轭梯度迭代法迭代求解法由给定初值通过若干迭代步骤而求得满足一定精度要求的近似解,由于无需进行三角分解,在计算过程中仅需进行矩阵与向量的乘积运算,所需的内存空间少,对于超

18、大规模稀疏线性方程组最为有效。经典的迭代法(Jacobi,Gauss-Seidel,SSOR收敛速度过慢或不稳定,故实际工程中应用较少。对称正定方程组的迭代求解目前最常用的是预处理共轭梯度法13-15,通过对矩阵进行预处理,将原先的方程组转换成一个等价的,但是系数矩阵条件数更小、更易于收敛的方程组,然后再用共轭梯度法(CG去求解。本文中采用基于稀疏矩阵的“不完全分解”的方法ICCG(0构造预优矩阵,为了提高不完全分解的稳定性,本文将非对角元素按一定比例缩小以保证整体刚度矩阵A“对角占优”。A=LDL T+R,下三角矩阵L的稀疏结构与相同A。4.3 数值算例:本文编制了预处理共轭梯度法和稀疏矩阵

19、直接解法的FORTRAN程序,对某地铁车辆承受500KN纵向拉力的工况进行了求解,并与ANSYS5.7的计算结果进行了对比,在表2中列出了采用三种方法求得结构的U X,U Y,U Z的最大值,和所耗费的计算时间,所有计算均在CPU2.8GHz,内存512M的Dell机器上进行:表2 不同方法得到的计算结果和计算时间时间/s U X/mm U Y/mm U Z/mmAnsys波前法450 -1.42687.9828 -8.1722稀疏直接法423 -1.41177.7775 -8.1607预优迭代法2104-1.43807.9617 -8.1840通过对结果的比较,无论是采用哪种方法,求出的位移

20、最大值非常接近,相差不超过2 %。由于稀疏直接法和预优迭代法均是基于本文中提出的“仅存贮非零元素”的组集法进行的,通过这个数值算例也说明该组集法是正确的。从三个方法耗费的时间来看,稀疏直接法与ansys5.7所用的波前法基本相当,预优迭代法的时间要明显长一些,这是由本文所采取的预优矩阵的特性决定的,如果能获得性能更好的预优矩阵,将会大大加快方程组迭代收敛的进程。但是对于超过机器内存限制的物理问题,由于迭代法最节省空间的优势,无疑是求解这类问题的有效算法。5 结 论基于“细胞元”索引存贮方案,本文提出了一种“仅存贮非零元素”的结构刚度矩阵的组集方法,该方法所需内存空间与有限元网格节点和单元的编号

21、模式无关,适于进行“自适应网格细化”有限元分析。本文针对刚度矩阵的“一维压缩存贮”格式,对稀疏矩阵直接解法和预处理共轭梯度法进行了探讨,采用这两种方法对某地铁车辆有限元模型进行了分析,通过与ANSYS5.7的计算结果比较,说明提出的存贮方案和求解方法是正确可靠的。参考文献1 Bathe KJ. Finite Element Procedures M. Engelwood Cliffs, NJ: Prentice Hall, 19962 王勖成. 有限单元法 M. 北京: 清华大学出版社, 2003Wang Xu-Cheng. Finite Element Method M. Beijing:

22、Tsinghua University Press, 20033 Zienkiewicz O C. The Finite Element M. Boston: Oxford University Press, 20004 I.M. Smith, D.V. Griffiths. Programming the Finite Element Method M. New Jersey: John Wiley & Sons Inc, 1998 5 J.Rajasankar, Nagesh R. Iyer. A procedure to assemble only non-zero coefficien

23、ts of global matrix in finite element analysis J. Computer & structures, 2000, 77(3: 595-599 6 Laura C. Dutto. Effect of the Storage Format of Sparse Linear Systems on Parallel CFD Computations J. Computer Methods Applied Mechanics and Engineering, 2000, 188(3: 441-453 7 Pu Chen, Dong Zheng. High Pe

24、rformance Sparse Static Solver in Finite Element Analyses with Loop-Unrolling J. Advances in Engineering Software, 2003, 34 (1: 203-215 8 Chen P, Runesha H. Sparse Algorithms for Indefinite Systems of Linear Equations J. Computer Mechanics, 2000, 25(1: 3342. 9 Felippa CA. Solution of Linear Equation

25、s with Skylined-Stored Symmetric Matrix J. Computer & Structures, 1975, 5(1: 1329 10 Duff IS, Erisman AM, Reid JK. Direct Methods for Sparse Matrices M. Oxford: Clarendon Press, 1989 11 Weinberg DJ. A Performance Assessment of NE/Nastrans New Sparse Direct and Iterative Solvers J. Advanced Engineeri

26、ng Software, 2000, 31(4: 547554 12 Ricardo D. Arantes. Code Orderings for the Solution of Sparse Positive Definite Linear Systems J. Applied Numerical Mathematics, 1999, 30 (2: 137140 13 Wen Yang Lin. Minimum Communication Cost Reordering for Parallel Sparse Cholesky Factorization J. Parallel Comput

27、ing, 1999, 25 (10: 943-967 14 J.W.-H. Liu. Computational Models and Task Scheduling for Parallel Sparse Cholesky Factorizations J. Parallel Computing, 1986,3(3: 327-343 15 Eurpides Montagne. An Optimal Storage Format for Sparse Matrices J. Information Processing Letters, 2004, 90 (1: 8792 A Compress

28、ed Storage Assembly Procedure of Global Stiffness Matrix in FEA And Fast Solution Yao song Tian hongqi Truck and Transportation Safety Key Laboratory of Ministry of Education, Central South University, Changsha, Hunan , 410075 Abstract A procedure to assemble only non-zero coefficients of a global s

29、tiffness matrix of finite element analysis in a one-dimensional array is proposed, which utilizes cell element index storage scheme. One of the important characteristics of the scheme is that the storage and solution requirements are independent of the numbering pattern of nodes and elements in fini

30、te element mesh, which is suitable for adaptively refined finite element analysis. For one-dimension compressed format of sparse structure matrix, two solution methods for large sparse matrix equations are discussed in this paper, including sparse direct solve method and preconditioned conjugate gra

31、dient method. The procedure has been implemented in a computer program. At last a finite element model of a subway car-body is carried out using the program. The results are in accordance with the analysis results obtained by commercial FEA codes ANSYS5.7. The difference doesnt exceed 2 percent. The

32、 fact shows that the storage scheme and solution method proposed in this paper are correct and reliable. Keywords: cell element；assemble；non-zero coefficients；sparse direct solver；preconditioned conjugate gradient 矩阵方程的求解问题白秀琴 杨宝玉（平顶山工业职业技术学院，基础部，河南 平顶山 467001）摘要：主要考察了矩阵方程的求解问题，给出了一般矩阵方程当系数矩阵满足不同条件时

33、的两种求解方法。关键词：矩阵；矩阵的逆；矩阵方程The square distance of matrix solves the problemBAI Xiu-qin,Yang Bao-yu(Pingdingshang Industrial College Of Technology,Pingdingshan,Henan,467001)Abstract: Mainly investigated the square distance of matrix to solve the problem, give the two kinds of general square distance of

34、matrix solve method. When the coefficient matrix to satisfy the different conditionsKey words: Matrix；Matrix of negative；Square distance of matrix矩阵是线性代数中的最重要的部分，它贯穿于线性代数的始终，可以说线性代数就是矩阵的代数，矩阵是处理高等数学很多问题的有力工具。矩阵方程是矩阵运算的一部分，这里我们主要讨论如何求解矩阵方程的问题。掌握简单的矩阵方程的求法，对于求解复杂的矩阵方程有很大帮助。 简单的矩阵方程有三种形式：如果这里的、都是可逆矩阵，则

35、求解时需要找出矩阵的逆，注意左乘和右乘的区别。它们的解分别为：例如，求解方程先考察是否可逆，如果可逆时，方程两边同时左乘，得即这里要注意只能左乘不能右乘，因为矩阵的乘法不满足交换律。同样，对于方程只能右乘，得即而对于方程只能是左乘而右乘，得即看下面解矩阵方程例题：例1： 解：先求出，则则 例2： 解：先求出，则则例3： 解：先求出，则则 例4：解矩阵方程其中，是三阶单位方阵。解：移项，将矩阵方程化为标准形式：由于可逆，两边同时左乘，得注：如果按计算，需要先求，再求，最后相乘，计算量大且易出错。因此应先尽量化简矩阵方程，再计算求解。当矩阵方程中的、不是方阵或者是不可逆的方阵时，前面的方法就不能用

36、了。这时，我们需要用待定元素法来求矩阵方程。设未知矩阵的元素为，即，然后由所给的矩阵方程列出所满足的线性方程组，通过解线性方程组求出所有元素，从而得到所求矩阵。例5：解矩阵方程 解：利用元素法，先确定的行数等于左边矩阵的行数的列数等于积矩阵的列数，则是的矩阵。 设，则 即，于是得方程组解得，所以，其中为任意实数。例6：解矩阵方程其中，由于，所以是不可逆矩阵，需要用元素法求解。设则，即比较第一列元素得，解得同样，比较第二、三列元素可得对应方程组，分别解得，所以可得，其中是任意实数。总之，对于矩阵方程，当系数矩阵是方阵时，先判断是否可逆。如果可逆，则可以利用左乘或右乘逆矩阵的方法求未知矩阵，如果方阵不可逆或是系数矩阵不是方阵，则需要用待定元素法通过解方程确定未知矩阵。作者简介：白秀琴（1965-），女，河南淇县人，平顶山工业职业技术学院副教授，主要从事基础数学的研究参考文献：1赵树塬。线性代数M。北京：中国人民大学出版社，19972李君文，线性代数理论与解题方法M。长沙 ：湖南大学出版社，2002