单晶塑性abaqus-umat 自学—续

6 篇文章 39 订阅
3 篇文章 33 订阅


本文是针对黄老师单晶umat程序中的一些细节问题。都是本人一下没看明白,通过思考和查资料弄明白后,写在这里。

1.刚度矩阵的旋转

C-----  Rotation matrix: ROTD to transform local elastic matrix DLOCAL 
C     to global elastic matrix D
C
      DO J=1,3
         J1=1+J/3
         J2=2+J/2

         DO I=1,3
            I1=1+I/3
            I2=2+I/2

            ROTD(I,J)=ROTATE(I,J)**2
            ROTD(I,J+3)=2.*ROTATE(I,J1)*ROTATE(I,J2)
            ROTD(I+3,J)=ROTATE(I1,J)*ROTATE(I2,J)
            ROTD(I+3,J+3)=ROTATE(I1,J1)*ROTATE(I2,J2)+
     2                    ROTATE(I1,J2)*ROTATE(I2,J1)

         END DO
      END DO

C-----  Elastic matrix in global system: D
C     {D} = {ROTD} * {DLOCAL} * {ROTD}transpose

旋转矩阵设为 R \bm{R} R ,向量旋转:
u ⃗ N e w = R ⋅ u ⃗ O l d u ⃗ i N e w = R i j u ⃗ j O l d \vec{u}^{New}=\bm{R}\cdot\vec{u}^{Old}\quad\vec{u}_{i}^{New}={R_{ij}}\vec{u}_{j}^{Old} u New=Ru Oldu iNew=Riju jOld

二阶张量旋转:
σ N e w = R ⋅ σ O l d ⋅ R T σ i j N e w = R i k σ k l O l d R j l \bm{\sigma}^{New}=\bm{R}\cdot\bm{\sigma}^{Old}\cdot \bm{R}^T\quad\bm{\sigma}_{ij}^{New}={R_{ik}}\bm{\sigma}_{kl}^{Old} R_{jl} σNew=RσOldRTσijNew=RikσklOldRjl

改写为: σ N e w = R : σ σ i j N e w = R i j k l σ k l O l d \bm{\sigma}^{New}=\mathbb{R}:\bm{\sigma}\quad \bm{\sigma}_{ij}^{New}=\mathbb{R}_{ijkl}\bm{\sigma}_{kl}^{Old} σNew=R:σσijNew=RijklσklOld
其中: R i j k l = R i k R j l \mathbb{R}_{ijkl}=R_{ik}R_{jl} Rijkl=RikRjl

对于四阶刚度张量:
D N e w = R : D O l d : R T D i j k l N e w = R i j n m D n m α β O l d R α β k l \bm{D}^{New}=\mathbb{R}:\bm{D}^{Old}:\mathbb{R}^T\quad \bm{D}_{ijkl}^{New}=\mathbb{R}_{ijnm}\bm{D}_{nm\alpha\beta}^{Old} \mathbb{R}_{\alpha\beta kl} DNew=R:DOld:RTDijklNew=RijnmDnmαβOldRαβkl

程序中ROTATE R R RROTD R \mathbb{R} R的voigt展开:

