傳遞矩陣法求簡支梁固有頻率的近似解 --matlab程序

%傳遞矩陣法求簡支梁固有頻率的近似解
clc
clear
syms p q
Sp = sym('[1 0 0 0;0 1 0 0;0 0 1 0;x 0 0 1]');        %點傳遞矩陣
Sf = sym('[1 1 1/2 1/6;0 1 1 1/2;0 0 1 1;0 0 0 1]');  %場傳遞矩陣
n = input('輸入划分單元數:');
S = ((Sf*Sp)^(n-1))*Sf;         %兩端支座之間的傳遞矩陣
%求固有頻率
xs = solve(S(1,2)*S(3,4)-S(1,4)*S(3,2));
xs = sort(double(xs));
xt = xs*n^4;
xt = sqrt(xt);
xe(1:n-1) = (pi*(1:n-1)).^2;    %精確解
xe = xe';
fprintf('傳遞矩陣法的結果：\n')
for i = 1:n-1
    fprintf('第%d階固有頻率：%8.4f(EI/ml^3)^(1/2)\n',i,xt(i))
end
%求模態
step = 1/n;
for i = 1:n-1
    f0 = -S(3,2)/S(3,4);
    f0 = subs(f0,'x',xs(i));
    xk(:,1) = [0 1 0 f0]';
    for j = 2:n+1
        xk(:,j) = Sf*Sp*xk(:,j-1);
        xk(:,j) = subs(xk(:,j),'x',xs(i));
    end
    xkk = xk(1,2:n);
    xkk = xkk/max(abs(xkk));
    xkk = double(xkk);
    xkk = real(xkk);
    if(xkk(1)<0)
        xkk = -xkk;
    end
    fprintf('第%d階模態為:',i)
    disp(xkk)
    figure()
    plot(0:step:1,[0 xkk 0].*abs(sin(i*pi*(0:step:1))),'ro')
    hold on
    xx = 0:pi/200:1;
    plot(xx,sin(i*pi*xx),'b')
end
fprintf('精確解的結果：\n')
for i = 1:n-1
    fprintf('第%d階固有頻率：%8.4f(EI/ml^3)^(1/2)\n',i,xe(i))
end

運行結果：

輸入划分單元數:5
傳遞矩陣法的結果：
第1階固有頻率：  9.8684(EI/ml^3)^(1/2)
第2階固有頻率： 39.3808(EI/ml^3)^(1/2)
第3階固有頻率： 87.1779(EI/ml^3)^(1/2)
第4階固有頻率：143.5557(EI/ml^3)^(1/2)
第1階模態為:    0.6180    1.0000    1.0000    0.6180

第2階模態為:    1.0000    0.6180   -0.6180   -1.0000

第3階模態為:    1.0000   -0.6180   -0.6180    1.0000

第4階模態為:    0.6180   -1.0000    1.0000   -0.6180

精確解的結果：
第1階固有頻率：  9.8696(EI/ml^3)^(1/2)
第2階固有頻率： 39.4784(EI/ml^3)^(1/2)
第3階固有頻率： 88.8264(EI/ml^3)^(1/2)
第4階固有頻率：157.9137(EI/ml^3)^(1/2)