求出的自振频率的平方是负数，非常郁闷特向大家请教
论坛上的朋友们大家好哦，我最近遇到一个问题，非常郁闷特向大家请教！大致是：我现在用fortran编程求出一个剪切型楼层的自振频率（此楼层的整个刚度矩阵是实对称的）。方法是雅克比法。很奇怪的是我求出的自振频率的平方是负数，也就是无法开平方求出自振频率了。现在很急，那位朋友能帮帮忙，讨论下是否也遇到过这个情况，以及如何去解决。先谢谢了。
本人QQ：40803237 .
首先应该考核一下所用的JACOBI程序是否正确；
在用自编程序计算一个简单的例题，看看结果是不是符合实际；
这样再去做其他规模交大一些的算例就不会出问题了，俗话说：“磨刀不误砍柴功”... .. :'( 这个我设置过断点，雅克比程序是没有问题的。而起按照数学理论实对称矩阵只能保证特征值是实数。不能保证为正数 .
实对称矩阵的确不会出现复特征值，但如果振动系统不考虑阻尼，且系统为线弹性，那么对应的特征值一定是正值，开放就是系统的圆频率形式的固有频率。 唉，那我的问题正好是考虑了阻尼和进入了塑性阶段!现在该怎么解决呢？好痛苦啊 .
考虑阻尼那要看用的是什么阻尼形式，如果是质量矩阵和刚度矩阵的线性组合，那么也还是正的特征值；
问题中的变形进入到塑性，就成了非线性问题了，非线性的模态概念还得相应去研究、研究吧.. ..
有什么进展也告诉大家呦 .. .. 呵呵 我用的阻尼就是瑞雷阻尼 （质量矩阵和刚度矩阵的线性组合）进入非线性特征值有部分变为负的，另外还有一个就是 我的程序编译 连接 都没有错误但是 运行 却输不出值来，这不知道是怎么回事？还望朋友指教啊， 这个具体解决需要看那些方面的书呢？:@Q .
非线性动力有限元... . 得到一个好消息：出现复数的话就不能用JACOBI了,用其它算法如实赫申伯格矩阵的QR算法
也有一个坏消息：这种计算方法我现在暂时找不到，朋友门帮帮忙找找啊
[ 本帖最后由 狼跃冲顶 于 2007-1-30 11:29 编辑 ] .
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----------------------------------------------------------------------- 中华 同仁 再次谢谢你 有没有关于实赫申伯格矩阵的QR算法的相关理论？？我想先看看理论才敢用这个算法 原帖由 狼跃冲顶 于 2007-1-28 12:27 发表
唉，那我的问题正好是考虑了阻尼和进入了塑性阶段!现在该怎么解决呢？好痛苦啊
如果是塑性问题，现有的算法是求不出模态的 ，因为这是一个非线性的模态分析问题，传统的不管是QR法还是Jaccobi法，乃至其它什么方法，都没戏。 个人以为，
如果矩阵阶数不大的话，直接用fortran的imsl库函数求解，
看看结果是否正确，
不要在求解特征值的方法上纠缠，
如果库函数求解没问题，说明编程没啥问题，
再考虑使用的特征值求解方法。
探讨 另外，
关于非线性模态，在kjbathe的书上略有介绍，
不知楼主可否分享一下心得，
谢谢。
页:
[1]
2