VUMAT的应力分量顺序为 { σ 11   σ 22   σ 33   σ 12   σ 23   σ 31 } \{\sigma_{11}\,\sigma_{22}\,\sigma_{33}\, \sigma_{12}\,\sigma_{23}\,\sigma_{31}\} {σ11σ22σ33σ12σ23σ31}
R i j k l = [ k l 11 22 33 12 23 31 21 32 13 i j 11 R 11 R 11 R 12 R 12 R 13 R 13 R 11 R 12 R 12 R 13 R 13 R 11 R 12 R 11 R 13 R 12 R 11 R 13 22 R 21 R 21 R 22 R 22 R 23 R 23 R 21 R 22 R 22 R 23 R 23 R 21 R 22 R 21 R 23 R 22 R 21 R 23 33 R 31 R 31 R 32 R 32 R 33 R 33 R 31 R 32 R 32 R 33 R 33 R 31 R 32 R 31 R 33 R 32 R 31 R 33 12 R 11 R 21 R 12 R 22 R 13 R 23 R 11 R 22 R 12 R 23 R 13 R 21 R 12 R 21 R 13 R 21 R 11 R 23 23 R 21 R 31 R 22 R 32 R 23 R 33 R 21 R 32 R 22 R 33 R 23 R 31 R 22 R 31 R 23 R 32 R 21 R 33 31 R 31 R 11 R 32 R 12 R 33 R 13 R 31 R 12 R 32 R 13 R 33 R 11 R 32 R 11 R 33 R 12 R 31 R 13 21 R 21 R 11 R 22 R 12 R 23 R 13 R 21 R 12 R 22 R 13 R 23 R 11 R 22 R 11 R 23 R 12 R 21 R 13 32 R 31 R 21 R 32 R 22 R 33 R 23 R 31 R 22 R 32 R 23 R 33 R 21 R 32 R 21 R 33 R 22 R 31 R 23 13 R 11 R 31 R 12 R 32 R 13 R 33 R 11 R 32 R 12 R 33 R 13 R 31 R 12 R 31 R 13 R 32 R 11 R 33 ] \mathbb{R}_{ijkl}=\left[\begin{array}{cc:ccc:ccc:} &kl&11&22&33&12&23&31&21&32&13\\ ij\\\hdashline 11&&R_{11}R_{11}&R_{12}R_{12}&R_{13}R_{13}&R_{11}R_{12}&R_{12}R_{13}&R_{13}R_{11}&R_{12}R_{11}&R_{13}R_{12}&R_{11}R_{13}\\ 22&&R_{21}R_{21}&R_{22}R_{22}&R_{23}R_{23}&R_{21}R_{22}&R_{22}R_{23}&R_{23}R_{21}&R_{22}R_{21}&R_{23}R_{22}&R_{21}R_{23}\\ 33&&R_{31}R_{31}&R_{32}R_{32}&R_{33}R_{33}&R_{31}R_{32}&R_{32}R_{33}&R_{33}R_{31}&R_{32}R_{31}&R_{33}R_{32}&R_{31}R_{33}\\\hdashline 12&&R_{11}R_{21}&R_{12}R_{22}&R_{13}R_{23}&R_{11}R_{22}&R_{12}R_{23}&R_{13}R_{21}&R_{12}R_{21}&R_{13}R_{21}&R_{11}R_{23}\\ 23&&R_{21}R_{31}&R_{22}R_{32}&R_{23}R_{33}&R_{21}R_{32}&R_{22}R_{33}&R_{23}R_{31}&R_{22}R_{31}&R_{23}R_{32}&R_{21}R_{33}\\ 31&&R_{31}R_{11}&R_{32}R_{12}&R_{33}R_{13}&R_{31}R_{12}&R_{32}R_{13}&R_{33}R_{11}&R_{32}R_{11}&R_{33}R_{12}&R_{31}R_{13}\\\hdashline 21&&R_{21}R_{11}&R_{22}R_{12}&R_{23}R_{13}&R_{21}R_{12}&R_{22}R_{13}&R_{23}R_{11}&R_{22}R_{11}&R_{23}R_{12}&R_{21}R_{13}\\ 32&&R_{31}R_{21}&R_{32}R_{22}&R_{33}R_{23}&R_{31}R_{22}&R_{32}R_{23}&R_{33}R_{21}&R_{32}R_{21}&R_{33}R_{22}&R_{31}R_{23}\\ 13&&R_{11}R_{31}&R_{12}R_{32}&R_{13}R_{33}&R_{11}R_{32}&R_{12}R_{33}&R_{13}R_{31}&R_{12}R_{31}&R_{13}R_{32}&R_{11}R_{33} \end{array}\right] Rijkl=ij112233122331213213kl11R11R11R21R21R31R31R11R21R21R31R31R11R21R11R31R21R11R3122R12R12R22R22R32R32R12R22R22R32R32R12R22R12R32R22R12R3233R13R13R23R23R33R33R13R23R23R33R33R13R23R13R33R23R13R3312R11R12R21R22R31R32R11R22R21R32R31R12R21R12R31R22R11R3223R12R13R22R23R32R33R12R23R22R33R32R13R22R13R32R23R12R3331R13R11R23R21R33R31R13R21R23R31R33R11R23R11R33R21R13R3121R12R11R22R21R32R31R12R21R22R31R32R11R22R11R32R21R12R3132R13R12R23R22R33R32R13R21R23R32R33R12R23R12R33R22R13R3213R11R13R21R23R31R33R11R23R21R33R31R13R21R13R31R23R11R33

R i j k l V o i g t = [ k l 11 22 33 12 23 31 i j 11 R 11 R 11 R 12 R 12 R 13 R 13 2 R 11 R 12 2 R 12 R 13 2 R 13 R 11 22 R 21 R 21 R 22 R 22 R 23 R 23 2 R 21 R 22 2 R 22 R 23 2 R 23 R 21 33 R 31 R 31 R 32 R 32 R 33 R 33 2 R 31 R 32 2 R 32 R 33 2 R 33 R 31 12 R 11 R 21 R 12 R 22 R 13 R 23 R 11 R 22 + R 12 R 21 R 12 R 23 + R 13 R 21 R 13 R 21 + R 11 R 23 23 R 21 R 31 R 22 R 32 R 23 R 33 R 21 R 32 + R 22 R 31 R 22 R 33 + R 23 R 32 R 23 R 31 + R 21 R 33 31 R 31 R 11 R 32 R 12 R 33 R 13 R 31 R 12 + R 32 R 11 R 32 R 13 + R 33 R 12 R 33 R 11 + R 31 R 13 ] \mathbb{R}_{ijkl}^{Voigt}=\left[\begin{array}{cc:ccc:ccc} &kl&11&22&33&12&23&31\\ ij\\\hdashline 11&&R_{11}R_{11}&R_{12}R_{12}&R_{13}R_{13}&2R_{11}R_{12}&2R_{12}R_{13}&2R_{13}R_{11}\\ 22&&R_{21}R_{21}&R_{22}R_{22}&R_{23}R_{23}&2R_{21}R_{22}&2R_{22}R_{23}&2R_{23}R_{21}\\ 33&&R_{31}R_{31}&R_{32}R_{32}&R_{33}R_{33}&2R_{31}R_{32}&2R_{32}R_{33}&2R_{33}R_{31}\\\hdashline 12&&R_{11}R_{21}&R_{12}R_{22}&R_{13}R_{23}&R_{11}R_{22}+R_{12}R_{21}&R_{12}R_{23}+R_{13}R_{21}&R_{13}R_{21}+R_{11}R_{23}\\ 23&&R_{21}R_{31}&R_{22}R_{32}&R_{23}R_{33}&R_{21}R_{32}+R_{22}R_{31}&R_{22}R_{33}+R_{23}R_{32}&R_{23}R_{31}+R_{21}R_{33}\\ 31&&R_{31}R_{11}&R_{32}R_{12}&R_{33}R_{13}&R_{31}R_{12}+R_{32}R_{11}&R_{32}R_{13}+R_{33}R_{12}&R_{33}R_{11}+R_{31}R_{13} \end{array}\right] RijklVoigt=ij112233122331kl11R11R11R21R21R31R31R11R21R21R31R31R1122R12R12R22R22R32R32R12R22R22R32R32R1233R13R13R23R23R33R33R13R23R23R33R33R13122R11R122R21R222R31R32R11R22+R12R21R21R32+R22R31R31R12+R32R11232R12R132R22R232R32R33R12R23+R13R21R22R33+R23R32R32R13+R33R12312R13R112R23R212R33R31R13R21+R11R23R23R31+R21R33R33R11+R31R13

