ひとり勉強会

ひとり楽しく勉強会

SICM 読書会 (1)

SICM

水曜日は Structure and Interpretation of Classical Mechanics (SICM, 古典力学の構造と解釈) の読書会をやります。

プログラミングの教科書として、SICP（計算機プログラムの構造と解釈）という有名な本があります。SICM は同じ著者による、古典力学についての教科書だそうです。SICP の文章中のプログラムが Scheme で書かれているのとおなじで、SICM の文章中のプログラムは Scheme で書かれているそうです。

・・・古典力学の本でプログラム？どうやら、方程式を記号的に記述したり変形したり計算したりするライブラリが用意されていて、それで実際に力学の式を動かしながらの解説が含まれているらしい。というわけで、ちょっと面白そうなので読んでみることにしました。

資料

SICM
//mitpress.mit.edu/SICM/">Structure and Interpretation of Classical Mechanics:本家本元です。HTMLで全文読めます。トップページの Looking for Scheme? のリンクからは、本で使われている力学計算のためのSchemeライブラリが入手できます。MIT/GNU Schemeの処理系ごとまとめて入ってます。
//groups.google.com.vn/group/sicm/">SICM の Google グループ:MIT/GNU Scheme on Linux 以外の環境で試すときの手がかりになりそうな情報。
//www.cs.rochester.edu/~gildea/guile-scmutils/">guile-scmutils:それで、オフィシャルな版をいろいろ頑張ったつもりなんですが全然動かせなかったので、Google グループで紹介されていたGuileへの移植版を使うことにしました。この読書会ではこちらを基本的に使っていくことにします。

ではスタート！

1 Lagrangian Mechanics

SICM

ラグランジュ力学。

物体がたくさんあって、それぞれ相互作用したりしなかったりしてるシステムを記述したい。

Configuration
ある時点、時間軸上での一点での、システムの「状態」のこと。システム内の全物体の位置、みたいな情報のことだと思う。具体的にどういう変数で記述するとかは今はまだ考えない。
Configuration path
時間からConfigurationへの関数 (time->configuration)。時間経過にともなうシステムの動きの記述。

このシステムでは物体がこー動いてあー動きます、という全てを記述したのがconfiguration path、と。さて、configuration pathはただの関数なので、どんな関数を考えることもできる。けれど実際には、「現実にありうる動き」と「現実にありえない動き」があるでしょう。

そこで、現実にありえるかどうかを判定するための、configuration path から何かへの関数が欲しい。例えば、configuration path から真偽値への関数 (configuration path->bool)。ありえるなら true、ありえないなら false を返すとか。ニュートンの運動方程式はこの形とみなすことができる。configuration path が ma = F を満たしてたら true、なければ false を返すみたいに思えば。

でも、「現実にありうる／えない」を区別する手段として、もっと強力なものがあるらしい。

Path-distinguishing function
Configuration pathを引数にとって実数を返す(?)関数で、「ありうる動き」のところで最小（極小？）値をとるようなもののこと。
Variational strategy
システムが与えられたら、まず path-distinguishing function を見つけ出して、それを極小にするものを探す、という形で現実の動きを計算する方式。

ここからはこの方式について勉強です。

1.1 The Principle of Stationary Action

SICM

停留作用の原理。

Path-distinguishing function て、どういう関数になるのか考えてみよう。いやその前に、「ありうる動き」の特徴について考えてみよう。

  1. 連続。滑らか。いきなり物体が瞬間移動したりはしない。
  2. 履歴に依存してない。"現在のconfiguration"が完全に同じなら、そこに至るまでの過程が違ってもその後は同じ。
  3. 決定的。現在のconfigurationが決まれば未来が確定的に決まる。
  4. 局所的。「ありうる動き」の一部分(path segment)は「ありうる動き」だし、逆にすべての部分が「ありうる動き」なら全体も「ありうる動き」だろう。path segmentの「ありうる度」は、path segment内の点の情報のみで定まるし、逆にpath segment内の点全ての情報を使って定まるはず。

量子レベルの話は考えないよーとのこと。

こう考えると、path-distinguishing function は、configuration path segment の各点各点での情報を平等に計算して足し合わせたものになるのではないか。つまり積分。各点各点での情報の計算の仕方 F はまだ不明だが、ともかく、Path-distinguishing function S はこういう形をしてるのではないか。
　　　S(\gamma)(t_1,t_2) = \int_{t_1}^{t_2} F(\gamma)
開始時刻 t_1 と終了時刻 t_2 を固定すると、S は \gamma (configuration path) から実数への関数になる。この値が極小になるような\gammaが、その時刻間で一番「ありうる動き」を与える。

