前記事では、幾何学的な特性と座標変換を組み合わせて、大圏コースをメルカトル図法座標上に変換することで最短経路を求めた。今度はメルカトル図法上で直接最短経路を導く方法を考えよう。
ある非ユークリッドな座標系に対する最短経路を満たす曲線を、媒介変数$t$を用いて以下で定義しよう。
$$X(t)=(x(t),y(t))$$
なんのことはなく、$x,y$は媒介変数$t$によって決まります、とだけ言っているのみである。
この関数を使って経路の長さを与える式をたてる。
$$I=\int_{t_1}^{t_2} L(X)dt$$
$$L(X)=\sqrt{\sum_i g_i(X)\frac{dX_i^2}{dt}}$$
ここで$g_i$は着目点$X$における距離の重みであり、目盛りの刻み(縮尺の逆数)に相当する。
結論から先にいえば、メルカトル図法に対する目盛りの刻みは座標によって異なり、以下の式で与えられる。
$$g_i(X)=\sin^2 y$$
ちなみに$X=(\theta,\phi)$とし、$g_i$が$i$に依存しないので縦横同じ刻みである。
リーマン多様体と計量、反変、共変、縮約
メルカトル図法をはじめ、各図法はリーマン多様体上の座標であると考えた方が便利である。
ここで、メルカトル図法を一旦忘れてリーマン多様体とその付近の用語を整理しよう。
刻みとは、たとえば1mを1mmの刻みの定規で表現しようとすると1000目盛りである。一方、1目盛りあたりの長さは1mの1/1000であるので、
$$g_mm=\frac{1}{1000}[\mathrm m],g_m=1[\mathrm m]$$
となる。
このとき、刻み1mm定規で1000目盛りの距離は
$$g_mm\cdot 1000=1[\mathrm m]$$
一方、刻み1mの定規で1目盛りの距離も
$$g_m\cdot 1=1[\mathrm m]$$
となり、使った定規に依存しない値として求まる。
この定規のことを計量と呼ぶ。
ここで、ある定規で定義された距離を刻みの小さい定規で表したとき、値は刻みの小ささと反比例することになることに着目しよう。座標は多次元であればベクトルで表現できるので、このような性質のベクトルを反変ベクトル、テンソルならば反変テンソルと呼ぶ。
同様に、目盛りの刻み自身は細かく刻むほど小さくなる。このような性質をもつベクトルやテンソルを共変ベクトル、共変テンソル呼ぶ。
定規(座標系)は縦横で同じである必要はなく、もっといえば二軸が直交する必要もない。
ユークリッド幾何では、多次元空間上の距離は三平方の定理から定義できるのであった。
二次元の場合、
$$r_2^2=x^2+y^2$$
多次元の場合は斜辺$r_2$を使って
$$r_3^2=r_2^2+z^2$$
と拡張していけば、
$$r^2=\sum_i x_i^2$$
刻みが縦横で違う場合には、それぞれ$g_i$として
$$r^2=\sum_i g_ix_i^2$$
軸が直交していないときは、三平方の定理の代わりに余弦定理を使えばいい。
$$r_2^2=g_xx^2+g_yy^2-2\sqrt{g_xg_y}xy\cos\theta$$
ちょっと書き換えると、
$$r_2^2=g_{xx}xx+g_{yy}yy+g_{xy}xy+g_{yx}yx$$
このとき、
$$g_{xy}=g_{yx}=-\sqrt{g_{xx}g_{yy}}\cos\theta$$
である。
多次元に拡張すると結局
$$r^2=\sum_{i,j}g_{ij}X_iX_j$$
となる。$g_{ij}$のことをとくに計量テンソルと呼称する。
今度は平坦(刻みが一様)な空間ではなく、歪み (刻みが場所によって違う)をもった空間を考える。このとき$g_{ij}$は座標依存で変化するため、距離はもはや上記で定義できない。
そこで、微小区域に着目すれぼ平坦だとみなして
$$ds^2=\sum_{ij}g_{ij}dX^idX^j$$
と定義する(線素とよび、この微小座標のことを局所座標系とよぶ)。ここから$X^i$のように、**反変ベクトル(テンソル)**であれば上付き添字で表すものとする。下付添字は共変テンソルをあらわす。下付と上付で同じ添字が使われた場合は和をとる(アインシュタインの縮約とよぶ)ルールで書き直すと、線素の定義式は
$$ds^2=g_{ij}dX^idX^j$$
となる。
大局的な距離は、線素を経路に沿って積分すればいいので、
$$I=\int ds=\int \sqrt{g_{ij}\dot X^i \dot X^j}dt$$
と表現できる。なお、$\dot X^i=\frac{dX^i}{dt}$である。
このように、距離が積分で定義可能だが一様でない計量で定義された空間をリーマン多様体とよぶ。
メルカトル座標系の計量
メルカトル図法は、経度直交かつ等角性を保つために、局所的に以下を満たすのであった。
$$du=d\phi$$
$$dv=\frac{d\theta}{\sin\theta}$$
いま$u^i=(u,v)$、$\theta^i=(\phi,\theta)$とすると線素は、
$$ds^2=g_{ij}'d\theta^id\theta^j=g_{ij}du^idu^j$$
とできる。$g'$は球座標上の微小変化が球面上でどの程度の変化につながるかを考えればよく、
$$g_{ij}'=\begin{pmatrix}\sin^2\theta &0\
0&1\end{pmatrix}$$
で表現できる。
メルカトル図法上の計量は、この線素の関係式を満たすように係数比較すれば求まる。
$$g_{11}'d\phi d\phi=\sin^2\theta d\phi d\phi=\sin^2\theta dudu$$
$$g_{22}'d\theta d\theta=\sin^2\theta dvdv$$
$$g_{12}d\phi d\theta=g_{21}d\theta d\phi=0$$
$$\therefore g_{ij}=\sin^2\theta\delta_{ij}$$
極小値計算
求めた計量を距離計算式に代入して、
$$I=\int \sin\theta \sqrt{\dot u^2+\dot v^2}dt$$
この式の最小値を求めよう。
まず、右辺に$\theta$と$u,v$が共存していると見通しが悪いので、なんとかしよう。
幸い、これは比較的変形が容易いことがわかっている。
$$v=\int \frac{1}{\sin\theta}d\theta=\ln\tan\frac{\theta}{2}$$
$$\sin\theta=\sin(2\arctan(\exp(v)))$$
$p=\tan t$のとき、
$$\sin t=\frac{p}{\sqrt{1+p^2}},\quad\cos t=\frac{1}{\sqrt{1+p^2}}$$
$$\sin2t=\frac{2p}{1+p^2}=\frac{2\tan t}{1+\tan^2 t}$$
であるので、
$$\sin\theta=\sin(2\arctan(\exp (v)))=\frac{2\exp(v)}{1+\exp(2v)}=\frac{1}{\cosh(v)}$$
つまり、
$$I=\int \frac{\sqrt{\dot u^2+\dot v^2}}{\cosh v}$$
$$\hat u(t),\hat v(t)=\arg\min_{u,v} I, \quad (\hat u(t_0),\hat v(t_0))=(u_0,v_0)\quad (\hat u(t_1),\hat v(t_1))=(u_1,v_1)$$
を解けばいい。何通りか解き方があるので、全部試してみよう。
変分問題として解く
またリーマン幾何の話に戻ろう。
ある計量で定義された空間の、二点間の最短経路を求める式は、結論を先にいえば測地線方程式である。やってみよう。
$$I=\int \sqrt{g_{ij}\dot X^i\dot X^j}dt=\int L(X,\dot X)dt$$
$$\delta I=\int \sum_{q=X,\dot X}\frac{\partial L}{\partial q}\delta q dt=0$$
端点の変分は固定なのでEuler-Lagrange方程式に代入。
$$\frac{d}{dt}\left(\frac{\partial L}{\partial \dot X}\right)-\frac{\partial L}{\partial X}=0$$
ここで、計量は座標のみに依存する。
※ここから先はwikiとだいたい同じ変形たが、wikiの変形は若干不親切に感じるので、敢えて導出を続ける。
第一項から変形していこう。
$$Nom_1=\frac{\partial L}{\partial \dot X^\mu}=\frac{g_{\mu j}X^j+g_{i\mu}X^i}{2\sqrt{g_{ij}\dot X^i\dot X^j}}$$
ここで、線素
$$ds^2=g_{ij}dX^idX^j$$
をつかえば、これは時間にしか依存しないので
$$\left(\frac{ds}{dt}\right)^2=g_{ij}\dot X^i\dot X^j$$
である。ゆえに
$$Nom_1=\frac{1}{2}\frac{dt}{ds}g_{\mu i} g_{i\mu}\dot X^i=\frac{g_{\mu i}+g_{j \mu}}{2}\frac{dX^i}{ds}=g_{\mu i}\frac{dX^i}{ds}$$
ここで、縮約の添字は入れ替え可能であることと計量は対称行列で定義できることを利用している。ゆえに、
$$\frac{d}{dt}\left(Nom_1\right)=\frac{ds}{dt}\frac{d}{ds}g_{\mu i}\frac{dX^i}{ds}=\frac{ds}{dt}\frac{\partial g_{\mu i}}{\partial X^j}\frac{dX^i}{ds}\frac{dX^j}{ds}+\frac{ds}{dt}g_{\mu i}\frac{d^2X^i}{ds^2}$$
と変形できる。第二項も変形してみよう。
$$Nom_2=\frac{dL}{dX^\mu}=\frac{1}{2}\frac{dt}{ds}\frac{\partial g_{ij}}{\partial X^\mu}=\frac{1}{2}\frac{ds}{dt}\frac{\partial g_{ij}}{\partial X^\mu}\frac{dX^i}{ds}\frac{dX^j}{ds}$$
さて、それぞれ代入すると以下になる。
$$\frac{ds}{dt}\left(g_{\mu \nu}\frac{d^2X^\nu}{ds^2}+g_{\mu \nu}\gamma^\nu_{ij}\frac{dX^i}{ds}\frac{dX^j}{ds}\right)==\frac{ds}{dt}g_{\mu\nu}\left(... \right)=0$$
$$\gamma^\nu_{ij}=\frac{1}{2}g^{\nu \alpha}\left(\frac{2\partial g_{\alpha i}}{\partial X^j}-\frac{\partial g_{ij}}{\partial X^\alpha}\right)$$
ここで、計量の反変テンソルは計量の逆行列であり、クロネッカーのデルタを用いて
$$g_{\mu\nu}g^{\nu \alpha}=\delta_\mu^\alpha$$
と定義される。
計量テンソルはゼロ因子行列ではないので共通因数を払うと、
$$\frac{d^2X^\mu}{ds^2}+\gamma^\mu_{ij}\frac{dX^i}{ds}\frac{dX^j}{ds}=0$$
計量が対称行列であり$g_{i\alpha}=g_{\alpha i}$であることから、$\gamma$の代わりに以下のクリストッフェル記号を使っても上記方程式は等価である。
$$\left(\because \gamma^\mu_{ij}\frac{\partial X^i}{\partial s}\frac{\partial X^j}{\partial s}=\Gamma^\mu_{ij}\frac{\partial X^i}{\partial s}\frac{\partial X^j}{\partial s}\right)$$
$$\Gamma^\mu_{ij}=\frac{1}{2}g^{\mu\alpha}\left(\frac{\partial g_{\alpha j}}{\partial X^i}+
\frac{\partial g_{i\alpha}}{\partial X^j}-\frac{\partial g_{ij}}{\partial X^\alpha}\right)$$
$$\frac{d^2X^\mu}{ds^2}+\Gamma^\mu_{ij}\frac{dX^i}{ds}\frac{dX^j}{ds}=0$$
これを測地線方程式とよぶ。
メルカトル図法の測地線方程式
メルカトル図法の計量は
$$g_{ij}=\sin^2\theta\delta_{ij}$$
であることが分かっているので、
$$g_{ij}=\frac{\delta_{ij}}{\cosh^2 v}$$
ただし、$v=e_{2;\nu}X^\nu$つまり、$v$方向単位ベクトルとの内積である。
$$\frac{\partial g_{ij}}{\partial X^k}=\frac{-2\sinh v}{\cosh^3 v}\delta_{ij}e_{2;k}=-2\tanh(v)g_{ij}e_{2;k}$$
$$g^{ij}g_{kl}=\delta_{kl}\delta^{ij}$$
これを組みあわせるとクリストッフェル記号は
$$\Gamma^\mu_{ij}=-\tanh v\delta^{\mu\alpha}(\delta_{\alpha j}e_{2;i}+\delta_{i\alpha}e_{2;j}-\delta_{ij}e_{2;\alpha})=-\tanh v(\delta^\mu_je_{2;i}+\delta^\mu_ie_{2;j}-\delta_{ij}e^\mu_{2;})$$
測地線方程式に代入すると、
$$\frac{d^2X^\mu}{ds^2}-\tanh v\left(2\frac{dv}{ds}\frac{dX^\mu}{ds}-e^\mu_{2;}\delta_{ij}\frac{dX^i}{ds}\frac{dX^j}{ds}\right)=0$$
これにそれぞれの単位ベクトルと内積をとれば、
$$\ddot u-2\tanh(v)\dot v\dot u=0$$
$$\ddot v-\tanh
(v)(\dot v^2-\dot u^2)=0$$
この方程式はそれぞれ2階微分について解けるので、別記事で触れた2解微分用のRunge-Kutta法を用いて数値解を計算できる。(あるいは、頑張ればこの式の一般解は計算できるのかもしれないが、そこは頑張らないことにする)
前記事でシミュレーションした開始位置と開始進路を初期条件として数値解を求めてみよう。
前記事で定義したRung-Kutta関数で方程式を解くには、二階微分方程式を変形して二階微分係数が出力されるような関数を定義してやればいいのであった。そういえば、微分方程式の記事で作ったRunge-Kutta関数は方程式がベクトル量で定義されていても成立する形状になっていた。やってみよう。
def func(t, x, dx):
ddu = 2*np.tanh(x[1])*dx[0]*dx[1]
ddv = np.tanh(x[1])*(dx[1]**2-dx[0]**2)
return np.stack((ddu,ddv))
さあ、数値計算で曲線計算だ。
x=np.array((u[0],v[0]))
dx=np.array(((u[1]-u[0]), (v[1]-v[0])))
h=0.003
xs=[]
for _ in range(10000):
xs.append((x, dx))
x,dx=solver_rungeKutta(t, x, dx, func, h)
t+=h
xs=np.stack(xs).T
ここで、u,vは前回記事で求めたメルカトル図法上の最短経路線を使う(ただし、スケールは描画用に係数をかける前の値を使う)。
前記事で求めた経路線と、上記計算で求めた線を重ねよう。
橙線が前回記事の線で、青線が今回求めた数値解である。というわけで、測地線方程式を使えばメルカトル図法の地図上で最短経路を直接求めることは可能ということが分かったといえる。
前記事とどっちが楽かといわれれば、圧倒的に前記事だろうけど。
それと、この解き方にはちょっとズルがあるので、ズルしない解き方についても次記事で触れよう。