UMAT的应力分量顺序为 { σ 11   σ 22   σ 33   σ 12   σ 13   σ 23 } \{\sigma_{11}\,\sigma_{22}\,\sigma_{33}\, \sigma_{12}\,\sigma_{13}\,\sigma_{23}\} {σ11σ22σ33σ12σ13σ23}

R i j k l V o i g t = [ k l 11 22 33 12 13 23 i j 11 R 11 R 11 R 12 R 12 R 13 R 13 2 R 11 R 12 2 R 11 R 13 2 R 12 R 13 22 R 21 R 21 R 22 R 22 R 23 R 23 2 R 21 R 22 2 R 21 R 23 2 R 22 R 23 33 R 31 R 31 R 32 R 32 R 33 R 33 2 R 31 R 32 2 R 31 R 33 2 R 32 R 33 12 R 11 R 21 R 12 R 22 R 13 R 23 R 11 R 22 + R 12 R 21 R 11 R 23 + R 13 R 21 R 12 R 23 + R 13 R 22 13 R 11 R 31 R 12 R 32 R 13 R 33 R 11 R 32 + R 12 R 31 R 11 R 33 + R 13 R 31 R 12 R 33 + R 13 R 32 23 R 21 R 31 R 22 R 32 R 23 R 33 R 21 R 32 + R 22 R 31 R 21 R 33 + R 23 R 31 R 22 R 33 + R 23 R 32 ] \mathbb{R}_{ijkl}^{Voigt}=\left[\begin{array}{cc:ccc:ccc} &kl&11&22&33&12&13&23\\ ij\\\hdashline 11&&R_{11}R_{11}&R_{12}R_{12}&R_{13}R_{13}&2R_{11}R_{12}&2R_{11}R_{13}&2R_{12}R_{13}\\ 22&&R_{21}R_{21}&R_{22}R_{22}&R_{23}R_{23}&2R_{21}R_{22}&2R_{21}R_{23}&2R_{22}R_{23}\\ 33&&R_{31}R_{31}&R_{32}R_{32}&R_{33}R_{33}&2R_{31}R_{32}&2R_{31}R_{33}&2R_{32}R_{33}\\\hdashline 12&&R_{11}R_{21}&R_{12}R_{22}&R_{13}R_{23}&R_{11}R_{22}+R_{12}R_{21}&R_{11}R_{23}+R_{13}R_{21}&R_{12}R_{23}+R_{13}R_{22}\\ 13&&R_{11}R_{31}&R_{12}R_{32}&R_{13}R_{33}&R_{11}R_{32}+R_{12}R_{31}&R_{11}R_{33}+R_{13}R_{31}&R_{12}R_{33}+R_{13}R_{32}\\ 23&&R_{21}R_{31}&R_{22}R_{32}&R_{23}R_{33}&R_{21}R_{32}+R_{22}R_{31}&R_{21}R_{33}+R_{23}R_{31}&R_{22}R_{33}+R_{23}R_{32} \end{array}\right] RijklVoigt=ij112233121323kl11R11R11R21R21R31R31R11R21R11R31R21R3122R12R12R22R22R32R32R12R22R12R32R22R3233R13R13R23R23R33R33R13R23R13R33R23R33122R11R122R21R222R31R32R11R22+R12R21R11R32+R12R31R21R32+R22R31132R11R132R21R232R31R33R11R23+R13R21R11R33+R13R31R21R33+R23R31232R12R132R22R232R32R33R12R23+R13R22R12R33+R13R32R22R33+R23R32

2. LU分解

      SUBROUTINE LUDCMP (A, N, NP, INDX, D)

C-----  LU decomposition

C-----  Use single precision on cray
C
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NMAX=200, TINY=1.0E-20)
      DIMENSION A(NP,NP), INDX(N), VV(NMAX)

      D=1.
      DO I=1,N
         AAMAX=0.

         DO J=1,N
            IF (ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J))
         END DO

         IF (AAMAX.EQ.0.) PAUSE 'Singular matrix.'
         VV(I)=1./AAMAX
      END DO

      DO J=1,N
         DO I=1,J-1
            SUM=A(I,J)

            DO K=1,I-1
               SUM=SUM-A(I,K)*A(K,J)
            END DO

            A(I,J)=SUM
         END DO
         AAMAX=0.

         DO I=J,N
            SUM=A(I,J)

            DO K=1,J-1
               SUM=SUM-A(I,K)*A(K,J)
            END DO

            A(I,J)=SUM
            DUM=VV(I)*ABS(SUM)
            IF (DUM.GE.AAMAX) THEN
               IMAX=I
               AAMAX=DUM
            END IF
         END DO

         IF (J.NE.IMAX) THEN
            DO K=1,N
               DUM=A(IMAX,K)
               A(IMAX,K)=A(J,K)
               A(J,K)=DUM
            END DO

            D=-D
            VV(IMAX)=VV(J)
         END IF

         INDX(J)=IMAX
         IF (A(J,J).EQ.0.) A(J,J)=TINY
         IF (J.NE.N) THEN
            DUM=1./A(J,J)
            DO I=J+1,N
               A(I,J)=A(I,J)*DUM
            END DO
         END IF

      END DO

      RETURN
      END

