meteor 发表于 2008-3-29 21:33

乘幂法求频率和振型

哪个大侠给我讲解一下乘幂法求频率和振型这种方法 有源程序也行!谢谢了

meteor 发表于 2008-3-29 23:11

***********************************************************************
*   程序说明:N   --自由度                                                                                        *
*                NS--需求的振型数目                                                                        *
*                          RTOL--控制精度                                                                                *
*                          SM--质量矩阵                                                                                *
*                          MS12--质量矩阵控制系数:1--对角阵                                                *
*                                                  2--满阵                                                *
*                          KD12--刚度柔度控制系数:1--柔度矩阵                                        *
*                                                  2--刚度矩阵                                        *
*                          W(N)--频率向量                                                                                *
*                          T(N)--周期向量                                                                                *
***********************************************************************
        REAL A(1000)
        OPEN(1,FILE='INPUT.DAT')
        OPEN(2,FILE='OUTPUT.DAT')
        READ(1,*)N,NS,MS12,KD12,RTOL
        KA=1
        KYO=KA+N*N
        KSM=KYO+N
        KU=KSM+N*N
        KYOP=KU+N*NS
        KR=KYOP+N*N
        KBT=KR+N
        KSX=KBT+N
        KYK=KSX+N
        KW=KYK+N
        KT=KW+N
        CALL ITER(A(KA),A(KYO),A(KSM),A(KU),A(KYOP),A(KR),
   !         A(KBT),A(KSX),A(KYK),A(KW),A(KT),N,NS,MS12,KD12,RTOL)
      STOP
        END
*
*-----------------------------------------------------------------------
        SUBROUTINE ITER(A,YO,SM,U,YOP,R,BT,SX,YK,W,T,N,MJ,MS12,KD12,RTOL)
*-----------------------------------------------------------------------
        REAL A(N,N),YO(N),SM(N,N),U(N,MJ),YOP(N),R(N),BT(N),SX(N)
   !   ,YK(N),W(N),T(N)
      WRITE(2,*)'THE SYSTEM STIFF MATRIX'
*
        DO 5 I=1,N
           READ(1,*)(A(I,J),J=1,N)
           WRITE(2,*)(A(I,J),J=1,N)
5        CONTINUE
*
      WRITE(2,*)'THE SYSTEM MASS MATRIX'
        IF (MS12.EQ.2) THEN
          DO 6 I=1,N
             READ(1,*)(SM(I,J),J=1,N)
             WRITE(2,*)(SM(I,J),J=1,N)
6          CONTINUE
*
      ELSE
          READ(1,*)(SM(I,I),I=1,N)
          WRITE(2,*)(SM(I,I),I=1,N)
        END IF
*
        IF (KD12.EQ.2) THEN
          CALL UTDU3(A,SX,N)
        END IF
*
        DO 300 J=1,MJ
          J1=J-1
          NITER=0
          W1=0
*          
          DO 10 I=1,N
          U(I,J)=1.0
10          CONTINUE
*
          DO 30 I=1,N
             SUM=0.0
             DO 20 K=1,N
                SUM=SUM+SM(I,K)*U(K,J)
20               CONTINUE
               YO(I)=SUM
30          CONTINUE
*
40          NITER=NITER+1
      IF(NITER.GE.40) STOP
*
          DO 50 I=1,N
             YOP(I)=YO(I)
50          CONTINUE
*
      IF(J.EQ.1)GOTO 80
          DO 70 L=1,J1
             RJ=0.0
             DO 60 I=1,N
                RJ=RJ+U(I,J)*YO(I)
60               CONTINUE
         BT(L)=R(L)*RJ

             write(*,*)R(L),rj

70          CONTINUE
*
80          IF (KD12.EQ.2)THEN
         CALL DECOMP(A,YO,N)
          ELSE

         END IF
*
           DO 120 I=1,N
              U(I,J)=YO(I)
120      CONTINUE
*
         IF(J.EQ.1) GOTO 150
           DO 140 I=1,N
              ST=0.0
              DO 130 L=1,J1
                 ST=ST+BT(L)*U(I,L)
130              CONTINUE
            YO(I)=YO(I)-ST
140           CONTINUE
*
150      DO 170 I=1,N
            SUM=0.0
                  DO 160 K=1,N
                     SUM=SUM+SM(I,K)*YO(K)
160              CONTINUE
            YK(I)=SUM
170           CONTINUE
*
         ST=0.0
           RJ=0.0
           DO 180 I=1,N
                  ST=ST+YO(I)*YOP(I)
                  RJ=RJ+YO(I)*YK(I)
180           CONTINUE
      
*   
           W2=ST/RJ
           ST=1.0/SQRT(RJ)
           DO 190 I=1,N
                  YO(I)=YK(I)*ST
                  U(I,J)=U(I,J)*ST
190           CONTINUE
*
           EL=ABS(W2-W1)
           W1=W2
           EL=EL/W2
           IF(EL.GT.RTOL) GOTO 40
*
           DO 210 I=1,N
              SUM=0.0
              DO 200 K=1,N
                 SUM=SUM+SM(I,K)*U(K,J)
200                  CONTINUE
              YK(I)=SUM
210           CONTINUE
*
           RJ=0.0
           DO 220 I=1,N
              RJ=RJ+YK(I)*U(I,J)
220           CONTINUE
*
           R(J)=1/(W2*RJ)
           W(J)=SQRT(W2)
           WRITE(2,230)J,W(J)
           WRITE(2,240)J,(U(I,J)/U(1,J),I=1,N)
           WRITE(*,*)NITER
230      FORMAT('W(',I2,')=',F15.7,5X)
240      FORMAT('A(',I2,')=',5F15.7)
300   CONTINUE
      RETURN
        END
*
*----------------------------------------------------------------------
        SUBROUTINE UTDU3(A,C,N)
*----------------------------------------------------------------------
        REAL A(N,N),C(N)
        IF (N.EQ.1) GOTO 40
        L=N-1
*
        DO 20 I=1,L
          IF(A(I,I).LE.0.0) GOTO 30
          M=I+1
          DO 10 J=M,N
          C(J)=A(I,J)
          A(I,J)=A(I,J)/A(I,I)
          DO 10 K=M,J
              A(K,J)=A(K,J)-A(I,J)*C(K)
10          CONTINUE
20    CONTINUE
*
      IF(A(N,N).GT.0) GOTO 40
        I=N
30        WRITE(*,*)I,A(I,J)
        STOP
40        RETURN
        END

*
*-------------------------------------------------------------------------
        SUBROUTINE DECOMP(A,B,N)
*-------------------------------------------------------------------------
        REAL A(N,N),B(N)
        IF (N.GT.1)GOTO 10
        B(1)=B(1)/A(1,1)
        GOTO 50
10    L=N-1
*
      DO 20 I=2,N
          K=I-1
          DO 20 J=1,K
          B(I)=B(I)-B(J)*A(J,I)
20        CONTINUE
*
      DO 30 I=1,N
          B(I)=B(I)/A(I,I)
30        CONTINUE
*
        DO 40 I=1,L
          J1=L-I+1
          IP=J1+1
          DO 40 J=IP,N
          B(J1)=B(J1)-B(J)*A(J1,J)
40        CONTINUE
*
50    RETURN
        END

meteor 发表于 2008-3-29 23:13

拜托各位大侠 帮我看看这个程序有点问题不能运行 下面是输入的值
3
3
1
2
5E-12
3 -2 0
-2 5 -3
0 -3 3
1
2
3
页: [1]
查看完整版本: 乘幂法求频率和振型