さらに、Fは局所性があるはず。path \gammaそのものじゃなくて、\gammaから計算できる各時刻tでの情報のみを使って定義されてるべきだ。T(\gamma)(t) = (t,\gamma(t), (D\gamma)(t), ...) というタプル（Dは微分演算子）の関数として定義されるべき。
　　　F(\gamma)(t) = L(T(\gamma)(t))
L : (time,configration,Dconfigration,...)->real は系に応じて定義される、各点各点での情報を計算する関数。まとめると、
　　　S(\gamma)(t_1,t_2) \mapsto \int_{t_1}^{t_2} L \circ T(\gamma)

Local tuple
T(\gamma)(t)のこと
Action
path-distinguish functionのことは、この本以外ではactionと呼ぶらしい
Lagrangian
L のこと
Lagrangian action
上のようにLを使って定義されたactionのこと

練習問題1.1

フェルマー曰く、光は２点間を最短時間で結ぶ道を通る。この法則から、光の反射の法則（入射角==反射角）と、屈折の法則（sin(屈折角)/sin(入射角)==光の速度@入射前/光の速度@屈折後）を説明せよ。

  • 始点→反射面→終点 という経路を最短にするには、終点を反射面に関して面対称の位置に持ってきた点（終点'）と始点を結ぶ最短経路を考えればよくて、それは始点→終点'の直線になる。で、その直線と反射面のなす角が入射角と反射角だけど、それらは対頂角なので等しい。
  • 屈折面がx軸、始点が(0,Ya)、終点が(X,-Yb)とする。光が(x,0)で屈折するとすると、進むのにかかる時間は以下の通りなので微分して0になる(=最短時間になる)xを求めるとその時ちょうど問題の比が等しくなっている。

\frac{\sqrt{x^2+{Ya}^2}}{Va}+\frac{\sqrt{(X-x)^2+{Yb}^2}}{Vb}
もっと微分とか出さないで説明できそうな気がするんだけど挫折しましたorz

1.2 Configuration Spaces

SICM
Configuration space
システムのとりうるConfiguration全ての集合
Dimension
システムのConfigurationを一意に特定するのに最小で何個のパラメタが要るか。自由度(degrees of freedom)ともいう。

例えば粒子が一個自由に飛んでるシステムの自由度は３。二個自由に飛んでたら６。二個飛んでるんだけど、その二個の距離は常にある定数、みたいな制約がかかってたら、自由度は５。

練習問題1.2

いろんなシステムの自由度を考えてみよう

  • a. 各pinの自由度6（位置3、向き3）*3 で 18。
  • b. 球面上を自由に動けるので緯度と経度みたいな表現で、2。重力は関係ないよね??
  • c. bが2個つながってるので、4。
  • d. ワイヤー上の位置は変数1個で表せるので、1。
  • e. 傾きと、軸が傾いてる向きで、2。
  • f. 軸対象じゃないと、さらにコマ自体の向きが加わって、3、かな。

1.3 Generalized Coordinates

SICM

一般座標。

自由度がnの系のconfigurationは、n個(以上)の変数で記述できる。簡単のためにしばらくは、常にピタリn個で記述する場合のみを考える。その方が変数間の制約などを考えなくていいので楽。

変数は、「各粒子のxyz座標」のようなお決まりの取り方じゃなくて、システムを記述できるような変数ならどうとっても自由。そういうように自由にとった座標系を一般座標という。いろいろ記号を導入。

  • \chi^i : configuration->第i座標
  • \chi(m) = (\chi^0(m), \chi^1(m), \dots, \chi^{n-1}(m))
  • q^i = \chi^i \circ \gamma : time->第i座標
  • q = \chi \circ \gamma
  • local tuple からその一般座標表示を返す関数 "chart"_\chi。TeXでどう出すのかわからない。。。
  • \Gamma(q) は T(\gamma) の座標表示版。時間tを受け取るとlocal tupleの座標表示を返す。
  • L_\chi : (real, real^n, real^n, ...) -> real はラグランジアン L の座標表示版

まとめると、path-distinguishing function を一般座標表示すると
S_\chi(q)(t_1,t_2) = \int_{t_1}^{t_2}L_\chi\circ\Gamma(q)
こうなるわけです。

練習問題1.3

1.2のシステムに一般座標を与えてみよう。

