VIB_Benaroya_11.21

目录

• C.1引言
• C.2单自由度无阻尼系统
• C.3单自由度阻尼系统
• C.4单自由度过阻尼系统
• C.5单自由度无阻尼系统受简谐激励
• C.6单自由度阻尼系统受简谐激励
• C.7单自由度阻尼系统受基础激励（续）
• C.8单自由度阻尼系统受旋转不平衡激励
• C.9单自由度阻尼系统受脉冲输入
• C.10单自由度阻尼系统受阶跃输入
• C.11单自由度阻尼系统受方波脉冲输入
• C.12单自由度阻尼系统受斜坡输入
• C.13单自由度系统受任意周期输入（续）
• C.14多自由度无阻尼系统
• C.15多自由度阻尼系统（续）
• C.16通用振动求解器（续）
• C.17范德波尔振子
• C.18随机振动
• C.19受随机激励的杜芬振子（续）
• C.20随机系统的蒙特卡洛模拟

Mechanical vibration analysis, uncertainties, and control (Benaroya, Haym Han, Seon Mi Nagurka, Mark L)

文档中提及的MATLAB程序在附录C中，主要用于解决书中涉及的振动问题，具体内容如下：

C.1引言

介绍了MATLAB在解决振动问题中的应用，提到附录C中的MATLAB程序可用于生成书中的一些结果。

C.2单自由度无阻尼系统

1. 计算自由振动响应
   • 给出了计算单自由度无阻尼系统自由振动响应的MATLAB代码，包括系统参数（质量、刚度、初始位移、初始速度）的设定，以及根据公式计算自然频率、周期、振幅和相位角，最后绘制位移响应曲线。
   • 示例代码如下：

% 单自由度无阻尼系统自由振动响应计算
m = 1; % 质量
k = 1; % 刚度
x0 = 1; % 初始位移
v0 = 0; % 初始速度

wn = sqrt(k/m); % 自然频率
T = 2*pi/wn; % 周期
A = sqrt(x0^2 + (v0/wn)^2); % 振幅
phi = atan2(v0,x0*wn); % 相位角

t = 0:0.01:T; % 时间向量
x = A*cos(wn*t - phi); % 位移响应

plot(t,x);
xlabel('时间 (s)');
ylabel('位移 (m)');
title('单自由度无阻尼系统自由振动响应');

C.3单自由度阻尼系统

1. 计算阻尼振动响应
   • 提供了计算单自由度阻尼系统振动响应的MATLAB程序，涵盖了阻尼比、自然频率、阻尼自然频率等参数的计算，以及根据阻尼情况（欠阻尼、临界阻尼、过阻尼）计算响应并绘制位移响应曲线。
   • 例如，欠阻尼情况下的代码如下：

% 单自由度阻尼系统振动响应计算（欠阻尼情况）
m = 1; % 质量
c = 0.2; % 阻尼系数
k = 1; % 刚度
x0 = 1; % 初始位移
v0 = 0; % 初始速度

wn = sqrt(k/m); % 自然频率
zeta = c/(2*m*wn); % 阻尼比
wd = wn*sqrt(1 - zeta^2); % 阻尼自然频率

if zeta < 1 % 欠阻尼情况
    C = sqrt((x0^2)*(wd^2) + (v0 + zeta*wn*x0)^2)/wd;
    phi = atan2((v0 + zeta*wn*x0),x0*wd);
    t = 0:0.01:10; % 时间向量
    x = C*exp(-zeta*wn*t).*cos(wd*t - phi); % 位移响应
    plot(t,x);
    xlabel('时间 (s)');
    ylabel('位移 (m)');
    title('单自由度欠阻尼系统振动响应');
end

C.4单自由度过阻尼系统

1. 计算过阻尼振动响应
   • 展示了计算单自由度过阻尼系统振动响应的MATLAB代码，涉及过阻尼系统的特征根计算，以及根据初始条件计算响应并绘制位移响应曲线。
   • 示例代码如下：

% 单自由度过阻尼系统振动响应计算
m = 1; % 质量
c = 3; % 阻尼系数
k = 1; % 刚度
x0 = 1; % 初始位移
v0 = 0; % 初始速度

wn = sqrt(k/m); % 自然频率
zeta = c/(2*m*wn); % 阻尼比
r1 = (-zeta + sqrt(zeta^2 - 1))*wn;
r2 = (-zeta - sqrt(zeta^2 - 1))*wn;

