求解四自由度，其中有一个转角自由度的有阻尼受迫振动系统的响应，固有模态和固有振型matlab源代码

以下是一个四自由度系统的阻尼受迫振动的matlab代码，其中包括了计算固有模态和固有振型的部分：

clear all;
clc;

% 定义系统参数
m1 = 1; m2 = 2; m3 = 3; m4 = 4;
k1 = 10; k2 = 20; k3 = 30;
c1 = 1; c2 = 2; c3 = 3;
f = 1;

% 定义系统矩阵
M = diag([m1, m2, m3, m4]);
K = [k1+k2, -k2, 0, 0;
     -k2, k2+k3, -k3, 0;
     0, -k3, k3, 0;
     0, 0, 0, 0];
C = [c1+c2, -c2, 0, 0;
     -c2, c2+c3, -c3, 0;
     0, -c3, c3, 0;
     0, 0, 0, 0];

% 计算固有频率和固有振型
[V, D] = eig(K, M);
wn = sqrt(diag(D));
phi = V;

% 定义初始条件
x0 = [0.1; 0.2; 0.3; 0];

% 求解系统响应
tspan = [0 10];
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);
[t, x] = ode45(@(t, x) sys(t, x, M, C, K, f), tspan, x0, options);

% 绘制系统响应
figure;
plot(t, x(:, 1), 'r-', t, x(:, 2), 'b-', t, x(:, 3), 'g-', t, x(:, 4), 'k-');
xlabel('时间');
ylabel('位移');
legend('x1', 'x2', 'x3', 'x4');

% 绘制固有振型
figure;
subplot(2, 2, 1);
plot(phi(:, 1), 'r-', 'LineWidth', 2);
title(['$\omega_n$=', num2str(wn(1))], 'Interpreter', 'latex');
subplot(2, 2, 2);
plot(phi(:, 2), 'b-', 'LineWidth', 2);
title(['$\omega_n$=', num2str(wn(2))], 'Interpreter', 'latex');
subplot(2, 2, 3);
plot(phi(:, 3), 'g-', 'LineWidth', 2);
title(['$\omega_n$=', num2str(wn(3))], 'Interpreter', 'latex');
subplot(2, 2, 4);
plot(phi(:, 4), 'k-', 'LineWidth', 2);
title(['$\omega_n$=', num2str(wn(4))], 'Interpreter', 'latex');

function dx = sys(t, x, M, C, K, f)
    dx = zeros(4, 1);
    dx(1) = x(2);
    dx(2) = -inv(M)*(K*x(1:4)+C*x(2:5))+inv(M)*[0; 0; 0; f*sin(2*pi*2*t)];
    dx(3) = x(4);
    dx(4) = -inv(M)*(K*x(1:4)+C*x(2:5))+inv(M)*[0; 0; 0; 0];
end

在这个代码中，固有频率和固有振型的计算使用了MATLAB内置的eig函数。计算系统响应使用了ODE45函数，该函数是MATLAB中用于求解常微分方程组的函数，详细使用方法可以参考MATLAB官方文档。