狼跃冲顶 发表于 2007-1-25 18:46

求出的自振频率的平方是负数,非常郁闷特向大家请教

论坛上的朋友们大家好哦,我最近遇到一个问题,非常郁闷特向大家请教!
大致是:我现在用fortran编程求出一个剪切型楼层的自振频率(此楼层的整个刚度矩阵是实对称的)。方法是雅克比法。很奇怪的是我求出的自振频率的平方是负数,也就是无法开平方求出自振频率了。现在很急,那位朋友能帮帮忙,讨论下是否也遇到过这个情况,以及如何去解决。先谢谢了。
本人QQ:40803237

欧阳中华 发表于 2007-1-25 19:51

.
    首先应该考核一下所用的JACOBI程序是否正确;

    在用自编程序计算一个简单的例题,看看结果是不是符合实际;

    这样再去做其他规模交大一些的算例就不会出问题了,俗话说:“磨刀不误砍柴功”... ..

狼跃冲顶 发表于 2007-1-26 15:35

:'( 这个我设置过断点,雅克比程序是没有问题的。而起按照数学理论实对称矩阵只能保证特征值是实数。不能保证为正数

欧阳中华 发表于 2007-1-26 20:03

.
   实对称矩阵的确不会出现复特征值,但如果振动系统不考虑阻尼,且系统为线弹性,那么对应的特征值一定是正值,开放就是系统的圆频率形式的固有频率。

狼跃冲顶 发表于 2007-1-28 12:27

唉,那我的问题正好是考虑了阻尼和进入了塑性阶段!现在该怎么解决呢?好痛苦啊

欧阳中华 发表于 2007-1-28 13:24

.
   考虑阻尼那要看用的是什么阻尼形式,如果是质量矩阵和刚度矩阵的线性组合,那么也还是正的特征值;

    问题中的变形进入到塑性,就成了非线性问题了,非线性的模态概念还得相应去研究、研究吧.. ..

    有什么进展也告诉大家呦 .. ..

狼跃冲顶 发表于 2007-1-28 15:20

呵呵 我用的阻尼就是瑞雷阻尼 (质量矩阵和刚度矩阵的线性组合)进入非线性特征值有部分变为负的,另外还有一个就是 我的程序编译 连接 都没有错误但是 运行 却输不出值来,这不知道是怎么回事?还望朋友指教啊,

狼跃冲顶 发表于 2007-1-29 17:01

这个具体解决需要看那些方面的书呢?:@Q

欧阳中华 发表于 2007-1-29 18:14

.
   非线性动力有限元... .

狼跃冲顶 发表于 2007-1-30 11:24

得到一个好消息:出现复数的话就不能用JACOBI了,用其它算法如实赫申伯格矩阵的QR算法
也有一个坏消息:这种计算方法我现在暂时找不到,朋友门帮帮忙找找啊

[ 本帖最后由 狼跃冲顶 于 2007-1-30 11:29 编辑 ]

欧阳中华 发表于 2007-1-30 11:40

.
   JACOBI方法是通过矩阵转换方法得到特征值的,前提是矩阵为实对称矩阵,书上都写的呀。

   QR法的程序到是可以找到,但还是应该研究一下再用。

C-----------------------------------------------------------------------
C
C-----------------------------------------------------------------------
      DIMENSION A(5,5),WR(5),WI(5),Z(5,5),SCALE(5),ORT(5)
      OPEN(10,FILE='I')
      READ(10,*) A
      CLOSE (10)
      N=5
      CALL QR1(N,A,WR,WI,Z,SCALE,ORT)
      OPEN(20,FILE='O')
      DO 300 I=1,N
      WRITE(20,100) I,WR(I),WI(I)