A1 = (v0 - r2*x0)/(r1 - r2);
A2 = x0 - A1;

t = 0:0.01:10; % 时间向量
x = A1*exp(r1*t) + A2*exp(r2*t); % 位移响应

plot(t,x);
xlabel('时间 (s)');
ylabel('位移 (m)');
title('单自由度过阻尼系统振动响应');

C.5单自由度无阻尼系统受简谐激励

1. 计算受简谐激励的响应
   • 给出了单自由度无阻尼系统受简谐激励时计算响应的MATLAB程序，包括设定激励频率、幅值、系统参数，计算静态位移、动态放大因子，根据公式计算响应并绘制位移响应曲线。
   • 示例代码如下：

% 单自由度无阻尼系统受简谐激励响应计算
m = 1; % 质量
k = 1; % 刚度
F0 = 1; % 激励幅值
omega = 0.5; % 激励频率
xst = F0/k; % 静态位移

wn = sqrt(k/m); % 自然频率
beta = omega/wn; % 频率比
D = xst/(1 - beta^2); % 动态放大因子

t = 0:0.01:10; % 时间向量
x = D*cos(omega*t); % 位移响应

plot(t,x);
xlabel('时间 (s)');
ylabel('位移 (m)');
title('单自由度无阻尼系统受简谐激励响应');

C.6单自由度阻尼系统受简谐激励

1. 计算受简谐激励的响应（考虑阻尼）
   • 提供了单自由度阻尼系统受简谐激励时计算响应的MATLAB代码，涵盖了阻尼比、自然频率、阻尼自然频率等参数的计算，以及根据公式计算响应（包括稳态响应和瞬态响应）并绘制位移响应曲线。
   • 例如，计算稳态响应的代码如下：

% 单自由度阻尼系统受简谐激励响应计算（考虑阻尼，稳态响应部分）
m = 1; % 质量
c = 0.2; % 阻尼系数
k = 1; % 刚度
F0 = 1; % 激励幅值
omega = 0.5; % 激励频率

wn = sqrt(k/m); % 自然频率
zeta = c/(2*m*wn); % 阻尼比
wd = wn*sqrt(1 - zeta^2); % 阻尼自然频率

beta = omega/wn; % 频率比
C = F0/k/sqrt((1 - beta^2)^2 + (2*zeta*beta)^2);
phi = atan2(2*zeta*beta,1 - beta^2);

t = 0:0.01:10; % 时间向量
x = C*cos(omega*t - phi); % 稳态位移响应

plot(t,x);
xlabel('时间 (s)');
ylabel('位移 (m)');
title('单自由度阻尼系统受简谐激励稳态响应');

C.7单自由度阻尼系统受基础激励（续）

1. 计算受基础激励的响应（完整代码）
   • 以下是完整的单自由度阻尼系统受基础激励时计算响应的MATLAB代码，包括计算传递函数、响应等步骤，并绘制位移响应曲线。

% 单自由度阻尼系统受基础激励响应计算
m = 1; % 质量
c = 0.2; % 阻尼系数
k = 1; % 刚度
omega = 0.5; % 激励频率（假设基础激励频率为0.5）
t = 0:0.01:10; % 时间向量

% 基础激励（假设为正弦位移激励）
y = 0.1*sin(omega*t); 
dy = omega*0.1*cos(omega*t); 
ddy = -omega^2*0.1*sin(omega*t);

wn = sqrt(k/m); % 自然频率
zeta = c/(2*m*wn); % 阻尼比
wd = wn*sqrt(1 - zeta^2); % 阻尼自然频率

% 计算传递函数
num = [k m];
den = [m 2*zeta*wn*m wn^2*m];
sys = tf(num,den);

% 计算响应
[yout,xout] = lsim(sys,ddy,t);

plot(t,xout);
xlabel('时间 (s)');
ylabel('位移 (m)');
title('单自由度阻尼系统受基础激励响应');

C.8单自由度阻尼系统受旋转不平衡激励

1. 计算受旋转不平衡激励的响应
   • 给出了单自由度阻尼系统受旋转不平衡激励时计算响应的MATLAB程序，包括设定旋转不平衡质量、偏心距、旋转速度等参数，系统参数（质量、刚度、阻尼系数），计算响应并绘制位移响应曲线。
   • 示例代码如下：

% 单自由度阻尼系统受旋转不平衡激励响应计算
m = 1; % 质量
c = 0.2; % 阻尼系数
k = 1; % 刚度
me = 0.1; % 旋转不平衡质量
e = 0.05; % 偏心距
omega = 0.5; % 旋转速度

