在地球上挖一条何种走向的隧道能使小球仅受引力作用不接触隧道壁从一端运动到另一端？

Emmm

其实我没看懂大家在闹腾啥，这么简单的一个问题偏偏要做数值解…

我们直接把其他答案里列出的方程抄过来

$\begin{aligned} \ddot x&=-kx-2\omega\dot y+\omega^2 x\\ \ddot y&=-ky+2\omega\dot x+\omega^2 y\\ \ddot z&=-kz \end{aligned}$

关于 z 的式子是可以直接肉眼看出解的 $z = z_0\cos\sqrt{k}t$ ；而对于前两个耦合的方程，标准套路是引入向量 $\bm v = (x,y,\dot x,\dot y)^T$ ，将它们变成一个四阶的矩阵方程

$\dot{\bm v} = A\bm v$

再寻求线性变换P将矩阵A对角化 $P^{-1}AP=\text{diag}\{\lambda_1,\lambda_2,\lambda_3,\lambda_4\}$ ，这样就可以拆成四个独立的简单方程直接写出解。剩下的事情就是拼回来根据初始条件待定系数了。

当然如果你稍微有点经验，比如说上力学课的时候学过怎么解傅科摆，那么你应该知道这里有种骚操作。我们引入复变量 $u=x+\mathrm iy$ ，把前面说的两个式子一个乘上 i 加上另外一个，得到

$\ddot u=-ku+2\mathrm i\omega u+\omega^2 u$

这种其实就类似于一个阻尼谐振子的方程，它的通解可以写作 $u=u_1\mathrm e^{\mathrm i\Omega_1 t}+u_2\mathrm e^{\mathrm i\Omega_2 t}$ ，而 $\Omega$ 是下面这个方程的两个根

$-\Omega^2=\omega^2-k-2\omega\Omega$

解出来的结果是 $\Omega_{1,2}=\omega\pm \sqrt{k}$ ，然后x,y就可以通过 $x={\rm Re}u,y={\rm Im}u$ 得到。

我们现在来看根据初始条件确定系数。假设人从纬度角 $\theta$ 处静止出发，那么初始条件可以取作

$x(0)=R\cos\theta, y(0)=\dot x(0)=\dot y(0)=0$

这样确定了系数之后，运动方程就可以写作

$u=\mathrm e^{\mathrm i\omega t}(\cos kt-\frac{\omega}{k}\sin kt)R\cos\theta$

或者写成 x, y 下的形式

$x=\cos{\omega t}(\cos \sqrt{k}t-\frac{\omega}{\sqrt{k}}\sin \sqrt{k}t)R\cos\theta\\ y=\sin{\omega t}(\cos \sqrt{k}t-\frac{\omega}{\sqrt{k}}\sin \sqrt{k}t)R\cos\theta$

从这个式子去看，事情就很明显了：物体在惯性系中作一个闭合的椭圆轨道运动，再附加上地球的自转。说到底，其实从一开始就可以换到惯性系中去作这个计算，只是不少人非得要在非惯性系里去折腾这些麻烦的东西，我也没办法。

图我就不画了，大家看看别的答案就好。值得指出的一点是，这类运动轨迹被称为旋轮线，可以理解为一个圆绕着另一个固定的圆滚动时，动圆上一点形成的轨迹。大家如果和我年龄相当的话，小时候应该玩过类似的玩具

下面回答 @王赟 Maigo 在评论区提出的问题：为什么轨迹是17角星？

事实上，这和 @山醒最先给出，然后被其他人借用的参数有关。或者说，就是个巧合：

前面我们已经给出了，在惯性系里椭圆运动的角频率为 $\sqrt{k}=\sqrt{\frac{GM}{R^3}}$ ，用这个角频率除以地球自转的角频率，就是两种运动周期之比，结果是

$\frac{1}{7.2921\times 10^{-5}}\sqrt{\frac{6.674\times10^{-11} 5.9722 \times10^{24}}{(6.3781\times 10^6)^3}} = 16.9969$

emmm

在千分之一的精度上等于17

所以对于地心轨道列车一个更简明的结论是，每天往返17趟（诶？

当然，同许多地方已经讨论过的一样，这个计算未曾考虑过地球内部密度的变化，因此具体数值肯定是不能采信的。

我们再来看看 @姜小白71 同志在数值上给出的讨论，这个和问题就稍微偏的有些远，不关心的读者可以点叉了。

首先可以确定，我们不能单凭下面这张图就断定数值结果有误差积累

因为很明显，这对应着两种周期运动的周期比并不完全是个整数，而是和17差了一点点，所产生的进动。要确认这是不是反映了误差，我们要去计算解析结果预测的进动角

$\delta=(\omega-\sqrt{k}/17)t=0.75^\circ$

远小于图中反映的数度的进动角。从这里，我们才能确认，他使用的第一种数值方法，相对后面使用的辛算法（欧拉中点格式），存在明显的累积误差。

下一个问题是，在面对这类问题，即准确计算长时演化的进动角，换言之准确估计传播子的本征值时，辛算法相对其他算法存在优势吗？

答案是否定的。事实上，辛算法的优势是在于长时演化时保守恒量，例如可以保证系统原有的二次守恒量在数值传播时守恒，保证存在一个守恒的准哈密顿量，保相空间体积元等等诸如此类。而对于估计本征值并没有特别的优势。例如，对于 @姜小白71 所使用的欧拉中点格式，其相当于用一个一阶的Cayley展开去近似传播子

$\mathrm e^{\tau A}= \frac{1+\tau A/2}{1-\tau A/2} + O(\tau^3)$

具有二阶精度。而对于各种数学软件常用的默认格式，四阶龙格库塔，则具有四阶精度。如果拿这两个算法，在共同的收敛域内取同样的时间步长，去比较谁算进动角算得更准，那显然是四阶龙格库塔完成

（未完待续）

这一块是毫无意义的脑洞：

根据 Christopher Wren 于1660 年向英国皇家学会提出的建议，一米长的单摆振动周期可以定为两秒，即

$2\pi\sqrt{\frac{1\,\rm m}{g}} = 2\,\rm s$

根据 1791 年法国科学院的定义，子午线的长度是一千万米，即

$\pi R/2 = 10^7\,\rm m$

根据我们的计算，地心列车的运行频次是每天17次，即

$2\pi\sqrt{\frac{R^3}{GM}} = \frac{1}{17} \,\text{day}$

综合以上三式，每天的时长应该是

$1\,\text{day} = 34\sqrt{2\times 10^7/\pi}\approx 85787\,\rm s$

不足传统认知上的86400 s.

来源：知乎 www.zhihu.com

作者：梁昊

【知乎日报】千万用户的选择，做朋友圈里的新鲜事分享大牛。
点击下载

此问题还有 46 个回答，查看全部。
延伸阅读：
在地球上挖一条何种走向的隧道能使小球仅受引力作用不接触隧道壁从一端运动到另一端？
以现在人类的能力，可以挖掘出一条贯穿地球的隧道吗？另外，地球中心的引力是什么样的？