重力列車の最速降下問題について、あれこれ考えてみました。

表題について

「重力列車」とは、重力のみを推進力として移動する列車です。地表の２点（始点、終点）を結ぶ地下トンネルを掘り、レールを敷いて列車を転がします。地球の密度分布が地球中心から見て等方的で、空気抵抗やレールとの摩擦などの非保存力を無視できる場合、重力列車は地下トンネルの形状に依らず、速度ゼロで地表の終点に到着します。

en.wikipedia.org

「最速降下問題」とは、任意の２点を結ぶ曲線のうち、質点が最短時間で転がりきるものを求める問題です。一様重力の場合、答えはサイクロイド（擺線）となることが知られています。

en.wikipedia.org

さて、重力列車の最速降下問題を考えると、重力の大きさと向きが場所場所で異なるため、答えはもはやサイクロイドとはなりません。ただし、始点と終点がごく近傍にある場合に限っては、重力がほぼ一様と見做せるため、答えはサイクロイドに近い曲線になろうと思われます。すると、この重力列車の最速降下曲線について次の２点が気になります：

サイクロイドのように簡単な式で表すことができるか
終点に到達するまでの所要時間について、サイクロイドとの差はいかほどか

この記事の目的は、上記 1., 2. について考えることです。

一様重力における最速降下問題

まずは、一様重力における最速降下問題を考えます。この段階を踏む目的は２つあります：

単位と座標を設定する
極限 $1/R \to 0$ （ $R$ は地球半径）の特別な場合として、あとで整合性のチェックに使う

問題設定

図のように、水平方向に $x$ 軸、鉛直上向きに $y$ 軸をとります。最速降下曲線は、パラメタ $\theta$ （始点は $\theta=-\pi$ 、終点は $\theta=+\pi$ ）を使って媒介変数表示の形で求めます。座標原点は、次の式が満たされるようにとります。

$\displaystyle x(+\pi) = -x(-\pi) \gt 0 \\ y(+\pi) = y(-\pi) = 0$

力学的エネルギー保存則と、運動エネルギーが非負であることから、列車は始点の高さ $y(-\pi)$ より上には登れません。よって

$\displaystyle y(\theta) \leq 0$

です。

最速降下曲線の方程式

列車の質量を $m$ 、重力加速度を $1$ とすると、力学的エネルギー保存則から

$\displaystyle \frac{m}{2}v^2 + my = 0$

これを速さ $v$ について解くと $v = \sqrt{-2y}$ です。所要時間 $T$ は次の通り書けます。

