求解压电梁的自由振动的matlab程序，需要将梁简化为欧拉伯努利梁模型

以下是使用欧拉伯努利梁模型求解压电梁自由振动的 MATLAB 程序：

% 参数设置
L = 1;           % 梁的长度
b = 0.01;        % 梁的宽度
h = 0.01;        % 梁的厚度
rho = 7800;      % 梁的密度
E = 2.1e11;      % 梁的杨氏模量
vp = 0.3;        % 梁的泊松比
d33 = 15.3e-12;  % 压电系数
V = 100;         % 激励电压
omega = 1000;    % 激励频率
n = 10;          % 梁的分段数

% 计算梁的惯性矩和刚度矩阵
I = b * h^3 / 12;
K = [12 6*L -12 6*L;
    6*L 4*L^2 -6*L 2*L^2;
    -12 -6*L 12 -6*L;
    6*L 2*L^2 -6*L 4*L^2];
K = E*I/L^3*K;

% 计算压电效应矩阵
P = [0 0 0 0;
    0 0 0 0;
    0 0 0 0;
    0 d33*h/2 0 -d33*h/2];
P = P*V;

% 拆分梁的分段数
x = linspace(0,L,n+1);
dx = L/n;

% 组装整个系统的刚度矩阵和质量矩阵
K_sys = zeros(2*n+2);
M_sys = zeros(2*n+2);
for i = 1:n
    K_sys(i:i+3,i:i+3) = K_sys(i:i+3,i:i+3) + K;
    M_sys(i:i+3,i:i+3) = M_sys(i:i+3,i:i+3) + rho*b*h*dx/6*[2 0 1 0;
        0 2 0 1;
        1 0 2 0;
        0 1 0 2];
    if i == n
        M_sys(i+2:end,i+2:end) = M_sys(i+2:end,i+2:end) + rho*b*h*dx/6*[1 0 2 0;
            0 1 0 2;
            2 0 4 0;
            0 2 0 4];
    end
    P_sys = zeros(2*n+2,1);
    P_sys(i+1:i+2) = P_sys(i+1:i+2) + P*[1; 1];
end
K_sys = K_sys + [zeros(2,n+2); P_sys'; zeros(2,n+2)];
M_sys = M_sys + [rho*b*h*dx/3*[1 0 0 0;
    0 1 0 0;
    0 0 2 0;
    0 0 0 2], zeros(4,2*n-4);
    zeros(2*n-4,4), rho*b*h*dx/3*eye(2*n-4)];

% 求解特征值和特征向量
[eigvec,eigval] = eig(K_sys,M_sys);

% 取前10个自由振动频率和模态形状
[omega_sort,ind] = sort(sqrt(diag(eigval)));
omega_sort = omega_sort(1:10);
eigvec_sort = eigvec(:,ind(1:10));

% 绘制模态形状
for i = 1:10
    figure;
    plot([0; x'; L],[0; eigvec_sort(2*i-1:2*i,1); 0]);
    xlabel('x');
    ylabel('y');
    title(['Mode Shape ' num2str(i) ', \omega = ' num2str(omega_sort(i)) ' Hz']);
end

该程序首先计算了梁的惯性矩和刚度矩阵，然后根据欧拉伯努利梁模型将梁分段，并组装整个系统的刚度矩阵和质量矩阵。接下来，程序计算了压电效应矩阵，并将其加入到整个系统的刚度矩阵中。最后，程序求解了整个系统的特征值和特征向量，并取前10个自由振动频率和模态形状进行绘制。

注意：该程序仅供参考，实际使用时需要根据具体情况进行修改和调整。