用Matlab微分方程组求解二自由度振动系统

x1''(t) + 2ξ1ω1x1'(t) + ω1^2x1(t) + k1(x1(t) - x2(t)) = 0
x2''(t) + 2ξ2ω2x2'(t) + ω2^2x2(t) + k1(x2(t) - x1(t)) = 0

function xdot = vibrationSystem(t, x)
    % 系统参数
    ξ1 = 0.1;
    ξ2 = 0.2;
    ω1 = 1;
    ω2 = 2;
    k1 = 1;

    % 位移和速度
    x1 = x(1);
    x2 = x(2);
    v1 = x(3);
    v2 = x(4);

    % 计算加速度
    a1 = -2*ξ1*ω1*v1 - ω1^2*x1 + k1*(x1 - x2);
    a2 = -2*ξ2*ω2*v2 - ω2^2*x2 + k1*(x2 - x1);

    % 返回速度向量
    xdot = [v1; v2; a1; a2];
end

% 初始条件
x0 = [0; 0; 1; 0];  % 初始位移和速度

% 时间范围
tspan = [0 10];  % 时间范围从0到10

% 求解微分方程组
[t, x] = ode45(@vibrationSystem, tspan, x0);

% 绘制位移随时间的变化曲线
plot(t, x(:, 1), 'r', t, x(:, 2), 'b');
xlabel('时间');
ylabel('位移');
legend('质点1', '质点2');