matlab求解双自由度系统频率以及振型的实振型分析代码

以下是matlab求解双自由度系统频率以及振型的实振型分析代码：

% 定义系统参数
m1 = 1; % 质量1
m2 = 2; % 质量2
k1 = 3; % 劲度系数1
k2 = 4; % 劲度系数2
c1 = 0.1; % 阻尼系数1
c2 = 0.2; % 阻尼系数2

% 定义状态空间矩阵A,B,C,D
A = [0 1 0 0; -k1/m1 -c1/m1 k1/m1 c1/m1; 0 0 0 1; k1/m2 c1/m2 -(k1+k2)/m2 -(c1+c2)/m2];
B = [0; 1/m1; 0; 0];
C = [1 0 0 0; 0 0 1 0];
D = [0; 0];

% 求解系统特征值与特征向量
[V, D] = eig(A);

% 取特征值中实部为零的部分
omega = sqrt(diag(D).^2);
omega = omega(imag(omega)==0);

% 求解振型
phi = C*V;
phi = phi./sqrt(sum(phi.^2)); % 归一化处理

% 绘制实振型分析图
figure;
subplot(2,1,1);
plot(phi(:,1),'r','LineWidth',2);
hold on;
plot(phi(:,2),'b','LineWidth',2);
title('振型');
ylabel('位移');
legend('模态1','模态2');
grid on;

subplot(2,1,2);
plot(omega,'ro','LineWidth',2);
title('频率');
ylabel('频率(rad/s)');
grid on;

在上述代码中，首先定义了双自由度系统的参数，然后根据状态空间法建立系统的状态空间矩阵A、B、C、D。接着使用eig函数求解系统的特征值与特征向量，取特征值中实部为零的部分计算自振频率。然后，使用C*V计算系统的振型，并进行归一化处理。最后，绘制实振型分析图，包括振型和频率两个子图。