wn = sqrt(k/m); % 自然频率
zeta = c/(2*m*wn); % 阻尼比
wd = wn*sqrt(1 - zeta^2); % 阻尼自然频率

A = me*e/m/sqrt((1 - (omega/wn)^2)^2 + (2*zeta*omega/wn)^2);
phi = atan2(2*zeta*omega/wn,1 - (omega/wn)^2);

t = 0:0.01:10; % 时间向量
x = A*sin(omega*t - phi); % 位移响应

plot(t,x);
xlabel('时间 (s)');
ylabel('位移 (m)');
title('单自由度阻尼系统受旋转不平衡激励响应');

C.9单自由度阻尼系统受脉冲输入

1. 计算受脉冲输入的响应
   • 提供了单自由度阻尼系统受脉冲输入时计算响应的MATLAB代码，包括设定脉冲力的参数（幅值、作用时间），系统参数（质量、刚度、阻尼系数），计算响应并绘制位移响应曲线。
   • 例如，以下是计算受脉冲输入响应的代码：

% 单自由度阻尼系统受脉冲输入响应计算
m = 1; % 质量
c = 0.2; % 阻尼系数
k = 1; % 刚度
F0 = 1; % 脉冲力幅值
t0 = 1; % 脉冲作用时间
t = 0:0.01:10; % 时间向量

wn = sqrt(k/m); % 自然频率
zeta = c/(2*m*wn); % 阻尼比
wd = wn*sqrt(1 - zeta^2); % 阻尼自然频率

u = zeros(size(t));
u(t == t0) = F0; % 脉冲输入

% 计算响应
num = [1];
den = [m c k];
sys = tf(num,den);
[yout,xout] = lsim(sys,u,t);

plot(t,xout);
xlabel('时间 (s)');
ylabel('位移 (m)');
title('单自由度阻尼系统受脉冲输入响应');

C.10单自由度阻尼系统受阶跃输入

1. 计算受阶跃输入的响应
   • 给出了单自由度阻尼系统受阶跃输入时计算响应的MATLAB程序，包括设定阶跃输入的幅值，系统参数（质量、刚度、阻尼系数），计算响应并绘制位移响应曲线。
   • 示例代码如下：

% 单自由度阻尼系统受阶跃输入响应计算
m = 1; % 质量
c = 0.2; % 阻尼系数
k = 1; % 刚度
F0 = 1; % 阶跃力幅值
t = 0:0.01:10; % 时间向量

wn = sqrt(k/m); % 自然频率
zeta = c/(2*m*wn); % 阻尼比
wd = wn*sqrt(1 - zeta^2); % 阻尼自然频率

u = F0*ones(size(t)); % 阶跃输入

% 计算响应
num = [1];
den = [m c k];
sys = tf(num,den);
[yout,xout] = lsim(sys,u,t);

plot(t,xout);
xlabel('时间 (s)');
ylabel('位移 (m)');
title('单自由度阻尼系统受阶跃输入响应');

C.11单自由度阻尼系统受方波脉冲输入

1. 计算受方波脉冲输入的响应
   • 展示了单自由度阻尼系统受方波脉冲输入时计算响应的MATLAB程序，包括设定方波脉冲的参数（幅值、周期、占空比），系统参数（质量、刚度、阻尼系数），计算响应并绘制位移响应曲线。
   • 示例代码如下：

% 单自由度阻尼系统受方波脉冲输入响应计算
m = 1; % 质量
c = 0.2; % 阻尼系数
k = 1; % 刚度
A = 1; % 方波脉冲幅值
T = 2; % 方波周期
duty = 0.5; % 占空比
t = 0:0.01:10; % 时间向量

wn = sqrt(k/m); % 自然频率
zeta = c/(2*m*wn); % 阻尼比
wd = wn*sqrt(1 - zeta^2); % 阻尼自然频率

y = A*square(2*pi*t/T,duty); % 方波脉冲输入

% 计算响应
num = [1];
den = [m c k];
sys = tf(num,den);
[yout,xout] = lsim(sys,y,t);

plot(t,xout);
xlabel('时间 (s)');
ylabel('位移 (m)');
title('单自由度阻尼系统受方波脉冲输入响应');

C.12单自由度阻尼系统受斜坡输入

