简谐振子受迫运动的简单数值计算

预备知识　简谐振子受迫运动， Matlab 编程基础

　　以下以弹簧振子的受迫运动为例，介绍一种解 $n$ 阶微分方程的简单方法，以便了解数值解微分方程的基本思想．但是这种方法误差较大，需要大量计算才能获得较精确的数值解．在实际运用中已有更复杂更成熟的算法（参考 MATLAB 常微分方程（组）数值解简介）．

　　在简谐振子受迫运动中，列出的二阶微分方程为

\begin{equation} m\ddot y = \alpha \dot y - ky + f(t) \end{equation}

若已知初值条件（可代入任意具体数值） $\dot y(0) = v_0$， $y(0) = y_0$，且已知驱动力 $f(t)$，此时可以把初值条件代入式 1 求出 $t = 0$ 时的加速度．

\begin{equation} \ddot y(0) = [- \alpha \dot y(0) - ky(0) + f(0)]/m \end{equation}

接下来的一小段微小的时间 $\Delta t$ 内（ $\Delta t$ 称为步长，步长越小误差越小），根据微分近似，可以算出 $t = \Delta t$ 时刻的状态．

\begin{gather} y(\Delta t) = y(0) + \Delta y \approx y(0) + \dot y(0) \Delta t\\ \dot y(\Delta t) = \dot y(0) + \Delta \dot y \approx \dot y(0) + \ddot y(0) \Delta t \end{gather}

微分近似在这里的物理意义是在 $\Delta t$ 内速度和加速度都近似为常数．把 $y(\Delta t)$ 和 $\dot y(\Delta t)$ 再次代入式 1

\begin{equation} \ddot y(\Delta t) = [- \alpha \dot y(\Delta t) - ky(\Delta t) + f(\Delta t)]/m \end{equation}

再次使用微分近似有

\begin{gather} y(2\Delta t) = y(\Delta t) + \Delta y \approx y(\Delta t) + \dot y(\Delta t) \Delta t\\ \dot y(2\Delta t) = \dot y(\Delta t) + \Delta \dot y \approx \dot y(\Delta t) + \ddot y(\Delta t) \Delta t \end{gather}

重复以上各步骤，就可以继续得到 $y(3\Delta t)$， $y(4\Delta t)$ 等的近似值．在 $y$-$t$ 图中把这些散点连接起来，就得到了 $y(t)$ 的函数图．

SHOf.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

% 参数设定
m = 0.1; % 质量
k = 1; % 劲度系数
a = 0.03; % 阻尼系数
T = 20; % 停止时间
Nstep = 10000; % 步数
A = 2; w = 3; % f(t)=A*sin(w*t);

dt = T/Nstep; % 计算步长
y2 = zeros(step,1); y1 = y2; y = y2; % 矩阵预赋值
y(1) = 0; y1(1) = 0; % 初值, y1 是 y 的一阶导数

% 迭代循环
for ii = 2:step
    y2(ii) = (-a*y1(ii)-k*y(ii)+2*sin(w*(ii*dt)))/m; % 代入微分方程求出 y''.
    y(ii) = y(ii-1) + y1(ii-1)*dt; % y 的微分近似
    y1(ii) = y1(ii-1) + y2(ii-1)*dt; % y' 的微分近似
end

% 画图
t=(0:step-1)*dt;
plot(t,y);

　　运行结果如图 1 ，可见开始时驱动力不断给弹簧振子补充能量，振幅变大．当补充的功率等于消耗的功率时，弹簧做稳定振动．

图1：运行结果

致读者：小时物理百科一直以来坚持所有内容免费且不做广告，这导致我们处于日渐严重的亏损状态。长此以往很可能会最终导致我们不得不选择商业化，例如大量广告，内容付费，会员制，甚至被收购。因此，我们鼓起勇气在此请求广大读者热心捐款，使网站得以健康发展。如果看到这条信息的每位读者能慷慨捐助 10 元，我们几天内就能脱离亏损状态，并保证网站能在接下来的一整年里向所有读者继续免费提供优质内容。感谢您的支持。