黄老师子程序中的LU分解,实际是改进后的LUP分解,LUP分解中每一步是找的是这一列的最大值并通过P矩阵交换的对角线上,黄老师程序并不是找最大值,而是相对最大值,当前值要除以此行的最大值,再比较大小DUM=VV(I)*ABS(SUM) IF (DUM.GE.AAMAX) THEN,找出相对最大值,且下三角矩阵 L L L 存在变化后矩阵的下三角部分(下三角矩阵主对角线上的值均为1),上三角矩阵 U U U 存在变化后矩阵的上三角部分。置换关系P,通过INDX(N)记录
以一个4×4的矩阵举例说明(数值任意取的):
A = [ A 11 A 12 A 13 A 14 A 21 A 22 A 23 A 24 A 31 A 32 A 33 A 34 A 41 A 42 A 43 A 44 ] = [ 1 2 3 4 2 3 4 1 3 3 6 7 5 3 2 1 ] A=\left[ \begin{matrix} A_{11}&A_{12}&A_{13}&A_{14}\\ A_{21}&A_{22}&A_{23}&A_{24}\\ A_{31}&A_{32}&A_{33}&A_{34}\\ A_{41}&A_{42}&A_{43}&A_{44}\\ \end{matrix}\right]=\left[ \begin{matrix} 1&2&3&4\\2&3&4&1\\3&3&6&7\\5&3&2&1 \end{matrix}\right] A=A11A21A31A41A12A22A32A42A13A23A33A43A14A24A34A44=1235233334624171
首先,计算每一行的最大值,并求倒数记为 vv(I),具体值vv(1)=1/4,vv(2)=1/4,vv(3)=1/7,vv(1)=1/5
其次,列循环:
第一列,当J=1I=1,J-1不存在,I=J,N不存在,
找出一列的相对最大值,IMAX=4
J ≠ I M A X \rm{J\ne IMAX} J=IMAXJ.NE.IMAX) 将第IMAX行和第J行互换,交换后第IMAX行的最大值倒数为 vv(IMAX)=vv(J),于此同时INDX(1)=4,即第1行和第4行互换
A = P 1 ⋅ A = [ 0 0 0 1 0 1 0 0 0 0 1 0 1 0 0 0 ] [ 1 2 3 4 2 3 4 1 3 3 6 7 5 3 2 1 ] = [ 5 3 2 1 2 3 4 1 3 3 6 7 1 2 3 4 ] A=P_1\cdot A=\left[ \begin{matrix} 0&0&0&1\\0&1&0&0\\0&0&1&0\\1&0&0&0\\ \end{matrix}\right]\left[ \begin{matrix} 1&2&3&4\\2&3&4&1\\3&3&6&7\\5&3&2&1 \end{matrix}\right]=\left[ \begin{matrix} 5&3&2&1\\2&3&4&1\\3&3&6&7\\1&2&3&4 \end{matrix}\right] A=P1A=00010100001010001235233334624171=5231333224631174
求第1列下对角阵
A = L 1 ⋅ A = [ 5 3 2 1 2 3 4 1 3 3 6 7 1 2 3 4 ] = [ 1 0 0 0 2 / 5 1 0 0 3 / 5 0 1 0 1 / 5 0 0 1 ] [ 5 3 2 1 0 3 − 3 ( 2 / 5 ) 4 − 2 ( 2 / 5 ) 1 − ( 2 / 5 ) 0 3 − 3 ( 3 / 5 ) 6 − 2 ( 3 / 5 ) 7 − ( 3 / 5 ) 0 2 − 3 ( 1 / 5 ) 3 − 2 ( 1 / 5 ) 4 − ( 1 / 5 ) ] A=L_1\cdot A=\left[ \begin{matrix} 5&3&2&1\\2&3&4&1\\3&3&6&7\\1&2&3&4 \end{matrix}\right]=\left[ \begin{matrix} 1&0&0&0\\{2}/{5}&1&0&0\\{3}/{5}&0&1&0\\1/5&0&0&1\\ \end{matrix}\right]\left[ \begin{matrix} 5&3&2&1\\0&3-3(2/5)&4-2(2/5)&1-(2/5)\\0&3-3(3/5)&6-2(3/5)&7-(3/5)\\0&2-3(1/5)&3-2(1/5)&4-(1/5) \end{matrix}\right] A=L1A=5231333224631174=12/53/51/50100001000015000333(2/5)33(3/5)23(1/5)242(2/5)62(3/5)32(1/5)11(2/5)7(3/5)4(1/5)
第1列循环完
A = [ 5 3 2 1 2 / 5 3 4 1 3 / 5 5 6 7 1 / 5 2 3 4 ] A=\left[ \begin{matrix} 5&3&2&1\\2/5&3&4&1\\3/5&5&6&7\\1/5&2&3&4 \end{matrix}\right] A=52/53/51/5335224631174
循环继续,第2列J=2
I=1,J-1 A 12 = A 12 \quad A_{12}=A_{12} A12=A12
I=J,N
DO K=1,J-1 SUM=SUM-A(I,K)*A(K,J) END DO
A = [ 5 3 2 1 2 / 5 3 − 3 ( 2 / 5 ) 4 1 3 / 5 5 − 3 ( 3 / 5 ) 6 7 1 / 5 2 − 3 ( 1 / 5 ) 3 4 ] = [ 5 3 2 1 2 / 5 9 / 5 4 1 3 / 5 6 / 5 6 7 1 / 5 7 / 5 3 4 ] A=\left[ \begin{matrix} 5&3&2&1\\2/5&3-3(2/5)&4&1\\3/5&5-3(3/5)&6&7\\1/5&2-3(1/5)&3&4\end{matrix}\right]=\left[ \begin{matrix} 5&3&2&1\\2/5&9/5&4&1\\3/5&6/5&6&7\\1/5&7/5&3&4\end{matrix}\right] A=52/53/51/5333(2/5)53(3/5)23(1/5)24631174=52/53/51/539/56/57/524631174
IMAX=2 \quad IMAX=J \quad 于此同时INDX(2)=2 \quad 求下三角:
A = [ 5 3 2 1 2 / 5 3 − 3 ( 2 / 5 ) 4 1 3 / 5 5 − 3 ( 3 / 5 ) 6 7 1 / 5 2 − 3 ( 1 / 5 ) 3 4 ] = [ 5 3 2 1 2 / 5 9 / 5 4 1 3 / 5 6 / 9 6 7 1 / 5 7 / 9 3 4 ] A=\left[ \begin{matrix} 5&3&2&1\\2/5&3-3(2/5)&4&1\\3/5&5-3(3/5)&6&7\\1/5&2-3(1/5)&3&4\end{matrix}\right]=\left[ \begin{matrix} 5&3&2&1\\2/5&9/5&4&1\\3/5&6/9&6&7\\1/5&7/9&3&4\end{matrix}\right] A=52/53/51/5333(2/5)53(3/5)23(1/5)24631174=52/53/51/539/56/97/924631174
循环继续,第3列J=3
I=1,J-1
DO K=1,I-1 SUM=SUM-A(I,K)*A(K,J) END DO A 23 = A 23 − A 21 A 13 \quad A_{23}=A_{23}-A_{21}A_{13} A23=A23A21A13
A = [ 5 3 2 1 2 / 5 9 / 5 4 − 2 ( 2 / 5 ) 1 3 / 5 6 / 9 6 7 1 / 5 7 / 9 3 4 ] A=\left[ \begin{matrix} 5&3&2&1\\2/5&9/5&4-2(2/5)&1\\3/5&6/9&6&7\\1/5&7/9&3&4\end{matrix}\right] A=52/53/51/539/56/97/9242(2/5)631174

