$
\newcommand{\I}{\mathrm{i}}
\newcommand{\E}{\mathrm{e}}
\newcommand{\bvec}[1]{\boldsymbol{\mathbf{#1}}}
\newcommand{\mat}[1]{\boldsymbol{\mathbf{#1}}}
\newcommand{\ten}[1]{\boldsymbol{\mathbf{#1}}}
\newcommand{\Nabla}{\boldsymbol{\nabla}}
\renewcommand{\Tr}{^{\mathrm{T}}}
\newcommand{\uvec}[1]{\hat{\boldsymbol{\mathbf{#1}}}}
\renewcommand{\pmat}[1]{\begin{pmatrix}#1\end{pmatrix}}
\newcommand{\ali}[1]{\begin{aligned}#1\end{aligned}}
\newcommand{\leftgroup}[1]{\left\{\begin{aligned}#1\end{aligned}\right.}
\newcommand{\vmat}[1]{\begin{vmatrix}#1\end{vmatrix}}
\newcommand{\Cj}{^*}
\newcommand{\Her}{^\dagger}
\newcommand{\Q}[1]{\hat #1}
\newcommand{\Qv}[1]{\uvec #1}
\newcommand{\sinc}{\operatorname{sinc}}
\newcommand{\Arctan}{\operatorname{Arctan}}
\newcommand{\erfi}{\operatorname{erfi}}
\newcommand{\Arctan}{\operatorname{Arctan}}
\newcommand{\Si}[1]{\,\mathrm{#1}}
\newcommand{\les}{\leqslant}
\newcommand{\ges}{\geqslant}
\newcommand{\qty}[1]{\left\{{#1}\right\}}
\newcommand{\qtyRound}[1]{\left({#1}\right)}
\newcommand{\qtySquare}[1]{\left[{#1}\right]}
\newcommand{\abs}[1]{\left\lvert{#1}\right\rvert}
\newcommand{\eval}[1]{\left.{#1}\right\rvert}
\newcommand{\comm}[2]{\left[{#1},{#2}\right]}
\newcommand{\commStar}[2]{[{#1},{#2}]}
\newcommand{\pb}[2]{\left\{{#1},{#2}\right\}}
\newcommand{\pbStar}[2]{\{{#1},{#2}\}}
\newcommand{\vdot}{\boldsymbol\cdot}
\newcommand{\cross}{\boldsymbol\times}
\newcommand{\grad}{\boldsymbol\nabla}
\newcommand{\div}{\boldsymbol{\nabla}\boldsymbol{\cdot}}
\newcommand{\curl}{\boldsymbol{\nabla}\boldsymbol{\times}}
\newcommand{\laplacian}{\boldsymbol{\nabla}^2}
\newcommand{\sinRound}[2][{}]{\sin^{#1}\left(#2\right)}
\newcommand{\cosRound}[2][{}]{\cos^{#1}\left(#2\right)}
\newcommand{\tanRound}[2][{}]{\tan^{#1}\left(#2\right)}
\newcommand{\cscRound}[2][{}]{\csc^{#1}\left(#2\right)}
\newcommand{\secRound}[2][{}]{\sec^{#1}\left(#2\right)}
\newcommand{\cotRound}[2][{}]{\cot^{#1}\left(#2\right)}
\newcommand{\sinhRound}[2][{}]{\sinh^{#1}\left(#2\right)}
\newcommand{\coshRound}[2][{}]{\cosh^{#1}\left(#2\right)}
\newcommand{\tanhRound}[2][{}]{\tanh^{#1}\left(#2\right)}
\newcommand{\arcsinRound}[2][{}]{\arcsin^{#1}\left(#2\right)}
\newcommand{\arccosRound}[2][{}]{\arccos^{#1}\left(#2\right)}
\newcommand{\arctanRound}[2][{}]{\arctan^{#1}\left(#2\right)}
\newcommand{\expRound}[1]{\exp\left(#1\right)}
\newcommand{\logRound}[2][{}]{\log^{#1}\left(#2\right)}
\newcommand{\lnRound}[2][{}]{\ln^{#1}\left(#2\right)}
\renewcommand{\Re}{\mathrm{Re}}
\renewcommand{\Im}{\mathrm{Im}}
\newcommand{\dd}[1][]{\,\mathrm{d}^{#1}}
\newcommand{\dv}[2][{}]{\frac{\mathrm{d}^{#1}}{\mathrm{d}{#2}^{#1}}}
\newcommand{\dvStar}[2][{}]{\mathrm{d}^{#1}/\mathrm{d}{#2}^{#1}}
\newcommand{\dvTwo}[3][{}]{\frac{\mathrm{d}^{#1}{#2}}{\mathrm{d}{#3}^{#1}}}
\newcommand{\dvStarTwo}[3][{}]{\mathrm{d}^{#1}{#2}/\mathrm{d}{#3}^{#1}}
\newcommand{\pdv}[2][{}]{\frac{\partial^{#1}}{\partial{#2}^{#1}}}
\newcommand{\pdvStar}[2][{}]{\partial^{#1}/\partial{#2}^{#1}}
\newcommand{\pdvTwo}[3][{}]{\frac{\partial^{#1}{#2}}{\partial{#3}^{#1}}}
\newcommand{\pdvStarTwo}[3][{}]{\partial^{#1}{#2}/\partial{#3}^{#1}}
\newcommand{\pdvThree}[3]{\frac{\partial^2{#1}}{\partial{#2}\partial{#3}}}
\newcommand{\pdvStarThree}[3]{\partial^2{#1}/\partial{#2}\partial{#3}}
\newcommand{\bra}[1]{\left\langle{#1}\right\rvert}
\newcommand{\braStar}[1]{\langle{#1}\rvert}
\newcommand{\ket}[1]{\left\lvert{#1}\right\rangle}
\newcommand{\ketStar}[1]{\lvert{#1}\rangle}
\newcommand{\braket}[1]{\left\langle{#1}\middle|{#1}\right\rangle}
\newcommand{\braketStar}[1]{\langle{#1}|{#1}\rangle}
\newcommand{\braketTwo}[2]{\left\langle{#1}\middle|{#2}\right\rangle}
\newcommand{\braketStarTwo}[2]{\langle{#1}|{#2}\rangle}
\newcommand{\ev}[1]{\left\langle{#1}\right\rangle}
\newcommand{\evStar}[1]{\langle{#1}\rangle}
\newcommand{\evTwo}[2]{\left\langle{#2}\middle|{#1}\middle|{#2}\right\rangle}
\newcommand{\evStarTwo}[2]{\langle{#2}|{#1}|{#2}\rangle}
\newcommand{\mel}[3]{\left\langle{#1}\middle|{#2}\middle|{#3}\right\rangle}
\newcommand{\melStar}[3]{\langle{#1}|{#2}|{#3}\rangle}
\newcommand{\order}[1]{\mathcal{O}\left(#1\right)}
\newcommand{\bmat}[1]{\begin{bmatrix}#1\end{bmatrix}}
\newcommand{\Bmat}[1]{\left\{\begin{matrix}#1\end{matrix}\right\}}
\newcommand{\sumint}[1]{\int\kern-1.4em\sum}
\newcommand{\Q}[1]{\hat{#1}}
\newcommand{\opn}{\operatorname}
\newcommand{\norm}[1]{\left\lVert{#1}\right\rVert}
$
简谐振子受迫运动的简单数值计算
预备知识 简谐振子受迫运动
， 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;
dt = T/Nstep;
y2 = zeros(step,1); y1 = y2; y = y2;
y(1) = 0; y1(1) = 0;
for ii = 2:step
y2(ii) = (-a*y1(ii)-k*y(ii)+2*sin(w*(ii*dt)))/m;
y(ii) = y(ii-1) + y1(ii-1)*dt;
y1(ii) = y1(ii-1) + y2(ii-1)*dt;
end
t=(0:step-1)*dt;
plot(t,y);
|
运行结果如图 1 ，可见开始时驱动力不断给弹簧振子补充能量，振幅变大． 当补充的功率等于消耗的功率时，弹簧做稳定振动．
图1：运行结果
致读者： 小时物理百科一直以来坚持所有内容免费且不做广告，这导致我们处于日渐严重的亏损状态。长此以往很可能会最终导致我们不得不选择商业化，例如大量广告，内容付费，会员制，甚至被收购。因此，我们鼓起勇气在此请求广大读者
热心捐款，使网站得以健康发展。如果看到这条信息的每位读者能慷慨捐助 10 元，我们几天内就能脱离亏损状态，并保证网站能在接下来的一整年里向所有读者继续免费提供优质内容。感谢您的支持。
—— 小时（项目创始人）
