$
\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 元,我们几天内就能脱离亏损状态,并保证网站能在接下来的一整年里向所有读者继续免费提供优质内容。感谢您的支持。
—— 小时(项目创始人)