双自由度有阻尼系统不同路面激励下的响应matlab求解

双自由度有阻尼系统的运动方程可以表示为:

$$\begin{Bmatrix} m_1 & 0 \\ 0 & m_2 \end{Bmatrix} \begin{Bmatrix} \ddot{x_1} \\ \ddot{x_2} \end{Bmatrix} + \begin{Bmatrix} c_1+c_2 & -c_2 \\ -c_2 & c_2 \end{Bmatrix} \begin{Bmatrix} \dot{x_1} \\ \dot{x_2} \end{Bmatrix} + \begin{Bmatrix} k_1+k_2 & -k_2 \\ -k_2 & k_2 \end{Bmatrix} \begin{Bmatrix} x_1 \\ x_2 \end{Bmatrix} = \begin{Bmatrix} F_1 \\ F_2 \end{Bmatrix}$$

其中，$m_1$和$m_2$分别是两个质点的质量，$c_1$和$c_2$分别是两个阻尼系数，$k_1$和$k_2$分别是两个弹簧系数，$F_1$和$F_2$分别是两个外力。

我们可以利用MATLAB中的ode45函数求解双自由度有阻尼系统的响应。具体步骤如下：

1. 定义运动方程

将双自由度有阻尼系统的运动方程转换为一阶微分方程组，即：

$$\begin{aligned} \dot{x_1} &= v_1 \\ \dot{v_1} &= \frac{1}{m_1}(F_1-c_1v_1-c_2(v_1-v_2)-k_1x_1-k_2(x_1-x_2)) \\ \dot{x_2} &= v_2 \\ \dot{v_2} &= \frac{1}{m_2}(F_2-c_2(v_2-v_1)-k_2(x_2-x_1)) \end{aligned}$$

2. 定义ODE函数

利用MATLAB的函数句柄定义ODE函数，代码如下：

function dy = double_pendulum(t,y,m1,m2,c1,c2,k1,k2,F1,F2)
    dy = zeros(4,1);
    dy(1) = y(2);
    dy(2) = (F1-c1*y(2)-c2*(y(2)-y(4))-k1*y(1)-k2*(y(1)-y(3)))/m1;
    dy(3) = y(4);
    dy(4) = (F2-c2*(y(4)-y(2))-k2*(y(3)-y(1)))/m2;
end

3. 定义初始条件和时间范围

定义初始条件和时间范围，代码如下：

tspan = [0 10];
y0 = [0; 0; 0; 0];

4. 调用ode45函数求解

调用ode45函数求解双自由度有阻尼系统的响应，代码如下：

[t,y] = ode45(@(t,y) double_pendulum(t,y,m1,m2,c1,c2,k1,k2,F1,F2),tspan,y0);

其中，@(t,y) double_pendulum(t,y,m1,m2,c1,c2,k1,k2,F1,F2)表示将double_pendulum函数作为输入参数传递给ode45函数。

5. 绘制响应曲线

利用MATLAB的plot函数绘制双自由度有阻尼系统的响应曲线，代码如下：

figure;
plot(t,y(:,1),'r',t,y(:,3),'b');
legend('x1','x2');
xlabel('time');
ylabel('displacement');

其中，y(:,1)和y(:,3)分别表示质点1和质点2的位移响应。

完整代码如下：

function dy = double_pendulum(t,y,m1,m2,c1,c2,k1,k2,F1,F2)
    dy = zeros(4,1);
    dy(1) = y(2);
    dy(2) = (F1-c1*y(2)-c2*(y(2)-y(4))-k1*y(1)-k2*(y(1)-y(3)))/m1;
    dy(3) = y(4);
    dy(4) = (F2-c2*(y(4)-y(2))-k2*(y(3)-y(1)))/m2;
end

m1 = 1;
m2 = 2;
c1 = 0.1;
c2 = 0.2;
k1 = 1;
k2 = 2;
F1 = 0;
F2 = 0;

tspan = [0 10];
y0 = [0; 0; 0; 0];

[t,y] = ode45(@(t,y) double_pendulum(t,y,m1,m2,c1,c2,k1,k2,F1,F2),tspan,y0);

figure;
plot(t,y(:,1),'r',t,y(:,3),'b');
legend('x1','x2');
xlabel('time');
ylabel('displacement');