I=J,N
DO K=1,J-1 SUM=SUM-A(I,K)*A(K,J) END DO
A = [ 5 3 2 1 2 / 5 9 / 5 16 / 5 1 3 / 5 6 / 9 6 − 2 ( 3 / 5 ) − ( 16 / 5 ) ( 6 / 9 ) 7 1 / 5 7 / 9 3 − 2 ( 1 / 5 ) − ( 16 / 5 ) ( 7 / 9 ) 4 ] = [ 5 3 2 1 0.4 1.8 3.2 1 0.6 6 / 9 2.667 7 0.2 7 / 9 0.1111 4 ] A=\left[ \begin{matrix} 5&3&2&1\\2/5&9/5&16/5&1\\3/5&6/9&6-2(3/5)-(16/5)(6/9)&7\\1/5&7/9&3-2(1/5)-(16/5)(7/9)&4\end{matrix}\right]=\left[ \begin{matrix} 5&3&2&1\\0.4&1.8&3.2&1\\0.6&6/9&2.667&7\\0.2&7/9&0.1111&4\end{matrix}\right] A=52/53/51/539/56/97/9216/562(3/5)(16/5)(6/9)32(1/5)(16/5)(7/9)1174=50.40.60.231.86/97/923.22.6670.11111174
IMAX=3 \quad IMAX=J \quad \quad 于此同时INDX(3)=3 \quad 求下三角:
A = [ 5 3 2 1 0.4 1.8 3.2 1 0.6 6 / 9 2.667 7 0.2 7 / 9 1 / 24 4 ] A=\left[ \begin{matrix} 5&3&2&1\\0.4&1.8&3.2&1\\0.6&6/9&2.667&7\\0.2&7/9&1/24&4\end{matrix}\right] A=50.40.60.231.86/97/923.22.6671/241174

循环继续,第4列J=4
I=1,J-1
DO K=1,I-1 SUM=SUM-A(I,K)*A(K,J) END DO
A = [ 5 3 2 1 0.4 1.8 3.2 1 − 0.4 0.6 6 / 9 2.667 7 − 0.6 − ( 6 / 9 ) ( 1 − 0.4 ) 0.2 7 / 9 1 / 24 4 ] = [ 5 3 2 1 0.4 1.8 3.2 0.6 0.6 6 / 9 2.667 6 0.2 7 / 9 1 / 24 4 ] A=\left[ \begin{matrix} 5&3&2&1\\0.4&1.8&3.2&1-0.4\\0.6&6/9&2.667&7-0.6-(6/9)(1-0.4) \\0.2&7/9&1/24&4\end{matrix}\right]=\left[ \begin{matrix} 5&3&2&1\\0.4&1.8&3.2&0.6\\0.6&6/9&2.667&6 \\0.2&7/9&1/24&4\end{matrix}\right] A=50.40.60.231.86/97/923.22.6671/24110.470.6(6/9)(10.4)4=50.40.60.231.86/97/923.22.6671/2410.664
I=J,N
DO K=1,J-1 SUM=SUM-A(I,K)*A(K,J) END DO
A = [ 5 3 2 1 0.4 1.8 3.2 0.6 0.6 6 / 9 2.667 6 0.2 7 / 9 1 / 24 4 − 0.2 − 0.6 ( 7 / 9 ) − 6 ( 1 / 24 ) ] = [ 5 3 2 1 0.4 1.8 3.2 0.6 0.6 6 / 9 2.667 6 0.2 7 / 9 1 / 24 3.083 ] A=\left[ \begin{matrix} 5&3&2&1\\0.4&1.8&3.2&0.6\\0.6&6/9&2.667&6 \\0.2&7/9&1/24&4-0.2-0.6(7/9)-6(1/24)\end{matrix}\right]=\left[ \begin{matrix} 5&3&2&1\\0.4&1.8&3.2&0.6\\0.6&6/9&2.667&6 \\0.2&7/9&1/24&3.083\end{matrix}\right] A=50.40.60.231.86/97/923.22.6671/2410.6640.20.6(7/9)6(1/24)=50.40.60.231.86/97/923.22.6671/2410.663.083
IMAX=4 \quad IMAX=J \quad \quad 于此同时INDX(4)=4 \quad J=N \quad 不用求下三角:
最终: A = [ 5 3 2 1 0.4 1.8 3.2 0.6 0.6 6 / 9 2.667 6 0.2 7 / 9 1 / 24 3.083 ] A=\left[ \begin{matrix} 5&3&2&1\\0.4&1.8&3.2&0.6\\0.6&6/9&2.667&6 \\0.2&7/9&1/24&3.083\end{matrix}\right] A=50.40.60.231.86/97/923.22.6671/2410.663.083

