season's quarterly

数学/物理/プログラミング

座標不変プログラミング

計算機における座標不変性

ディスプレイのピクセルは縦横に並んでいるため、コンピュータ上の位置の指定は最終的にはデカルト座標で行う必要がある。物理シミュレーションのライブラリでも物体の運動はデカルト座標で指定するのが普通だ。しかしながら高レイヤ側では座標系をハードウェアに依存せず自由に取れる方が好ましい。デカルト座標への変換を与えたらプログラムが自動でデカルト座標での値を求めるようにしたい。

ハミルトン力学

ハミルトン形式の解析力学によれば、一般化座標 q_i、一般化運動量 p_i、ハミルトニアン H(p_i, q_i)に対して、正準方程式

\frac{dq_i}{dt} = \frac{\partial H}{\partial p_i}

\frac{dp_i}{dt} = -\frac{\partial H}{\partial q_i}

が成り立つ。ハミルトニアンの偏微分が求められれば座標と運動量の時間微分が分かるので、適当に前進オイラー法で一般化座標と一般化運動量を求めた上でデカルト座標に変換する。このようにすることで座標に依存しないプログラムが可能となる。

つまり、簡単のため2次元で考えると、予めユーザーが与える情報は

  1. 座標変換 x(q_1, q_2), y(q_1, q_2)
  2. デカルト座標系におけるハミルトニアン H(x, y, p_x, p_y)

であり、プログラムが実行する内容は

  1. 座標変換によりデカルト座標のハミルトニアン H(x, y, p_x, p_y)を一般化座標のハミルトニアン H(q_1, q_2, p_1, p_2)に変換。
  2. ハミルトニアンの偏微分を計算することで、一般化座標と一般化運動量の時間微分を計算。
  3. 前進オイラー法で各時刻の q_1, q_2, p_1, p_2を求める。
  4. 座標変換から各時刻の x, yを求め、描画ライブラリに渡す。

である。

2, 3, 4については、変数を計算グラフのノードで表し、記号微分と評価を行えば良い。1については、 x, yは自明に q_1, q_2で表現できるので、 p_x, p_yが q_1, q_2, p_1, p_2で表現可能か考える。

p_x, p_yと q_1, q_2, p_1, p_2の関係を導く。一般化運動量はラグランジアン L(q_i, \dot{q_i})を用いて

p_i = \frac{\partial L}{\partial \dot{q_i}}

と定義される。デカルト座標における運動量は

p_x = \frac{\partial L}{\partial \dot{x}} = \frac{\partial L}{\partial q_1}\frac{\partial q_1}{\partial \dot{x}} + \frac{\partial L}{\partial q_2}\frac{\partial q_2}{\partial \dot{x}} + \frac{\partial L}{\partial \dot{q_1}}\frac{\partial \dot{q_1}}{\partial \dot{x}} + \frac{\partial L}{\partial \dot{q_2}}\frac{\partial \dot{q_2}}{\partial \dot{x}} = p_1\frac{\partial \dot{q_1}}{\partial \dot{x}} + p_2\frac{\partial \dot{q_2}}{\partial \dot{x}}

ここで

\frac{dq_1}{dt} = \frac{\partial q_1}{\partial x}\frac{dx}{dt} + \frac{\partial q_1}{\partial y}\frac{dy}{dt}

なので

\frac{\partial \dot{q_1}}{\partial \dot{x}} = \frac{\partial q_1}{\partial x}

つまり

p_x = p_1 \frac{\partial q_1}{\partial x} + p_2 \frac{\partial q_2}{\partial x}

p_y = p_1 \frac{\partial q_1}{\partial y} + p_2 \frac{\partial q_2}{\partial y}

(p_x p_y) = (p_1 p_2) \frac{\partial (q_1, q_2)}{\partial (x, y)} = (p_1 p_2) \frac{\partial (x, y)}{\partial (q_1, q_2)}^{-1}