$\displaystyle T=\int_0^T{\textrm{d}t} = \int_{\theta=-\pi}^{\theta=+\pi} \frac{\textrm{d}l}{v} = \int_{\theta=-\pi}^{\theta=+\pi} f(y,y';x)\textrm{d}x$

$\textrm{d}l=\sqrt{\textrm{d}x^2+\textrm{d}y^2}$ は曲線の微小長さです。 $x$ についての微分をプライム ( $'$ ) で表します。被積分関数

$\displaystyle f(y,y';x) = \sqrt{\frac{1+y'^2}{-2y}}$

は $x$ に陽に依存しないので、 $T$ が勝手な微小変形 $\delta y$ に依らない条件はベルトラミの恒等式

$\displaystyle f-y'\frac{\partial f}{\partial y'} = \frac{1}{\sqrt{-2y(1+y'^2)}} = \textrm{const.}$

すなわち

$\displaystyle y(1+y'^2)=\textrm{const.}$

で与えられます。

en.wikipedia.org

この微分方程式は変数分離法で解くことができますが、ここでは少し趣向を変え、アタリをつけて考えてみます。

最速降下曲線の概形について考えてみる（少し脱線）

列車にかかる重力をレールに平行／垂直な成分に分解すると、列車の加速に関わるのは平行成分のみです。特に、レールが鉛直のとき加速度最大（自由落下）、水平のとき加速度ゼロ（等速直線運動）になります。最短時間を目指すとき、次の２つがトレードオフの関係になっています：

レールが急なほど、大きい加速が得られる
レールが緩やかなほど、大きい 水平方向の 速度が得られる

これを踏まえると、低速な始点付近ではレールを急にして加速し、スピードがのったらレールを緩やかにして終点直下に向かって進み、最後にレールを急にして終点まで駆け上ると良さそうです（例えば、せっかちな人の車の運転：急アクセル・急ブレーキを想起してください）。すると、最速降下曲線の傾きの候補として

$\displaystyle y'=\tan{\frac{\theta}{2}}$

が思いつきます。

最速降下曲線

実際に試してみます。 $y'=\tan(\theta/2)$ より

$\displaystyle 1+y'^2=\frac{1}{\cos^2\frac{\theta}{2}}=\frac{2}{1+\cos\theta}$

これをベルトラミの恒等式に代入して

$\displaystyle y \propto 1+\cos\theta$

を得ます。長さの単位を決めてませんでしたので、 $y$ の最小値を $-2$ として

$\displaystyle y = -(1+\cos\theta)$

とします。すると

$\displaystyle \frac{\textrm{d}x}{\textrm{d}\theta}=\frac{1}{y'}\frac{\textrm{d}y}{\textrm{d}\theta}=\frac{1}{\tan\frac{\theta}{2}}\sin\theta=\frac{\cos\frac{\theta}{2}}{\sin\frac{\theta}{2}}2\sin\frac{\theta}{2}\cos\frac{\theta}{2}=2\cos^2\frac{\theta}{2}=1+\cos\theta$

となりますから、 $\theta$ について積分して

$\displaystyle x=\theta+\sin\theta$

ただし、問題設定の $x(\pi) = -x(-\pi)$ を満たすように積分定数を決めました。まとめると

$\displaystyle x=\theta+\sin\theta \\ y = -(1+\cos\theta)$

です。この曲線、サイクロイドは、確かにベルトラミの恒等式を満たしています。

en.wikipedia.org

傾斜の変化率

あとで整合性チェックに用いるため、傾斜の変化率 $y''$ も求めておきます。

$\displaystyle y'' = \frac{\textrm{d}}{\textrm{d}x}\left(y'\right)=\frac{1}{\textrm{d}x/\textrm{d}\theta}\frac{\textrm{d}}{\textrm{d}\theta}\left(\tan\frac{\theta}{2}\right) = \frac{1}{1+\cos\theta}\frac{\frac{1}{2}}{\cos^2\frac{\theta}{2}}= \frac{1}{(1+\cos\theta)^2}$

$y''>0$ より、最速降下曲線は下に凸です。

所要時間 $T$

以上より、所要時間 $T$ は

$\displaystyle T= \int_{\theta=-\pi}^{\theta=+\pi} \sqrt{\frac{1+y'^2}{-2y}}\textrm{d}x = 2\pi$

です。

せっかくなので、１つ具体的に計算してみます。人類が到達した最深地点が「コラ半島超深度掘削坑」の約 $12 \textrm{km}$ だそうなので、深さ $12.0 \textrm{km}$ のサイクロイド地下トンネルを考えます。

en.wikipedia.org

まず、サイクロイドの深さを $2$ としたので、長さの $1$ 単位は $6.0 \textrm{km}$ です。次に、重力加速度を $1$ とおいたことから時間の $1$ 単位は $\sqrt{6.0\textrm{km} / (9.8\textrm{m}/\textrm{s}^2)}\sim 25\textrm{s}$ となります。よって所要時間は $2\pi \sim 160\textrm{s}$ 、およそ $2$ 分 $40$ 秒です。始点と終点との距離は $2\pi \sim 38\textrm{km}$ ですが、東京～八王子間の直線距離がおよそ $40 \textrm{km}$ だそうなので、すごい速さですね。

重力列車の最速降下問題

問題設定

座標原点で $x$ 軸が地表に接するよう、地球中心の座標を $(0, -R)$ とします。 $y$ 軸は、次の式が満たされるようにとります。

$\displaystyle x(+\pi) = -x(-\pi) > 0 \\ y(+\pi) = y(-\pi) > -R$

極限 $1/R \to 0$ で、前節の問題設定と一致します。

最速降下曲線の方程式

地球の質量を $1$ とします。重力加速度を $1$ とおいたため、質量 $m$ の質点が地表で受ける重力は

$\displaystyle m = G\frac{m}{R^2}$

です。よって、万有引力定数は $G=R^2$ となります。

地球を密度が一様な球体と考えると、中心から距離 $r (\leq R)$ 以内の部分の質量は $(r/R)^3$ となります。列車が地球の中心から距離 $r$ にあるとき、距離 $r$ より外側の部分から受ける重力はすべて相殺されることが知られていますので、このとき列車が受ける重力 $F(r)$ は

$\displaystyle F(r) = G\frac{m(r/R)^3}{r^2} = m\frac{r}{R}$

です。

en.wikipedia.org

ゆえに、地球内部の $r=\rho$ における位置エネルギーは、地球中心を基準として

$\displaystyle -\int_\rho^0 F(r) \textrm{d}r= \frac{m}{2}\frac{\rho^2}{R}$

すると、力学的エネルギー保存則から

$\displaystyle \frac{m}{2}v^2 + \frac{m}{2}\frac{\rho^2}{R} = 0 + \frac{m}{2}R$

これを $v$ について解くと

$\displaystyle v = \sqrt{\frac{R^2-\rho^2}{R}} = \sqrt{\frac{R^2-x^2-(y+R)^2}{R}} = \sqrt{-2y-\frac{x^2+y^2}{R}}$

$1/R\to 0$ で $v=\sqrt{-2y}$ となり、整合性がとれています。所要時間 $T$ は次の通り書けます。

$\displaystyle T=\int_0^T{\textrm{d}t} = \int_{\theta=-\pi}^{\theta=+\pi} \frac{\textrm{d}l}{v} = \sqrt{R}\int_{\theta=-\pi}^{\theta=+\pi} h(y,y',x)\textrm{d}x$

ここで、被積分関数

$\displaystyle h(y,y',x) = \sqrt{\frac{1+y'^2}{R^2-x^2-(y+R)^2}}$

は残念ながら $x$ に陽に依存しているため、オイラーの方程式

$\displaystyle \frac{\partial h}{\partial y} - \frac{\textrm{d}}{\textrm{d}x}\left(\frac{\partial h}{\partial y'}\right)=0$

を解かなければなりません。

en.wikipedia.org

部分ごとに分けて計算しますと、まず左辺第１項は

$\displaystyle \frac{\partial h}{\partial y} = \frac{\sqrt{1+y'^2}}{\sqrt{R^2-x^2-(y+R)^2}^3}(y+R)$

となります。次に、第２項で

$\displaystyle \frac{\partial h}{\partial y'} = \frac{y'}{\sqrt{(1+y'^2)\left[R^2-x^2-(y+R)^2\right]}}(y+R)$

なので、

$\displaystyle \frac{\textrm{d}}{\textrm{d}x}\left(\frac{\partial h}{\partial y'}\right) = \frac{y''\left[R^2-x^2-(y+R)^2\right]+y'(1+y'^2)[x+(y+R)y']}{\sqrt{(1+y'^2)\left[R^2-x^2-(y+R)^2\right]}^3}$

となります（しんどくなってきました）。以上をオイラーの方程式に代入すると

$\displaystyle (1+y'^2)^2(y+R)-y''\left[R^2-x^2-(y+R)^2\right]-y'(1+y')[x+(y+R)y']=0$

を得ます。これを展開して整理すると、次の式が得られます。

$\displaystyle (1+y'^2)\left(y+R-xy'\right)+y''\left(x^2+y^2+2yR\right)=0$

この微分方程式の解について、またアタリをつけて考えてみます。

最速降下曲線の概形について考えてみる（少し脱線）

最速降下曲線の方程式を $R$ で割り

$\displaystyle (1+y'^2)\left(1+\frac{y-xy'}{R}\right)+2yy''\left(1+\frac{x^2+y^2}{2yR}\right)=0$

極限 $1/R \to 0$ をとると

$\displaystyle 1+y'^2+2yy''=\frac{1}{y'}\frac{\textrm{d}}{\textrm{d}x}\left(y(1+y'^2)\right)=0$

となり、一様重力の場合と整合しています。すると、式中の

$\displaystyle 1+\frac{y-xy'}{R},$

$\displaystyle 1+\frac{x^2+y^2}{2yR}$

の２つの因子が、一様重力の場合との違いそのものであることが分かります。
ここで、第２の因子は地表

$\displaystyle x^2+(y+R)^2=R^2$

においてゼロになりますが、このとき最速降下曲線の方程式が満たされるには第１項がゼロでなければならず、 $1+y'^2 > 0$ より第１の因子がゼロとなります。 $y$ について解けば

$\displaystyle y=y'x-R$

です。これは、地表における最速降下曲線の接線が地球の中心 $(0, -R)$ を通ることを示しています。つまり、一様重力の場合と同様に、始点と終点においてレールが重力に平行になります。このことを踏まえると、レールの傾きは

$\displaystyle y'=\tan\left[\left(\frac{1}{2}-\frac{1}{R}\right)\theta\right]$

となるのでは？と予想できます。
逆に、第１の因子は地球中心を通る直線 $y=y'x-R$ 上でゼロとなります（直線なので $y'=\textrm{const.}$ ）が、このとき $y''=0$ なので方程式が満たされています。当たり前のようですが、対蹠地（地球の反対側）を結ぶ最速降下曲線は、地球中心を通る直線です。整理すると、最速降下曲線は、

極限 $1/R \to 0$ でサイクロイドになる
始点と終点で、その接線が地球の中心を通る
対蹠地を結ぶとき、地球の中心を通る直線になる

ような曲線です。最も自然に思いつくのはハイポサイクロイド（内擺線）です。

en.wikipedia.org

最速降下曲線

実際に試してみます。まずハイポサイクロイドの式を求め、次に最速降下曲線の方程式に代入し、果たして１つの解になっているか確かめます。

極限 $1/R\to 0$ で前節のサイクロイドと一致するように、転がす円の半径を $1$ とします。パラメタ $\theta$ は、図の $\angle ABX$ です。すると、その対頂角を中心角とする円弧（黄線）の長さが $\theta$ になります。この弧長は、転がす円が地表に接する点と $O$ とを結ぶ円弧（緑線）の長さと等しいはずですから、この円弧に対する中心角は $\angle OAB = \theta/R$ となります。また、直線 $BX$ が $y$ 軸となす角は $\angle ABX - \angle OAB = (1-1/R)\theta$ となります。
さて、 $\vec{OX}$ を求めます。 $\vec{OX}=\vec{OA}+\vec{AB}+\vec{BX}$ と分解すると、それぞれのベクトルの向き（角度）は今みた通りですから、

$\displaystyle \vec{OA}=-R\hat{y}$

$\displaystyle \vec{AB}= (R-1)\left(\sin\frac{\theta}{R}\hat{x}+\cos\frac{\theta}{R}\hat{y}\right) \\$

$\displaystyle \vec{BX}=\sin\left[\left(1-\frac{1}{R}\right)\theta\right]\hat{x}-\cos\left[\left(1-\frac{1}{R}\right)\theta\right]\hat{y}$

です。ただし、 $\hat{x}, \hat{y}$ は、それぞれ $x, y$ 軸正方向の単位ベクトルです。結局、ハイポサイクロイドの媒介変数表示は

$\displaystyle x=(R-1)\sin\frac{\theta}{R}+\sin\left[\left(1-\frac{1}{R}\right)\theta\right]$

$\displaystyle y=(R-1)\cos\frac{\theta}{R}-\cos\left[\left(1-\frac{1}{R}\right)\theta\right]-R$

となります。すると、傾き $y'$ は

$\displaystyle y'=\frac{\textrm{d}y/\textrm{d}\theta}{\textrm{d}x/\textrm{d}\theta}=\frac{\left(1-\frac{1}{R}\right)\left(\sin\frac{\theta}{R}+\sin\left[\left(1-\frac{1}{R}\right)\theta\right]\right)}{\left(1-\frac{1}{R}\right)\left(\cos\frac{\theta}{R}+\cos\left[\left(1-\frac{1}{R}\right)\theta\right]\right)}=\frac{2\cos\frac{\theta}{2}\sin\left[\left(\frac{1}{2}-\frac{1}{R}\right)\theta\right]}{2\cos\frac{\theta}{2}\cos\left[\left(\frac{1}{2}-\frac{1}{R}\right)\theta\right]}=\tan\left[\left(\frac{1}{2}-\frac{1}{R}\right)\theta\right]$

となり、まさに予想通りの結果となりました。
それでは、最後に最速降下曲線の方程式にハイポサイクロイドの式を代入します。方程式を再掲します。

$\displaystyle (1+y'^2)\left(y+R-xy'\right)+y''\left(x^2+y^2+2yR\right)=0$

これを、４つのパートに分けて計算していきます。まず、左辺第１項について

$\displaystyle 1+y'^2=\frac{1}{\cos^2\left[\left(\frac{1}{2}-\frac{1}{R}\right)\theta\right]}$

また

$\displaystyle y+R-xy' \\ =(R-1)\cos\frac{\theta}{R}-\cos\left[\left(1-\frac{1}{R}\right)\theta\right]-\left((R-1)\sin\frac{\theta}{R}+\sin\left[\left(1-\frac{1}{R}\right)\theta\right]\right)\tan\left[\left(\frac{1}{2}-\frac{1}{R}\right)\theta\right] \\ =\frac{(R-1)\left(\cos\frac{\theta}{R}\cos\left[\left(\frac{1}{2}-\frac{1}{R}\right)\theta\right]-\sin\frac{\theta}{R}\sin\left[\left(\frac{1}{2}-\frac{1}{R}\right)\theta\right]\right)-\left(\cos\left[\left(1-\frac{1}{R}\right)\theta\right]\cos\left[\left(\frac{1}{2}-\frac{1}{R}\right)\theta\right]+\sin\left[\left(1-\frac{1}{R}\right)\theta\right]\sin\left[\left(\frac{1}{2}-\frac{1}{R}\right)\theta\right]\right)}{\cos\left[\left(\frac{1}{2}-\frac{1}{R}\right)\theta\right]} \\ =(R-2)\frac{\cos\frac{\theta}{2}}{\cos\left[\left(\frac{1}{2}-\frac{1}{R}\right)\theta\right]}$

次に、第２項について

$\displaystyle y''=\frac{\textrm{d}}{\textrm{d}x}y'=\frac{1}{\frac{\textrm{d}x}{\textrm{d}\theta}}\frac{\textrm{d}}{\textrm{d}\theta}\left(\tan\left[\left(\frac{1}{2}-\frac{1}{R}\right)\theta\right]\right) \\ =\frac{1}{\left(1-\frac{1}{R}\right)\left(\cos\frac{\theta}{R}+\cos\left[\left(1-\frac{1}{R}\right)\theta\right]\right)}\frac{\frac{1}{2}-\frac{1}{R}}{\cos^2\left[\left(\frac{1}{2}-\frac{1}{R}\right)\theta\right]} \\ =\frac{1-\frac{2}{R}}{1-\frac{1}{R}}\frac{1}{\left(\cos\frac{\theta}{R}+\cos\left[\left(1-\frac{1}{R}\right)\theta\right]\right)\left(1+\cos\left[\left(1-\frac{2}{R}\right)\theta\right]\right)}$

これは、確かに極限 $1/R\to 0$ で一様重力の場合と整合しています。もう少し変形を進めて

$\displaystyle y''=\frac{R-2}{R-1}\frac{1}{4\cos\frac{\theta}{2}\cos^3\left[\left(\frac{1}{2}-\frac{1}{R}\right)\theta\right]}$

次に、

$\displaystyle x^2+y^2+2yR \\ =x^2+(y+R)^2-R^2 \\ =\left((R-1)\sin\frac{\theta}{R}+\sin\left[\left(1-\frac{1}{R}\right)\theta\right]\right)^2+\left((R-1)\cos\frac{\theta}{R}+\cos\left[\left(1-\frac{1}{R}\right)\theta\right]\right)^2-R^2 \\ =(R-1)^2+1-2(R-1)\left(\cos\frac{\theta}{R}\cos\left[\left(1-\frac{1}{R}\right)\theta\right]-\sin\frac{\theta}{R}\sin\left[\left(1-\frac{1}{R}\right)\theta\right]\right)-R^2\\ =-2(R-1)(1+\cos\theta)\\ =-4(R-1)\cos^2\frac{\theta}{2}$

です。以上より、左辺は

$\displaystyle (1+y'^2)\left(y+R-xy'\right)+y''\left(x^2+y^2+2yR\right)\\ =\frac{1}{\cos^2\left[\left(\frac{1}{2}-\frac{1}{R}\right)\theta\right]}(R-2)\frac{\cos\frac{\theta}{2}}{\cos\left[\left(\frac{1}{2}-\frac{1}{R}\right)\theta\right]}+\frac{R-2}{R-1}\frac{1}{4\cos\frac{\theta}{2}\cos^3\left[\left(\frac{1}{2}-\frac{1}{R}\right)\theta\right]}\left[-4(R-1)\cos^2\frac{\theta}{2}\right]\\ =0$

となり、ハイポサイクロイドは最速降下曲線の方程式の１つの解であることが分かりました。

所要時間 $T$

以上より、所要時間 $T$ は

$\displaystyle T= \sqrt{R}\int_{\theta=-\pi}^{\theta=+\pi} \sqrt{\frac{1+y'^2}{R^2-x^2-(y+R)^2}}\textrm{d}x = 2\pi\sqrt{1-\frac{1}{R}}$

となり、 $1/R\to 0$ の極限において一様重力の場合と一致しています。 $R\sim 6400\textrm{km}$ に対して深さ $12 \textrm{km}$ 程度では $0.073$ 秒しか変わりません。

一方、 $R=2$ （最小の $R$ ）は対蹠地を結ぶ直線の場合で、このとき $T=\sqrt{2}\pi$ となります。長さの $1$ 単位が $3200 \textrm{km}$ 、時間の $1$ 単位が $\sqrt{3200\textrm{km}/(9.8\textrm{m}/\textrm{s}^2)}\sim 570\textrm{s}$ となるため、所要時間は $\sqrt{2}\pi\sim 2500\textrm{s}$ 、およそ $42$ 分となります。