3. 反向代入法

      SUBROUTINE LUBKSB (A, N, NP, INDX, B)

C-----  Linear equation solver based on LU decomposition

C-----  Use single precision on cray
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(NP,NP), INDX(N), B(N)

      II=0
      DO I=1,N
         LL=INDX(I)
         SUM=B(LL)
         B(LL)=B(I)

         IF (II.NE.0) THEN
            DO J=II,I-1
               SUM=SUM-A(I,J)*B(J)
            END DO
         ELSE IF (SUM.NE.0.) THEN
            II=I
         END IF

         B(I)=SUM
      END DO

      DO I=N,1,-1
         SUM=B(I)

         IF (I.LT.N) THEN
            DO J=I+1,N
               SUM=SUM-A(I,J)*B(J)
            END DO
         END IF

         B(I)=SUM/A(I,I)
      END DO

      RETURN
      END

继续用上节4×4矩阵演示:
求解方程组:
A x = y Ax=y Ax=y

[ 1 2 3 4 2 3 4 1 3 3 6 7 5 3 2 1 ] { x 1 x 2 x 3 x 4 } = { y 1 y 2 y 3 y 4 } \left[ \begin{matrix} 1&2&3&4\\2&3&4&1\\3&3&6&7\\5&3&2&1\end{matrix}\right]\left\{ \begin{matrix} x_1\\x_2\\x_3\\x_4 \end{matrix}\right\}=\left\{ \begin{matrix} y_1\\y_2\\y_3\\y_4 \end{matrix}\right\} 1235233334624171x1x2x3x4=y1y2y3y4

LUP分解后的矩阵 A L U P A^{LUP} ALUP
A L U P = [ 5 3 2 1 0.4 1.8 3.2 0.6 0.6 6 / 9 2.667 6 0.2 7 / 9 1 / 24 3.083 ] A^{LUP}=\left[ \begin{matrix} 5&3&2&1\\0.4&1.8&3.2&0.6\\0.6&6/9&2.667&6 \\0.2&7/9&1/24&3.083\end{matrix}\right] ALUP=50.40.60.231.86/97/923.22.6671/2410.663.083

置换关系INDX=4,2,3,4
y y y已知求解 x x x

