单摆模拟计算物理python 单摆的数学模型

求解了非小角单摆（数学摆）的周期，破除中学物理5度神话。是不是可以看成单摆，不是看是不是摆角小于5度，而是看计时仪器的精度。

参考：

• Moshe Gitterman，THE CHAOTIC PENDULUM，Chaper 1
• 非线性振动初步

体系

如图所示，体系由一个轻质细杆和末端带有的质点组成，细杆长\(l\)，质点质量为\(m\)。

图为数学摆

设摆角为\(\phi\)，忽略所有阻力和摩擦力，有

$\begin{equation*} ml^2\frac{d^2\phi}{dt^2}=-mgl\sin\phi \end{equation*}$

整理得

$\begin{equation*} \frac{d^2\phi}{dt^2}+\frac{g}{l}\sin\phi=0 \end{equation*}$

令\(\omega_0=\sqrt{g/l}\)，上式即为

$\begin{equation} \frac{d^2\phi}{dt^2}+\omega_0^2\sin\phi=0 \label{newtondyn} \end{equation}$

上式两边乘以\(\frac{d\phi}{dt}\)，

$\begin{equation*} \begin{split} &\frac{d\phi}{dt}\cdot\frac{d^2\phi}{dt^2}+\omega_0^2\sin\phi\frac{d\phi}{dt}=\frac{d\phi}{dt}\cdot\frac{d}{dt}\left (\frac{d\phi}{dt}\right )-\omega_0^2\frac{d\cos\phi}{dt}\\ &\frac{1}{2}\frac{d}{dt}\left (\frac{d\phi}{dt}\right )^2-\omega_0^2\frac{d\cos\phi}{dt}=0 \end{split} \end{equation*}$

上式对\(t\)积分，得

$\begin{equation} E=\frac{1}{2}\left (\frac{d\phi}{dt}\right )^2+\omega_0^2(1-\cos\phi) \label{energy} \end{equation}$

积分常数\(E\)的物理意义是摆的能量。\(E\)选得摆在最低点时势能为0。

在最大角位移\(\phi_0\)处，角速度\(\frac{d\phi}{dt}=0\)，代入\eqref{energy}式，得

$\begin{equation} E= \omega_0^2(1-\cos\phi_0) \label{energyexpl} \end{equation}$

将上式代入\eqref{energy}式，得

$\begin{equation} \frac{d\phi}{dt} =\omega_0\sqrt{2(\cos\phi-\cos\phi_0)} \label{dphidt} \end{equation}$

设\(t=0\)时，\(\phi=0\)，经过1/4周期，\(t=T/4\)时，摆到达最大位移\(\phi=\phi_0\)。对\eqref{dphidt}式积分，有

$\begin{equation*} \int_0^{T/4}\omega_0 dt=\int_0^{\phi_0}\frac{d\phi}{\sqrt{2(\cos\phi-\cos\phi_0)}} \end{equation*}$

即

$\begin{equation} \begin{split} \frac{\omega_0 T}{4}=&\int_0^{\phi_0}\frac{d\phi}{\sqrt{2(\cos\phi-\cos\phi_0)}}\\ =&\int_0^{\phi_0}\frac{d\phi}{2\sqrt{\sin^2(\phi_0/2)-\sin^2(\phi/2)}} \end{split} \label{intphi} \end{equation}$

上式用到了\(\cos\phi=1-2\sin^2(\phi/2)\)。

令

$\begin{equation} \sin(\phi/2)=\sin(\phi_0/2)\sin\psi \label{newvar} \end{equation}$

\(\phi\)的变化范围是从0到\(\phi_0\)，那么\(\psi\)的变化范围是从0到\(\pi/2\)。再考虑到\(\cos(\phi/2)/2d\phi=\sin(\phi_0/2)\cos\psi d\psi\)，代入\eqref{intphi}有

$\begin{equation} \begin{split} \frac{\omega_0 T}{4}=&\int_0^{\phi_0}\frac{d\phi}{2\sqrt{\sin^2(\phi_0/2)-\sin^2(\phi/2)}}\\ =&\int_0^{\phi}\frac{d\phi}{2\sin(\phi_0/2)\cos\psi }\\ =&\int_0^{\pi/2}\frac{d\psi}{\cos(\phi/2)}\\ =&\int_0^{\pi/2}\frac{d\psi}{\sqrt{1-\sin^2(\phi_0/2)\sin^2\psi}} \end{split} \label{intpsi} \end{equation}$

将上式中的被积函数展开成级数，

$\begin{equation*} \begin{split} \left [1-\sin^2(\phi_0/2)\sin^2\psi \right ]^{-0.5}=&1+\frac{1}{2}\sin^2(\phi_0/2)\sin^2\psi + \\ &\frac{3}{8}\sin^4(\phi_0/2)\sin^4\psi +\cdots \end{split} \end{equation*}$

代入\eqref{intpsi}，右边

$\begin{equation} \begin{split} & \int_0^{\pi/2}\frac{d\psi}{\sqrt{1-\sin^2(\phi_0/2)\sin^2\psi}}=\\ & \int_0^{\pi/2}\left [1+\frac{1}{2}\sin^2(\phi_0/2)\sin^2\psi +\frac{3}{8}\sin^4(\phi_0/2)\sin^4\psi \right ] d\psi = \\ & \frac{\pi}{2}+\frac{\pi}{8}\sin^2(\phi_0/2)+\frac{9}{128}\sin^4(\phi_0/2)+\cdots \end{split} \label{right} \end{equation}$

代入\eqref{intpsi}，考虑到小角摆动周期为\(T_0=2\pi/\omega_0\)，得数学摆的周期为

$\begin{equation} \begin{split} T=&\frac{2\pi}{\omega_0}\left [1+\frac{1}{4}\sin^2(\phi_0/2)+\frac{9}{64}\sin^4(\phi_0/2)+\cdots \right ]\\ \approx & T_0\left [1+\frac{1}{4}\sin^2(\phi_0/2) \right ] \end{split} \label{period} \end{equation}$

可见，周期与摆幅\(\phi_0\)有关。只有当摆幅\(\phi_0\)非常小的时候，周期才与摆幅无关。下表是周期与摆幅的对应的关系：

\(\phi_0\)	\(0^{\circ}\)	\(3^{\circ}\)	\(5^{\circ}\)	\(10^{\circ}\)	\(30^{\circ}\)	\(45^{\circ}\)
\(T/T_0\)	1	0.0002	1.0005	1.0019	1.0174	1.0369

即便摆幅很小，\(\phi_0\lt 5^{\circ}\)，周期也是随摆幅而变的。即便摆幅比较大，如\(\phi_0 = 45^{\circ}\)，周期也只是比0摆幅极限下的周期\(T_0\)多了3.7%而已。

中学物理讲单摆有个\(5^{\circ}\)神话，即要求单摆摆幅要小于\(5^{\circ}\)。应该是因为中学物理实验室的计时仪器最高能精确到毫秒的原因吧。如果计时仪器只能精确到秒，摆幅\(45^{\circ}\)也是可以的。