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

这问题是在七年前提出来的，在知乎上可以说是一个古老的问题了。

然而遗憾的是，这么多年过去了，此问题依然没有得到很好地解决。

那么现在我来解决这个问题。

如上图所示，在地球的旋转参考系下，建立质点 $m$ 的运动方程。

质点受地球的引力，在地球内部时，引力表现为一谐振力

$\vec{F}=-\frac{GMm}{R^3}\vec r=-k\vec r$

加上科里奥利力 $-2m\vec\omega \times \dot{\vec r}$ 和惯性离心力 $-m\vec\omega \times ( \vec\omega \times \vec r )$ ，可写出控制方程

$m\ddot{ \vec r}=-k \vec r-2m\vec\omega \times \dot{\vec r}-m\vec\omega \times ( \vec\omega \times \vec r )$

其分量形式为

$m\ddot{x}=-kx+2m\omega\dot{ y}+m\omega^2x$

$m\ddot{y}=-ky-2m\omega \dot{x}+m\omega^2y$

$m\ddot{z}=-kz$

给定初值条件,例如在北纬30°处释放，

$x(0)=R\cos(\pi/6),y(0)=0,z(0)=R\sin(\pi/6)$

$x'(0)=y'(0)=z'(0)=0$

通过数值求解，我们发现，经过约40分钟（2550秒）后，质点从地球另一端出来了。

从 $z$ 轴上方俯视，轨迹是这样子，

可以看出，轨迹近似为一条直线，不经过地球中心，但这仅仅是参数使然，曲线的

解析形式已有答主给出。解析解的好处之一是可以用于分析下文所述的轨道“”花瓣“”

瓣数，但唯一不足的是存在解析解的复杂问题很少，因此适用性不够广泛。

继续分析轨道的特征，我们知道， 中国上海市的经纬度约为北纬30°，东经120°，

如果挖掘一条隧道穿过地球，那么这条轨道的出口将会在智利的

拉塞雷纳市（南纬30°，西经71°），距离智利首都圣地亚哥直线距离约为400公里。

而中国上海正对着的是阿根廷圣达菲省。现在看来，两者差的还是挺远的。

总结下来，如果科技进步到那一天，我们想去阿根廷首都布宜诺斯艾利斯看

足球赛的话，那么最佳的交通选择是先花40分钟穿过地球，然后转机飞往

布宜诺斯艾利斯。哈哈，值得期待吧。

然而事情远没有结束，接下来才是真正的波诡云谲。在讨论完上述问题之后，

很多人的第一想法就是延长计算时间，看看轨迹将会如何发展。

这个问题下的大部分答案都是如此。为了对比，我们也忽略惯性离心力，

沿用高票答案的初值。计算结果（计算时间10^5s）如下，

俯视图也给上

计算到这里都没什么问题，但是把时间延长到10^6s后，俯视图就变成了这般模样，

从参差不齐的轨迹边缘可以看出，误差的积累已经十分明显了，这在一个守恒

系统中显然是不准确的？为什么会这样呢？

其实在考虑哈密顿系统的长时间演化时，数值计算的格式非常重要。我国数学家

冯康院士提出，一切解哈密顿方程的“正确”离散算法都应当是辛变换的

（冯康，1997年国家自然科学一等奖，《哈密顿系统的辛几何算法》）。

在求解常微分方程时，普通的欧拉格式以及龙格-库塔格式等算法都是不保辛的。

因此在计算长时间演化的问题时，上述格式无一例外都失效了。

在常见数学软件Matlab,Mathematica中，默认算法都是龙格-库塔格式等

不保辛的传统算法，因此使用软件求解的结果基本都是失效的。那么应该怎样做呢？

这里根据上面的问题，我简单介绍一下。首先写出系统的哈密顿函数，

$H(p_x,p_y,x,y)=\frac{p_x^2+p_y^2+p_z^2}{2m}+\frac{1}{2}k(x^2+y^2+z^2)+\omega p_x y-\omega x p_y$

记 $z=\left\{p_x,p_y,x,y \right\}$ ，据此我们将哈密顿正则方程写成如下形式，

$z'=J^{-1}\nabla H(z)$

其中 $J$ 表示辛矩阵(Symplectic Matrix)。针对上述方程，最简单的一个二阶

辛差分格式是中点格式，

$\frac{z_{n+1}-z_{n}}{\Delta t}=J^{-1}\nabla H(\frac{z_{n+1}+z_{n}}{2})$

容易证明这个格式是保辛的。这是一个隐式格式，但方程是线性的，

因此求解并无困难。我们来看看计算结果，

结果很nice,未出现传统算法的误差积累。其实从误差分析中可以看出，辛差分格式

（Symplectic Difference Scheme）的误差始终是有界的，并不会无限增长，

这里不多赘述。具体分析，可以阅读文末的参考文献。值得一提的是，

mathematica中已经内置了哈密顿系统的辛几何算法，主要是

辛分区龙格-库塔算法（SymplecticPartitionedRungeKutta），

用来求解可分离的哈密顿系统

$H(p,q)=T(p)+V(q)$ 。常见的哈密顿系统都是可分离的，

因此使用mathematica能解决绝大部分类似的问题。最后哔哔几句吧，

《哈密顿系统的辛几何算法》是冯康院士世界级的科研成果，更为重要的是，

它是基础而实用的科研成果，了解并应用相关的知识并不需要多么

深厚的数学基础，希望能在中文互联网上见到更多的相关内容。

【参看文献】

【1】Feng K, Qin M. Symplectic geometric algorithms for Hamiltonian systems[M]. Berlin: Springer, 2010.

【2】Hairer E, Lubich C, Wanner G. Geometric numerical integration: structure-preserving algorithms for ordinary differential equations[M]. Springer Science & Business Media, 2006.

【3】“SymplecticPartitionedRungeKutta” Method for NDSolve

【4】冯康, 秦孟兆. 哈密顿系统的辛几何算法[M]. 杭州：浙江科学技术出版社，2003.

【哈密顿正则方程推导Mathematica代码】

Clear["Global`*"];
\[CapitalOmega] = {0, 0, \[Omega]};
r = {x[t], y[t], z[t]};
v = {x'[t], y'[t], z'[t]};
Vcor = -m Cross[\[CapitalOmega], r].v;
Vcen = -m Cross[\[CapitalOmega], r].Cross[\[CapitalOmega], r]/2;
Vela = k (r.r)/2;
T = m (v.v)/2;
L = T - Vcen - Vcor - Vela;
rulex = Solve[px[t] == D[L, x'[t]], x'[t]] // First;
ruley = Solve[py[t] == D[L, y'[t]], y'[t]] // First;
rulez = Solve[pz[t] == D[L, z'[t]], z'[t]] // First;
rule = {rulex // First, ruley // First, rulez // First};
H = ((-L + px[t] x'[t] + py[t] y'[t] + pz[t] z'[t]) /. rule) // 
  Simplify
eqs = {
   px'[t] == (-D[H, x[t]] // Simplify),
   py'[t] == (-D[H, y[t]] // Simplify),
   pz'[t] == (-D[H, z[t]] // Simplify),
   x'[t] == D[H, px[t]],
   y'[t] == D[H, py[t]],
   z'[t] == D[H, pz[t]]};
eqs // Column

来源：知乎 www.zhihu.com

作者：姜小白71

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

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