安徽建设工程学院计算结构力学6.ppt

第六章 线性代数方程组的求解 　　因为消元分解法是按行进行消元的，故一般采用一维变带宽上三角按行存贮(LDLT法采用)。 　　这时只需141+30=171个存贮单元，对于大型复杂结构来说，节约的存贮空间则非常可观。 　　KD数组的程序设计见讲义所附结构静动力综合分析程序PSTDY中的SUB.QKD(P106)，或参考教材P196。 3. 带状矩阵的解法 (1)形成一维变带宽上三角按行存贮的结构刚度矩阵ZK(NZY)，NZY:存贮的刚度系数的总数。其框图设计同SUB. KJX1，只不过增加: 　　①调用KD数组； 　　②由KD数组寻址。 　　程序见后附结构静动力综合分析程序PSTDY中子程序SUB.KJX(P110)，或参考教材P199。 * 首页 上页 返回 下页 * 计算结构力学 第六章 线性代数方程组的求解 结构刚度方程 　　　　　　[K]{Δ}={P} 　　求解方法有直接法和迭代法两大类。直接法从Gauss消元法衍生而来，在程序设计中又称为消元分解法。 解:由5—3中例可知刚度方程为 Gauss消元法解线性方程组的步骤： 第一步： (令[K](1)=[K]) 　　除K11外，使第1列的元素为零，即(1)式中第二式+第一式 6-1 Gauss消元法 例:消元法求解 结构刚度方程 　　由于 ，故第三行类似的运算就不必进行，式中上标(1)、(2)表示消元序号。 经过第一步消元后得： 将第1列主元以下的元素变为0 第二步，除上三角外，使第二列元素均变为 零，即在(2)中的第三式+第二式 经过第二步消元后得： 　　将第2列主元以下的元素变为0，当该矩阵的下三角为零时，便完成了Gauss消元法的消元过程。 注意： 　　每步消元后的[K]仍是正定对称的。 由(3)式的最后一式定得: 将(4)代入(3)的第二式得： 将(4)式与(5)式代入(3)的第一式得： 以上的运算称为Gauss消元法的回代过程。 对于n阶线性方程组： 第r步消元，就是要把第r列主元以下的元素消为零，并同时处理荷载项，消元公式为： 回代过程:由(10)式倒数第一式，得: 消元完毕，(8)式成为： 再代入(10)式倒数第二式，得: 故有 　　通常将(10)式中的[K(n)]记成[S]，若从主元开始对各行进行规格化，即用主元除各行元素,则得到对角线为1的单位上三角阵，记为: 。 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 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相比，仅区别于第r列的负号！由式(4)可知： [K]=[L][S] 　　　　(5) 在上式中由初等阵的性质可知： 代入(5)式: 由于[D]是对角阵,由[K]的对称性可知: 由(8)式可得: [K]=[L][D][L]T (10) 我们把(10)式称为对[K]的正消(分解)，这样通过消元公式自然可把[K]分解为三个矩阵的乘积，而不必再找[Lr]做乘法了。(10)式亦称为消元分解法的消元过程。现利用矩阵运算求解刚度方程: [K]{Δ}={P} 由(10)式: [L][D][L]T{Δ}={P} (11) 设: [L][Y]={P} (12) 则有: