用MATLAB求二自由度非线性系统的强迫响应

二自由度非线性系统的强迫响应需要先建立系统的数学模型，然后使用MATLAB进行求解。以下是一个示例模型和求解过程：

模型：
考虑一个二自由度非线性系统，其动力学模型可以表示为：

$m_1\ddot{x}_1 + k_1x_1 + c_1(\dot{x}_1-\dot{x}_2) + g_1(x_1) = f_1(t)$

$m_2\ddot{x}_2 + k_2x_2 + c_2(\dot{x}_2-\dot{x}_1) + g_2(x_2) = 0$

其中，$m_1$和$m_2$分别为两个质点的质量，$k_1$和$k_2$分别为两个质点的弹性系数，$c_1$和$c_2$分别为两个质点的阻尼系数，$g_1(x_1)$和$g_2(x_2)$分别为两个质点的非线性阻尼力，$f_1(t)$为外力。

求解：
将模型转化为状态空间形式：

$\dot{\textbf{x}} = \textbf{Ax} + \textbf{Bu}$

$\textbf{y} = \textbf{Cx}$

其中，

$\textbf{x} = \begin{bmatrix} x_1 \\ \dot{x}_1 \\ x_2 \\ \dot{x}_2 \end{bmatrix}$

$\textbf{y} = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}$

$\textbf{u} = \begin{bmatrix} f_1(t) \\ 0 \end{bmatrix}$

$\textbf{A} = \begin{bmatrix} 0 & 1 & 0 & 0 \\ -\frac{k_1}{m_1} & -\frac{c_1}{m_1} & \frac{k_1}{m_1} & \frac{c_1}{m_1} \\ 0 & 0 & 0 & 1 \\ \frac{k_1}{m_2} & \frac{c_1}{m_2} & -\frac{k_1+k_2}{m_2} & -\frac{c_1+c_2}{m_2} \end{bmatrix}$

$\textbf{B} = \begin{bmatrix} 1/m_1 \\ 0 \\ 0 \\ 0 \end{bmatrix}$

$\textbf{C} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{bmatrix}$

根据以上模型和参数，可以使用MATLAB中的ode45函数求解系统的强迫响应。以下是一个示例代码：

% 定义系统参数
m1 = 1;
m2 = 2;
k1 = 10;
k2 = 20;
c1 = 0.1;
c2 = 0.2;
g1 = @(x) -0.1*x^3;
g2 = @(x) -0.2*x^3;

% 定义系统动力学模型
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 = [1/m1 0 0 0]';
C = [1 0 0 0;
    0 0 1 0];
f = @(t) 5*sin(2*pi*t); % 外力信号

% 定义强迫响应求解函数
function [t, y] = forced_response(f, tspan, x0, A, B, C, g1, g2)
    function dx = state_eqn(t, x)
        dx = A*x + B*f(t);
        dx(1) = dx(1) + g1(x(1));
        dx(3) = dx(3) + g2(x(3));
    end
    [t, y] = ode45(@state_eqn, tspan, x0);
    y = C*y';
end

% 求解强迫响应
tspan = [0 10];
x0 = [0 0 0 0]';
[t, y] = forced_response(f, tspan, x0, A, B, C, g1, g2);

% 绘制结果
subplot(2,1,1);
plot(t, y(1,:), 'b', t, y(2,:), 'r');
legend('x1', 'x2');
xlabel('Time (s)');
ylabel('Position (m)');
subplot(2,1,2);
plot(y(1,:), y(2,:));
xlabel('x1');
ylabel('x2');

运行以上代码，得到的结果如下所示：


从图中可以看出，系统的两个质点的位置随时间变化呈现出相位差，且存在一定的非线性响应。