进入程序,现根据置换关系和下三角矩阵调整 y y y:
I=1 y 1 y_1 y1 y 4 y_4 y4互换 且II=1 { y 4 y 2 y 3 y 1 } \quad \left\{ \begin{matrix} y_4&y_2&y_3&y_1 \end{matrix}\right\} {y4y2y3y1}
I=2利用下三角矩阵调整 y ( 2 ) y(2) y(2) 的值 { y 4 y 2 − 0.4 ∗ y ( 1 ) y 3 y 1 } \quad\left\{ \begin{matrix} y_4&y_2-0.4*y(1)&y_3&y_1 \end{matrix}\right\} {y4y20.4y(1)y3y1}
即: { y 4 y 2 − 0.4 ∗ y 4 y 3 y 1 } \quad\left\{ \begin{matrix} y_4&y_2-0.4*y_4&y_3&y_1 \end{matrix}\right\} {y4y20.4y4y3y1}
I=3利用下三角矩阵调整 y ( 2 ) y(2) y(2) 的值 { y ( 1 ) y ( 2 ) y ( 3 ) − 0.4 y ( 1 ) − 0.6 y ( 2 ) y ( 4 ) } \quad\left\{ \begin{matrix} y(1)&y(2)&y(3)-0.4y(1)-0.6y(2)&y(4) \end{matrix}\right\} {y(1)y(2)y(3)0.4y(1)0.6y(2)y(4)}
I=4利用下三角矩阵调整 y ( 2 ) y(2) y(2) 的值 { y ( 1 ) y ( 2 ) y ( 3 ) y ( 4 ) − 0.4 y ( 1 ) − 0.6 y ( 2 ) − 0.2 y ( 3 ) } \quad\left\{ \begin{matrix} y(1)&y(2)&y(3)&y(4)-0.4y(1)-0.6y(2)-0.2y(3) \end{matrix}\right\} {y(1)y(2)y(3)y(4)0.4y(1)0.6y(2)0.2y(3)}
此时已经转换成求解:
A U x = [ 5 3 2 1 0 1.8 3.2 0.6 0 0 2.667 6 0 0 0 3.083 ] { x 1 x 2 x 3 x 4 } = { y ( 1 ) y ( 2 ) y ( 3 ) y ( 4 ) } A^{U} x=\left[ \begin{matrix} 5&3&2&1\\0&1.8&3.2&0.6\\0&0&2.667&6 \\0&0&0&3.083\end{matrix}\right]\left\{ \begin{matrix} x_1\\x_2\\x_3\\x_4 \end{matrix}\right\}=\left\{ \begin{matrix} y(1)\\y(2)\\y(3)\\y(4) \end{matrix}\right\} AUx=500031.80023.22.667010.663.083x1x2x3x4=y(1)y(2)y(3)y(4)
反向代入:
x ( 4 ) = [ y ( 4 ) ] / A ( 4 , 4 ) x(4)=[y(4)]/A(4,4) x(4)=[y(4)]/A(4,4)
x ( 3 ) = [ y ( 3 ) − A ( 3 , 4 ) ∗ x ( 4 ) ] / A ( 3 , 3 ) x(3)=[y(3)-A(3,4)*x(4)]/A(3,3) x(3)=[y(3)A(3,4)x(4)]/A(3,3)
x ( 2 ) = [ y ( 2 ) − A ( 2 , 4 ) ∗ x ( 4 ) − A ( 2 , 3 ) ∗ x ( 3 ) ] / A ( 2 , 2 ) x(2)=[y(2)-A(2,4)*x(4)-A(2,3)*x(3)]/A(2,2) x(2)=[y(2)A(2,4)x(4)A(2,3)x(3)]/A(2,2)
x ( 1 ) = [ y ( 1 ) − A ( 1 , 4 ) ∗ x ( 4 ) − A ( 1 , 3 ) ∗ x ( 3 ) − A ( 1 , 2 ) ∗ x ( 2 ) ] / A ( 1 , 1 ) x(1)=[y(1)-A(1,4)*x(4)-A(1,3)*x(3)-A(1,2)*x(2)]/A(1,1) x(1)=[y(1)A(1,4)x(4)A(1,3)x(3)A(1,2)x(2)]/A(1,1)

4.利用LUP分解和反向代入法求逆矩阵

如果 y y y 的初始值为 y 1 = { 1 0 0 0 } T y_1=\left\{\begin{matrix} 1&0&0&0\end{matrix}\right\}^T y1={1000}T y 4 = { 0 0 0 1 } T y_4=\left\{\begin{matrix} 0&0&0&1 \end{matrix}\right\}^T y4={0001}T,求出的 [ x 1 x 2 x 3 x 4 ] \left[x_1\quad x_2\quad x_3\quad x_4\right] [x1x2x3x4] A A A的逆矩阵。
L 1 L 2 . . . L N − 1 U N × N = P N − 1 . . . P 2 P 1 A N × N L_{1}L_{2}...L_{N-1}U_{N\times N}=P_{N-1}...P_2P_{1}A_{N\times N} L1L2...LN1UN×N=PN1...P2P1AN×N

求出 U U U之后,通过 U U U的主对角线乘积算出 det ⁡ ( A ) = ( − 1 ) k det ⁡ ( U ) \det(A)=(-1)^{k}\det(U) det(A)=(1)kdet(U)( k为行交换的次数),判断逆矩阵是否存在。
由于 P i − 1 = P i P_i^{-1}=P_{i} Pi1=Pi L i − 1 = 2 I − L i L_i^{-1}=2I-L_{i} Li1=2ILi ,因此:
U − 1 L N − 1 − 1 . . . L 2 − 1 L 1 − 1 = A N × N − 1 P 1 P 2 . . . P N − 1 U^{-1}L_{N-1}^{-1}...L_2^{-1}L_1^{-1}=A_{N\times N}^{-1}P_{1}P_2...P_{N-1} U1LN11...L21L11=AN×N1P1P2...PN1

A N × N − 1 = U − 1 L N − 1 − 1 . . . L 2 − 1 L 1 − 1 P N − 1 . . . P 2 P 1 A_{N\times N}^{-1}=U^{-1}L_{N-1}^{-1}...L_2^{-1}L_1^{-1}P_{N-1}...P_2P_1 AN×N1=U1LN11...L21L11PN1...P2P1
U − 1 U^{-1} U1通过反向代入法求出。

程序:

      Program Main
      Implicit none
      Integer,parameter:: N=4
      Integer I,J
      Real A(N,N),B(N,N),C(N,N),Det
      DATA A/1,2,3,4,34,15,2,1,3,22,14,7,7,81,2,13/ !任意输入矩阵
      Do I=1,N   !备份一个相同的矩阵,方便后续验证
          B(I,:)=A(I,:)
          write(*,*) B(I,:)
      End Do
      write(*,*) "*****"
      Call matInv(A,N,Det)  !求逆
      write(*,*) "Det=",Det
      write(*,*) "*****"
      Do I=1,N
          write(*,*) A(I,:)
      End Do
      write(*,*) "*****"
      C=matmul(A,B)  !验证
      Do I=1,N
          write(*,*) C(I,:)
      End Do
      write(*,*) "*****"
      End Program Main
      