1. 计算受斜坡输入的响应
   • 提供了单自由度阻尼系统受斜坡输入时计算响应的MATLAB代码，包括设定斜坡输入的斜率，系统参数（质量、刚度、阻尼系数），计算响应并绘制位移响应曲线。
   • 例如，以下是计算受斜坡输入响应的代码：

% 单自由度阻尼系统受斜坡输入响应计算
m = 1; % 质量
c = 0.2; % 阻尼系数
k = 1; % 刚度
k1 = 0.1; % 斜坡输入斜率
t = 0:0.01:10; % 时间向量

wn = sqrt(k/m); % 自然频率
zeta = c/(2*m*wn); % 阻尼比
wd = wn*sqrt(1 - zeta^2); % 阻尼自然频率

y = k1*t; % 斜坡输入

% 计算响应
num = [1];
den = [m c k];
sys = tf(num,den);
[yout,xout] = lsim(sys,y,t);

plot(t,xout);
xlabel('时间 (s)');
ylabel('位移 (m)');
title('单自由度阻尼系统受斜坡输入响应');

C.13单自由度系统受任意周期输入（续）

1. 计算受任意周期输入的响应（傅里叶级数法）
   • 介绍了使用傅里叶级数法计算单自由度系统受任意周期输入响应的MATLAB程序，包括设定周期函数的参数（周期、傅里叶级数项数），系统参数（质量、刚度、阻尼系数），计算傅里叶系数，构建输入函数，计算响应并绘制位移响应曲线。
   • 示例代码如下：

% 单自由度系统受任意周期输入响应计算（傅里叶级数法）
m = 1; % 质量
c = 0.2; % 阻尼系数
k = 1; % 刚度
T = 2*pi; % 周期函数的周期
N = 10; % 傅里叶级数项数
t = 0:0.01:10; % 时间向量

% 计算傅里叶系数（假设周期函数为f(t)）
a0 = 1/T * integral(@(t) f(t), 0, T);
an = zeros(1, N);
bn = zeros(1, N);
for n = 1:N
    an(n) = 2/T * integral(@(t) f(t).*cos(n*2*pi*t/T), 0, T);
    bn(n) = 2/T * integral(@(t) f(t).*sin(n*2*pi*t/T), 0, T);
end

% 构建输入函数
y = zeros(size(t));
for n = 1:N
    y = y + an(n)*cos(n*2*pi*t/T) + bn(n)*sin(n*2*pi*t/T);
end

wn = sqrt(k/m); % 自然频率
zeta = c/(2*m*wn); % 阻尼比
wd = wn*sqrt(1 - zeta^2); % 阻尼自然频率

% 计算响应
num = [1];
den = [m c k];
sys = tf(num,den);
[yout,xout] = lsim(sys,y,t);

plot(t,xout);
xlabel('时间 (s)');
ylabel('位移 (m)');
title('单自由度系统受任意周期输入响应');

C.14多自由度无阻尼系统

1. 计算无阻尼多自由度系统的响应
   • 给出了计算多自由度无阻尼系统响应的MATLAB程序，包括设定质量矩阵、刚度矩阵，计算固有频率和振型，根据初始条件计算响应并绘制位移响应曲线（以两自由度系统为例）。
   • 示例代码如下：

% 多自由度无阻尼系统响应计算（以两自由度系统为例）
m = [1 0; 0 1]; % 质量矩阵
k = [2 -1; -1 2]; % 刚度矩阵

% 计算固有频率和振型
[V,D] = eig(k,m);
wn = sqrt(diag(D));

% 初始条件
x0 = [1; 0];
v0 = [0; 0];

% 计算响应
t = 0:0.01:10;
phi = zeros(2,2);
for i = 1:2
    phi(:,i) = [V(1,i); V(2,i)] / sqrt(x0'*m*x0 + (v0'*m*V(:,i)/(wn(i)*V(1,i)))^2);
end
x = zeros(2,length(t));
for i = 1:2
    x(:,i) = phi(:,i).*cos(wn(i)*t - atan2(v0'*m*V(:,i)/(wn(i)*V(1,i)),x0'*m*x0));
end
x = x';

plot(t,x);
xlabel('时间 (s)');
ylabel('位移 (m)');
title('多自由度无阻尼系统响应');
legend('质量1位移','质量2位移');

C.15多自由度阻尼系统（续）

1. 计算阻尼多自由度系统的响应（完整代码）
   • 以下是完整的计算多自由度阻尼系统响应的MATLAB代码，包括设定质量矩阵、阻尼矩阵、刚度矩阵，计算相关参数，根据初始条件计算响应并绘制位移响应曲线（以两自由度系统为例）。