100FORMAT(//,5X,I3,2F10.4//)
      WRITE(20,200) (Z(I,J),J=1,N)
200FORMAT(5X,5F10.4//)
300CONTINUE
      CLOSE (20)
      STOP
      END
C-----------------------------------------------------------------------
C
C-----------------------------------------------------------------------
      SUBROUTINE QR1(N,A,WR,WI,Z,SCALE,ORT)
      DIMENSION A(N,N),WR(N),WI(N),Z(N,N),SCALE(N),ORT(N)
C
200FORMAT(//5(5X,5F10.1/))
      CALL BALAN(N,A,LOW,IGH,SCALE)
      WRITE(*,200) ((A(I,J),J=1,N),I=1,N)
      WRITE(*,*) LOW,IGH
      CALL ORTH(N,LOW,IGH,A,ORT)
      WRITE(*,200) ((A(I,J),J=1,N),I=1,N)
      CALL ORTR(N,LOW,IGH,A,ORT,Z)
      WRITE(*,200) ((A(I,J),J=1,N),I=1,N)
      CALL HQR2(N,LOW,IGH,A,WR,WI,Z,IERR)
C
      WRITE(*,100) IERR
100FORMAT(/,5X,'IERR OF QR1=',I5,/)
      IF(IERR.NE.0) RETURN
C
      CALL BALBA(N,LOW,IGH,SCALE,N,Z)
C
      RETURN
      END
C-----------------------------------------------------------------------
C
C-----------------------------------------------------------------------
      SUBROUTINE BALAN(N,A,LOW,IGH,SCALE)
      DIMENSION A(N,N),SCALE(N)
      RADIX=4.0
      B2=RADIX*RADIX
      K=1
      L=N
      GOTO 100
20    SCALE(M)=J
      WRITE(*,*)M
      IF(J.EQ.M) GO TO 50
      DO 30 I=1,L
      F=A(I,J)
      A(I,J)=A(I,M)
      A(I,M)=F
30    CONTINUE
      DO 40 I=K,N
      F=A(J,I)
      A(J,I)=A(M,I)
      A(M,I)=F
40    CONTINUE
50    GO TO(80,130),IEXC
80    IF(L.EQ.1) GO TO 280
      L=L-1
100   DO 120 JJ=1,L
      J=L+1-JJ
      DO 110 I=1,L
      IF(I.EQ.J) GO TO 110
      IF(A(J,I).NE.0.0) GO TO 120
110   CONTINUE
      M=L
      IEXC=1
      GO TO 20
120   CONTINUE
      GO TO 140
130   K=K+1
140   DO 170 J=K,L
      DO 150 I=K,L
      IF(I.EQ.J) GO TO 150
      IF(A(I,J).NE.0.0) GO TO 170
150   CONTINUE
      M=K
      IEXC=2
      GO TO 20
170   CONTINUE
      DO 180 I=K,L
180   SCALE(I)=1.0
190   NOCONV=0
      DO 270 I=K,L
      C=0.0
      R=0.0
      DO 200 J=K,L
      IF(J.EQ.I) GO TO 200
      C=C+ABS(A(J,I))
      R=R+ABS(A(I,J))
200   CONTINUE
      IF(C.EQ.0.0) GO TO 270
      IF(R.EQ.0.0) GO TO 270
      G=R/RADIX
      F=1.0
      S=C+R
210   IF(C.GE.G) GO TO 220
      F=F*RADIX
      C=C*B2
      GO TO 210
220   G=R*RADIX
230   IF(C.LT.G) GO TO 240
      F=F/RADIX
      C=C/B2
      GO TO 230
240   IF((C+R)/F.GE.0.95*S) GO TO 270
      G=1.0/F
      SCALE(I)=SCALE(I)*F
      NOCONV=1
      DO 250 J=K,N
250   A(I,J)=A(I,J)*G
      DO 260 J=1,L
260   A(J,I)=A(J,I)*F
270   CONTINUE
      IF(NOCONV.EQ.1) GO TO 190
280   LOW=K
      IGH=L
      RETURN
      END
C-----------------------------------------------------------------------
C
C-----------------------------------------------------------------------
      SUBROUTINE ORTR(N,LOW,IGH,A,ORT,Z)
      DIMENSION A(N,N),ORT(N),Z(N,N)
      DO 80 I=1,N
      DO 60 J=1,N
60    Z(I,J)=0.0
      Z(I,I)=1.0
80    CONTINUE
      KL=IGH-LOW-1
      IF(KL.LT.1) GO TO 200
      DO 140 MM=1,KL
      MP=IGH-MM
      IF(A(MP,MP-1).EQ.0.0) GO TO 140
      MP1=MP+1
      DO 100 I=MP1,IGH
100   ORT(I)=A(I,MP-1)
      DO 130 J=MP,IGH
      G=0.0
      DO 110 I=MP,IGH
110   G=G+ORT(I)*Z(I,J)
      G=(G/ORT(MP))/A(MP,MP-1)
      DO 120 I=MP,IGH
120   Z(I,J)=Z(I,J)+G*ORT(I)
130   CONTINUE
140   CONTINUE
200   RETURN
      END
C-----------------------------------------------------------------------
C
C-----------------------------------------------------------------------
      SUBROUTINE ORTH(N,LOW,IGH,A,ORT)
      DIMENSION A(N,N),ORT(N)
      LA=IGH-1
      KP1=LOW+1
      IF(LA.LT.KP1) GO TO 200
      DO 180 M=KP1,LA
      H=0.0
      ORT(M)=0.0
      SCALE=0.0
      DO 90 I=M,IGH
90    SCALE=SCALE+ABS(A(I,M-1))
      IF(SCALE.EQ.0.0) GO TO 180
      MP=M+IGH
      DO 100 II=M,IGH
      I=MP-II
      ORT(I)=A(I,M-1)/SCALE
      H=H+ORT(I)*ORT(I)
100   CONTINUE
      G=-SIGN(SQRT(H),ORT(M))
      H=H-ORT(M)*G
      ORT(M)=ORT(M)-G
      DO 130 J=M,N
      F=0.
      DO 110 II=M,IGH
      I=MP-II
      F=F+ORT(I)*A(I,J)
110   CONTINUE
      F=F/H
      DO 120 I=M,IGH
120   A(I,J)=A(I,J)-F*ORT(I)
130   CONTINUE
      DO 160 I=1,IGH
      F=0.
      DO 140 JJ=M,IGH
      J=MP-JJ
      F=F+ORT(J)*A(I,J)
140   CONTINUE
      F=F/H
      DO 150 J=M,IGH
150   A(I,J)=A(I,J)-F*ORT(J)
160   CONTINUE
      ORT(M)=SCALE*ORT(M)
      A(M,M-1)=SCALE*G
180   CONTINUE
200   RETURN
      END
C----------------------------------------------------------------------
C
C----------------------------------------------------------------------
      SUBROUTINE HQR2(N,LOW,IGH,H,WR,WI,Z,IERR)
      INTEGER EN,ENM2
      REAL NORM,MACHEP
      DIMENSION H(N,N),WR(N),WI(N),Z(N,N)
      MACHEP=1.E-8
      IERR=0
      NORM=0.0
      K=1
      DO 50 I=1,N
      DO 40 J=K,N
40    NORM=NORM+ABS(H(I,J))
      K=I
      IF(I.LT.LOW) GO TO 45
      IF(I.LE.IGH) GO TO 50
45   WR(I)=H(I,I)
      WI(I)=0.0
50   CONTINUE
      EN=IGH
      T=0.0
60   IF(EN.LT.LOW) GO TO 340
      ITS=0
      NA=EN-1
      ENM2=NA-1
70   N11=N-1
      DO 80 LL=LOW,EN
      L=EN+LOW-LL
      IF(L.EQ.LOW) GO TO 100
      S=ABS(H(L-1,L-1))+ABS(H(L,L))
      IF(S.EQ.0.0) S=NORM
      IF(ABS(H(L,L-1)).LE.MACHEP*S) GO TO 100
80    CONTINUE
100   X=H(EN,EN)
      IF(L.EQ.EN) GO TO 270
      Y=H(NA,NA)
      W=H(EN,NA)*H(NA,EN)
      IF(L.EQ.NA) GO TO 280
      IF(ITS.EQ.30) GO TO 1000
      IF(ITS.EQ.10) GO TO 110
      IF(ITS.NE.20) GO TO 130
110   T=T+X
      DO 120 I=LOW,EN
120   H(I,I)=H(I,I)-X
      S=ABS(H(EN,NA))+ABS(H(NA,ENM2))
      X=0.75*S
      Y=X
      W=-0.4375*S*S
130   ITS=ITS+1
      DO 140 MM=L,ENM2
      M=ENM2+L-MM
      ZZ=H(M,M)
      R=X-ZZ
      S=Y-ZZ
      P=(R*S-W)/H(M+1,M)+H(M,M+1)
      Q=H(M+1,M+1)-ZZ-R-S
      R=H(M+2,M+1)
      S=ABS(P)+ABS(Q)+ABS(R)
      P=P/S
      Q=Q/S
      R=R/S
      IF(M.EQ.L) GO TO 150
      IF(ABS(H(M,M-1))*(ABS(Q)+ABS(R)).LE.MACHEP*ABS(P)
   1*(ABS(H(M-1,M-1))+ABS(ZZ)+ABS(H(M+1,M+1)))) GO TO 150
140   CONTINUE
150   MP2=M+2
      DO 160 I=MP2,EN
      H(I,I-2)=0.0
      IF(I.EQ.MP2) GO TO 160
      H(I,I-3)=0.0
160   CONTINUE
      DO 260 K=M,NA
      NOTLAS=1
      IF(K.EQ.NA) NOTLAS=0
      IF(K.EQ.M) GO TO 170
      P=H(K,K-1)
      Q=H(K+1,K-1)
      R=0.0
      IF(NOTLAS.EQ.1) R=H(K+2,K-1)
      X=ABS(P)+ABS(Q)+ABS(R)
      IF(X.EQ.0.0) GO TO 260
      P=P/X
      Q=Q/X
      R=R/X
170   WOR=P*P+Q*Q+R*R
      SW=SQRT(WOR)
      S=SIGN(SW,P)
      IF(K.EQ.M) GO TO 180
      H(K,K-1)=-S*X
      GO TO 190
180   IF(L.NE.M) H(K,K-1)=-H(K,K-1)
190   P=P+S
      X=P/S
      Y=Q/S
      ZZ=R/S
      Q=Q/P
      R=R/P
      DO 210 J=K,N
      P=H(K,J)+Q*H(K+1,J)
      IF(NOTLAS.EQ.0) GO TO 200
      P=P+R*H(K+2,J)
      H(K+2,J)=H(K+2,J)-P*ZZ
200   H(K+1,J)=H(K+1,J)-P*Y
      H(K,J)=H(K,J)-P*X
210   CONTINUE
      J=K+3
      IF(EN.LT.J) J=EN
      DO 230 I=1,J
      P=X*H(I,K)+Y*H(I,K+1)
      IF(NOTLAS.EQ.0) GO TO 220
      P=P+ZZ*H(I,K+2)
      H(I,K+2)=H(I,K+2)-P*R
220   H(I,K+1)=H(I,K+1)-P*Q
      H(I,K)=H(I,K)-P
230   CONTINUE
      DO 250 I=LOW,IGH
      P=X*Z(I,K)+Y*Z(I,K+1)
      IF(NOTLAS.EQ.0) GO TO 240
      P=P+ZZ*Z(I,K+2)
      Z(I,K+2)=Z(I,K+2)-P*R
240   Z(I,K+1)=Z(I,K+1)-P*Q
      Z(I,K)=Z(I,K)-P
250   CONTINUE
260   CONTINUE
      GO TO 70
270   H(EN,EN)=X+T
      WR(EN)=H(EN,EN)
      WI(EN)=0.0
      EN=NA
      GO TO 60
280   P=(Y-X)/2.0
      Q=P*P+W
      ZZ=SQRT(ABS(Q))
      H(EN,EN)=X+T
      X=H(EN,EN)
      H(NA,NA)=Y+T
      IF(Q.LT.0.0) GO TO 320
      ZZ=P+SIGN(ZZ,P)
      WR(NA)=X+ZZ
      WR(EN)=WR(NA)
      IF(ZZ.NE.0.0) WR(EN)=X-W/ZZ
      WI(NA)=0.0
      WI(EN)=0.0
      X=H(EN,NA)
      S=ABS(X)+ABS(ZZ)
      P=X/S
      Q=ZZ/S
      WOR=P*P+Q*Q
      R=SQRT(WOR)
      P=P/R
      Q=Q/R
      DO 290 J=NA,N
      ZZ=H(NA,J)
      H(NA,J)=Q*ZZ+P*H(EN,J)
      H(EN,J)=Q*H(EN,J)-P*ZZ
290   CONTINUE
      DO 300 I=1,EN
      ZZ=H(I,NA)
      H(I,NA)=Q*ZZ+P*H(I,EN)
      H(I,EN)=Q*H(I,EN)-P*ZZ
300   CONTINUE
      DO 310 I=LOW,IGH
      ZZ=Z(I,NA)
      Z(I,NA)=Q*ZZ+P*Z(I,EN)
      Z(I,EN)=Q*Z(I,EN)-P*ZZ
310   CONTINUE
      GO TO 330
320   WR(NA)=X+P
      WR(EN)=X+P
      WI(NA)=ZZ
      WI(EN)=-ZZ
330   EN=ENM2
      GO TO 60
340   IF(NORM.EQ.0.0) GO TO 1001
      DO 800 NN=1,N
      EN=N+1-NN
      P=WR(EN)
      Q=WI(EN)
      NA=EN-1
      IF(Q) 710,600,800
600   M=EN
      H(EN,EN)=1.0
      IF(NA.EQ.0) GO TO 800
      DO 700 II=1,NA
      I=EN-II
      W=H(I,I)-P
      R=H(I,EN)
      IF(M.GT.NA) GO TO 620
      DO 610 J=M,NA
610   R=R+H(I,J)*H(J,EN)
620   IF(WI(I).GE.0.0) GO TO 630
      ZZ=W
      S=R
      GO TO 700
630   M=I
      IF(WI(I).NE.0.0) GO TO 640
      T=W
      IF(W.EQ.0.0) T=MACHEP*NORM
      H(I,EN)=-R/T
      GO TO 700
640   X=H(I,I+1)
      Y=H(I+1,I)
      Q=(WR(I)-P)*(WR(I)-P)+WI(I)*WI(I)
      T=(X*S-ZZ*R)/Q
      H(I,EN)=T
      IF(ABS(X).LE.ABS(ZZ)) GO TO 650
      H(I+1,EN)=(-R-W*T)/X
      GO TO 700
650   H(I+1,EN)=(-S-Y*T)/ZZ
700   CONTINUE
      GO TO 800
710   M=NA
      IF(ABS(H(EN,NA)).LE.ABS(H(NA,EN))) GO TO 720
      H(NA,NA)=-(H(EN,EN)-P)/H(EN,NA)
      H(NA,EN)=-Q/H(EN,NA)
      GO TO 730
720Z3=H(NA,NA)-P
      WOR=-H(NA,EN)
C
          CALL CDIV(WOR,0.0,Z3,Q,RES1,RES2)
C
      H(NA,NA)=RES1
      H(NA,EN)=RES2
730   H(EN,NA)=1.0
      H(EN,EN)=0.0
      ENM2=NA-1
      IF(ENM2.EQ.0.0) GO TO 800
      DO 790 II=1,ENM2
      I=NA-II
      W=H(I,I)-P
      RA=H(I,EN)
      SA=0.0
      DO 760 J=M,NA
      RA=RA+H(I,J)*H(J,NA)
      SA=SA+H(I,J)*H(J,EN)
760   CONTINUE
      IF(WI(I).GE.0.0) GO TO 770
      ZZ=W
      R=RA
      S=SA
      GO TO 790
770   M=I
      IF(WI(I).NE.0.0) GO TO 780
C
          CALL CDIV(-RA,-SA,W,Q,RES1,RES2)
C
      H(I,NA)=RES1
      H(I,EN)=RES2
      GO TO 790
780   X=H(I,I+1)
      Y=H(I+1,I)
      VR=(WR(I)-P)*(WR(I)-P)+WI(I)*WI(I)-Q*Q
      VI=(WR(I)-P)*2.0*Q
      IF(VR.NE.0.0) GO TO 781
      IF(VI.EQ.0.0) VR=MACHEP*NORM*(ABS(W)+ABS(Q)+ABS(X)+ABS(Y)
   1+ABS(ZZ))
C
781       CALL CDIV(X*R-ZZ*RA+Q*SA,X*S-ZZ*SA-Q*RA,VR,VI,RES1,RES2)
C
      H(I,NA)=RES1
      H(I,EN)=RES2
      IF(ABS(X).LE.ABS(ZZ)+ABS(Q)) GO TO 785
      H(I+1,NA)=(-RA-W*H(I,NA)+Q*H(I,EN))/X
      H(I+1,EN)=(-SA-W*H(I,EN)-Q*H(I,NA))/X
      GO TO 790
785   Z3=-R-Y*H(I,NA)
      Z4=-S-Y*H(I,EN)
C
          CALL CDIV(Z3,Z4,ZZ,Q,RES1,RES2)
C
      H(I+1,NA)=RES1
      H(I+1,EN)=RES2
790   CONTINUE
800   CONTINUE
      DO 840 I=1,N
      IF(I.LT.LOW) GO TO 810
      IF(I.LE.IGH) GO TO 840
810   DO 820 J=I,N
820   Z(I,J)=H(I,J)
840   CONTINUE
      DO 910 JJ=LOW,N
      J=N+LOW-JJ
      M=J
      IF(J.GT.IGH) M=IGH
      L=J-1
      IF(WI(J).GE.0.0) GO TO 900
      DO 880 I=LOW,IGH
      Y=0.0
      ZZ=0.0
      DO 860 K=LOW,M
      Y=Y+Z(I,K)*H(K,L)
860   ZZ=ZZ+Z(I,K)*H(K,J)
      Z(I,J)=ZZ
      Z(I,L)=Y
880   CONTINUE
900   IF(WI(J).NE.0.0) GO TO 910
      DO 920 I=LOW,IGH
      ZZ=0.0
      DO 930 K=LOW,M
930   ZZ=ZZ+Z(I,K)*H(K,J)
920   Z(I,J)=ZZ
910   CONTINUE
      GO TO 1001
1000 IERR=EN
1001 RETURN
      END
C-----------------------------------------------------------------------
C
C-----------------------------------------------------------------------
      SUBROUTINE CDIV(A,B,C,D,E,F)
      IF(ABS(C).LT.ABS(D)) GO TO 1
      R=D/C
      DEN=C+R*D
      E=(A+B*R)/DEN
      F=(B-A*R)/DEN
      GO TO 2
1       R=C/D
      DEN=D+R*C
      E=(A*R+B)/DEN
      F=(B*R-A)/DEN
2      RETURN
      END
C-----------------------------------------------------------------------
C
C-----------------------------------------------------------------------
      SUBROUTINE BALBA(N,LOW,IGH,SCALE,M,Z)
      DIMENSION SCALE(N),Z(N,N)
      IF(M.EQ.0) GO TO 200
      IF(IGH.EQ.LOW) GO TO 120
      DO 110 I=LOW,IGH
      S=SCALE(I)
      DO 100 J=1,M
100   Z(I,J)=Z(I,J)*S
110   CONTINUE
120   DO 140 II=1,N
      I=II
      IF(I.LT.LOW) GO TO 125
      IF(I.LE.IGH) GO TO 140
125   IF(I.LT.LOW) I=LOW-II
      K=SCALE(I)
      IF(K.EQ.I) GO TO 140
      DO 130 J=1,M
      S=Z(I,J)
      Z(I,J)=Z(K,J)
      Z(K,J)=S
130   CONTINUE
140   CONTINUE
200   RETURN
      END
C-----------------------------------------------------------------------

狼跃冲顶 发表于 2007-1-30 16:59

中华 同仁 再次谢谢你 有没有关于实赫申伯格矩阵的QR算法的相关理论??我想先看看理论才敢用这个算法

wsjxaut 发表于 2007-1-30 23:54

原帖由 狼跃冲顶 于 2007-1-28 12:27 发表
唉,那我的问题正好是考虑了阻尼和进入了塑性阶段!现在该怎么解决呢?好痛苦啊

如果是塑性问题,现有的算法是求不出模态的 ,因为这是一个非线性的模态分析问题,传统的不管是QR法还是Jaccobi法,乃至其它什么方法,都没戏。

linqus 发表于 2007-2-13 08:40

个人以为,
如果矩阵阶数不大的话,直接用fortran的imsl库函数求解,
看看结果是否正确,
不要在求解特征值的方法上纠缠,

如果库函数求解没问题,说明编程没啥问题,
再考虑使用的特征值求解方法。

探讨

linqus 发表于 2007-2-13 08:42

另外,
关于非线性模态,在kjbathe的书上略有介绍,
不知楼主可否分享一下心得,
谢谢。
页: [1] 2
查看完整版本: 求出的自振频率的平方是负数,非常郁闷特向大家请教