C *********************************

      Subroutine matInv(A,N,Det)
      Integer I,J,K,P(N)
      Real A(N,N),B(N,N),U(N,N),L(N,N),Tempt(N),Det
     
      
      Call LUP(A,N,P,Det)  !LUP分解

      Do I=1,N
          Det=Det*A(I,I)
      End Do

      Call Back(A,N)  !求上对角矩阵的逆
      Do I=1,N
          Do J=1,N
              If (I .GT. J) Then
                  U(I,J)=0.0
              Else
                  U(I,J)=A(I,J)
              End If
          End Do
      End Do
    
      Do J=N-1,1,-1
                    
          Do K=1,N
              Do I=1,N
                  If (I .EQ. K) Then
                      L(K,K)=1.0
                  Else
                      L(I,K)=0.0
                  End If
              End Do
          End Do
          Do K=J+1,N 
              L(K,J)=-A(K,J)
          End Do
          
          U=matmul(U,L)
      End Do
      Do J=N-1,1,-1
          If(P(J).NE.J) Then
              Do K=1,N
                  Tempt(K)=U(K,J)
                  U(K,J)=U(K,P(J))
                  U(K,P(J))=Tempt(K)
              End Do          
          End If
      End Do
      Do I=1,N
      Do J=1,N
          A(I,J)=U(I,J)
      End Do
      End Do
      
      End Subroutine matInv
      
C *********************************
      Subroutine LUP(A,N,P,Det)
      Integer I,J,K,N,P(N)
      Real A(N,N),Temp(N),Det,Jmax
      
      Det=1.0
      Do J=1,N
          P(J)=J
          Jmax=Abs(A(J,J))  !先假设每一行的最大值是其对角线值,这样就不用了换行了
          Do I=J,N
              If (Abs(A(I,J)).GT.Jmax) Then  !寻找此列最大值
                  Jmax=Abs(A(I,J))
                  P(J)=I
              End If
          End Do
          If (J .NE. P(J)) Then !换行
              Det=-Det
              Do K=1,N
                  Temp(K)=A(J,K)
                  A(J,K)=A(P(J),K)
                  A(P(J),K)=Temp(K)
              End Do
          End If
          
          Do I=J+1,N  !计算下三角矩阵,并更新相应的上三角矩阵
              A(I,J)=A(I,J)/A(J,J)
              Do K=J+1,N
                  A(I,K)=A(I,K)-A(I,J)*A(J,K)
              End Do
          End Do
      End Do
      Return
      End Subroutine LUP
      
C *********************************
      Subroutine Back(A,N)
      Integer I,J,K,N
      Real A(N,N),Temp
      
      Do J=N,1,-1
          Do I=J,1,-1
              If (I.EQ.J) Then
                  A(I,J)=1.0/A(I,J)
              Else
                  Temp=0.0
                  Do K=J,I+1,-1
                      Temp=Temp+A(K,J)*A(I,K)
                  End Do
                  A(I,J)=-Temp/A(I,I)
              End If          
          End Do
      End Do
     
      Return
      End Subroutine Back

5. 行列式知识点补充

5.1 初等变换

定义:如果B可以由A经过一系列初等变换得到,则称矩阵A与B称等价。
初等行(列)变换:
定义:所谓数域P上矩阵的初等行变换是指下列3种变换:
1)以P中一个非零的数乘矩阵的某一行
2)把矩阵的某一行的c倍加到另一行,这里c是P中的任意一个数
3)互换矩阵中两行的位置
一般来说,一个矩阵经过初等行变换后就变成了另一个矩阵,当矩阵A经过初等行变换变成矩阵B时,一般写作
可以证明:任意一个矩阵经过一系列初等行变换总能变成阶梯型矩阵。

5.2 行列式初等变换的相关性质

性质1:行列互换,行列式不变
性质2:一数乘行列式的一行就相当于这个数乘此行列式
性质3:如果行列式中有两行相同,那么行列式为0,所谓两行相同,即两行对应的元素都相等
性质4:如果行列式中,两行成比例,那么该行列式为0
性质5:把一行的倍数加到另一行,行列式不变
性质6:对换行列式中两行的位置,行列式反号

5.3 矩阵乘积的行列式

引理:一个n阶矩阵A总可以通过第二种行和列的初等变换化为一个对角矩阵
A = [ A 11 A 12 . . . A 1 N A 21 A 22 . . . A 2 N . . . . . . . . . . . . A N 1 A N 2 . . A N N ] → A ˉ [ A ˉ 11 0 . . . 0 0 A ˉ 22 . . . 0 . . . . . . . . . . . . 0 0 . . A ˉ N N ] A=\left[ \begin{matrix} A_{11}&A_{12}&...&A_{1N}\\A_{21}&A_{22}&...&A_{2N}\\...&...&...&...\\A_{N1}&A_{N2}&..&A_{NN}\\ \end{matrix}\right] \to \bar{A}\left[ \begin{matrix} \bar{A}_{11}&0&...&0\\0&\bar{A}_{22}&...&0\\...&...&...&...\\0&0&..&\bar{A}_{NN}\\ \end{matrix}\right] A=A11A21...AN1A12A22...AN2...........A1NA2N...ANNAˉAˉ110...00Aˉ22...0...........00...AˉNN
第二种行和列的初等变换化,不改变行列式大小,因此 det ⁡ ( A ) = det ⁡ ( A ˉ ) \det(A)=\det(\bar A) det(A)=det(Aˉ)

定理:设A,B是任意两个n阶矩阵,那么 det ⁡ ( A B ) = det ⁡ ( A ) ⋅ det ⁡ ( B ) \det(AB)=\det(A)\cdot \det(B) det(AB)=det(A)det(B)

由此可知: det ⁡ ( A ) = ( − 1 ) k det ⁡ ( U ) \det(A)=(-1)^{k}\det(U) det(A)=(1)kdet(U) k为行交换的次数。

评论 21
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值