1.2のとこで書いたパラメタを座標にすればよさそう。

1.4 Computing Actions

SICM

さて、ここまでの知識を実際にSchemeで動かしてみようタイムです。guile版を準備。

guile> (load "/usr/local/scim/load.scm")
guile> (set-current-module generic-environment)
guile> (define D derivative)

一番簡単な例として、１個の粒子が自由に飛んでいるケースを考えてみよう！とのこと。座標は、普通にx,y,zにしましょう。q(t)=(x(t),y(t),z(t))

(define x (literal-function 'x))
(define y (literal-function 'y))
(define z (literal-function 'z))
(define (q t) (up (x t) (y t) (z t)))

昔の偉いひとによると、このシステムのラグランジアンは、物体の速さをv、質量をmとすると\frac{1}{2}mv^2だそうです。

(define ((L m) q)
  (let ((v (velocity q)))
    (* 1/2 m (dot-product v v))))

で、actionの積分の中身は L \circ \Gamma(q) でした。

(define sekibun-nakami (compose (L 'm) (Gamma q)))

いきなり変数名が酷くなりましたが気にしないことにしましょう。あ、ちなみに、compose, Gamma, velocity, dot-product などは全てライブラリ側で提供されている関数です。

guile> (pe (sakibun-nakami 't))
(+ (* 0.5 m (expt ((D x) t) 2))
   (* 0.5 m (expt ((D y) t) 2))
   (* 0.5 m (expt ((D z) t) 2)))

積分の中身を表示してみました。pe = print-expression。
このライブラリは、こういう風に記号的に式を式として扱うこともできますし、実際に計算させることもできます。

guile> (define (fp-action q t1 t2)
         (definite-integral (compose (L 3) (Gamma q)) t1 t2))

configuraion path と開始終了時刻を与えると、重さ 3 の粒子の action を計算してくれる関数。直線状にぶっとんで行くpathを与えると…

guile> (fp-action
         (lambda (t) (up (+ (* 4 t) 7)
                         (+ (* 3 t) 5)
                         (+ (* 2 t) 1))) 0 10)
435.0

最小作用のpath

ここまでの知識によれば、fp-actionは、「ありえる動き」のところで極小値をとるはずです。何も他に要因がなければ、一定速度でどこまでも飛んでいくのが「ありえる動き」で、それ以外の、微妙に速度が変わったりするpathは「ありえない動き」でしょう。というわけで、そういうpathを与えると、値は435.0より大きくなるはず。

guile> (fp-action
         (lambda (t) (up (+ (* 4 t) 7)
                         (+ (* 3 t) 5)
                         (+ (* 2 t t) 1))) 0 10)
8375.00000000001

何故かz軸方向にものすごい加速してみるとか。。。ちょっと極小とか以前にこれは動きが違いすぎますかね。

guile> (define ((v-action q t1 t2) eps)
         (fp-action
           (lambda (t) (+ (q t) (up 0 0 (* eps (- t t1) (- t t2)))))
           t1 t2))
guile> ((v-action
         (lambda (t) (up (+ (* 4 t) 7)
                         (+ (* 3 t) 5)
                         (+ (* 2 t) 1))) 0 10) 0.01)
435.05

qとepsを与えると、微妙に揺れてるpath q'(t) = (x(t), y(t), z(t)+{eps}(t-t1)(t-t2)) に関するactionを計算してくれる関数です。時刻t1とt2では揺れが0になるようにしてあるので、始点終点がちゃんと同じconfigurationになります。微妙にずらすと微妙にactionが増えてますね。

guile> (minimize (v-action
...          (lambda (t) (up (+ (* 4 t) 7)
...                          (+ (* 3 t) 5)
...                          (+ (* 2 t) 1))) 0 10) -2 2)
(-1.11022302462516e-16 435.0 5)

上のようにminimize関数を使うと、-2〜2の範囲で、v-actionを最小にするepsを探すことができます。実際、epsがほぼ0、つまり直線のところでactionが最小になってるようです。

まとめ

SICM

力尽きたのでこの辺りで今回はおしまいにします。

次回は1.4の途中「最小作用のpathを探す」節から再開します。上では、直線で動くのが「ありうる動き」だとわかってたので、あー確かにその付近では最小だねーという「検算」ができました。では、元々よくわかってない時はどうするか？というお話。

この調子だと本の翻訳というか写しになってしまいますね。どうしたものか。もっと難しくなってくれば状況も違うかなあ。。。