强迫振动matlab

### 使用MATLAB实现强迫振动的模拟与分析

对于强迫振动的研究，通常涉及外部周期力作用下的单自由度系统。这类系统的运动可以通过带有外加激励项的二阶线性常微分方程来描述[^3]。

#### 数学建模
假设存在一个由质量 \(m\)、弹簧刚度 \(k\) 和阻尼器组成的一维单自由度系统，在受到谐波形式的外界驱动力 \(F_0\sin(\omega t)\) 影响下工作，则该系统的动态行为遵循下面给出的动力学方程式：

\[ m \ddot{x} + c \dot{x} + kx = F_0\sin(\omega t), \]

其中 \(c\) 是粘性阻尼系数，\(x(t)\) 表示位移响应，而 \(\omega\) 则指代激振频率。为了简化表达并便于后续处理，引入无量纲变量定义如下：

- 自然角频率：\[\Omega_n=\sqrt{\frac{k}{m}};\]
- 阻尼比：\[\zeta=\frac{c}{2\sqrt{mk}}.\]

因此原方程可以重写成标准形式：

\[ \ddot{x}(t)+2\zeta\Omega_n\dot{x}(t)+\Omega^2_nx(t)=f_{ext}\sin(\omega t).\]

这里 \(f_{ext}=F_0/m\) 代表单位质量所受的有效外力幅值。

#### MATLAB 实现方法
在MATLAB环境中执行上述过程的一个简单例子展示如下：

% 参数设定
m = 1; % 质量 kg
k = 4*pi^2; % 弹簧刚度 N/m, 设定使得自然频率接近于1Hz
c = 0.5*2*sqrt(m*k); % 中等程度的阻尼
Fo = 1; % 外部力幅度 N
w_ext = sqrt(k/m)*1.2; % 不同于固有频率的驱动频率 rad/s

% 时间向量设置
Tmax = 20;
dt = Tmax / 1e3;
time = linspace(0,Tmax,floor(Tmax/dt));

% 初始化状态空间矩阵A,B,C,D用于ode求解器
A=[0 1;-k/m -c/m];
B=[0 Fo/m]';
C=eye(size(A));
D=zeros(length(B),length(C));

sys = ss(A,B,C,D);

% 计算稳态响应以及瞬态部分
[y,t,x]=lsim(sys,sin(w_ext*time)',time,[0;0]);

figure();
plot(time,y,'LineWidth',1.5);
xlabel('Time(s)');
ylabel('Displacement');
title(['Forced Vibration Response at Driving Frequency ',num2str(w_ext/(2*pi)),' Hz']);
grid on;

% 绘制相轨迹图（Lissajous图形）
figure();
plot(y(:,end-1),y(:,end),'r','LineWidth',1.5);
axis equal;
xlabel('Position');
ylabel('Velocity');
title('Phase Trajectory of Forced Oscillator');

这段脚本首先设置了物理参数，并构建了一个状态空间模型 `sys` 来表征给定条件下的机械系统。接着调用了内置函数 `lsim()` 对此非齐次ODE进行了数值积分运算，从而获得了随着时间演化的位移历史记录 `y` 及其对应的时刻序列 `t` 。最后还额外制作了一张相平面内的轨迹图表以便更直观地观察系统的长期演化趋势[^1]。