matlab求解多自由度体系强迫振动采用振型叠加法正则坐标求解代码实现

以下是MATLAB代码实现多自由度体系强迫振动采用振型叠加法正则坐标求解：

%定义系统参数
m = [1, 2, 3]; %质量矩阵
k = [10, 20, 30]; %刚度矩阵
c = [0.1, 0.2, 0.3]; %阻尼矩阵
f = @(t) [5*sin(t), 10*sin(2*t), 15*sin(3*t)]; %外力函数

%定义时间范围和步长
tspan = [0, 10];
dt = 0.01;
t = tspan(1):dt:tspan(2);

%定义初始状态和初始速度
q0 = [0, 0, 0];
v0 = [0, 0, 0];

%求解系统响应
options = odeset('RelTol',1e-6,'AbsTol',1e-6);
[t, x] = ode45(@(t, x) dynamic_equation(t, x, m, k, c, f(t)), t, [q0, v0], options);

%计算振型系数
U = zeros(length(t), length(m));
for i = 1:length(t)
for j = 1:length(m)
U(i,j) = x(i,j)/sqrt(m(j));
end
end
U = U';

%计算振型频率和阻尼比
[omega, zeta] = mode_parameters(U, t);

%计算正则坐标
q = zeros(length(t), length(m));
for i = 1:length(t)
for j = 1:length(m)
q(i,j) = dot(U(:,j), x(i,1:length(m))');
end
end

%绘制正则坐标随时间变化的图像
figure
plot(t, q)
xlabel('Time (s)')
ylabel('Displacement (m)')

%定义动力学方程
function dx = dynamic_equation(t, x, m, k, c, f)
n = length(m);
q = x(1:n);
v = x(n+1:2*n);
dx = zeros(2*n, 1);
dx(1:n) = v;
dx(n+1:2*n) = f(t)/m - k/m*q - c/m*v;
end

%计算振型频率和阻尼比
function [omega, zeta] = mode_parameters(U, t)
n = size(U, 1);
omega = zeros(size(U, 2), 1);
zeta = zeros(size(U, 2), 1);
for i = 1:size(U, 2)
[pks, locs] = findpeaks(U(i,:));
T = diff(t(locs));
omega(i) = 2*pi/mean(T);
zeta(i) = -log(pks(2)/pks(1))/sqrt((2*pi/mean(T))^2+log(pks(2)/pks(1))^2);
end
end