ヤコビアンは座標変換から記号微分によって計算可能であり、 p_x, p_yを q_1, q_2, p_1, p_2で表すことができる。

実装

例として単振り子の物理シミュレーションをHaskellで実装する。

まず計算グラフのノードを定義して、値を評価する関数evalと微分を計算する関数diffを実装。

data Node = Num Float
          | Var String
          | Add Node Node
          | Sub Node Node
          | Mul Node Node
          | Div Node Node
          | Pow Node Int
          | Neg Node
          | Exp Node
          | Log Node
          | Sin Node
          | Cos Node
          | Tan Node

eval :: Map.Map String Float -> Node -> Float
diff :: Node -> String -> Node

evalの第一引数は系の状態（一般化座標/運動量）を変数名と値の対応として表す。

次に上で示したように、座標変換を定義し、一般化座標/運動量からデカルト座標における運動量を計算。

    -- coodinate transformation
    let x = q1 * cos q2
        y = q1 * sin q2

    let a11 = diff x "q1"
        a12 = diff x "q2"
        a21 = diff y "q1"
        a22 = diff y "q2"

    let det = a11 * a22 - a12 * a21

    let b11 = a22 / det
        b12 = -a12 / det
        b21 = -a21 / det
        b22 = a11 / det

    let px = p1 * b11 + p2 * b21
        py = p1 * b12 + p2 * b22

単振り子のハミルトニアンを定義し、正準方程式から状態の時間微分を求める。

    -- Hamiltonian in Cartesian coodinate
    let h = (px /\ 2 + py /\ 2) / (2 * m) + m * g * y

    -- Euler method
    let dotq1 = diff h "p1"
        dotq2 = diff h "p2"
        dotp1 = - diff h "q1"
        dotp2 = - diff h "q2"

前進オイラー法で現在の状態から次の状態を計算する関数を定義。単振り子の場合は張力によって半径が一定になる拘束条件があるので、q1の更新をコメントアウトしている。

type Model = Map.Map String Float

eulerMethod :: (Node, Node, Node, Node) -> Float -> Model -> Model
eulerMethod (dotq1, dotq2, dotp1, dotp2) dt (state) = do
    let q1 = state Map.! "q1"
        q2 = state Map.! "q2"
        p1 = state Map.! "p1"
        p2 = state Map.! "p2"

    let q1Next = q1 -- + eval state dotq1 * dt
        q2Next = q2 + eval state dotq2 * dt
        p1Next = p1 + eval state dotp1 * dt
        p2Next = p2 + eval state dotp2 * dt

    let stateNext = Map.insert "q1" q1Next . Map.insert "q2" q2Next . Map.insert "p1" p1Next . Map.insert "p2" p2Next $ state

    stateNext

最後に初期条件と描画関数を与えてglossでアニメーションさせる。

-- initial model
    let initModel = Map.fromList [("q1", 150), ("q2", -pi / 3), ("p1", 0), ("p2", 0)]

    simulate window white 24 initModel (draw (x, y)) (\_ -> eulerMethod (dotq1, dotq2, dotp1, dotp2))

動いた（精度が悪いのか振れ幅が少しずつ大きくなっている）。

先行研究

あった...

手順は

  1. 一般化座標とそのデカルト座標への変換を与える。
  2. デカルト座標でのポテンシャルエネルギーの表示を与える。

で全く同じですね。。。大部分完成した段階で見つけてしまったのでめっちゃ萎えた。そう...ラグランジュ形式だと無理なのでハミルトン形式でやる必要があって...Haskellで書きたくなるよね...いやーその通りです。基本的に作りたいものとかない人間なんですが、今回は物理×情報で双方の知見がバランス良く役に立つ題材（そんな題材めったに思い付かないだろう）で結構テンション上がってたのであ～って感じですね。まあそりゃあるか...って感じだけど。

やる気なくなったのであんまりリファクタリングしてません笑。

github.com