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

                     

贡献者: addis; yuxuanwzy

预备知识 简谐振子受迫运动,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)$ 的函数图。

代码 1:SHOf.m
% 参数设定
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(Nstep,1); y1 = y2; y = y2; % 矩阵预赋值
y(1) = 0; y1(1) = 0; % 初值, y1 是 y 的一阶导数

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

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

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

图
图 1:运行结果

致读者: 小时百科一直以来坚持所有内容免费,这导致我们处于严重的亏损状态。 长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。 因此,我们请求广大读者热心打赏 ,使网站得以健康发展。 如果看到这条信息的每位读者能慷慨打赏 10 元,我们一个星期内就能脱离亏损, 并保证在接下来的一整年里向所有读者继续免费提供优质内容。 但遗憾的是只有不到 1% 的读者愿意捐款, 他们的付出帮助了 99% 的读者免费获取知识, 我们在此表示感谢。

                     

友情链接: 超理论坛 | ©小时科技 保留一切权利