计算物理第八次作业

物基一班 李云龙 2013301020065

背景

振动是自然界最常见的物理现象之一，在所有振动现象中，简谐振动(simple harmonic motion)是最简单的一类，我们可以给简谐振动施加适当的阻尼和驱动力，使简谐振动模式多样化，其中阻尼振动根据阻尼的大小，又可以分为过阻尼、欠阻尼和临界阻尼振动。

摘要

本文根据课本思路，分别使用欧拉方法和欧拉-克拉默方法分析了简谐振动过程，简单比较了这两种方法的不同，同时分析了单摆、阻尼摆、驱动摆的物理运动过程。

正文

对单摆受力分析如下图所示:

单摆仅受重力，受力关系式为：

$$F_\theta=-mgsin\theta\tag 1$$
假设单摆进行小角度振动，则$sin\theta\approx\theta$,
我们可以得到运动方程： $$\frac{d^2\theta}{dt^2}=-\frac{g}{l}\tag 2$$
该方程的解为: $$\theta=\theta_0sin(\Omega t+\phi)\tag 3$$
其中：$\Omega=\sqrt{g/l}$, $\theta_0$和$\phi$是与初始状态有关的常数。

• 欧拉方法：

由（2）式我们可以得到单摆运动的一阶方程组：

$$\frac{d\omega}{dt}=-\frac{g}{l}\theta\tag 4$$
$$\frac{d\theta}{dt}=\omega\tag 5$$
于是我们很容易得到欧拉分布计算方法：
$$\omega_{i+1}=\omega_i-(g/l)\theta_i\Delta t\tag 6$$
$$\theta_{i+1}=\theta_i+\omega_i\Delta t\tag 7$$
$$t_{i+1}=t_i+\Delta t\tag 8$$

• 欧拉-克拉默方法：

考虑到系统的总能量保持不变：

$$E=\frac{1}{2}ml^2\omega^2+mgl(1-cos\theta)\tag 9$$
根据小角近似，同时将欧拉方法的（6）（7）式带入到$\omega$,$\theta$中，得到：
$$E_{i+1}=E_i+\frac{1}{2}mgl(\omega_i^2+\frac{g}{l}\theta_i^2)(\Delta t)^2\tag{10}$$
可以得到欧拉-克拉默方法的的分步计算方法为：
$$\omega_{i+1}=\omega_i-(g/l)\theta_i\Delta t\tag{11}$$
$$\theta_{i+1}=\theta_i+\omega_{i+1}\Delta t\tag{12}$$
$$t_{i+1}=t_i+\Delta t\tag{13}$$
可见，它们的区别仅仅在于用欧拉法时我们是用上一个$\omega$和上一个$\theta$来得到新的$\omega$和$\theta$；而用欧拉-克拉默法则是用上一个$\omega$和上一个$\theta$来得到新的$\omega$，再利用新的$\omega$来得到新的$\theta$。

利用两种方法来获得单摆的摆角随时间的变化曲线。
实验代码如下：（simple_pendulum_tz.py）
我们假设摆长为l=1m，初始角速度$\omega_0=0rad/s$，初始角度为$\theta_0=10^\circ$,
时间间隔为$\Delta=0.04s$,得到曲线图：

时间间隔为$\Delta=0.01s$,得到曲线图：

由图我们可以看出，利用欧拉法得出的曲线与实际情况并不符合，摆的振幅及能量都随着时间增加而增大，这显然不正确,并且欧拉方法随着计算的时间间隔即步长越大，偏差会越明显。而利用欧拉-克拉默法得到的曲线则不存在这样的问题，可以看到，欧拉-克拉默方法计算的结果比较稳定，系统的能量始终保持不变，为稳定的正弦图像，并且不随步长的变化有明显的偏差。因此下面对摆的问题的处理我们都采用欧拉-克拉默方法。

• 阻尼摆：

对阻尼摆，运动方程可以写为：

$$\frac{d^2\theta}{dt^2}=-\frac{g}{l}-q\frac{d\theta}{dt}\tag{14}$$
欠阻尼情况下，解析解为：
$$\theta(t)=\theta_0e^{-qt/2}sin(\sqrt{\Omega^2-q^2/4}t+\Omega)\tag{15}$$
过阻尼情况下，解析解为：
$$\theta(t)=\theta_0e^{-（q/2\pm\sqrt{q^2/4-\Omega^2)}t}\tag{16}$$
临界阻尼情况下，解析解为：
$$\theta(t)=(\theta_0+Ct)e^{-qt/2}\tag{17}$$
利用欧拉-克拉默法可以得到三种情况下的曲线，
实验代码如下：（damped_pendulum_tz.py）
假设摆长为l=1m，初始角速度$\omega_0$=0rad/s，初始角度为$\theta_0=10^\circ$,时间间隔为$\Delta t$=0.01s，如图所示：

由实验结果结合理论推导，可以得到：
欠阻尼情况下，单摆以频率$\sqrt{\Omega^2-q^2/4}$摆动，并且振幅随时间增加逐渐减小。
过阻尼情况下，单摆摆角随时间增加而以指数形式单调递减。
临界阻尼情况下，单摆恰好不起振。

• 驱动摆：

在阻尼摆的基础上加上驱动项得到驱动摆，驱动摆运动方程形如：

$$\frac{d^2\theta}{dt^2}=-\frac{g}{l}-q\frac{d\theta}{dt}+F_Dsin(\Omega_D t)\tag{18}$$
其中$F_Dsin(\Omega_D t)$为驱动项。
下面研究驱动项对摆的影响。假设摆长为l=1m，初始角速度$\omega_0$=0rad/s，初始角度为$\theta_0=10^\circ$,时间间隔为$\Delta t$=0.01s，阻尼系数q=1.0。
先使$\Omega_D=2.0$保持不变，分别使$F_D$等于0.2,0.6,1.0，
实验代码如下：（driven_pendulum_tz）
得到曲线如下：

然后使$F_D=0.2$保持不变，分别使$\Omega_D$等于1.0,2.0,3.0，
得到曲线如下：

上两种情况有一共同点，摆在经历一段不稳定的情况后均以稳定的振幅及频率摆动，相当于单摆的运动。当我们保持$\Omega_D$不变改变$F_D$时，不同情况下摆在稳定摆动时的频率不变，振幅有明显差别；当我们保持$F_D$不变改变$\Omega_D$时，不同情况下摆在稳定摆动时的频率有明显区别，振幅大致相同。

总结

本文主要研究了简谐运动——单摆的相关问题，我们讨论了欧拉法和欧拉-克拉默两种方法，后者在分析时更具有优势。同时研究了单摆的阻尼振动，并在此基础上加上驱动项，得到单摆在一段时间后会趋于稳定。

致谢

本文参考了舟舟同学的作业和计算物理教材，特此感谢。