function [] = String_vib_2_L3 % vibration of a string %...................................................... % x_1=0 x_2 x_3 x_4 x_5= L % T <==*-----(1)----*-----(2)-----*==> T % Fixed U_1 U_2 U_3 U_4 U_5 Fixed % Connectivity : e i j k % 1 1 2 3 % L3 % 2 3 4 5 % L3 %...................................................... n_e = 2 ; n_g = 1 ; n_n = 3 ; n_m = 5 ; % constants n_i = n_n*n_g ; n_d = n_g*n_m ; % elem & system DOFs L = 2. ; T = 6e5 ; rho = 0.0234 ; % given m, N, kg/m X = [0 0.25 0.5 0.75 1]*L ; % node coordinates L_e = L / n_e ; % element length Free = [2,3,4] ; % free displacements % Constant quadratic L3 element square matrices K_e = T * [ 7, -8, 1, ; ... % stiffness L3 -8, 16, -8, ; ... % 1, -8, 7 ] / (3*L_e) ; M_e = rho *L_e*[4, 2, -1 ; ... % mass L3 2, 16, 2 ; ... % -1, 2, 4 ] / 30 ; % Assemble two 3x3 element sq matrices into a 5x5 K=zeros (n_d, n_d); M=zeros (n_d, n_d) ; % allocate for k = 1:n_e; % loop over elements ===> ===> ===> ===> rows = [1:n_n] + (k - 1)*(n_n - 1) ; % connectivity K (rows, rows) = K (rows, rows) + K_e ; % add to stiff M (rows, rows) = M (rows, rows) + M_e ; % add to mass end ; % for element k <== <=== <=== <=== <=== <=== <=== % Solve for eigenvalues and eigenvectors, general form [Vec, Diag] = eig(K(Free,Free), M(Free,Free)) ; % part % Vec = eigenvector cols, Diag = eigenvalue^2 on diag Omega = sqrt(real(abs(diag(Diag)))) ; % diagonal freqs fprintf ('String natural frequencies \n') % heading exact = pi() * sqrt(T/(rho*L^2)) ; % exact constant V_abs = max (abs (Vec)) ; % relative amplitude for k = 1:n_d-2;% loop over DOFs - two EBCs ==> ==> ==> omega = Omega(k) ; % natural freq true = k * exact ; % exact natural freq fprintf ('%i, FEA = %8.3e, Exact = % 8.3e \n', ... k, omega, true) % pretty print Vec (:, k) = Vec (:, k) / V_abs (k) ; % scale abs = 1 end ; % for all active DOF <== <== <== <== <== <== <== % Pad mode shape with end zero BC displacements Pad = [0, 0, 0; Vec; 0, 0, 0] ; % insert rows of zero fprintf (' Mode Shapes \n') ; % heading fprintf ('Node 1 2 3 \n');% heading for k = 1:n_m ; % loop over nodes ---> ---> ---> ---> fprintf ('%i %9.3e %9.3e %9.3e \n', ... k, Pad(k,:)) ; % relative shape end ; % for node k relative shape <--- <--- <--- <--- % end of String_vib_2_L3 ============================== % Running gives: % String natural frequencies % 1, FEA = 7.984e+03, Exact = 7.954e+03 % 2, FEA = 1.601e+04, Exact = 1.591e+04 % 3, FEA = 2.873e+04, Exact = 2.386e+04 % Mode Shapes % Node 1 2 3 % 1 0.000e+00 0.000e+00 0.000e+00 % 2 7.068e-01 -1.000e+00 4.068e-01 % 3 1.000e+00 1.473e-16 -1.000e+00 % 4 7.068e-01 1.000e+00 4.068e-01 % 5 0.000e+00 0.000e+00 0.000e+00