《最新安徽建筑工程学院计算结构力学6PPT课件》由会员分享，可在线阅读，更多相关《最新安徽建筑工程学院计算结构力学6PPT课件（38页珍藏版）》请在金锄头文库上搜索。
1、安徽建筑工程学院计算结构安徽建筑工程学院计算结构力学力学6首页首页上页上页返回返回下页下页第六章第六章 线性代数方程组的求解线性代数方程组的求解结构刚度方程结构刚度方程K=P求求解解方方法法有有直直接接法法和和迭迭代代法法两两大大类类。直直接接法法从从Gauss消消元元法法衍衍生生而而来来，在在程程序设计中又称为序设计中又称为消元分解法消元分解法。首页首页上页上页返回返回下页下页首页首页上页上页返回返回下页下页首页首页上页上页返回返回下页下页首页首页上页上页返回返回下页下页首页首页上页上页返回返回下页下页首页首页上页上页返回返回下页下页首页首页上页上页返回返回下页下页对于对于n阶线性方程组：阶
2、线性方程组：第第r步消元，就是要把第步消元，就是要把第r列主元以下的元素消列主元以下的元素消为零，并同时处理荷载项，为零，并同时处理荷载项，消元公式消元公式为：为：首页首页上页上页返回返回下页下页回代过程回代过程:由由(10)式倒数第一式，得式倒数第一式，得:消元完毕，消元完毕，(8)式成为：式成为：再代入再代入(10)式倒数第二式，得式倒数第二式，得:首页首页上页上页返回返回下页下页故有故有通通常常将将(10)式式中中的的K(n)记记成成S，若若从从主主元元开开始始对对各各行行进进行行规规格格化化，即即用用主主元元除除各各行行元元素素,则则得得到到对对角角线线为为1的的单单位位上三角阵上三角
3、阵，记为，记为:。首页首页上页上页返回返回下页下页 SUBROUTINE SOLV(ZK,P,N) REAL*8 ZK(50,50),P(50),C DO 20 K=1,N-1 DO 20 I=K+1,N C=ZK(K,I)/ZK(K,K) DO 10 J=1,N 10 ZK(I,J)=ZK(I,J)-C*ZK(K,J) 20 P(I)=P(I)-C*P(K) P(N)=P(N)/ZK(N,N)Gauss消元法的程序设计为：消元法的程序设计为：首页首页上页上页返回返回下页下页 DO 40 K=1,N-1 I=N-K DO 30 J=I+1,N P(I)=P(I)-ZK(I,J)*P(J) 30
4、 CONTINUE 40 P(I)=P(I)/ZK(I,I) WRITE(*,*)THIS IS DISPLACEMENT * OF THE STRUCTURE WRITE(*,100)(P(I),I=1,N) 100 FORMAT(/,2X,6F12.8) RETURN END首页首页上页上页返回返回下页下页6-2消元分解法的矩阵表示消元分解法的矩阵表示由线性代数可知，第由线性代数可知，第r步消元就相当于步消元就相当于K(r)左乘初等阵左乘初等阵Lr-1式中式中:首页首页上页上页返回返回下页下页首页首页上页上页返回返回下页下页对于对于n阶情形阶情形:且且Lr与与Lr-1相比，仅区别于第相比，
5、仅区别于第r列列的负号！由式的负号！由式(4)可知：可知：K=LS (5)首页首页上页上页返回返回下页下页在上式中由初等阵的性质可知：在上式中由初等阵的性质可知：代入代入(5)式式:首页首页上页上页返回返回下页下页由于由于D是对角阵是对角阵,由由K的对称性可知的对称性可知:由由(8)式可得式可得: K=LDLT (10)我们把我们把(10)式称为对式称为对K的正消的正消(分解分解)，这样通，这样通过消元公式自然可把过消元公式自然可把K分解为三个矩阵的乘分解为三个矩阵的乘积，而不必再找积，而不必再找Lr做乘法了。做乘法了。(10)式亦称为式亦称为消元分解法的消元过程。现利用矩阵运算求消元分解法的
6、消元过程。现利用矩阵运算求解刚度方程解刚度方程:K=P首页首页上页上页返回返回下页下页由由(10)式式: LDLT=P (11)设设: LY=P (12)则有则有: DLT=Y (13)由由(12)式对式对P进行进行正消正消，则可得，则可得Y，这说明对这说明对P的正消与对的正消与对K的消元可以分开的消元可以分开进行，而不必进行，而不必“同时同时”进行。进行。首页首页上页上页返回返回下页下页 L=L1L2Ln-1Ln-2 (14) Y=Ln-1 -1Ln-2 -1 L2 -1L1 -1P (15)由此可知由此可知:Y=P1(1) P 2(2) Pn(n) (16)对对(13)式回代就得到式回代就
7、得到:LT=D-1 Y (17)(10)(17)式是式是消元分解法消元分解法的全过程，又称的全过程，又称为为LDLT法。法。首页首页上页上页返回返回下页下页=本本质质上上仍仍是是Gauss消消元元法法，但但消消元元法法要要同同时时处处理理K和和P，进行正消与回代。进行正消与回代。=LDLT法法可可分分别别处处理理K和和P，见见(10)式式、(12)式式，当当得得到到S即即DLT阵阵后后再再处处理理P，亦亦即即可可将将K 、 P分分开开处处理理，这这样样对对同同一一结结构构的的K只只要要进进行行一一次次分分解解，便便可可对对多多组组荷荷载载进进行行计计算，实现了算，实现了方程的再解功能方程的再解
8、功能。=程程序序存存贮贮的的是是K的的上上三三角角，消消元元后后存存贮贮的的是是D-1 与与 LT，这样方便回代，见这样方便回代，见(17)式。式。 LDLT法的法的程序设计如下：程序设计如下：LDLT法与法与Gauss消元法的比较消元法的比较首页首页上页上页返回返回下页下页 SUBROUTINE SOLV(IND,ZK,P,N) ！IND：入口参数：入口参数 REAL*8 ZK(50,50),P(50),C ！ IND0P的正消回代的正消回代 15 DO 20 K=1,N-1 DO 20 I=K+1,N C=ZK(K,I)/ZK(K,K) DO 10 J=1,N 10 ZK(I,J)=ZK(
9、I,J)-C*ZK(K,J) 20 CONTINUE GO TO 200 50 DO 60 K=1,N-1 DO 60 I=K+1,N C=ZK(K,I)/ZK(K,K) 首页首页上页上页返回返回下页下页60 P(I)=P(I)-C*P(K) P(N)=P(N)/ZK(N,N) DO 40 K=1,N-1 I=N-K DO 30 J=I+1,N P(I)=P(I)-ZK(I,J)*P(J) 30 CONTINUE 40 P(I)=P(I)/ZK(I,I) WRITE(*,*)THIS IS DISPLACEMENT * OF THE STRUCTURE WRITE(*,100)(P(I),I=
10、1,N)100 FORMAT(/,2X,6F12.8)200 RETURN END首页首页上页上页返回返回下页下页6-3带状矩阵的存贮与求解带状矩阵的存贮与求解例例:用用图图示示的的十十层层刚刚架架研研究究结结构构刚刚阵阵的的带带状状特特性性，不不考考虑虑轴轴向向变变形形，每每层层3个个变变量量，共共30个个变量。变量。首页首页上页上页返回返回下页下页首页首页上页上页返回返回下页下页根据根据相关变量相关变量的概念的概念(若若i与与j不相关不相关,则则Kij=0)可知可知:=结构刚度矩阵是相当稀疏的结构刚度矩阵是相当稀疏的,非零元素分布在主非零元素分布在主对角线附近而成带状。对角线附近而成带状。
11、=根据矩阵的对称性可只存贮上三角部分。根据矩阵的对称性可只存贮上三角部分。=图中红色粗线条轮廓线以外的零元素在消元分图中红色粗线条轮廓线以外的零元素在消元分解后仍然为零，因此不需要存贮。解后仍然为零，因此不需要存贮。在在不不考考虑虑轴轴向向变变形形情情况况下下，每每一一层层只只有有3个个结结点点位位移移未未知知量量:两两个个转转角角和和一一个个线线位位移移。考考察察半半刚刚阵阵的的每每一一行行，发发现现一一行行内内最最多多有有6个个元元素素(主主元元的的最最大大相相关关变变量量差差值值+1)，最最少少为为4个个元元素素，故故有两种存贮方式有两种存贮方式:首页首页上页上页返回返回下页下页每每行行
12、都都按按6个个元元素素存存贮贮，不不足足的的补补零零，这这样样把把总总刚刚K看看成成是是一一个个306的的矩矩阵阵，而而不不再再看看成成是是3030的方阵，可以省很多存贮单元。的方阵，可以省很多存贮单元。等带宽存贮时的寻址工作等带宽存贮时的寻址工作： (1)带宽带宽NDK的计算的计算 NDK =各单元的最大未知量编号差值各单元的最大未知量编号差值+1。 本例本例NDK=6。带宽带宽:从主元起到该行最后一个非零元素的个数。从主元起到该行最后一个非零元素的个数。1. 等带宽存贮等带宽存贮首页首页上页上页返回返回下页下页(2)寻址寻址(即求即求Kij的存贮位置的存贮位置)等带宽存贮的行号不变等带宽存
13、贮的行号不变;列列号号:新新列列号号从从主主元元开开始始为为第第一一列列，由由于于主主元元的的旧旧列列号号=主主元元的的行行号号，故故新新列列号号与与旧旧列列号号有有下下列关系：列关系： J=J-I+1 (寻址公式寻址公式) (1)这样，即可将这样，即可将首页首页上页上页返回返回下页下页2. 变带宽存贮变带宽存贮由由上上个个例例子子可可知知，等等带带宽宽存存贮贮是是按按最最大大带带宽宽来来计计算算存存贮贮量量的的, 而而根根据据消消元元分分解解法法的的要要求求，按按图图中中阶阶梯梯形形粗粗轮轮廓廓线线存存贮贮元元素素即即可可，而而因因此此浪浪费费了了一一些些存存贮贮单单元元,把把由由图图中中阶
14、阶梯梯形形粗粗轮轮廓廓线线所所勾勾勒勒的的称称为为实实际际带带宽宽(变变带带宽宽)，由由此此可可按按实实际际带带宽宽来来存存贮贮刚刚度度系系数数Kij，这这就就是是变变带带宽宽存存贮贮的的由来。由来。由由于于每每行行的的带带宽宽不不相相同同，所所以以就就不不能能在在二二维维矩矩阵阵中中存存贮贮，只只能能是是一一行行接接一一行行(或或一一列列接接一一列列)地存贮，即采用地存贮，即采用一维存贮方式一维存贮方式。本例中，若按行存贮，本例中，若按行存贮，30行存成一行；若按行存成一行；若按列存贮，列存贮，30列存成一列，列存成一列，如何寻址如何寻址?首页首页上页上页返回返回下页下页一维变带宽存贮的寻址
15、工作一维变带宽存贮的寻址工作仔仔细细观观察察变变带带宽宽存存贮贮，发发现现只只要要知知道道主主元元在在一一维维数数组组的的位位置置,那那么么还还是是能能搞搞清清副副元元地地址址码码的，即二维下标与一维下标的变换关系。的，即二维下标与一维下标的变换关系。设设主主元元的的地地址址码码(主主元元在在一一维维数数组组中中的的位位置置)存存放放在在KD(N)数数组组中中，并并确确定定是是按按行行还还是是按按列列存存贮贮，即即可可确确定定各各副副元元的的地地址址码码，这这项项工工作作称称为寻址为寻址。本例若按行存贮。本例若按行存贮(参考前图参考前图)，这时：，这时：KD(30)=1 7 12 16 22
16、27139 141 (2)可可推推出出Kij在在一一维维存存贮贮时时的的序序号号(KijKIJ的的地地址码，址码，IJ为双下标为双下标)为：为： IJ=KD(I)+J-I (3)首页首页上页上页返回返回下页下页因因为为消消元元分分解解法法是是按按行行进进行行消消元元的的，故故一一般般采采用用一一维维变变带带宽宽上上三三角角按按行行存存贮贮(LDLT法采用法采用)。这这时时只只需需141+30=171个个存存贮贮单单元元，对对于于大大型型复复杂杂结结构构来来说说，节节约约的的存存贮贮空空间间则则非非常常可可观。观。KD数数组组的的程程序序设设计计见见讲讲义义所所附附结结构构静静动动力力综综合合分
17、分析析程程序序PSTDY中中的的SUB.QKD(P106)，或或参考教材参考教材P196。首页首页上页上页返回返回下页下页3. 带状矩阵的解法带状矩阵的解法(1)形形成成一一维维变变带带宽宽上上三三角角按按行行存存贮贮的的结结构构刚刚度度矩矩阵阵ZK(NZY)，NZY:存存贮贮的的刚刚度度系系数数的的总总数数。其框图设计同其框图设计同SUB. KJX1，只不过增加，只不过增加:调用调用KD数组；数组；由由KD数组寻址。数组寻址。程程序序见见后后附附结结构构静静动动力力综综合合分分析析程程序序PSTDY中中子子程程序序SUB.KJX(P110)，或或参参考考教教材材P199。首页首页上页上页返回
18、返回下页下页 (2)消元分解法求解的子程序消元分解法求解的子程序程序设计思想程序设计思想:消消元元分分解解使使用用公公式式(6-1-9)，但但将将K的的正正消消与与P的正消分别进行的正消分别进行;K消消元元分分解解后后进进行行规规格格化化处处理理, KSS， S= DLT，程序保存，程序保存D-1LT，方便回代；，方便回代；按按公公式式(6-1-13)回回代代，求求出出的的存存贮贮在在P中中。程程序序设设计计见见后后附附结结构构静静动动力力综综合合分分析析程程序序PSTDY中中子子程程序序SUB.NXFJ，其其中中IND为为入入口口参数参数: IND0 进进行行消消元元分分解解；IND0 进进
19、行行回回代代求求解。解。首页首页上页上页返回返回下页下页6-4 解的误差分析解的误差分析目的目的：了解误差，尽量避免或减小误差。：了解误差，尽量避免或减小误差。误差误差：和真解的差距。：和真解的差距。引起误差的原因很多，很复杂，大致可分为：引起误差的原因很多，很复杂，大致可分为：1. 由于结构的计算模型所引起的误差由于结构的计算模型所引起的误差这种误差不是改进算法所能解决，需在确定这种误差不是改进算法所能解决，需在确定结构计算模型时考虑，结构工程师应予注意并做结构计算模型时考虑，结构工程师应予注意并做到胸中有数。到胸中有数。2. 计算误差计算误差截截断断误误差差：计计算算机机有有效效数数字字，
20、单单精精度度7位位，双双精精度度14位。位。舍入误差：运算截断时的四舍五入。舍入误差：运算截断时的四舍五入。首页首页上页上页返回返回下页下页3.若若K和和P有有微微小小的的变变化化便便会会引引起起有有很很大大的的变变化化，或或K有有很很大大的的变变化化，就就称称K为为病病态态矩矩阵阵或或不不稳稳定定矩矩阵。阵。例例1: 这里这里则会解出：则会解出：方方程程右右端端项项变变化化很很小小，而而解解的的变变化化却却很很大大，称称这这类方程是类方程是病态方程病态方程，或称其系数矩阵为病态矩阵。，或称其系数矩阵为病态矩阵。首页首页上页上页返回返回下页下页对称矩阵的病态程度可用矩阵的对称矩阵的病态程度可用
21、矩阵的条件数条件数来描述：来描述： 分别是分别是K的特征值的模的最大的特征值的模的最大与最小值。与最小值。最最小小的的条条件件数数是是1，条条件件数数越越大大，矩矩阵阵病病态程度越厉害。对于非对称矩阵，条件数为：态程度越厉害。对于非对称矩阵，条件数为：分别是分别是K的最大特征值与最小特征值。的最大特征值与最小特征值。首页首页上页上页返回返回下页下页运算误差对条件数的计算有重大影响。运算误差对条件数的计算有重大影响。推荐：尽可能采用双精度变量和数组。推荐：尽可能采用双精度变量和数组。避免刚度系数相互变化很大，方法：主从关系避免刚度系数相互变化很大，方法：主从关系减小截断误差的影响，尽量减小截断误差的影响，尽量采用双精度采用双精度对对称称正正定定且且主主角角占占优优的的矩矩阵阵误误差差小小，当当主主角角不不占优时在程序设计中要重新遴选主元占优时在程序设计中要重新遴选主元带宽越窄，误差越小，故应优化结点编号带宽越窄，误差越小，故应优化结点编号为为避避免免刚刚度度系系数数相相差差太太大大，最最好好采采用用无无量量纲纲的的刚度系数。如对于梁，刚度系数的量纲为刚度系数。如对于梁，刚度系数的量纲为结束语结束语谢谢大家聆听！谢谢大家聆听！38
