如下图所示的单自由度系统: m=20kg, c=30 Ns/m, and k=3000 N/m.
1、应用数学、机械等理论知识推导出以下系统的数学模型：
其中：
x坐标如上图所示，取其稳态点为零点。
M为重物质量
K为弹簧弹性系数
C为阻尼系数
F(t)为施加在重物上的外激励力量
2、计算当外激励为周期性激励：F(t)=40cos(10t)(N),而所有的初值条件为零，即 x(0)=dx(0)/dt=0, 应用MATLAB计算系统的动态响应x(t)=?并绘制图形。（计算时间域为[0 5]秒）。c=30 Ns/m, and k=3000 N/m
% 定义系统参数
m = 20.0; % 质量，单位：kg
c = 30.0; % 阻尼系数，单位：Ns/m
k = 3000.0; % 弹簧刚度，单位：N/m
% 定义外部激励函数
F = @(t) 40*cos(10*t);
% 定义一阶方程系统的匿名函数
system_equation = @(t, y) [y(2); (F(t) - c*y(2) - k*y(1)) / m];
% 初始条件
y0 = [0; 0]; % 初始位移和初始速度
% 时间范围
tspan = [0 5];
% 使用ode45求解
[t, y] = ode45(system_equation, tspan, y0);
% 提取位移
x = y(:, 1);
figure;
plot(t, x);
xlabel('T');
title('动态响应');
3、当外激励为脉冲激励时, 计算系统的稳态响应x(t)=?并绘制图形。（计算时间域为[0 5]秒）。 IMPLUSE
c=30 Ns/m, and k=3000 N/m
即：考虑具有粘性阻尼的弹簧 - 质量系统在 t = 0 时受到一个单位冲量作用:
% 定义系统参数
m = 20.0; % 质量，单位：kg
c = 30.0; % 阻尼系数，单位：Ns/m
k = 3000.0; % 弹簧刚度，单位：N/m % 定义ODE的匿名函数
ode_system = @(t, y) [y(2); (-c*y(2) - k*y(1)) / m]; % 初始条件（假设脉冲导致了一个初始速度）
y0 = [0; 1]; % 初始位移为0，初始速度为1 m/s % 时间范围
tspan = [0 5]; % 使用ode45求解
[t, y] = ode45(ode_system, tspan, y0);
% 提取位移
x = y(:, 1);
% 绘制图形
figure;
plot(t, x);
xlabel('T');
title('脉冲激励响应');
4、当外激励为阶跃激励时, 计算系统的稳态响应x(t)=?并绘制图形。（计算时间域为[0 5]秒）。 STEP
即：F(t)=20 N
c=30 Ns/m, and k=3000 N/m
% 参数
m = 20.0; % 质量，单位：kg
c = 30.0; % 阻尼系数，单位：Ns/m
k = 3000.0; % 弹簧刚度，单位：N/m
F = 20; % 阶跃激励的大小，单位：N
% 定义ODE的匿名函数，考虑阶跃激励
ode_system = @(t, y) [y(2); (F * (t > 0) - c*y(2) - k*y(1)) / m];
% 初始条件
y0 = [0; 0]; % 初始位移和速度都为0
% 时间范围
tspan = [0 5];
% 使用ode45求解
[t, y] = ode45(ode_system, tspan, y0); % 提取位移
x = y(:, 1); % 绘制图形
figure;
plot(t, x);
xlabel('T');
title('阶跃激励响应');
5、当外激励为阶跃激励（即：F(t)=20 N）时, 改变不同弹簧弹性系数k、阻尼系数c的值（见下述条件），计算系统的稳态响应x(t)=?并绘制图形，分析结果。（计算时间域为[0 5]秒）。
条件：
c=10；15 ；20 ；25；30；35；40；Ns/m,
k=1000；2000；3000；4000 N/m
% 参数定义
F = 20; % 阶跃激励的大小，单位：N
m = 20.0; % 质量，单位：kg
tspan = [0 5]; % 时间范围
y0 = [0; 0]; % 初始条件：位移和速度都为0 % 定义ODE函数（不需要单独的文件，直接作为匿名函数）
ode_system = @(t, y, F, m, c, k) [y(2); (F - c*y(2) - k*y(1)) / m]; % 初始化绘图参数
figure; % 创建一个新的图形窗口
hold on; % 保持当前图形，以便在同一窗口上绘制多个曲线
grid on; % 打开网格
xlabel('Time (s)'); % x轴标签
ylabel('Displacement (m)'); % y轴标签
title('Response of the System to a Step Input'); % 图形标题 % 绘制不同k和c值的响应
colors = 'rbgcykm'; % 颜色向量，用于区分不同曲线
k_values = [1000, 2000, 3000, 4000]; % 弹簧弹性系数
c_values = [10, 15, 20, 25, 30, 35, 40]; % 阻尼系数 subplot(221)
k = 1000;
for i_c = 1:length(c_values)c = c_values(i_c);% 使用ode45求解ODE[t, y] = ode45(@(t, y) ode_system(t, y, F, m, c, k), tspan, y0);x = y(:, 1); % 提取位移% 绘制位移响应曲线plot(t, x, 'LineWidth', 2, 'Color', colors(mod(i_k + i_c, length(colors)) + 1)); % 使用不同的颜色和线宽hold on;
end
title("K=1000")
legend('c=10','c=15','c=20','c=25','c=30','c=35','c=40')subplot(222)
k = 2000;
for i_c = 1:length(c_values)c = c_values(i_c);% 使用ode45求解ODE[t, y] = ode45(@(t, y) ode_system(t, y, F, m, c, k), tspan, y0);x = y(:, 1); % 提取位移% 绘制位移响应曲线plot(t, x, 'LineWidth', 2, 'Color', colors(mod(i_k + i_c, length(colors)) + 1)); % 使用不同的颜色和线宽hold on;
end
title("K=2000")
legend('c=10','c=15','c=20','c=25','c=30','c=35','c=40')subplot(223)
k = 3000;
for i_c = 1:length(c_values)c = c_values(i_c);% 使用ode45求解ODE[t, y] = ode45(@(t, y) ode_system(t, y, F, m, c, k), tspan, y0);x = y(:, 1); % 提取位移% 绘制位移响应曲线plot(t, x, 'LineWidth', 2, 'Color', colors(mod(i_k + i_c, length(colors)) + 1)); % 使用不同的颜色和线宽hold on;
end
title("K=3000")
legend('c=10','c=15','c=20','c=25','c=30','c=35','c=40')subplot(224)
k = 4000;
for i_c = 1:length(c_values)c = c_values(i_c);% 使用ode45求解ODE[t, y] = ode45(@(t, y) ode_system(t, y, F, m, c, k), tspan, y0);x = y(:, 1); % 提取位移% 绘制位移响应曲线plot(t, x, 'LineWidth', 2, 'Color', colors(mod(i_k + i_c, length(colors)) + 1)); % 使用不同的颜色和线宽hold on;
end
title("K=4000")
legend('c=10','c=15','c=20','c=25','c=30','c=35','c=40')
