给我一个matlab 悬臂梁有限元分析代码，要求计算结构的自然频率，总刚度，总质量矩阵，采用newmark 法进行动力学计算结构的动力响应。

好的，以下是一个简单的 Matlab 悬臂梁有限元分析代码，用于计算结构的自然频率、总刚度、总质量矩阵，并采用 Newmark 法进行动力学计算结构的动力响应。请注意，此代码仅供参考，可能需要根据您的具体问题进行修改和优化。

% 悬臂梁有限元分析代码

% 定义结构参数
L = 1; % 梁长度
H = 0.1; % 梁高度
B = 0.05; % 梁宽度
E = 2e11; % 弹性模量
rho = 7800; % 密度
A = B * H; % 梁截面积
I = B * H^3 / 12; % 梁转动惯量

% 定义有限元参数
ne = 10; % 元素个数
nx = ne + 1; % 节点个数
dx = L / ne; % 元素长度
m = rho * A * dx / 2 * [2, 1; 1, 2]; % 质量矩阵
k = E * I / dx^3 * [12, 6*dx, -12, 6*dx; 6*dx, 4*dx^2, -6*dx, 2*dx^2;...
    -12, -6*dx, 12, -6*dx; 6*dx, 2*dx^2, -6*dx, 4*dx^2]; % 刚度矩阵

% 组装总质量矩阵和总刚度矩阵
M = zeros(nx*2);
K = zeros(nx*2);
for i = 1:ne
    M((i-1)*2+1:i*2+2, (i-1)*2+1:i*2+2) = M((i-1)*2+1:i*2+2, (i-1)*2+1:i*2+2) + m;
    K((i-1)*2+1:i*2+2, (i-1)*2+1:i*2+2) = K((i-1)*2+1:i*2+2, (i-1)*2+1:i*2+2) + k;
end

% 求解结构的自然频率和振型
[V, D] = eig(K, M);
omega = sqrt(diag(D));
[freq, index] = sort(omega / (2*pi));
modes = V(:, index);

% 输出结构的自然频率和振型
fprintf('The first 3 natural frequencies of the beam are: %.2f, %.2f, %.2f Hz.\n', freq(1:3));
figure;
for i = 1:3
    subplot(3, 1, i);
    plot(linspace(0, L, nx), [0; modes(1:2:end, i); 0]);
    xlabel('x (m)');
    ylabel('Displacement (m)');
    title(sprintf('Mode %d: %.2f Hz', i, freq(i)));
end

% 定义求解动力响应的参数
dt = 0.001; % 时间步长
tmax = 1; % 最大计算时间
beta = 0.25; % Newmark 法中的参数
gamma = 0.5; % Newmark 法中的参数
C = beta * M + gamma * dt * K; % Newmark 法中的阻尼矩阵
u = zeros(nx*2, 1); % 初始位移
v = zeros(nx*2, 1); % 初始速度
a = zeros(nx*2, 1); % 初始加速度
f = zeros(nx*2, 1); % 外力向量

% 定义外力函数
function y = force(t)
    if t <= 0.1
        y = 1000 * t;
    elseif t <= 0.3
        y = 100;
    elseif t <= 0.5
        y = -1000 * t + 600;
    else
        y = 0;
    end
end

% 求解动力响应
t = 0:dt:tmax;
n = length(t);
u_history = zeros(nx*2, n);
for i = 1:n
    f(end-1) = force(t(i));
    a = inv(C) * (f - K*u - M*v - (1-beta)*M*a);
    v = v + (1-gamma)*dt*a;
    u = u + dt*v + dt^2/2 * ((1-2*beta)*a + 2*beta*a_prev);
    u_history(:, i) = u;
    a_prev = a;
end

% 输出动力响应结果
figure;
plot(t, u_history(end-1, :));
xlabel('Time (s)');
ylabel('Displacement (m)');
title('Dynamic response of the beam');

希望这个代码对您有所帮助！