/******************************************************************** * 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) = 0.105173 * * L(2) = 0.808612 * * L(3) = 3.919454 * * * * Pulsations: * * * * W(1) = 0.324304 * * W(2) = 0.899229 * * W(3) = 1.979761 * * * * Frequencies: * * * * F(1) = 0.051615 * * F(2) = 0.143117 * * F(3) = 0.315089 * * * * Eigenvectors: * * * * E.V. 1 * * 1.000000 * * 1.394827 * * 2.037782 * * * * E.V. 2 * * 1.000000 * * 0.691383 * * -0.484901 * * * * E.V. 3 * * 1.000000 * * -2.419593 * * 0.224903 * * * * 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]. * * * * C++ Release By J-P Moreau, Paris. * * (www.jpmoreau.fr) * ********************************************************************/ #include #include #define MACH_EPS 2e-16 #define NMAX 25 #define REAL double #define ABS fabs //Labels: e800,e820,e970,e1070,e1160,e1290,e1490. typedef double MAT[NMAX][NMAX]; typedef double VEC[NMAX]; MAT M, K; MAT A, B, C, D, K1, S, R, R1, R2, S1, T, T1, V; VEC L, X, F, F1; double A0,A1,A2,D1,DET,M1,P,PI,SUM,TMP,W2; int I,J,J1,K0,K8,K9,N,N1,N2,N3; /****************************************** * SOLVING A LINEAR MATRIX SYSTEM AX = B * * with Gauss-Jordan method using full * * pivoting at each step. During the pro- * * cess, original AA and BB matrices are * * destroyed to spare storage location. * * --------------------------------------- * * INPUTS: AA MATRIX N*N * * BB MATRIX N*M * * --------------------------------------- * * OUTPUS: AA INVERSE OF AA N*N * * DET DETERMINANT OF AA * * BB SOLUTION MATRIX N*M * * --------------------------------------- * * NOTA - If M=0 inversion of AA matrix * * only. * ******************************************/ void MATINV(int N,int M,MAT AA,MAT BB,REAL *DET) { int PC[NMAX],PL[NMAX]; REAL CS[NMAX]; REAL PV,PAV,temp,TT; int I,IK,J,JK,K; // Initializations : *DET= 1.0; for (I=1; I<=N; I++) { PC[I]= 0; PL[I]= 0; CS[I]= 0.0; } // Main loop for (K=1; K<=N; K++) { // Searching greatest pivot : PV=AA[K][K]; IK=K; JK=K; PAV= ABS(PV); for (I=K; I<=N; I++) for (J=K; J<=N; J++) { temp = ABS(AA[I][J]); if (temp > PAV) { PV=AA[I][J]; PAV= ABS(PV); IK=I; JK=J; } } // Search terminated, the pivot is in location I=IK, J=JK. // Memorizing pivot location: PC[K]=JK; PL[K]=IK; // DETERMINANT DET is actualised // If DET=0, ERROR MESSAGE and STOP // Machine dependant EPSMACH equals here 1e-20 if (IK!=K) *DET=-*DET; if (JK!=K) *DET=-*DET; *DET=*DET*PV; temp= ABS(*DET); if (temp < MACH_EPS) { printf("\n The determinant equals ZERO !!!\n"); return; } // POSITIONNING PIVOT IN K,K: if(IK!=K) for (I=1; I<=N; I++) { // EXCHANGE LINES IK and K of matrix AA: TT=AA[IK][I]; AA[IK][I]=AA[K][I]; AA[K][I]=TT; } if(M!=0) for (I=1; I<=M; I++) { TT=BB[IK][I]; BB[IK][I]=BB[K][I]; BB[K][I]=TT; } // PIVOT is at correct line if(JK!=K) for (I=1; I<=N; I++) { // EXCHANGE COLUMNS JK and K of matrix AA: TT=AA[I][JK]; AA[I][JK]=AA[I][K]; AA[I][K]=TT; } // The PIVOT is at correct column. // and is located in K,K. // Column K of matrix AA is stored in CS vector // then column K is set to zero. for (I=1; I<=N; I++) { CS[I]=AA[I][K]; AA[I][K]= 0.0; } CS[K]= 0.0; AA[K][K]= 1.0; // Line K of matrix AA is modified: temp= ABS(PV); if(temp < MACH_EPS) { printf("\n PIVOT TOO SMALL - STOP\n"); return; } for (I=1; I<=N; I++) AA[K][I]=AA[K][I]/PV; if (M!=0) for (I=1; I<=M; I++) BB[K][I]=BB[K][I]/PV; // Other lines of matrix AA are modified: for (J=1; J<=N; J++) { if (J==K) J++; for (I=1; I<=N; I++) // Line J of matrix AA is modified: AA[J][I]=AA[J][I]-CS[J]*AA[K][I]; if (M!=0) for (I=1; I<=M; I++) BB[J][I]=BB[J][I]-CS[J]*BB[K][I]; } // Line K is ready } //End of K loop // MATRIX AA INVERSION IS DONE - REARRANGEMENT OF MATRIX AA // EXCHANGE LINES for (I=N; I>0; I--) { IK=PC[I]; if (IK==I) goto fin1; // EXCHANGE LINES I and PC(I) of matrix AA: for (J=1; J<=N; J++) { TT=AA[I][J]; AA[I][J]=AA[IK][J]; AA[IK][J]=TT; } if (M!=0) for (J=1; J<=M; J++) { TT=BB[I][J]; BB[I][J]=BB[IK][J]; BB[IK][J]=TT; } // NO MORE EXCHANGE is NECESSARY // Go to next line fin1: ; } // for i // EXCHANGE COLUMNS: for (J=N; J>0; J--) { JK=PL[J]; if (JK==J) goto fin2; // EXCHANGE COLUMNS J ET PL(J) of matrix AA: for (I=1; I<=N; I++) { TT=AA[I][J]; AA[I][J]=AA[I][JK]; AA[I][JK]=TT; } fin2: ; // NO MORE EXCHANGE is NECESSARY // Go to next column. } // REARRANGEMENT TERMINATED. return; } //Matinv() /****************************************** * 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 * * * ******************************************/ void MATMUL(MAT A,MAT B, MAT C, int N) { REAL SUM; int I,J,K; for (I=1; I<=N; I++) for (J=1; J<=N; J++) { SUM= 0.0; for (K=1; K<=N; K++) SUM=SUM+A[I][K]*B[K][J]; C[I][J]=SUM; } return; } void main() { PI = 4*atan(1); printf("\n"); printf(" Order of system: "); scanf("%d",&N); printf("\n"); for (I=1; I<=N; I++) for (J=I; J<=N; J++) { printf(" M(%d,%d) = ", I,J); scanf("%lf", &TMP); M[I][J]=TMP; } printf("\n"); for (I=1; I<=N; I++) for (J=I; J<=N; J++) { printf(" K(%d,%d) = ", I,J); scanf("%lf", &TMP); K[I][J]=TMP; } A1 = 0.0; A2 = 0.0; for (I=1; I<=N; I++) for (J=I; J<=N; J++) { A1 = A1 + M[I][I]; A2 = A2 + K[I][I]; M[J][I] = M[I][J]; K[J][I] = K[I][J]; } printf("\n"); printf(" Number of eigenvectors asked for: "); scanf("%d", &N1); K8 = 1; printf("\n"); printf(" What precision for eigenvalues: "); scanf("%lf", &P); P = P * P; printf("\n"); printf(" Maximum number of iterations: "); scanf("%d", &N3); printf("\n"); A0 = A2 / A1 / 10.0; for (I=1; I<=N; I++) for (J=1; J<=N; J++) K1[I][J] = K[I][J] + A0 * M[I][J]; //K1=INV(K1) MATINV(N,0,K1,S,&DET); //D=K1 MPY M MATMUL(K1, M, D, N); //V=0 for (I=1; I<=N; I++) for (J=1; J<=N1; J++) V[I][J] = 0.0; //S1=D for (I=1; I<=N; I++) for (J=1; J<=N; J++) S1[I][J] = D[I][J]; e800: for (I=1; I<=N; I++) X[I] = 1.0; K9 = 0; e820: K9++; //F=S1 MPY X for (I=1; I<=N; I++) { SUM = 0.0; for (J=1; J<=N; J++) SUM += S1[I][J] * X[J]; F[I] = SUM; } W2 = 1.0 / F[1]; L[K8] = W2 - A0; for (I=1; I<=N; I++) F[I] = F[I] * W2; for (I=1; I<=N; I++) F1[I] = F[I] - X[I]; D1 = 0.0; M1 = 0.0; for (I=1; I<=N; I++) { D1 += F1[I] * F1[I]; M1 += F[I] * F[I]; X[I] = F[I]; } if (K9 - N3 <= 0) goto e1070; printf(" No convergence after %d iterations.\n, K9"); return; e1070: if (P - D1 / M1 < 0.0) goto e820; printf(" Eigenvector #%d, Convergence after %d iterations.\n", K8, K9); for (I=1; I<=N; I++) V[I][K8] = X[I]; if (K8 >= N1) goto e1490; //print results N2 = N - K8; //T=TRN(V) for (I=1; I<=N; I++) for (J=1; J<=N1; J++) T[J][I] = V[I][J]; //T1=T MPY M MATMUL(T, M, T1, N); for (I=1; I<=K8; I++) for (J=1; J<=N; J++) { if (J <= K8) { R[I][J] = T1[I][J]; goto e1290; } J1 = J - K8; R1[I][J1] = T1[I][J]; e1290:;} //R=INV(R) MATINV(K8,0,R,S,&DET); //R2(N,N2)=R(N,K8) MPY R1(N,N2) for (I=1; I<=K8; I++) for (J=1; J<=N2; J++) { SUM = 0.0; for (K0=1; K0<=K8; K0++) { SUM += R[I][K0] * R1[K0][J]; R2[I][J] = SUM; } } //S=Identity Matrix for (I=1; I<=N; I++) for (J=1; J<=N; J++) if (J == I) S[I][J] = 1.0; else S[I][J] = 0.0; for (I=1; I<=K8; I++) { S[I][I] = 0.0; for (J=K8+1; J<=N; J++) { J1 = J - K8; S[I][J] = -R2[I][J1]; } } //S1=D MPY S MATMUL(D, S, S1, N); K8++; goto e800; //go to next eigenvalue/eigenvector //results section e1490: printf("\n"); printf(" Eigenvalues:\n\n"); for (I=1; I<=N1; I++) printf(" L(%d) = %f\n", I, L[I]); printf("\n"); printf(" Pulsations:\n\n"); for (I=1; I<=N1; I++) printf(" W(%d) = %f\n", I, sqrt(L[I])); printf("\n"); printf(" Frequencies:\n\n"); for (I=1; I<=N1; I++) printf(" F(%d) = %f\n", I, sqrt(L[I])/2.0/PI); printf("\n"); printf("\n"); printf(" Eigenvectors:"); printf("\n"); for (J=1; J<=N1; J++) { printf("\n"); printf(" E.V. %d\n", J); for (I=1; I<=N; I++) printf(" %f\n", V[I][J]); } printf("\n"); //S1(N,N1)=M(N,N) MPY V(N,N1) for (I=1; I<=N; I++) for (J=1; J<=N1; J++) { SUM = 0.0; for (K0=1; K0<=N; K0++) { SUM += M[I][K0] * V[K0][J]; S1[I][J] = SUM; } } //T1=TRN(V) for (I=1; I<=N; I++) for (J=1; J<=N1; J++) T1[J][I] = V[I][J]; //R(N1,N1)=T1(N1,N) MPY S1(N,N1) for (I=1; I<=N1; I++) for (J=1; J<=N1; J++) { SUM = 0.0; for (K0=1; K0<=N; K0++) { SUM += T1[I][K0] * S1[K0][J]; R[I][J] = SUM; } } //print modal mass matrix printf("\n"); printf(" Modal Mass Matrix:\n\n"); for (I=1; I<=N1; I++) { for (J=1; J<=N1; J++) printf("%10.5f", R[I][J]); printf("\n"); } //S1(N,N1)=K(N,N) MPY V(N,N1) for (I=1; I<=N; I++) for (J=1; J<=N1; J++) { SUM = 0.0; for (K0=1; K0<=N; K0++) { SUM += K[I][K0] * V[K0][J]; S1[I][J] = SUM; } } //R(N1,N1)=T1(N1,N) MPY S1(N,N1) for (I=1; I<=N1; I++) for (J=1; J<=N1; J++) { SUM = 0.0; for (K0=1; K0<=N; K0++) { SUM += T1[I][K0] * S1[K0][J]; R[I][J] = SUM; } } //print modal stiffness matrix printf("\n"); printf(" Modal Stiffness Matrix:\n\n"); for (I=1; I<=N1; I++) { for (J=1; J<=N1; J++) printf("%10.5f", R[I][J]); printf("\n"); } printf("\n"); } // end of file ndof02.cpp