倒立单摆之分离频率
By 苏剑林 | 2014-03-18 | 33089位读者 |Mathieu方程
在文章《有质动力:倒立单摆的稳定性》中,我们分析了通过高频低幅振荡来使得倒立单摆稳定的可能性,并且得出了运动方程
$$l\ddot{\theta}+[h_0 \omega^2 \cos(\omega t)-g]\sin\theta=0$$
由此对单摆频率的下限提出了要求$\omega \gg \sqrt{\frac{g}{h_0}}$。然而,那个下限只不过是必要的,却不是充分的。如果要完整地分析该单摆的运动方程,最理想的方法当然是写出上述常微分方程的解析解。不过很遗憾,我们并没有办法做到这一点。我们只能够采取各种近似方法来求解。近似方法一般指数值计算方法,然后笔者偏爱的是解析方法,也就是说,即使是近似解,也希望能够求出近似的解析解。
首先不妨假设$\theta \ll 1$,那么可以使用近似$\sin \theta \approx \theta$,即
$$l\ddot{\theta}+[h_0 \omega^2 \cos(\omega t)-g]\theta=0$$
上述方程是一种特殊的Mathieu方程,与一般的Mathieu方程不同的是,参考书上的Mathieu方程的$g$是负数,而且$h_0 \omega^2$是小量。相反,这里的$h_0 \omega^2$是主要项,而$g$是正数。我们将采用一种特殊的技巧求得近似解析解。
数值分析
尽管我偏爱解析方法,但是数值分析发能够让我们先一瞧解的大概性质,以下是用Matlab求得的在$l=1,g=10,h=0.01,\omega=500$的情况下、初始条件为$\theta_0=1,\dot{\theta}_0=0$的解
在$l=1,g=10,h=0.01,\omega=500$的情况下的解,图像类似于简谐运动
图像表明周期解还是比较理想的,有点像简谐运动的解了,事实上也正是如此。但我们来仔细观察局部情况
上图中的“简谐运动”实际上叠加了一个高频振荡。
不难发现实际上的运动是在简谐运动基础上叠加了一个高频振荡。这也就不难理解为什么解析方法求解此题会遇到很多困难,主要是高频函数的存在。比如$\sin\omega t$是一个非常理想的函数,但是如果展开为泰勒级数就成了$\omega t-\frac{\omega^3 t^3}{6}+\dots$,如果只取前几项,收敛区域是非常小的。
现在关键的问题是,既然已经初步观察到是由双频率叠加的,我们最好能够把这两个频率分离开来,分别研究。这就是本文的核心。
分离频率
为了分离频率,我构思了以下方程组
$$\left\{ \begin{array}{l}
{\ddot{\phi}_1+\Omega^2 \phi_1=\varepsilon\left[\Omega^2 \phi_1+g\left(\phi_1+\phi_2\right)-\left(h_0 \omega^2 \cos\omega t\right)\phi_2\right]}\\
{\ddot{\phi}_2+\omega^2 \phi_2=\varepsilon\left[\omega^2 \phi_2-\left(h_0 \omega^2 \cos\omega t\right)\phi_1\right]}
\end{array} \right.$$
将两道方程叠加起来,取$\varepsilon=1$就得到了原来的Mathieu方程,其中$\theta=\phi_1+\phi_2$,而$\Omega$是待定的。初始条件是$\phi_1(0)=\theta_0,\dot{\phi}_1(0)=0;\phi_2(0)=h_0 \theta_0,\dot{\phi}_2(0)=0$。这样就相当于求解原Mathieu方程在$\theta(0)=(1+h_0)\theta_0,\dot{\theta}(0)=0$下的解。
把$\varepsilon$看作小参数,先求$\phi_2$,略去含有$\varepsilon$的项,即
$$\ddot{\phi}_2+\omega^2 \phi_2=0$$
结合初始条件得到
$$\phi_2=\theta_0 h_0 \cos\omega t$$
同理,对$\phi_1$进行同样的处理,我们求得
$$\phi_1=\theta_0 \cos\Omega t$$
要检验此解的合理性,主要是看略去的项是否合理。对于$\phi_2$,我们略去了$\varepsilon\left[\omega^2 \phi_2-\left(h_0 \omega^2 \cos\omega t\right)\phi_1\right]$,将所求得的$\phi_1$和$\phi_2$代进去,有
$$\begin{aligned}
\varepsilon\left[\omega^2 \theta_0 h_0 \cos\omega t-\left(h_0 \omega^2 \cos\omega t\right)\theta_0 \cos\Omega t\right]\\
=\varepsilon\omega^2 \theta_0 h_0 \cos\omega t(1-\cos\Omega t)
\end{aligned}$$
其中$\cos\Omega t$是一个低频项,在开始的长时间内(相对于高频项的周期)它都近似为1(这可以简单总结为在低频项中高频取0,在高频项中低频取常值),因此这一整项可以近似看作零。这表明,我们的略去是“自洽”的
对于$\phi_1$,我们略去了$\varepsilon\left[\Omega^2 \phi_1+g\left(\phi_1+\phi_2\right)-\left(h_0 \omega^2 \cos\omega t\right)\phi_2\right]$,代入即得
$$\begin{aligned}
\varepsilon \left[\Omega^2 \theta_0 \cos\Omega t+g\left(\theta_0 \cos\Omega t+\theta_0 h_0 \cos\omega t\right)-\left(h_0 \omega^2 \cos\omega t\right)\theta_0 h_0 \cos\omega t\right]\\
=\varepsilon\left[(\Omega^2+g)\theta_0 \cos\Omega t-\frac{1}{2}\theta_0 h_0^2 \omega^2+g\theta_0 h_0 \cos\omega t -\frac{1}{2}\theta_0 h_0^2 \omega^2 \cos 2\omega t\right]
\end{aligned}$$
继续我们的“高频取0,低频取1”,就得到
$$\varepsilon\left[(\Omega^2+g)\theta_0 -\frac{1}{2}\theta_0 h_0^2 \omega^2\right]$$
为了让项的略去更合理,这一项应该要近似为0,因此
$$\Omega^2=\frac{1}{2} h_0^2 \omega^2-g\geq 0$$
我们在上面的计算中默认了$l=1$,一般地有
$$\Omega^2=\frac{1}{2} h_0^2 \omega^2-gl\geq 0$$
这和朗道的《力学》中求得的结果是一样的,也和范翔的结果一致,可谓殊途同归,异曲同工!
所以最终我们求得近似解为
$$\theta\approx \theta_0 \cos\Omega t+h_0 \theta_0 \cos\omega t$$
为了求更高精度的近似解,我们可以把上面的方程组继续对$\varepsilon$展开。由于过程过于繁琐,在此从略。
转载到请包括本文地址:https://kexue.fm/archives/2471
更详细的转载事宜请参考:《科学空间FAQ》
如果您还有什么疑惑或建议,欢迎在下方评论区继续讨论。
如果您觉得本文还不错,欢迎分享/打赏本文。打赏并非要从中获得收益,而是希望知道科学空间获得了多少读者的真心关注。当然,如果你无视它,也不会影响你的阅读。再次表示欢迎和感谢!
如果您需要引用本文,请参考:
苏剑林. (Mar. 18, 2014). 《倒立单摆之分离频率 》[Blog post]. Retrieved from https://kexue.fm/archives/2471
@online{kexuefm-2471,
title={倒立单摆之分离频率},
author={苏剑林},
year={2014},
month={Mar},
url={\url{https://kexue.fm/archives/2471}},
}
March 18th, 2014
新皮肤看起来不错啊呵呵
最近用了下Mathematica,下次我也来仿真下~~
我就只会用mathematica解解方程,画画图,其他也就不会了...^_^
今天我想了一下这个问题,可以以那个移动的原点为参考系,单摆除了受重力外还受到一个非惯性力。当这个非惯性力是正弦形式时方程太复杂了我的数学水平根本搞不定,于是我就把这个非惯性力当成了方波形式,这样分析了一下。呵呵过几天我也写篇。
确实会是正弦形式,不过估计用方波也可以进行一定的近似。正弦的话就类似于我文章中的方程了,只能求助于各种近似方法。本文的主要介绍重点是“摄动中的分离频率”,我觉得这个值得进一步研究