'*********************************************************************
'* Frequencies and eigenmodes, masses and modal stiffnesses of a *
'* Mass-Spring undamped system represented by Motion Equation: *
'* [M] . [X]" + [K] . [X] = [0} *
'* ----------------------------------------------------------------- *
'* EXPLANATION: *
'* Let us take the 3 D.O.F. Mass-Spring System below: *
'* *
'* \| k ------ 2k ------ k ------ *
'* \|--/\/\/\--| 2m |--/\/\/\--| m |--/\/\/\--| 3m | *
'* \| ------ ------ ------ *
'* *
'* | 2 0 0 | *
'* The Mass Matrix M is | 0 1 0 | assuming m=1. *
'* | 0 0 3 | *
'* *
'* | 3 -2 0 | *
'* The Stiffness Matrix K is | -2 3 -1 | assuming k=1. *
'* | 0 -1 1 | *
'* *
'* The resonance fresuencies are: *
'* F1 = w1/2pi = 0.05161 sqr(k/m) *
'* F2 = w2/2pi = 0.1431 sqr(k/m) *
'* F3 = w3/2pi = 0.3151 sqr(k/m) *
'* with the corresponding modal vectors: *
'* | 1 | | 1 | | 1 | *
'* phi1 = | 1.395 | phi2 = | 0.6914 | phi3 = | -2.420 | *
'* | 2.038 | |-0.4849 | | 0.2249 | *
'* *
'* (See program ndof01.bas). *
'* ----------------------------------------------------------------- *
'* SAMPLE RUN: *
'* *
'* Order of system: 3 *
'* *
'* M(1,1) = 2 *
'* M(1,2) = 0 *
'* M(1,3) = 0 *
'* M(2,2) = 1 *
'* M(2,3) = 0 *
'* M(3,3) = 3 *
'* *
'* K(1,1) = 3 *
'* K(1,2) = -2 *
'* K(1,3) = 0 *
'* K(2,2) = 3 *
'* K(2,3) = -1 *
'* K(3,3) = 1 *
'* *
'* Number of eigenvectors asked for: 3 *
'* *
'* What precision for eigenvalues: 1e-5 *
'* *
'* Maximum number of iterations: 100 *
'* *
'* *
'* Eigenvector #1, Convergence after 10 iterations. *
'* Eigenvector #2, Convergence after 9 iterations. *
'* Eigenvector #3, Convergence after 2 iterations. *
'* *
'* *
'* Eigenvalues: *
'* *
'* L(1) = .1051733215666635 *
'* L(2) = .8086120631010683 *
'* L(3) = 3.919453983079086 *
'* *
'* Pulsations: *
'* *
'* W(1) = .3243043656299796 *
'* W(2) = .8992285933515839 *
'* W(3) = 1.979761092424812 *
'* *
'* Frequencies: *
'* *
'* F(1) = 5.161464285629262D-02 *
'* F(2) = .1431166756014764 *
'* F(3) = .3150887640004196 *
'* *
'* Eigenvectors: *
'* *
'* E.V. 1 *
'* 1 *
'* 1.394826678433337 *
'* 2.037781999496389 *
'* *
'* E.V. 2 *
'* 1 *
'* .691383466647192 *
'* -.4849013087124237 *
'* *
'* E.V. 3 *
'* 1 *
'* -2.41959287715014 *
'* .2249034647250999 *
'* *
'* Modal Mass Matrix: *
'* *
'* 16.40321 -0.00001 0.00000 *
'* -0.00001 3.18340 -0.00003 *
'* 0.00000 -0.00003 8.00617 *
'* *
'* Modal Stiffness Matrix: *
'* *
'* 1.72517 0.00000 -0.00000 *
'* 0.00000 2.57413 0.00001 *
'* -0.00000 0.00001 31.38059 *
'* *
'* ----------------------------------------------------------------- *
'* REFERENCE: "Mécanique des vibrations linéaires By M. Lalanne, *
'* P. Berthier, J. Der Hagopian, Masson, Paris 1980" *
'* [BIBLI 16]. *
'* *
'* Quick Basic Release By J-P Moreau, Paris. *
'* (www.jpmoreau.fr) *
'*********************************************************************
'Program NDOF02
DEFDBL A-H, O-Z
DEFINT I-N
PI = 4# * ATN(1#)
CLS
PRINT
INPUT " Order of system: ", N
PRINT
DIM M(N, N) AS DOUBLE
DIM K(N, N) AS DOUBLE
DIM A(N, N), B(N, N), C(N, N), S(N, N)
DIM R(N, N), R1(N, N), R2(N, N)
'auxiliary vectors for subroutine 1000
DIM PC(N), PL(N), CS(N)
DIM M1 AS DOUBLE
FOR I = 1 TO N
FOR J = I TO N
PRINT " M("; I; ","; J; ") = "; : INPUT "", M(I, J)
NEXT J
NEXT I
PRINT
FOR I = 1 TO N
FOR J = I TO N
PRINT " K("; I; ","; J; ") = "; : INPUT "", K(I, J)
NEXT J
NEXT I
A1 = 0#: A2 = 0#
FOR I = 1 TO N
FOR J = I TO N
A1 = A1 + M(I, I)
A2 = A2 + K(I, I)
M(J, I) = M(I, J)
K(J, I) = K(I, J)
NEXT J
NEXT I
PRINT
INPUT " Number of eigenvectors asked for: ", N1
K8 = 1
DIM K1(N, N) AS DOUBLE
DIM D(N, N)
PRINT
INPUT " What precision for eigenvalues: ", P
P = P * P
PRINT
INPUT " Maximum number of iterations: ", N3
PRINT
A0 = A2 / A1 / 10#
FOR I = 1 TO N
FOR J = 1 TO N
K1(I, J) = K(I, J) + A0 * M(I, J)
NEXT J
NEXT I
'K1=INV(K1)
FOR I = 1 TO N
FOR J = 1 TO N
A(I, J) = K1(I, J)
NEXT J
NEXT I
M = 0: GOSUB 1000'A=INV(A)
FOR I = 1 TO N
FOR J = 1 TO N
K1(I, J) = A(I, J)
NEXT J
NEXT I
'D=K1 MPY M
FOR I = 1 TO N
FOR J = 1 TO N
A(I, J) = K1(I, J)
B(I, J) = M(I, J)
NEXT J
NEXT I
GOSUB 2000 'C=A MPY B
FOR I = 1 TO N
FOR J = 1 TO N
D(I, J) = C(I, J)
NEXT J
NEXT I
DIM V(N, N1), S1(N, N)
DIM L(N1) AS DOUBLE
'V=0
FOR I = 1 TO N
FOR J = 1 TO N1
V(I, J) = 0#
NEXT J
NEXT I
'S1=D
FOR I = 1 TO N
FOR J = 1 TO N
S1(I, J) = D(I, J)
NEXT J
NEXT I
DIM X(N), F(N), F1(N)
DIM T(N, N), T1(N, N)
800 FOR I = 1 TO N
X(I) = 1#
NEXT I
K9 = 0
820 K9 = K9 + 1
'F=S1 MPY X
FOR I = 1 TO N
SUM = 0#
FOR J = 1 TO N
SUM = SUM + S1(I, J) * X(J)
NEXT J
F(I) = SUM
NEXT I
W2 = 1# / F(1)
L(K8) = W2 - A0
FOR I = 1 TO N
F(I) = F(I) * W2
NEXT I
FOR I = 1 TO N
F1(I) = F(I) - X(I)
NEXT I
D1 = 0#: M1 = 0#
FOR I = 1 TO N
D1 = D1 + F1(I) * F1(I)
M1 = M1 + F(I) * F(I)
X(I) = F(I)
NEXT I
970 IF K9 - N3 <= 0 THEN GOTO 1070
PRINT " No convergence after"; K9; " Iterations."
STOP
1070 IF P - D1 / M1 < 0# THEN GOTO 820
PRINT " Eigenvector #"; K8; ", Convergence after"; K9; " iterations."
FOR I = 1 TO N
V(I, K8) = X(I)
NEXT I
IF K8 >= N1 THEN GOTO 1490 'print results
1160 N2 = N - K8
'T=TRN(V)
FOR I = 1 TO N
FOR J = 1 TO N1
T(J, I) = V(I, J)
NEXT J
NEXT I
'T1=T MPY M
FOR I = 1 TO N1
FOR J = 1 TO N
A(I, J) = T(I, J)
B(I, J) = M(I, J)
NEXT J
NEXT I
GOSUB 2000 'C=A MPY B
FOR I = 1 TO N1
FOR J = 1 TO N
T1(I, J) = C(I, J)
NEXT J
NEXT I
FOR I = 1 TO K8
FOR J = 1 TO N
IF J <= K8 THEN
R(I, J) = T1(I, J)
GOTO 1290
END IF
J1 = J - K8
R1(I, J1) = T1(I, J)
1290 NEXT J
NEXT I
'R=INV(R)
FOR I = 1 TO K8
FOR J = 1 TO K8
A(I, J) = R(I, J)
NEXT J
NEXT I
NOLD = N: N = K8: M = 0: GOSUB 1000'A=INV(A)
FOR I = 1 TO K8
FOR J = 1 TO K8
R(I, J) = A(I, J)
NEXT J
NEXT I
'Restore N
N = NOLD
'R2(N,N2)=R(N,K8) MPY R1(N,N2)
FOR I = 1 TO K8
FOR J = 1 TO N2
SUM = 0#
FOR K = 1 TO K8
SUM = SUM + R(I, K) * R1(K, J)
R2(I, J) = SUM
NEXT K
NEXT J
NEXT I
'S=Identity Matrix
FOR I = 1 TO N
FOR J = 1 TO N
IF J = I THEN
S(I, J) = 1#
ELSE
S(I, J) = 0#
END IF
NEXT J
NEXT I
FOR I = 1 TO K8
S(I, I) = 0
FOR J = K8 + 1 TO N
J1 = J - K8
S(I, J) = -R2(I, J1)
NEXT J
NEXT I
'S1=D MPY S
FOR I = 1 TO N
FOR J = 1 TO N
A(I, J) = D(I, J)
B(I, J) = S(I, J)
NEXT J
NEXT I
GOSUB 2000 'C=A MPY B
FOR I = 1 TO N
FOR J = 1 TO N
S1(I, J) = C(I, J)
NEXT J
NEXT I
K8 = K8 + 1
GOTO 800 'go to next eigenvalue/eigenvector
'results section
1490 PRINT
PRINT " Eigenvalues:"
PRINT
FOR I = 1 TO N1
PRINT " L("; I; ") = "; L(I)
NEXT I
PRINT
PRINT " Pulsations:"
PRINT
FOR I = 1 TO N1
PRINT " W("; I; ") = "; SQR(L(I))
NEXT I
PRINT
PRINT " Frequencies:"
PRINT
FOR I = 1 TO N1
PRINT " F("; I; ") = "; SQR(L(I)) / 2# / PI
NEXT I
PRINT
PRINT " Eigenvectors:"
PRINT
FOR J = 1 TO N1
PRINT
PRINT " E.V. "; J
FOR I = 1 TO N
PRINT " "; V(I, J)
NEXT I
NEXT J
'S1(N,N1)=M(N,N) MPY V(N,N1)
FOR I = 1 TO N
FOR J = 1 TO N1
SUM = 0#
FOR K = 1 TO N
SUM = SUM + M(I, K) * V(K, J)
S1(I, J) = SUM
NEXT K
NEXT J
NEXT I
'T1=TRN(V)
FOR I = 1 TO N
FOR J = 1 TO N1
T1(J, I) = V(I, J)
NEXT J
NEXT I
'R(N1,N1)=T1(N1,N) MPY S1(N,N1)
FOR I = 1 TO N1
FOR J = 1 TO N1
SUM = 0#
FOR K = 1 TO N
SUM = SUM + T1(I, K) * S1(K, J)
R(I, J) = SUM
NEXT K
NEXT J
NEXT I
'print modal mass matrix
PRINT
PRINT " Modal Mass Matrix:"
PRINT
FOR I = 1 TO N1
FOR J = 1 TO N1
PRINT USING "####.#####"; R(I, J);
NEXT J
PRINT
NEXT I
'S1(N,N1)=K(N,N) MPY V(N,N1)
FOR I = 1 TO N
FOR J = 1 TO N1
SUM = 0#
FOR K = 1 TO N
SUM = SUM + K(I, K) * V(K, J)
S1(I, J) = SUM
NEXT K
NEXT J
NEXT I
'R(N1,N1)=T1(N1,N) MPY S1(N,N1)
FOR I = 1 TO N1
FOR J = 1 TO N1
SUM = 0#
FOR K = 1 TO N
SUM = SUM + T1(I, K) * S1(K, J)
R(I, J) = SUM
NEXT K
NEXT J
NEXT I
'print modal stiffness matrix
PRINT
PRINT " Modal Stiffness Matrix:"
PRINT
FOR I = 1 TO N1
FOR J = 1 TO N1
PRINT USING "####.#####"; R(I, J);
NEXT J
PRINT
NEXT I
PRINT
END 'of main program
'*******************************************
'* SOLVING A LINEAR MATRIX SYSTEM AX = B *
'* with Gauss-Jordan method using full *
'* pivoting at each step. During the pro- *
'* cess, original A and B matrices are *
'* destroyed to spare storage location. *
'* --------------------------------------- *
'* INPUTS: A MATRIX N*N *
'* B MATRIX N*M *
'* --------------------------------------- *
'* OUTPUS: A INVERSE OF A N*N *
'* DET DETERMINANT OF A *
'* B SOLUTION MATRIX N*M *
'* --------------------------------------- *
'* NOTA - If M=0 inversion of A matrix *
'* only. *
'*******************************************
1000 'Subroutine MATINV
'Labels: 100, 200, 300
'Initializations
DET = 1#
FOR I = 1 TO N
PC(I) = 0#
PL(I) = 0#
CS(I) = 0#
NEXT I
'Main loop
FOR K = 1 TO N
'Searching greatest pivot
PV = A(K, K)
IK = K: JK = K
PAV = ABS(PV)
FOR I = K TO N
FOR J = K TO N
temp = ABS(A(I, J))
IF temp > PAV THEN
PV = A(I, J)
PAV = ABS(PV)
IK = I: JK = J
END IF
NEXT J
NEXT I
'Search terminated, the pivot is in location I=IK, J=JK.
'Memorizing pivot location:
PC(K) = JK
PL(K) = IK
'DETERMINANT DET is updated
'If DET=0, ERROR MESSAGE and STOP
'Machine dependant EPSMACH equals here 2E-16
IF IK <> K THEN
DET = -DET
END IF
IF JK <> K THEN
DET = -DET
END IF
DET = DET * PV
temp = ABS(DET)
IF temp < EPSMACH THEN
PRINT #2,
PRINT #2, " The determinant equals ZERO !!!"
END 'stop program
END IF
'POSITIONNING PIVOT IN K,K:}
IF IK <> K THEN
FOR I = 1 TO N
'EXCHANGE LINES IK and K of matrix A:}
TT = A(IK, I)
A(IK, I) = A(K, I)
A(K, I) = TT
NEXT I
END IF
IF M <> 0 THEN
FOR I = 1 TO M
TT = B(IK, I)
B(IK, I) = B(K, I)
B(K, I) = TT
NEXT I
END IF
'PIVOT is at correct line
IF JK <> K THEN
FOR I = 1 TO N
'EXCHANGE COLUMNS JK and K of matrix A:
TT = A(I, JK)
A(I, JK) = A(I, K)
A(I, K) = TT
NEXT I
END IF
'The PIVOT is at correct column.
'and is located in K,K.
'Column K of matrix A is stored in CS vector
'then column K is set to zero.
FOR I = 1 TO N
CS(I) = A(I, K)
A(I, K) = 0#
NEXT I
CS(K) = 0#
A(K, K) = 1#
'Line K of matrix A is modified:
temp = ABS(PV)
IF temp < EPSMACH THEN
PRINT #2,
PRINT #2, " PIVOT TOO SMALL - STOP"
END 'stop program
END IF
FOR I = 1 TO N
A(K, I) = A(K, I) / PV
NEXT I
IF M <> 0 THEN
FOR I = 1 TO M
B(K, I) = B(K, I) / PV
NEXT I
END IF
'Other lines of matrix A are modified:
FOR J = 1 TO N
IF J = K THEN
GOTO 100
END IF
FOR I = 1 TO N
'Line J of matrix A is modified:
A(J, I) = A(J, I) - CS(J) * A(K, I)
NEXT I
IF M <> 0 THEN
FOR I = 1 TO M
B(J, I) = B(J, I) - CS(J) * B(K, I)
NEXT I
END IF
100 NEXT J
'Line K is ready
NEXT K 'end K loop
'MATRIX A INVERSION IS DONE - REARRANGEMENT OF MATRIX A
'EXCHANGE LINES
FOR I = N TO 1 STEP -1
IK = INT(PC(I))
IF IK = I THEN
GOTO 200
END IF
'EXCHANGE LINES I and PC(I) of matrix A:
FOR J = 1 TO N
TT = A(I, J)
A(I, J) = A(IK, J)
A(IK, J) = TT
NEXT J
IF M <> 0 THEN
FOR J = 1 TO M
TT = B(I, J)
B(I, J) = B(IK, J)
B(IK, J) = TT
NEXT J
END IF
'NO MORE EXCHANGE is NECESSARY
'Go to next line.
200 NEXT I 'End of loop i=N to 1 step -1
'EXCHANGE COLUMNS
FOR J = N TO 1 STEP -1
JK = INT(PL(J))
IF JK = J THEN
GOTO 300
END IF
'EXCHANGE COLUMNS J ET PL(J) of matrix A:
FOR I = 1 TO N
TT = A(I, J)
A(I, J) = A(I, JK)
A(I, JK) = TT
NEXT I
'NO MORE EXCHANGE is NECESSARY
'Go to next column.}
300 NEXT J
'REARRANGEMENT TERMINATED
RETURN 'End of subroutine 1000
'*******************************************
'* MULTIPLICATION OF TWO SQUARE REAL *
'* MATRICES *
'* --------------------------------------- *
'* INPUTS: A MATRIX N*N *
'* B MATRIX N*N *
'* N INTEGER *
'* --------------------------------------- *
'* OUTPUTS: C MATRIX N*N PRODUCT A*B *
'* *
'*******************************************
2000 'Start multiplication
FOR I = 1 TO N
FOR J = 1 TO N
SUM = 0#
FOR K = 1 TO N
SUM = SUM + A(I, K) * B(K, J)
C(I, J) = SUM
NEXT K
NEXT J
NEXT I
RETURN
'end of file ndof02.bas