構造物の固有周期、固有ベクトル計算に使う、プログラムと入力データを
塩井君に作っていただきました。
感謝！
固有値解析プログラム
C EIGEN VALUE AND EIGEN VECTOR
C IMPLICIT REAL*8(A-H,O-Z)
C EIGEN VALUE AND EIGEN VECTOR
COMMON /A/ M,MN,II,TINTAV,FAI(70,70),WEIGT(70),OMEGA(70),Z(600)
&,CO1(70),CO2(70),CO3(70),CO4,CO5,AA(600),BB(70,600),CC(70,600),
&EXPO(600),SINN(600),B(70,70),BMI(110,70),SFI(110,70),WM(70),
&MEM,AFI(110,70)
DIMENSION YY(70,70),T(70,70),TX(110),EIGV(70)
OPEN(5,FILE='dataspecial2.DAT',STATUS='UNKNOWN')
OPEN(6,FILE='ANSspecial2.ANS',STATUS='UNKNOWN')
READ(5,1001) MP,IPP
1001 FORMAT(I2,I6)
MN=MP
M=5
READ(5,1002) (WM(I),I=1,M)
1002 FORMAT(5F10.0)
WRITE(6,1022) (WM(I),I=1,M)
1022 FORMAT(5F10.0)
DO 10 I=1,M
10 WEIGT(I)=WM(I)/980.
DO 20 I=1,M
READ (5,100) (B(I,J),J=1,5)
100 FORMAT (5F15.7)
WRITE (6,110) (B(I,J),J=1,5)
110 FORMAT (5E15.7)
20 CONTINUE
READ(5,1003) EE
1003 FORMAT(F10.5)
210 DO 231 J=1,M
DO 231 I=1,M
231 B(I,J)=B(I,J)*WM(J)
C*** EIGEN (M,MP,IPP,EE,WEIGT,B,YY,EIGV,T)
CALL EBERLN (M,B,EE,YY,IPP,EIGV,WEIGT)
C*** EBERLN(N,A,EP,T,ITMX,EIGV,WEIGT)
DO 213 IV=1,MP
DO 213 I=1,M
213 FAI(I,IV)=YY(I,IV)
DO 214 IV=1,MP
OMEGA(IV)=EIGV(IV)
214 TX(IV)=2.*3.14159/EIGV(IV)
2000 WRITE(6,2001)
2001 FORMAT(28HEIGEN VALUE AND EIGEN VECTOR)
WRITE(6,2002)
2002 FORMAT(7X,17HNATURAL FREQUENCY,10X,14HNATURAL PERIOD)
WRITE(6,2003)
2003 FORMAT(9X,12H(RADIAN/SEC),16X,5H(SEC)/)
DO 215 IV=1,MP
WRITE(6,2004) IV,OMEGA(IV),TX(IV)
2004 FORMAT(2X,I2,6X,F11.4,16X,F7.3)
215 CONTINUE
WRITE(6,2005)
2005 FORMAT(6X,15HNORMALIZED MODE//)
DO 216 IV=1,MP
WRITE(6,2006) (FAI(I,IV),I=1,M)
2006 FORMAT(8F13.5)
WRITE(6,2007)
2007 FORMAT(//)
216 CONTINUE
END
C
SUBROUTINE EBERLN(N,A,EP,T,ITMX,EIGV,WEIGT)
C
C 改良された固有値のサブルーチン ヤコビ法
C 1991.5.22
C
C EIGENPROBLEM SOLVER FOR GENERAL MATRIX BY JACOBI-TYPE TR.S.SL***
C
C IMPLICIT REAL*8(A-H,O-Z)
DIMENSION A(70,70),T(70,70),EIGV(70),CCC(70,70),YYK(70)
DIMENSION WEIGT(70)
C
C
N=5
DO 1 I=1,N
DO 1 J=1,N
1 CCC(I,J)=A(I,J)
MARK =0
LEFT =0
IRIGHT = 0
IF (ITMX) 10,20,10
10 IF (ITMX) 30,40,40
30 LEFT = 1
GO TO 50
40 IRIGHT = 1
50 DO 60 I=1,N
T(I,I) =1.0
IP1= I+1
DO 60 J=IP1,N
T(I,J)= 0.0
60 T(J,I)= 0.0
20 EPS =SQRT(EP)
NM1=N-1
DO 70 IT=1,50
IF (MARK) 80,80,90
90 ITMX =1-IT
GO TO 410
80 DO 100 I=1,NM1
AII = CCC(I,I)
IP1 = I+1
DO 110 J=IP1,N
AIJ = CCC(I,J)
AJI = CCC(J,I)
IF (ABS(AIJ+AJI) - EPS) 120,120,130
120 IF (ABS(AIJ-AJI) - EPS) 110,110,150
150 IF (ABS(AII-CCC(J,J))-EPS) 110,110,130
130 GO TO 140
110 CONTINUE
100 CONTINUE
ITMX =IT -1
GO TO 410
140 MARK = 1
DO 160 K=1,NM1
KP1 = K+1
DO 160 M=KP1,N
H = 0.0
G = 0.0
HJ = 0.0
YH = 0.0
DO 180 I=1,N
AIK = CCC(I,K)
AIM = CCC(I,M)
TE = AIK*AIK
TEE = AIM*AIM
YH = YH + TE - TEE
IF (I-K) 190,200,190
190 IF (I-M) 210,200,210
210 AKI = CCC(K,I)
AMI = CCC(M,I)
H = H + AKI * AMI - AIK * AIM
TEP = TE + AMI * AMI
TEM = TEE+ AKI * AKI
G = G + TEP + TEM
HJ = HJ - TEP + TEM
200 CONTINUE
180 CONTINUE
H = H + H
D = CCC(K,K) - CCC(M,M)
AKM = CCC(K,M)
AMK = CCC(M,K)
C = AKM + AMK
E = AKM - AMK
IF (ABS(C) - EP) 220,220,230
220 CX = 1.0
SX = 0.0
GO TO 240
230 COT2X = D/C
IF (COT2X) 250,260,260
250 SIG = -1.0
GO TO 270
260 SIG = 1.0
270 COTX = COT2X + (SIG * SQRT(1.0 + COT2X*COT2X))
SX = SIG / SQRT(1.0+COTX*COTX)
CX = SX * COTX
240 IF (YH) 280,290,290
280 TEM = CX
CX = SX
SX = -TEM
290 COS2X = CX * CX - SX * SX
SIN2X = 2.0 * SX * CX
D = D * COS2X + C * SIN2X
H = H * COS2X - HJ* SIN2X
DEN = G + 2.0 *( E*E + D*D)
TANHY = (E * D - H/2.0) / DEN
IF (ABS(TANHY) - EP) 300,300,310
300 CHY = 1.0
SHY = 0.0
GO TO 320
310 CHY = 1.0/SQRT(1.0 - TANHY * TANHY)
SHY = CHY * TANHY
320 C1 = CHY * CX - SHY * SX
C2 = CHY * CX + SHY * SX
S1 = CHY * SX + SHY * CX
S2 =-CHY * SX + SHY * CX
IF (ABS(S1) - EP) 340,340,350
340 IF (ABS(S2) - EP) 160,160,350
350 MARK = 0
DO 360 I=1,N
AKI = CCC(K,I)
AMI = CCC(M,I)
CCC(K,I) = C1 * AKI + S1 * AMI
CCC(M,I) = S2 * AKI + C2 * AMI
IF (LEFT) 360,360,380
380 TKI = T(K,I)
TMI = T(M,I)
T(K,I) = C1 * TKI + S1 * TMI
T(M,I) = S2 * TKI + C2 * TMI
360 CONTINUE
DO 370 I=1,N
AIK = CCC(I,K)
AIM = CCC(I,M)
CCC(I,K) = C2 * AIK - S2 * AIM
CCC(I,M) =-S1 * AIK + C1 * AIM
IF (IRIGHT) 370,370,400
400 TIK = T(I,K)
TIM = T(I,M)
T(I,K) = C2 * TIK - S2 * TIM
T(I,M) =-S1 * TIK + C1 * TIM
370 CONTINUE
160 CONTINUE
70 CONTINUE
ITMX = 50
410 CONTINUE
DO 2 I=1,N
2 EIGV(I)=CCC(I,I)
DO 3 I=1,N
3 EIGV(I)=SQRT(980./EIGV(I))
DO 6 I=1,N
YYK(I)=0.
DO 4 J=1,N
4 YYK(I)=YYK(I)+T(J,I)**2*WEIGT(J)
YYK(I)=SQRT(YYK(I))
DO 5 J=1,N
5 T(J,I)=T(J,I)/YYK(I)
6 CONTINUE
RETURN
END
入力データの例
10900000
2383275. 720000. 2217000. 720000. 2383275.
0.1789414E-05 -0.2235023E-06 -0.8047154E-07 0.1900006E-07 0.2613326E-09
-0.2235026E-06 0.1092144E-03 0.5850431E-05 -0.1381326E-05 -0.1899919E-07
-0.8047149E-07 0.5850427E-05 0.2476859E-04 -0.5850350E-05 -0.8046669E-07
0.1900007E-07 -0.1381328E-05 -0.5850349E-05 0.1092140E-03 0.2234957E-06
0.2613326E-09 -0.1899919E-07 -0.8046673E-07 0.2234956E-06 0.1789389E-05
1.0E-06
Q 0