% 多自由度阻尼系统响应计算（以两自由度系统为例）
m = [1 0; 0 1]; % 质量矩阵
c = [0.2 0; 0 0.3]; % 阻尼矩阵
k = [2 -1; -1 2]; % 刚度矩阵

% 计算固有频率、阻尼比和振型
[V,D] = eig(k,m);
wn = sqrt(diag(D));
zeta = zeros(1,2);
for i = 1:2
    zeta(i) = c(i,i)/(2*m(i,i)*wn(i));
end

% 初始条件
x0 = [1; 0];
v0 = [0; 0];

% 计算响应
t = 0:0.01:10;
phi = zeros(2,2);
for i = 1:2
    phi(:,i) = [V(1,i); V(2,i)] / sqrt(x0'*m*x0 + (v0'*m*V(:,i)/(wn(i)*V(1,i)))^2);
end
x = zeros(2,length(t));
for i = 1:2
    if zeta(i) < 1 % 欠阻尼情况
        wd = wn(i)*sqrt(1 - zeta(i)^2);
        C = sqrt((x0'*m*x0)*(wd^2) + (v0'*m*V(:,i)/(wn(i)*V(1,i)) + zeta(i)*wn(i)*x0'*m*x0)^2)/wd;
        phi(:,i) = phi(:,i)*C;
        x(:,i) = phi(:,i).*exp(-zeta(i)*wn(i)*t).*cos(wd*t - atan2(v0'*m*V(:,i)/(wn(i)*V(1,i)) + zeta(i)*wn(i)*x0'*m*x0,x0'*m*x0));
    elseif zeta(i) == 1 % 临界阻尼情况
        x(:,i) = (x0'*m*x0 + (v0'*m*V(:,i)/(wn(i)*V(1,i)))*t).*exp(-wn(i)*t);
    else % 过阻尼情况
        r1 = (-zeta(i) + sqrt(zeta(i)^2 - 1))*wn(i);
        r2 = (-zeta(i) - sqrt(zeta(i)^2 - 1))*wn(i);
        A1 = (v0'*m*V(:,i)/(wn(i)*V(1,i)) - r2*x0'*m*x0)/(r1 - r2);
        A2 = x0'*m*x0 - A1;
        x(:,i) = (A1*exp(r1*t) + A2*exp(r2*t)).*[V(1,i); V(2,i)];
    end
end
x = x';

plot(t,x);
xlabel('时间 (s)');
ylabel('位移 (m)');
title('多自由度阻尼系统响应');
legend('质量1位移','质量2位移');

C.16通用振动求解器（续）

1. 介绍通用振动求解器（ode45函数）（完整代码）
   • 以下是使用MATLAB的ode45函数作为通用振动求解器的完整代码，以单自由度阻尼系统受简谐激励为例，包括设定系统参数、初始条件，定义微分方程函数，使用ode45求解并绘制位移响应曲线。

% 通用振动求解器（以单自由度阻尼系统受简谐激励为例）
m = 1; % 质量
c = 0.2; % 阻尼系数
k = 1; % 刚度
F0 = 1; % 激励幅值
omega = 0.5; % 激励频率

% 初始条件
x0 = 0;
v0 = 0;

% 定义微分方程函数
odefun = @(t,x) [x(2); (F0*cos(omega*t) - c*x(2) - k*x(1))/m];

% 时间向量
tspan = [0 10];

% 使用ode45求解
[t,x] = ode45(odefun,tspan,[x0;v0]);

% 绘制位移响应曲线
plot(t,x(:,1));
xlabel('时间 (s)');
ylabel('位移 (m)');
title('单自由度阻尼系统受简谐激励响应（ode45求解）');

C.17范德波尔振子

1. 计算范德波尔振子的响应
   • 给出了计算范德波尔振子响应的MATLAB程序，包括设定范德波尔振子的参数（(\mu)值），初始条件，定义微分方程函数，使用ode45求解并绘制位移和速度响应曲线。
   • 示例代码如下：

% 范德波尔振子响应计算
mu = 1; % 范德波尔振子参数

% 初始条件
x0 = 1;
v0 = 0;

% 定义微分方程函数
odefun = @(t,x) [x(2); mu*(1 - x(1)^2)*x(2) - x(1)];

% 时间向量
tspan = [0 20];

% 使用ode45求解
[t,x] = ode45(odefun,tspan,[x0;v0]);

% 绘制位移和速度响应曲线
subplot(2,1,1);
plot(t,x(:,1));
xlabel('时间 (s)');
ylabel('位移');
title('范德波尔振子响应');
subplot(2,1,2);
plot(t,x(:,2));
xlabel('时间 (s)');
ylabel('速度');

C.18随机振动

1. 生成随机振动信号并计算响应谱
   • 提供了生成随机振动信号并计算其响应谱的MATLAB代码，包括设定采样频率、信号长度，生成随机振动信号，计算自相关函数、功率谱密度，绘制信号波形、自相关函数曲线和功率谱密度曲线。
   • 例如，以下是生成随机振动信号并计算响应谱的部分代码：

% 随机振动信号生成与响应谱计算
fs = 1000; % 采样频率
N = 10000; % 信号长度
t = (0:N-1)/fs; % 时间向量

% 生成随机振动信号（假设为白噪声）
x = randn(1,N);

% 计算自相关函数
[Rxx, lags] = xcorr(x);

% 计算功率谱密度
[Pxx,f] = pwelch(x,[],[],[],fs);

% 绘制信号波形
subplot(3,1,1);
plot(t,x);
xlabel('时间 (s)');
ylabel('幅值');
title('随机振动信号');

% 绘制自相关函数曲线
subplot(3,1,2);
plot(lags,fliplr(Rxx));
xlabel('滞后');
ylabel('自相关函数');
title('自相关函数');

% 绘制功率谱密度曲线
subplot(3,1,3);
plot(f,Pxx);
xlabel('频率 (Hz)');
ylabel('功率谱密度');
title('功率谱密度');

C.19受随机激励的杜芬振子（续）

1. 计算受随机激励的杜芬振子的响应（完整代码）
   • 以下是计算受随机激励的杜芬振子响应的完整MATLAB代码，包括设定杜芬振子的参数、随机激励的参数，定义微分方程函数，使用ode45求解并绘制位移响应曲线。

% 受随机激励的杜芬振子响应计算
alpha = 1; % 杜芬振子参数
beta = 1;
delta = 0.2;
F0 = 1; % 激励幅值

% 初始条件
x0 = 0;
v0 = 0;

% 定义微分方程函数
odefun = @(t,x) [x(2); -delta*x(2) - alpha*x(1) - beta*x(1)^3 + F0*randn(1)];

% 时间向量
tspan = [0 100];

% 使用ode45求解
[t,x] = ode45(odefun,tspan,[x0;v0]);

% 绘制位移响应曲线
plot(t,x(:,1));
xlabel('时间 (s)');
ylabel('位移');
title('受随机激励的杜芬振子响应');

C.20随机系统的蒙特卡洛模拟

1. 进行随机系统的蒙特卡洛模拟
   • 给出了对随机系统进行蒙特卡洛模拟的MATLAB程序，包括设定模拟参数（样本数量、时间步长等），系统参数（质量、刚度、阻尼系数等），定义随机变量（如随机激励），循环进行模拟计算，统计响应的均值和标准差，绘制响应曲线（以单自由度系统受随机激励为例）。
   • 示例代码如下：

% 随机系统的蒙特卡洛模拟（以单自由度系统受随机激励为例）
N = 100; % 样本数量
dt = 0.01; % 时间步长
t = 0:dt:10; % 时间向量

m = 1; % 质量
c = 0.2; % 阻尼系数
k = 1; % 刚度

% 初始化响应矩阵
x_all = zeros(length(t),N);

for i = 1:N
    % 生成随机激励（假设为正态分布随机力）
    F = randn(size(t));
    
    % 初始条件
    x0 = 0;
    v0 = 0;
    
    % 定义微分方程函数
    odefun = @(t,x) [x(2); (F(t) - c*x(2) - k*x(1))/m];
    
    % 使用ode45求解
    [~,x] = ode45(odefun,t,[x0;v0]);
    
    % 存储响应
    x_all(:,i) = x(:,1);
end

% 计算响应的均值和标准差
x_mean = mean(x_all,2);
x_std = std(x_all,0,2);

% 绘制响应曲线
subplot(2,1,1);
plot(t,x_mean);
xlabel('时间 (s)');
ylabel('位移均值');
title('随机系统蒙特卡洛模拟结果');
subplot(2,1,2);
plot(t,x_std);
xlabel('时间 (s)');
ylabel('位移标准差');