求一个单自由度齿轮副的matlab程序

Question

用变步长四阶龙格库塔法计算位移x1随参数啮合频率w的动态响应分叉图。有偿，各位会的朋友请联系下我。（最后一张图片是结果图）
其动力学模型如下：

其动力学方程如下

结果类似图

Answer 1

引自微信小程序：皆我百晓生：

以下是一个用MATLAB编写的求解单自由度齿轮副动力学方程的程序：

%% 声明变量
syms x1 x2 t w m k c F theta omega_0 c_p Err;
x = [x1; x2];  % 位移向量

%% 动力学参数
m = 1;  % 质量
k = 100;  % 刚度
c = 0.5;  % 阻尼
F = 1;  % 外部激励

%% 动力学方程
dx1_dt = x2;  % x1的一阶导数
dx2_dt = (F*cos(w*t) - k*x1 - c*x2) / m;  % x2的一阶导数
dx_dt = [dx1_dt; dx2_dt];  % 转换成向量形式

%% 龙格库塔法求解
t_start = 0;  % 起始时间
t_end = 50;  % 结束时间
dt_initial = 0.01;  % 初始步长
tolerance = 1e-6;  % 容忍度

% 定义龙格库塔法的函数
rk4 = @(h, t, x) h / 6 * (k1(t, x) + 2*k2(t, x) + 2*k3(t, x) + k4(t, x));
k1 = @(t, x) dx_dt;
k2 = @(t, x) subs(dx_dt, [t, x], [t + h/2, x + h/2*rk4(h, t, x)]);
k3 = @(t, x) subs(dx_dt, [t, x], [t + h/2, x + h/2*rk4(h, t, x)]);
k4 = @(t, x) subs(dx_dt, [t, x], [t + h, x + h*rk4(h, t, x)]);

t_values = t_start;  % 存储时间的数组
x1_values = x(1);  % 存储x1的数组
x2_values = x(2);  % 存储x2的数组

while t_start < t_end
    % 根据当前步长求解下一个时间步长
    h = dt_initial;
    if t_start + h > t_end
        h = t_end - t_start;
    end
    
    % 龙格库塔法求解位移
    x_next = x + rk4(h, t_start, x);
    x = x_next;
    
    % 更新时间和位移数组
    t_start = t_start + h;
    t_values = [t_values, t_start];
    x1_values = [x1_values, x(1)];
    x2_values = [x2_values, x(2)];
end

%% 绘制动态响应分叉图
figure;
plot(t_values, x1_values, 'k');
xlabel('时间');
ylabel('位移x1');
title('动态响应分叉图');

请注意，此程序仅提供了一种使用变步长四阶龙格库塔法解决单自由度齿轮副动力学方程的方法。在实际使用中，您可能需要根据具体的问题修改参数和方程，以获得准确和合适的结果。