matlab用中心差分法求解振动力学多自由度位移响应代码

以下是一个简单的二自由度系统的matlab代码示例，使用中心差分法求解其位移响应：

% 定义系统参数
m1 = 1; m2 = 2; % 质量
k1 = 2; k2 = 4; % 刚度
f = [0; 1]; % 外力

% 定义时间参数
t = 0:0.01:10; % 时间向量
dt = t(2) - t(1); % 时间间隔

% 初始化位移和速度向量
u = zeros(2, length(t)); % 位移向量
v = zeros(2, length(t)); % 速度向量

% 中心差分法求解位移响应
for i = 2:length(t)
    % 求解当前时刻的加速度
    M = [m1, 0; 0, m2];
    K = [k1+k2, -k2; -k2, k2];
    a = (f - K*u(:,i-1)) / M;
    
    % 利用中心差分法求解位移和速度
    u(:,i) = u(:,i-1) + v(:,i-1)*dt + 0.5*a*dt^2;
    v(:,i) = v(:,i-1) + 0.5*(a + (f - K*u(:,i)) / M)*dt;
end

% 绘制位移响应曲线
plot(t, u(1,:), 'r', t, u(2,:), 'b');
xlabel('Time (s)');
ylabel('Displacement (m)');
legend('Node 1', 'Node 2');

需要注意的是，这只是一个简单的示例代码，对于更复杂的多自由度系统需要根据具体情况进行调整。同时，由于中心差分法的数值误差可能会导致计算结果不够精确，需要进行一定的误差分析和修正。