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

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

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

论坛上的朋友们大家好哦，我最近遇到一个问题，非常郁闷特向大家请教！
大致是：我现在用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)
100  FORMAT(//,5X,I3,2F10.4//)
      WRITE(20,200) (Z(I,J),J=1,N)
200  FORMAT(5X,5F10.4//)
300  CONTINUE
      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
200  FORMAT(//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
100  FORMAT(/,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
720  Z3=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的书上略有介绍，
不知楼主可否分享一下心得，
谢谢。

账号		自动登录	找回密码
密码			我要加入

[非线性振动] 求出的自振频率的平方是负数，非常郁闷特向大家请教

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

评分