でこぼこ面の測地線

三次元で z=Z(x,y)と表される面上の最短経路=測地線を考えてみる。

まず一般的に、n次元での運動の速度の二乗が速度の各成分viを使って V^2 = \sum_{ij}g_{ij}v_i v_jという二次式で表されるとする。gij は計量(measure)。
たとえばz=Z(x,y)面での速度はxとyの速度を使って(v_x, v_y, \partial Z/\partial x v_x + \partial Z/\partial y v_y)とかける。\partial Z/\partial x =Z_xなどとおくと、 V^2 = v_x^2 + v_y^2 + (v_x Z_x + v_y Z_y)^2 = (1+Z_x^2)v_x^2 + (1+Z_y^2)v_y^2 + (2 Z_x Z_y)v_x v_y
したがってこの場合、gxx=1+Zx^2, gyy=1+Zy^2, gxy=gyx=Zx Zy。まとめると gij=δij +Zi Zj。

計量gijの空間を一定速度で測地線にそって進む物体の運動方程式は、各成分の加速度が速度の二次式で表されるものになる。この二次式の係数がクリストッフェル記号。
a_i = -\sum_{jk} \Gamma^i_{jk} v_j v_k
 \Gamma^i_{jk} = 1/2 \sum_l (g^{-1})_{il} (g_{lj,k}+g_{lk,j}-g_{jk,l})
ここでg^{-1}_{ab}はgijの逆行列のab成分の意味。また
g_{ab,c} = \frac{\partial g_{ab}}{\partial x_c}
http://en.wikipedia.org/wiki/Geodesic

gij=δij +Zi Zjである場合、gij,k = Zik Zj+Zjk Zi。ただしZijはZをxiとxjで微分した二回微分で、Zij=Zji。 これをΓの式に代入すると
 \Gamma^i_{jk} = 1/2 \sum_l (g^{-1})_{il} (g_{lj,k}+g_{lk,j}-g{jk,l})= 1/2 \sum_l (g^{-1})_{il} (Z_{lk} Z_j + Z_{jk} Z_l + Z_{lj} Z_k + Z_{kj} Z_l - Z_{jl} Z_k - Z_{kl} Z_j)
 = \sum_l (g^{-1})_{il} Z_{jk} Z_l = Z_{jk} \sum_l (g^{-1})_{il} Z_l

二次元の場合、2x2行列の逆を計算して g^-1 ={1+ZxZx, ZxZy, ZxZy, 1+ZyZy}^-1 = {1+ZyZy, -ZxZy, -ZxZy, 1+ZxZx}/(1+ZxZx+ZyZy)
したかって \sum_l (g^{-1})_{il} Z_l は g^-1.(Zx, Zy) の第l成分だから計算すると (Zx,Zy)/(1+ZxZx+ZyZy)。
したがって

  • ax = - Zx/(1+Zx*Zx+Zy*Zy) * (Zxx * vx * vx + Zyy * vy * vy + 2 * Zxy * vx * vy)
  • ay = - Zy/(1+Zx*Zx+Zy*Zy) * (Zxx * vx * vx + Zyy * vy * vy + 2 * Zxy * vx * vy)

というわけで加速度は常に勾配の方向を向いている。

Z(x,y)= 1/(x*x+y*y) という発散するトゲの場合に、中心近くからいろんな角度へ発射した場合の軌跡

穴との角度が60度以下の場合は吸い込まれて何度も巻きつく。70,80度の場合は何回か巻いてから脱出する。90度以上の場合は巻かずに離れる。

レポート問題

シュワルツシルト計量の二次元極座標表示

  • gtt = -(1-1/r)
  • grr = 1/(1-1/r)
  • ghh = r^2/(1-1/r)

の場合に同様の計算を行い、自由落下の経路 (t(s),r(s),h(s))を計算しトラックバックを送信せよ。hはθのつもり。
rとtをペンローズ図で使う変数に変換してからのほうがいいのかなぁ。http://d.hatena.ne.jp/ita/20051106/p1
普通にやると、落ちる場合はr=1付近でいろいろ発散して、経路がt方向に向いてしまい事象の地平は越えられない。外から見てる分には正しいけど、落ちる人の主観では有限時間で越えてr=0まで行く。