かたちづくり

つれづれに、だらだらと、おきらくに

【解析入門Ⅰ 杉浦】第Ⅰ章 実数と連続（メモ）

備忘録として読書メモ。今まで無理に背伸びして多様体だのホモロジーだの齧ってみたものの、やはり基礎基本が分かっていないと身につかず表面だけをなぞって終わってしまうことを実感。では集合位相論を、とトライしたが抽象的すぎて迷子になり途中で有耶無耶に。やはりある程度具体的な実数体Rで位相的な性質を理解した上で抽象論に進まないといけないようだ。というわけで本書を購入。

解析入門 (1)

解析入門 (1)

※ これは自分用のメモです。他人に分かるように書きませんのであしからず。

1. 実数

  • 実数Rの公理が示される。
  • 四則演算の定義 → Rは体
  • 比較演算の定義 → Rは順序体
  • 連続の公理
    • \forall A \subset R, A \neq \emptyset \land Aは上に有界 \Rightarrow \sup A \in R
  • 有理数Qも順序体
  • 実数Rと有理数Qの違いは連続の公理にある
※

以降全体に言えることだが、実数の数直線による直感的イメージに頼りすぎると全てが自明なことに思えてしまって何をやっているのか分からなくなってしまう。数学ガールで言うところの「知らないふりゲーム」をやらなくてはならない。つまり、直感的イメージは捨てて（知らないふりをして）手元にある道具は公理系だけという視点で読まなくてはならない。数直線の直感的イメージを隙のない論理で厳密に表現するとどうなるのか、という視点で読まなくてはならない。

2. 実数列の極限

【定義】数列の収束

\lim_{n \to \infty} a_n = a \Leftrightarrow (\forall \epsilon > 0)(\exists n_0 \in N)(\forall n \in N)(n \ge n_0 \Rightarrow \left|a - a_n\right| < \epsilon)

【定義】数列の発散

\lim_{n \to \infty} a_n = +\infty \Leftrightarrow (\forall M \in R)(\exists n_0 \in N)(\forall n \in N)(n \ge n_0 \Rightarrow a_n > M)

【命題】はさみうちの原理

\left[(\forall n \in N)(a_n \le c_n \le b_n) \land \lim a_n = \lim b_n = a \right] \Rightarrow \lim c_n = a

3. 実数の連続性

【定理】上に有界な単調増加数列は上限に収束する

(a_n)_{n \in N}が上に有界な単調増加列 \Rightarrow a_n \to \sup \{a_n|n \in N\} (n \to \infty)

【定理】アルキメデスの原理

(\forall a, b > 0)(\exists n \in N)(na > b)
すなわち
単調増加列 (na)_{n \in N} は上に有界ではない
すなわち
a > 0 \Rightarrow \lim_{n \to \infty}(na) = +\infty

※
  • どんなに小さな a、どんなに大きな b を選んでも na > b となる n が存在する
  • 自明な定理に見えるが、この原理が成立しない順序体が存在する
  • 例えば有理式 R(t) に辞書順を入れると順序体になるが、アルキメデスの原理は成立しない

【定理】区間縮小法

二分探索アルゴリズムの基礎となる定理。

  • I_n = [a_n, b_n]（有界閉区間）とする
  • \forall n\in N(I_n \supset I_{n+1} \Rightarrow \cap_{n \in N}I_n \neq \emptyset)
  • 特に \lim_{n\to\infty}(b_n - a_n) = 0 \Rightarrow \cap_{n\in N}I_n=\{\exists a\}
証明
  • a_0 \le a_1 \le \ldots \le a_n \le b_n \le \ldots \le b_1 \le b_0
  • つまり a_n は上に有界な単調増加列、 b_n は下に有界な単調減少列
  • つまり a_n, b_n は収束する： \lim_{n\to\infty}a_n=a, \lim_{n\to\infty}b_n=b
  • \lim_{n\to\infty}(b_n - a_n) = 0 ならば a = b
  • c_n = (a_n + b_n)/2 とすれば、はさみうちの原理により \lim_{n\to\infty}c_n=a
※

とかく抽象的になりがちな議論の中で、この定理だけが（プログラマである自分にとっては）妙に具象的なものに感じられて分かりやすい。この区間縮小法を利用した証明を読んでいると、二分探索によるアルゴリズムを読んでいるような気分になる。

【定理】ボルツァーノ・ワイヤストラスの定理

任意の有界実数列は常に収束する部分列を持つ

証明
  • (a_n)_{n\in N} は有界 \Leftrightarrow \forall a_n\in I=[b, c] なる区間 I が存在する
  • I_0=Iとして数列の要素を無限に含む区間を二分探索すること（区間縮小法）により、収束する部分列を選び取れる

【定義】コーシー列

\lim_{m,n\to\infty}(a_m - a_n) = 0
丁寧に書くと
(\forall \epsilon>0)(\exists n_0\in N)(\forall m,n\in N)(m,n\ge n_0 \Rightarrow \left|a_m - a_n\right|<\epsilon)

【定理】コーシーの収束条件

実数列が収束する ⇔ 実数列はコーシー列である

証明
  • 十分条件（⇐）の証明にボルツァーノ・ワイヤストラスの定理を使う
  • コーシー列は有界である（命題3.5(1)）
  • したがって、コーシー列は収束する部分列を持つ（ボルツァーノ・ワイヤストラスの定理）
  • コーシー列のある部分列が収束すれば、元のコーシー列も収束する（命題3.5(2)）
  • したがって、任意のコーシー列は収束する

4. R^nとC

5. 級数

（あとで書く）

6. 極限と連続

（あとで書く）

7. コンパクト集合

以前に位相論の本でコンパクト集合の定義を読んだ時には、コンパクトの概念というかイメージがいまひとつ掴めなかった。
実数体でコンパクト集合のイメージを掴むことが今回の目的の一つ。

【定義】点列コンパクト

R^nの部分集合Kが点列コンパクト ＝ Kの任意の点列がKの点に収束する部分列を含む
※ これはボルツァーノ・ワイヤストラスの定理の概念を R^n に拡張したもの

【定義】コンパクト

  • 集合族 (U_\lambda)_{\lambda\in L}がKの被覆 ＝ K \subset \cup_{\lambda\in L}U_\lambda を満たす
  • すべての U_\lambda が開集合の場合は「開被覆」という
  • Kの任意の開被覆から有限個の開集合を選んでKを被覆できるとき、コンパクトという

【定理】点列コンパクト ⇔ コンパクト ⇔ 有界閉集合

KをR^nの部分集合とする。次の3つは同値。

  • K は点列コンパクトである
  • K はコンパクトである
  • K は有界閉集合である
※
  • コンパクトの定義においては、「任意の」開被覆から有限個の開集合を選べるところがミソである。ある特定の開被覆から有限個の開集合を選べるだけではコンパクトとはいえない。
  • 開集合はコンパクトではない。つまり開集合の場合には、開被覆の作り方によっては開集合をどう選んでも元の開集合を有限個では被覆しきれない。例えば開集合の境界に無限小のε近傍を無限個並んでいるような開被覆を考えると、そこから有限個の開集合をどう選んでも全体を被覆することは出来ないことが分かる。
  • 点列コンパクトは、点列の収束の定義において距離の演算に依存している。一方、コンパクトの定義は開集合による被覆の概念のみで構成されており、距離の概念からは独立している。つまり点列コンパクトや有界閉集合の概念を「距離が定義されない世界」に拡張した概念が「コンパクト」といえるのではないか。

8. 中間値の定理

【定義】連続写像

f は a で連続 ⇔ \lim_{x \to a}f(x)=f(a) が存在

【定理】中間値の定理

次の3つを考える。

  • 実数値関数 f : R \to R
  • 閉区間 I = [a, b]
  • 閉区間 J = [\min(f(a), f(b)), \max(f(a), f(b))]

このとき以下が成り立つ。
f が I で連続 \Rightarrow (\forall \gamma\in J)(\exists c\in I)(f(c)=\gamma)

（あとで書く）

某 C# コードを高速化してみた

aokomoriuta さんが面白いパフォーマンステストをやっていました。aokomoriuta.hateblo.jp

こちらが aokomoriuta さんによる計測結果。

SIMDとかは私は分からないので華麗にスルーして、C# もうちっと速くならんかなー、と思いまして。ゴニョゴニョやった結果、割と高速化出来ました。

結果

  • 高速化前：6秒くらい
  • 高速化後：1秒くらい 1.3秒くらい

※ 当初はX,Y,Zしか計算してなくてWを忘れていたので、修正しました＞＜

やったこと

  • Vector4 を class から struct に変えた。
  • Vector4 の中身も配列をやめて X, Y, Z, W の4変数を持たせるようにした。
  • unsafe コードによって Vector4 の配列をポインタで扱うようにした。
  • ちなみに Vector4 に定義した演算子のオーバーロードを利用するとパフォーマンスは劇的に落ちて6.5秒くらい。

感想

書いたコード

※ 当初は W の計算を忘れていたので修正しました。

  struct Vector4
  {
    public double X, Y, Z, W;

    public Vector4( double x, double y, double z, double w = 0.0 )
    {
      this.X = x;
      this.Y = y;
      this.Z = z;
      this.W = w;
    }

    public static Vector4 operator +( Vector4 v1, Vector4 v2 )
    {
      return new Vector4( v1.X + v2.X, v1.Y + v2.Y, v1.Z + v2.Z, v1.W + v2.W );
    }

    public static Vector4 operator *( double scalar, Vector4 v )
    {
      return new Vector4( scalar * v.X, scalar * v.Y, scalar * v.Z, scalar * v.W );
    }
  }

  static class PerformanceTest
  {
    // こっちは遅い
    static unsafe void RunWithOperatorOverload( Vector4* x, Vector4* v, Vector4* f, double m, double dt, int n )
    {
      double tmp = dt * dt / 2;
      double rm = 1.0 / m;
      for ( int i = 0; i < n; i++ ) {
        // a = f/m
        var a = rm * f[i];

        // x += v*dt + a*dt*dt/2
        x[i] += dt * v[i] + tmp * a;

        // v += a*dt
        v[i] += dt * a;
      }
    }

    // こっちは速い
    static unsafe void RunWithoutOperatorOverload( Vector4* x, Vector4* v, Vector4* f, double m, double dt, int n )
    {
      double tmp = dt * dt / 2;
      double rm = 1.0 / m;
      for ( int i = 0; i < n; i++ ) {
        // a = f/m
        var a = new Vector4( rm * f[i].X, rm * f[i].Y, rm * f[i].Z, rm * f[i].W );

        // x += v*dt + a*dt*dt/2
        x[i] = new Vector4(
          x[i].X + dt * v[i].X + tmp * a.X,
          x[i].Y + dt * v[i].Y + tmp * a.Y,
          x[i].Z + dt * v[i].Z + tmp * a.Z,
          x[i].W + dt * v[i].W + tmp * a.W );

        // v += a*dt
        v[i] = new Vector4(
          v[i].X + dt * a.X,
          v[i].Y + dt * a.Y,
          v[i].Z + dt * a.Z,
          v[i].W + dt * a.W );
    }

    static Vector4[] CreateRandomVectors( int n )
    {
      var rand = new System.Random();
      var v = new Vector4[n];
      for ( int i = 0; i < n; ++i ) v[i] = new Vector4( rand.NextDouble(), rand.NextDouble(), rand.NextDouble() );
      return v;
    }

    public static unsafe void Run()
    {
      const int n = 100000;
      const int loop = 1000;
      const double dt = 0.1;
      const double m = 2.5;

      {
        var f = CreateRandomVectors( n );
        var v = CreateRandomVectors( n );
        var x = CreateRandomVectors( n );
        fixed ( Vector4* pX = x, pV = v, pF = f ) {
          Console.Write( "RunWithOperatorOverload: " );
          var stopwatch = new System.Diagnostics.Stopwatch();
          stopwatch.Start();
          for ( int i = 0; i < loop; i++ ) {
            RunWithOperatorOverload( pX, pV, pF, m, dt, n );
          }
          Console.WriteLine( "{0} [ms]", stopwatch.ElapsedMilliseconds ); // 6.5秒くらい
        }
      }
      {
        var f = CreateRandomVectors( n );
        var v = CreateRandomVectors( n );
        var x = CreateRandomVectors( n );
        fixed ( Vector4* pX = x, pV = v, pF = f ) {
          Console.Write( "RunWithoutOperatorOverload: " );
          var stopwatch = new System.Diagnostics.Stopwatch();
          stopwatch.Start();
          for ( int i = 0; i < loop; i++ ) {
            RunWithoutOperatorOverload( pX, pV, pF, m, dt, n );
          }
          Console.WriteLine( "{0} [ms]", stopwatch.ElapsedMilliseconds ); // 1.3秒くらい
        }
      }

      System.Console.WriteLine( "Press Any key..." );
      System.Console.ReadKey();
    }
  }

マックスウェル方程式を整理してみた（自分用メモ）

最近、この本で電磁気学を復習しました。

電磁気学の考え方 (物理の考え方 2)

電磁気学の考え方 (物理の考え方 2)

軽妙な語り口でするすると頭に入ってくる感じで、大変分かりやすい素晴らしい本でした。おすすめ。

さて、本書では後半に電磁気学の集大成たるマックスウェル方程式に辿り着くわけですが、そのマックスウェル方程式を自分なりにゴニョゴニョいじってみたので、自分用のメモ書きとしてここに記すものであります。

マックスウェル方程式

まずは教科書に載っていたマックスウェル方程式をそのまま転記。
\mathrm{div} \vec{B} = 0 \qquad \mathrm{rot} \vec{E} + \frac{\partial\vec{B}}{\partial t} = \vec{0}
\mathrm{div} \vec{D} = \rho \qquad \mathrm{rot} \vec{H} - \frac{\partial\vec{D}}{\partial t} = \vec{i}
ただし
\vec{D} = \epsilon_0\vec{E}, \qquad \vec{B} = \mu_0\vec{H}

それぞれの記号の意味は下記の通り。

\epsilon_0 誘電率
\mu_0 透磁率
\vec{E} 電場の強さ
\vec{D} 電束密度
\vec{B} 磁束密度
\vec{H} 磁場の強さ
\vec{i} 電流密度

ちなみにローレンツ力は次式。
\vec{F} = e\left(\vec{E} + \vec{v}\times\vec{B}\right)
ただし e は荷電粒子の電荷、 \vec{v} は荷電粒子の運動速度。

さて、美しいと言われているマックスウェル方程式ですが。美しいような、そうでもないような、と言いますか。誘電率やら透磁率やら色々出てきて割とややこしいので、もうちっと整理できんのかなぁ、というのが式変形の動機。

\vec{E}, \vec{B} を消去

\vec{E}=\vec{D}/\epsilon_0, \vec{B}=\mu_0\vec{H}, \vec{i}=\rho\vec{v} を用いてマックスウェル方程式を式変形する。

\mathrm{div} \vec{H} = 0 \qquad \mathrm{rot} \vec{D} + \frac{1}{c^2}\frac{\partial\vec{H}}{\partial t} = \vec{0}
\mathrm{div} \vec{D} = \rho \qquad \mathrm{rot} \vec{H} - \frac{\partial\vec{D}}{\partial t} = \rho\vec{v}

ただし cは光速。

マックスウェル方程式から誘電率や透磁率が消えて、ちょっとスッキリ。代わりに光速 c が増えましたが、こちらのほうが割と本質的な物理定数だと思うので、こっちの表記のほうがいいんじゃないかな、と。

ミンコフスキー空間を意識して式変形

次の記号を導入します。
\tau = ct, \quad \vec{\tilde{H}}=\frac{\vec{H}}{c}, \quad \vec{\tilde{v}}=\frac{\vec{v}}{c}
するとマックスウェル方程式は次のように変形できます。

\mathrm{div} \vec{\tilde{H}} = 0 \qquad \mathrm{rot} \vec{D} + \frac{\partial\vec{\tilde{H}}}{\partial\tau} = \vec{0}
\mathrm{div} \vec{D} = \rho \qquad \mathrm{rot} \vec{\tilde{H}} - \frac{\partial\vec{D}}{\partial\tau} = \rho\vec{\tilde{v}}

おー。光速 c も式から消えました。電場と磁場の関係式の対称性がよりくっきりと見えるようになったと感じます。
さらに、次元解析すると単位もシンプルになっていることが分かります。

  • 電場 \vec{D} と磁場 \vec{\tilde{H}} の単位は共に [C/m^2]
  • 4つの微分方程式の単位はすべて [C/m^3]
  • 時間軸が \tau = ct により長さ次元に置き換わっているため、式に時間の次元が出てこない。

微分形式の導入

ベクトル解析は微分形式を用いるともっとシンプルに書けます。というわけで微分形式を導入。電場と磁場を次のように微分形式で表します。

D = D_x dy\wedge dz + D_y dz \wedge dx + D_z dx \wedge dy
H = \left(\tilde{H_x} dx + \tilde{H_y} dy + \tilde{H_z} dz\right)\wedge d\tau

ここで表記を簡潔にするために次の記号を導入します。

\vec{\sigma^1}=\left(\begin{array}{c} dx \\ dy \\ dz \end{array}\right), \quad \vec{\sigma^2}=\left(\begin{array}{c} dy \wedge dz \\ dz \wedge dx \\ dx \wedge dy \end{array}\right), \quad \sigma^3 = dx \wedge dy \wedge dz

するとD, Hは次のように書けます。

D = \vec{D}\cdot\vec{\sigma^2}, \quad H = \vec{\tilde{H}}\cdot\vec{\sigma^1}\wedge d\tau

これらにホッジ作用素を適用すると次のようになります。

*D = -\vec{D}\cdot\vec{\sigma^1}\wedge d\tau, \quad *H = \vec{\tilde{H}}\cdot\vec{\sigma^2}

*D で符号がマイナスになるのは、ミンコフスキー空間の計量が時間軸で -1 になるため。（たぶんこれでいいと思うんだけど、あんまり自信はない）

微分形式によるマックスウェル方程式

天下り的ですが、次の微分形式を作ります。

F = -*D + *H = \vec{D}\cdot\vec{\sigma^1}\wedge d\tau + \vec{\tilde{H}}\cdot\vec{\sigma^2}

これを全微分してみます。

\begin{array}{ccl} dF &=& \mathrm{rot}\vec{D}\cdot\vec{\sigma^2}\wedge d\tau + \left( \mathrm{div}\vec{\tilde{H}}\sigma^3 + \frac{\partial\vec{\tilde{H}}}{\partial \tau}\cdot\vec{\sigma^2}\wedge d\tau \right) \\ &=& \left(\mathrm{rot}\vec{D} + \frac{\partial\vec{\tilde{H}}}{\partial\tau}\right)\cdot\vec{\sigma^2}\wedge d\tau + \mathrm{div}\vec{\tilde{H}}\sigma^3 \end{array}

以上から、マックスウェル方程式の4つの式のうちの2つが、次のシンプルな1本の式で表せてしまうことが分かります。

dF = 0

残りの二つの式も微分形式で表せると素敵ですね。そのためにFにホッジ作用素を適用してみます。

*F = D - H = \vec{D}\cdot\vec{\sigma^2} + \vec{\tilde{H}}\cdot\vec{\sigma^1}\wedge d\tau

これを全微分してみます。

\begin{array}{ccl} d*F &=& \left( \mathrm{div}\vec{D}\sigma^3 + \frac{\partial\vec{D}}{\partial \tau}\cdot\vec{\sigma^2}\wedge d\tau \right) - \mathrm{rot}\vec{\tilde{H}}\cdot\vec{\sigma^2}\wedge d\tau \\ &=& \mathrm{div}\vec{D}\sigma^3 - \left(\mathrm{rot}\vec{\tilde{H}} - \frac{\partial\vec{D}}{\partial \tau} \right)\cdot\vec{\sigma^2}\wedge d\tau \end{array}

これをマックスウェル方程式と比べると、残りの2式は次のように表せることが分かります。

d*F = \rho\left(\sigma^3 - \vec{\tilde{v}}\cdot\vec{\sigma^2}\wedge d\tau\right)

以上でマックスウェル方程式を微分形式でシンプルに書き表すことに成功しました。

電磁場ポテンシャル

ここで「ポアンカレの補題（の逆）」ってヤツを使います。

ポアンカレの補題(の逆)：
k次微分形式ωが
d\omega = 0
を満たすとする（このようなωを「閉形式」という）。このとき、
\omega = d\eta
を満たすk-1次形式ηが存在する（このようなωを「完全形式」という）。

つまり「閉形式ならば完全形式である」というのがこの定理の主張。
えっと、うろ覚えなので怪しいですが、確か「完全形式ならば閉形式である」が「ポアンカレの補題」で、その逆の「閉形式ならば完全形式である」が「ポアンカレの補題の逆」だったはず。で、前者は常に真だが後者（逆）は必ずしも真ではなかったはず。ですが、逆が成立するには何か条件があって、ユークリッド空間ではその条件を満たしているから逆も成立するんじゃなかったかな。でもここではミンコフスキー空間を扱ってますからユークリッド空間じゃないですね。でもミンコフスキー空間でもポアンカレの補題の逆は成り立つってことなんでしょう、たぶん。あー、理解が怪しい…ゴニョゴニョ。

こまけーことはいいんだよ！（逆切れ）ということで、この定理を適用します。先ほどマックスウェル方程式では次式が成立すると書きました。

dF = 0

ということは、この定理によれば F は次のように書けるってことです。

F = d\Phi, \quad \Phi は1次微分形式

そこで \Phi を次のように置きましょう。

\begin{array}{ccl} \Phi &=& A_xdx + A_ydy + A_zdz - \phi d\tau \\ &=& \vec{A}\cdot\vec{\sigma^1} - \phi d\tau \end{array}

これを電磁場ポテンシャルと呼びます。
これをマックスウェル方程式の次の式に代入してみます。

d*F = \rho\left(\sigma^3 - \vec{\tilde{v}}\cdot\vec{\sigma^2}\wedge d\tau\right)

すると次式が得られます。

d*d\Phi = \rho\left(\sigma^3 - \vec{\tilde{v}}\cdot\vec{\sigma^2}\wedge d\tau\right)

おお、なんということでしょう、マックスウェル方程式が一つの式で表せてしまいました。

※ 計算間違いに気づいた方は教えてください。あんまり自信ないです。特に符号とか怪しい。

継続は力なり #FsAdvent

この記事は F# Advent Calendar 2014 - connpass の18日目の記事です。昨日は teramonagi さんの FsLab JournalでReproducible Research(レポート)を簡単に作りたい - My Life as a Mock Quant でした。

分をわきまえずに継続モナドとかいうヤツについて書きます。
ちょっと前のことですが、継続（continuation）の概念がストンと腑に落ちる瞬間を味わいまして割と気持ち良かったので、その感動を忘れないうちに書きます。まあ、経験上最初の「分かった！」はぬか喜びであることが多いので今回も誤解している可能性はありそうですが、多少の間違いはあっても感動の熱さが残っているうちに書かないと書く気が失せてしまうこともやはり経験から分かっているので、勇気を出してこの機会に書いてみようと思うのです。

なお、私はこの程度のオッサンですので過度な期待は禁物ですよ。

  • Haskell 分かりません。
  • 圏論知りません。
  • モナド則は見たことあるけど「ふーん？」という感じ。
  • 継続モナドで調べてると call/cc というのが出てくるのですが、これはまだ理解できてません。
  • 限定継続という言葉を見た気がしますが、何のことやらさっぱり。

書いたコードは gist にアップしてあります。
https://gist.github.com/u1roh/5e12b8819a5aa0e5eca5

課題として直線作図機能を考えてみよう

私はCAD系ツールの開発がお仕事なので、作図ツールの設計には関心があります。ですのでここでは、画面上の2点をクリックするとその2点を結ぶ直線を作図する機能を課題として考えていくことにします。

アプリケーション外観

具体的なイメージを持ってもらうため、サンプルとして作ったアプリケーションの外観を見せておきます。何やら怪しげなメニューが見えていますが今は気にしない気にしない。
まず直線作図メニューを選んで、次に画面上を2点クリックすると直線が作図されます。

準備

簡単な Forms アプリケーションを作ることにしましょう。準備として次のようなコードを用意しておきます。

[<EntryPoint>]
let main argv = 
  let form =
    let form = new Form ()

    let drawObjs = List<Graphics -> unit> ()
    form.Paint |> Event.add (fun e -> for draw in drawObjs do draw e.Graphics)

  Application.Run form
  0

drawObjs に描画処理を登録しておくと Paint イベントで描画が実行されるようになっています。drawObjs に直線を作図する関数を Add すれば直線が作図できるというわけです。

コンソールプログラムだったら…

突然ですが、コンソールプログラムでユーザー入力を受け付けるのは簡単でした。次のプログラムはユーザー入力を2行受け付けて、入力文字列を連結して出力します。

let line1 = Console.ReadLine ()
let line2 = Console.ReadLine ()
Console.WriteLine (line1 + line2)

これと似たような感じでマウスクリック入力を2回受け付けることができたらいいですね！

しかしマウスイベントの場合は…

マウスイベントはコンソールの入力のようには扱えません。イベントハンドラを登録しておいてクリックイベントが飛んでくるのを待ち構えるしかありません。ですので直線作図機能を作るためには例えば次のようなコードを書くことになります。

let firstPos = ref None
let subscription = ref Option<IDisposable>.None
subscription :=
  form.MouseClick |> Observable.subscribe (fun e ->
    match !firstPos with
    | None -> firstPos := Some e.Location
    | Some p ->
      drawObjs.Add (fun g -> g.DrawLine (Pens.Blue, p, e.Location))
      form.Invalidate ()
      !subscription |> Option.iter (fun d -> d.Dispose ())
      subscription := None
      firstPos := None
  ) |> Some)

飛んできたクリックイベントが1点目のクリックなのか2点目のクリックなのかを判定する必要があります。また、1点目のクリック位置は一旦どこかに保持しておかなくてはなりません。結果として出来上がったプログラムは、「2回クリックを受け付けて直線を作図する」というシンプルな脳内モデルとは似ても似つかない構造のプログラムとなってしまいます。

継続モナドを使うと！

結論から書きましょう。継続モナドを使うとこう書けるんです！

cont {
  let! click1 = form.MouseClick |> Cont.ofObservable
  let! click2 = form.MouseClick |> Cont.ofObservable
  drawObjs.Add (fun g -> g.DrawLine (Pens.Black, click1.Location, click2.Location))
  form.Invalidate ()
} |> Cont.run
  • コンソールプログラムと同じような流れでクリックイベントを取得しているところにご注目！
  • cont { ... } は継続モナドのためのコンピュテーション式です
  • Cont は継続モナドのための関数が定義されているモジュールです
  • コードには見えていませんが、Cont<'a> という継続モナドの型を定義して利用しています

少し不思議な印象を受けるかもしれません。このコードを実行すると、クリック待ち状態でいったんイベントループに処理が戻り、クリック入力を受け取ると再びコードに戻ってきてさっきの続きから処理が再開（つまり処理が継続！）されます。この不思議さは seq { yield ... } に似ていると思いますし、裏を返せば yield と同程度の不思議さでしかないとも言えます。
このコードを実現している仕組みを以下でひも解いていきましょう。

継続って何？

さて、いったん上記の直線作図の課題は忘れて、基礎的なところからコツコツと継続について説明していくことにします。例として次のような単純な処理を考えます。
f:id:u_1roh:20141217230441p:plain

  • 関数 f1 は a0 を受け取って a1 を返します
  • 関数 f2 は a1 を受け取って a2 を返します
  • 関数 f3 は a2 を受け取って a3 を返します

「普通」に関数の戻り値を使って実装すると…

let a1 = f1 a0
let a2 = f2 a1
let a3 = f3 a2

(* パイプライン演算子を使って次のように書くこともできる *)
a0 |> f1 |> f2 |> f3

何の変哲もない、普通の処理ですね。このコードは明らかに一気に最後まで実行されますから、上で見たような「処理がいったん中断して途中から再開する」ようなトリックが入り込む余地はどこにもありません。一応これらの関数の型を見ておきましょう。

val f1 : A0 -> A1
val f2 : A1 -> A2
val f3 : A2 -> A3
※ a0, a1, a2, a3 の型をそれぞれ A0, A1, A2, A3 としました

継続渡し（CPS; Continuation Passing Style) について

上記コードを「継続渡し」というスタイル（＝CPS）に書き換えて行きたいのですが、その前にここでCPSとは何か説明しておきたいと思います。
CPSでは関数の結果を戻り値として受け取るのではなく、コールバック関数に結果を入れるスタイルとなります。…と言われても分からないと思いますから、まずは f1 だけを例にとってCPSに書き換えてみましょう。
まず f1 を次のような形式に書き換えます。

(* CPS で呼び出せる f1。型は A0 -> (A1 -> unit) -> unit *)
let f1 (a0 : A0) (continuation : A1 -> unit) =
  ...

そして次のようにCPSで関数呼び出しをすることになります。

f1 a0 (fun a1 -> ...)

要点を箇条書きしてみます。

  • f1 の戻り値は unit になっています。つまり結果を戻り値として返していません。
  • f1 の引数が1つ増えていて、2つ目の引数にコールバック関数を渡すようになっています。
  • コールバック関数のシグネチャが A1 -> unit であることから分かるように、f1 の結果をコールバック関数で受け取るようになっています。

このコールバック関数のことを「継続」と呼び、このような継続を引数に渡すスタイルを「継続渡し」と呼ぶのです。
私見になりますが、「継続」などといういかめしい用語に惑わされている人って結構いるのではないでしょうかね。というか私は惑わされていましたね、とほほ。代わりに「続き」とか「続き渡し」と呼んでみれば一気に理解のハードルが下がる気がします。単に計算の「続き」をコールバック関数として渡しているだけです。

ではCPSは戻り値を使う方法と比べて何が優れているのでしょうか。何が本質的な違いなのでしょうか。
戻り値を使う方法は、f1 を呼び出すことで結果 a1 を pull するスタイルとなります。しかし CPS では f1 に「続き」を渡して結果を push してもらうスタイルとなります。つまり、CPSでは主導権を呼び出し先に委ねてしまっているのですね。ですから「続き」を呼び出すタイミングは呼び出し先である f1 が自由に決められるのです。いったんメッセージループに戻ってクリックイベントを検知してから「続き」を呼び出してもいいし、非同期計算の終了後に「続き」を呼び出してもいいのです。ここに継続モナドのトリックの種明かしがあるのです。

余談ですが、「続きはウェブで」っていうTVCMありますよね。あれも一種の継続渡しと言えるかもしれません。CMの最後に表示される検索キーワードが「継続」なんですね。CMの作り手側の視点に立つと分かります。TVCMは時間の制約が厳しく、TV放送という「スレッド」を長時間占有することは許されませんが、視聴者にはより長いCMを見てほしいので最後に継続（＝続き）を渡してくるわけです。もちろん、継続を受け取った視聴者のその後の行動の主導権は視聴者自身にあります。すぐにスマホで検索してもよし、後でPCでググってもよし、もちろん無視してもよし（大半は無視でしょうね）。

CPSに書き換えてみる

さて、f1, f2, f3 をすべてCPS用に書き換えるとこうなります。

let f1 (a0 : A0) (continuation : A1 -> unit) = ...
let f2 (a1 : A1) (continuation : A2 -> unit) = ...
let f3 (a2 : A2) (continuation : A3 -> unit) = ...

関数の型は下記の通りです。

val f1 : A0 -> (A1 -> unit) -> unit
val f2 : A1 -> (A2 -> unit) -> unit
val f3 : A2 -> (A3 -> unit) -> unit

これら f1, f2, f3 をCPSで呼び出すコードは下記のようになります。

f1 a0 (fun a1 ->
f2 a1 (fun a2 ->
f3 a2 (fun a3 -> ...)))

「継続（＝続き）」として渡しているコールバック関数が入れ子になって呼び出されている構造に注目してください。なお、入れ子になっていることをはっきりさせるためにはインデントを入れるべきですが、ここでは敢えてインデントを入れませんでした。その理由は f1, f2, f3 があたかも対等な立場でフラットに連続して処理されているように見せかけるためですが、好みに合わせて適宜インデントを入れて読んでください。
さて、次はこの入れ子構造を少し考察してみましょう。

ラッキョウとリストと継続

少し話が飛びますが、ここでリストのデータ構造を復習してみます。リスト型は次のように自己再帰的に定義されています。

type List<'a> =
  | Empty
  | List of Head:'a * Tail:List<'a>

Head が先頭要素で、Tail が「残り」です。Head を剥がしても剥がしても中からリストが出てくる構造です。つまり

1 :: 2 :: 3 :: 4 :: []

というリストは、各要素がフラットに並んでいるのではなく次のような自己再帰的な入れ子構造になっているわけです。

1 :: (2 :: (3 :: (4 :: [])))

まるでラッキョウのような構造ですね。ラッキョウは皮を剥いても剥いても中からラッキョウが出てくる構造をしています（実際に確かめたことはありませんが）。ラッキョウ型は次のように定義できるでしょうか。

type ラッキョウ =
  | 空
  | ラッキョウ of 皮 * ラッキョウ

そして「継続」にも、このようなラッキョウ型の自己再帰的な構造が見られるように思うのです。割と一般的な見方かと思いますが、プログラムというものを何らかの処理命令を連続して次々に処理していくものと捉えてることができます。つまり処理命令がフラットにシーケンシャルに並んでいる構造です（ループとか条件分岐は考えないとして）。しかし「継続」を理解するためには、プログラムを次のような自己再帰的な構造に捉える視点が必要なのではないか、と思うに至りました。

type プログラム =
  | 空
  | プログラム of 処理命令 * 継続（＝続き） : プログラム

つまり命令を処理しても処理しても中から続きのプログラムが出てくるという自己再帰的な構造です。ラッキョウのようにプログラムが入れ子になっていると捉えるわけです。前節に示した継続渡しのコードには、そういう入れ子構造が現れているのです。
次はいよいよ継続モナドの導出です。CPSのコードを蒸留して純粋な継続渡しの構造だけを抽出し、モナドという構造に昇華する作業です。

継続モナドの導出

さて、前節で書いた継続渡しのコードを書き換えて継続モナドを導出してみたいと思います。まず関数の型を再掲します。

val f1 : A0 -> (A1 -> unit) -> unit
val f2 : A1 -> (A2 -> unit) -> unit
val f3 : A2 -> (A3 -> unit) -> unit

型名にエイリアスをつけて分かりやすく

継続（＝続き）として渡すコールバック関数の型に別名をつけます。

type Callback<'a> = 'a -> unit

するとf1,f2,f3の型は次のように読み替えることができます。

val f1 : A0 -> Callback<A1> -> unit
val f2 : A1 -> Callback<A2> -> unit
val f3 : A2 -> Callback<A3> -> unit

さらに次の型を定義します。（これが継続モナド型！）

type Cont<'a> = Callback<'a> -> unit

するとf1,f2,f3の型は更に次のように読み替えることが出来ます。

val f1 : A0 -> Cont<A1>
val f2 : A1 -> Cont<A2>
val f3 : A2 -> Cont<A3>

ずいぶんスッキリして構造が分かりやすくなりました。重要なポイントは次の2つです。

  • Cont<'a> 型として「継続を受け取る構造」だけが抽出されていること
  • f1,f2,f3 が、元は戻り値なしの関数（unitを返す関数）であったにも関わらず、Cont<'a> という「値」を返す関数であるかのように書き換わっていること

合わせて f1, f2, f3 の定義も次のように書き換えておくと、「Cont<'a> という値を返している」ことがハッキリしますね。

let f1 (a0 : A0) : Cont<A1> = fun (continuation : Callback<A1>) -> ...
let f2 (a1 : A1) : Cont<A2> = fun (continuation : Callback<A2>) -> ...
let f3 (a2 : A2) : Cont<A3> = fun (continuation : Callback<A3>) -> ...

bind関数の導出

f1, f2, f3 を呼び出しているコードを再掲します。

f1 a0 (fun a1 ->
f2 a1 (fun a2 ->
f3 a2 (fun a3 -> ...)))

f1 の処理を f2 の処理を「くっつける」ことを考えていきます。つまり、f1 と f2 の処理をまとめて行う関数 f12 : A0 -> Cont を定義したいのです。

(* f1 と f2 をくっつけた関数 *)
let f12 (a0 : A0) : Cont<A2> =
  fun (continuation : Callback<A2>) ->
    f1 a0 (fun a1 -> f2 a1 continuation)

(* f12 を使うと呼び出し側のコードはこうなる *)
f12 a0 (fun a2 ->
f3  a2 (fun a3 -> ...))

さらに f12 と f3 もくっつけることができます。

(* f12 と f3 をくっつけた関数 *)
let f123 (a0 : A0) : Cont<A3> =
  fun (continuation : Callback<A3>) ->
    f12 a0 (fun a2 -> f3 a2 continuation)

(* f123 を使うと呼び出し側のコードはこうなる *)
f123 a0 (fun a3 -> ...)

このように次々と関数呼び出しをくっつけることができます。この「くっつける」操作だけを bind 関数として定義してみましょう。

(* 「くっつける」操作だけを取り出した関数 *)
let bind (f : 'a1 -> Cont<'a2>) (x : Cont<'a1>) : Cont<'a2> =
  fun (continuation : Callback<'a2>) ->
    x (fun a1 -> f a1 continuation) 

(* bind を使うと呼び出し側のコードはこうなる *)
let a1 = f1 a0
let a2 = bind f2 a1
let a3 = bind f3 a2
a3 (fun a3 -> ...)

(* パイプライン演算子を使うとこうなる *)
(f1 a0 |> bind f2 |> bind f3) (fun a3 -> ...)

おおお！？なんだか「らしく」なってきたと思いませんか？
ここで「継続って何？」の冒頭に戻って、普通の戻り値による関数の呼び出しとソックリであることを確認してください。

コードの整理

ここで一旦、コードを整理しなおしておきます。

(* ↓単純な型エイリアスはやめて、要素が1つだけの判別共用体として定義しなおした *)
type Callback<'a> = Callback  of ('a -> unit)
type Cont<'a>     = Cont      of (Callback<'a> -> unit)

(* Cont モナド用の関数をモジュールに整理 *)
module Cont =

  let inline bind binder (Cont invoke) =
    Cont <| fun callback ->
      invoke <| Callback (fun a -> let (Cont invoke) = binder a in invoke callback)

  (* Cont モナドを実行する *)
  let inline run (Cont invoke) = invoke (Callback (fun () -> ()))

  (* 単位元 *)
  let inline ofValue x =
    Cont <| fun (Callback callback) -> callback x

継続モナドを使って直線作図機能を作ろう

ようやく私たちは継続モナドを手に入れました（← ホントはモナド則を満たしているべきかきちんと確認するべきなのでしょうが私自身の理解があやふやなのでパス m(_ _)m）。これを使って、いよいよ本記事の冒頭の課題「画面上の2点をクリックするとその2点を結ぶ直線を作図する機能」に挑戦していきます。

ofObservable 関数

課題を解くには、クリックイベントを継続モナドとして扱えるようにしなくてはなりません。F# ではイベントは IObservable 型として扱うことができます。ですので、IObservable 型を「イベントを検知すると継続を呼び出す継続モナド」に変換する関数が定義できれば汎用性がありそうです。私は次のように定義してみました。

module Cont =
  let ofObservable o =
    Cont <| fun (Callback callback) ->
      let subscr  = ref Option<IDisposable>.None
      subscr := o |> Observable.subscribe (fun a ->
        callback a
        !subscr |> Option.iter (fun d -> d.Dispose (); subscr := None)) |> Some

ちょっと分かりにくいので簡単に解説します。
IObservable はイベントがストリームとして次々と流れてくるイメージで、そのイベントの「流れ」を subscribe することができるというものです。これに対して今回必要なのは、クリックイベントを1回だけ検知して継続に渡して終了するというモノですから、イベントを検知した時点でイベントの購読を終了させなくてはなりません。そのため subscribe に渡しているイベントハンドラの中で自身をイベントから外す処理を行っています。
これを使うことで、次のようにして MouseClick イベントを Cont に変換することができます。

let clickCont = Cont.ofObservable form.MouseClick

直線作図機能

準備が整いましたので直線作図機能を書いてみます。

Cont.ofObservable form.MouseClick |> Cont.bind (fun click1 ->
Cont.ofObservable form.MouseClick |> Cont.bind (fun click2 ->
  drawObjs.Add (fun g -> g.DrawLine (Pens.Black, click1.Location, click2.Location))
  form.Invalidate ()))
|> Cont.run

素晴らしい！

コンピュテーション式

モナドができたらコンピュテーション式も定義しておくと素敵になります。モナドがまるで言語自体に組み込まれたような感じで利用することができます。（コンピュテーション式については既に様々な方が素晴らしい説明を書いてくださっているので説明は割愛します。）
次のようなビルダーを定義しておくと、

type ContBuilder () =
  member __.Bind (x, f)  = x |> Cont.bind f
  member __.Return x  = Cont.ofValue x
  member __.Zero ()   = Cont.ofValue ()

let cont = ContBuilder ()

直線作図機能は次のようになります。（本記事の冒頭のコードの再掲）

cont {
  let! click1 = form.MouseClick |> Cont.ofObservable
  let! click2 = form.MouseClick |> Cont.ofObservable
  drawObjs.Add (fun g -> g.DrawLine (Pens.Black, click1.Location, click2.Location))
  form.Invalidate ()
} |> Cont.run

ますますもって素晴らしい！
これにてようやく冒頭のコードの種明かしが終了したわけですが、ご理解いただけましたでしょうか。とかく複雑で訳が分からなくなりがちなモナドですが、私なりにできる限り分かりやすい説明を試みたつもりです。私が感じた「腑に落ちた時の気持ちよさ」を少しでもお裾分け出来たなら嬉しいです。

（ここで記事としては一旦終わり。あとはゆるゆるダラダラと、応用とか、拡張とか）

応用

直線作図機能以外にも応用してみましょう。

モードレスダイアログ

モーダルダイアログはとても簡単ですが、モードレスダイアログになると途端に面倒になる気がしませんか。モーダルダイアログのような気軽さでモードレスダイアログが利用できたら良いと思いませんか。
下記はダイアログをモードレスで表示して、OKボタンが押されるとダイアログを閉じて MessageBox を表示するプログラムです。

cont {
  let dialog = new Form (Text = "dialog")
  let! _ =
    let btnOK = new Button (Text = "OK")
    dialog.Controls.Add btnOK
    dialog.Show ()
    btnOK.Click |> Cont.ofObservable
  dialog.Close ()
  MessageBox.Show "OK" |> ignore
  return ()
} |> Cont.run)

非同期計算

何らかの時間のかかる処理をバックグラウンドで実行し、終了したらUIスレッドに戻って続きを実行する、というのはよくあるパターンですね。もちろん async だけで十分簡単に出来るのですが、ここではあえて継続モナド化してみます。
まず準備として下記の ofAsync という関数を定義します。これは Async<'a> 型の計算オブジェクトを Cont<'a> に変換する関数です。

module Cont =
  ...
  let ofAsync (cts : CancellationTokenSource) computation =
    Cont <| fun (Callback callback) ->
      let con = SynchronizationContext.Current
      let cts = if cts = null then new CancellationTokenSource () else cts
      Async.StartWithContinuations (
        async {
          do! Async.SwitchToThreadPool ()
          return! computation cts.Token
        },
        (fun a -> con.Post ((fun _ -> callback a), null)), ignore, ignore, cts.Token)
  • 非同期計算をキャンセルすることも考えて CancellationToken が登場させましたが、結局このサンプルではキャンセル処理は扱っていないので例として少し良くなかったかもしれません。CancellationToken 周りは無視して読んでください。
  • 大事な点として、継続の callback がUIスレッドで実行されるように　SynchronizationContext に Post しています。

これを使って、1秒間のダミー処理をバックグラウンドで動かしてみます。プログレスバーの更新も行っています。

cont {
  let ui = SynchronizationContext.Current
  let! result = Cont.ofAsync null (fun token ->
    async {
      for i = 1 to 10 do
        Thread.Sleep 100
        printf "."
        ui.Post ((fun _ -> progressBar.Value <- progressBar.Maximum * i / 10), null)
      return "Async calculation result"
    })
  MessageBox.Show result |> ignore
  progressBar.Value <- 0
} |> Cont.run)

全部ちゃんぽん

全然意味ない処理ですが、全部ちゃんぽん出来ちゃいます。

cont {
  // クリックイベントを取得
  let! click = form.MouseClick |> Cont.ofObservable
  printfn "clicked: Location = %A" click.Location

  // モードレスダイアログを表示
  let dialog = new Form (Text = "dialog")
  let! _ =
    let btnOK = new Button (Text = "OK")
    dialog.Controls.Add btnOK
    dialog.Show ()
    btnOK.Click |> Cont.ofObservable
  dialog.Close ()

  // 非同期計算
  let ui = SynchronizationContext.Current
  do! Cont.ofAsync null (fun _ ->
    async {
      for i = 1 to 10 do
        Thread.Sleep 100
        ui.Post ((fun _ -> progressBar.Value <- progressBar.Maximum * i / 10), null)
    })

  // 最後にメッセージを表示
  MessageBox.Show "done" |> ignore
  progressBar.Value <- 0
} |> Cont.run)

コンピュテーション式の拡張

ここからは全くもって理解不足で、今後の私の課題になるのですが…。たぶん、コンピュテーション式のビルダーに Combine とか Using とか For とかも定義してあげたほうが色々と便利そうです。しかしまだコンピュテーション式の変換ルールが理解し切れておらず…という状態です。トライしてみたところまで書きます。

評価の遅延

例えば次のように書いたとき、

cont { printfn "hoge" }

これがすぐ評価されて "hoge" がプリントされてしまうよりは、Cont.run の呼び出しまで遅延されるほうが使いやすそうです。そのためにビルダーに Delay と Run を定義してみました。

    // Delay と Run はコンピュテーション式の評価を遅延させるための仕掛け
    member __.Delay f = f
    member __.Run   f = ofValue () |> bind f

Combine

これでいいのかな？

    // Combine は例えば if による条件分岐の後に続けて処理を書きたいときに必要になる
    // - x の後に「継続」して f を実行する。その際、x の結果の値は無視する。
    // - 第二引数の f が関数なのは Delay によって処理が遅延されているから（たぶん…）
    member __.Combine (x, f) = x |> bind (fun _ -> f ())

Using

概ねこんな感じだと思うのですが、例外が起きた場合に対応できてないです。

    // use 構文が使えるようにする。
    // f x で得られる計算処理の終了時に x.Dispose () を呼べば良いはず。
    // しかし例外が起きた時には対応できていない。
    member __.Using (x: #IDisposable, f : _ -> Cont<_>) =
      Cont <| fun (Callback callback) ->
        let (Cont invoke) = f x
        invoke <| Callback (fun a -> x.Dispose(); callback a)

というか、Using に限らず今回のサンプルコード全体にわたって例外のことは考慮できていないですね。例外に対応するには継続のコールバック関数が単純な関数ではダメで、 IObserver みたいに OnError メソッドがないとダメかな？

池袋物理学勉強会 演習問題 2-4 #ikeph

明日の池袋物理学勉強会では演習問題2-4が私の担当となっていますが、残念ながら都合により参加できませんので演習問題の解答のみブログで参加します。

演習問題2
4. 剛体の運動方程式を、変分原理から導け。

ああ、なんてシンプルな問題文なのでしょう。こういう時はいきなり巻末の解答を見ます（え？

演習問題略解
重心の座標を x_G、慣性テンソルを
{\displaystyle I_G=\int d^3 x \left(r^2\delta_{ij} - x_ix_j\right)\rho }
とすると、
{\displaystyle T=\frac{1}{2}M \dot{x}_G^2+\frac{1}{2}\sum_{i,j}I_{i,j}\omega_i\omega_j }
ただし、 \omega は重心の周りの角速度。外力がなければ, L = T である。

これは解答というよりヒントという感じですね。

問題のポイント

剛体について

質点ではなくて剛体なので、物体には大きさがあることを考慮しなくてはなりません。すなわち、物体の回転運動も考慮に入れる必要が生じます。ですから「略解」にも回転運動による運動エネルギーの項（慣性テンソルが含まれる項）が含まれています。

外力について

問題文には剛体に働く外力に関する記述が全くないのですが、ヒントに「外力がなければ L = T」と書いてあることですし、ここでは簡単のため外力はないものとして解きます。

導出までのストーリー

  1. 剛体の運動エネルギーの式を立てる（略解に記載済み）
  2. 運動エネルギーの作用積分の式を立てる
  3. 変分法により作用積分が最少となる条件を求める

運動方程式の導出

ここでは略解に記載されている運動エネルギーの式をそのまま用いて、変分法により運動方程式を導きます。
剛体の回転状態を回転ベクトル \vec{\theta}=(\theta_1, \theta_2, \theta_3) で表し、その時間微分が角速度ベクトルに一致するものとします。
{\displaystyle \dot{\vec{\theta}}=\vec{\omega} }
これを用いて、一般化座標 \vec{q} を次のように定めます。
{\displaystyle \vec{q}=\left(\vec{x_G}, \vec{\theta}\right)=\left(x_1, x_2, x_3, \theta_1, \theta_2, \theta_3\right) }
すると作用積分は次式となります。
{\displaystyle I=\int T(\vec{q}, \dot{\vec{q}}) dt = \int T(\vec{x_G}, \vec{\theta}, \dot{\vec{x_G}}, \dot{\vec{\theta}}) dt }
「ヒント」を見ると運動エネルギーには \vec{x_G} や \vec{\theta} の項は含まれていませんので、作用積分の式は次のように簡略化できます。
{\displaystyle I=\int T(\dot{\vec{q}}) dt = \int T(\dot{\vec{x_G}}, \dot{\vec{\theta}}) dt }
ここで次のような摂動を考え
{\displaystyle x_G \rightarrow x_G + \delta x }
{\displaystyle \theta \rightarrow \theta + \delta\theta }
作用積分の変分を計算します。
{\displaystyle \delta I = \int T(\dot{x_G} + \delta\dot{x}, \dot{\theta} + \delta\dot{\theta}) - T(\dot{x_G}, \dot{\theta}) dt = \int \left(\frac{\partial T}{\partial \dot{x_G}}\cdot\delta\dot{x} + \frac{\partial T}{\partial \dot{\theta}}\cdot\delta\dot{\theta}\right) dt }
運動エネルギーの式を代入すると、
{\displaystyle \delta I = \int \left(M\dot{x_G}\cdot\delta\dot{x_G} + I\dot{\theta}\cdot\delta\dot{\theta}\right)dt }
部分積分を用いると、
{\displaystyle \delta I = -\int \left(M\ddot{x_G}\cdot\delta x_G + I\ddot{\theta}\cdot\delta\theta\right)dt }
\delta I = 0 と置くことにより、運動方程式は次式となる。
{\displaystyle M\ddot{x_G}=0, I\ddot{\theta}=0 }

剛体の運動エネルギーについて

上では略解の運動エネルギーの式を所与のものとして使いましたが、ここでは運動エネルギーの式の導出を行ってみます。解析力学とはあまり関係はありませんが…。
剛体の運動エネルギー T は、並進運動による運動エネルギー T_1 と回転運動による運動エネルギー T_2 に分解できます。
{\displaystyle T = T_1 + T_2 }
並進運動のエネルギーは簡単です。
{\displaystyle T_1 = \frac{1}{2}M \dot{x}_G^2 }
ただし x_G は重心位置、 M は剛体全体の質量です。
問題は回転運動による運動エネルギーです。これを求めるために、まず角速度ベクトルについておさらいしておきます（←自分が忘れてた）。次のような角速度ベクトルを考えます。
{\displaystyle \vec{\omega}=\left(\begin{array}{c} \omega_x \\ \omega_y \\ \omega_z \end{array}\right) }
これが表す回転運動とはすなわち、 \vec{\omega} の方向を回転軸として、毎秒 ||\vec{\omega}|| = \sqrt{\omega_x^2+\omega_y^2+\omega_z^2} [rad] の速さで回転する運動です。これを用いて、重心位置からの相対座標 \vec{r} にある微小要素の運動エネルギーを求めます（下図参照）。
f:id:u_1roh:20140930211248p:plain
微小要素の速度 \vec{v} は次式で求まります。
{\displaystyle \vec{v}=\vec{\omega} \times \vec{r} = \left(\begin{array}{ccc} 0 & r_3 & -r_2 \\ -r_3 & 0 & r_1 \\ r_2 & -r_1 & 0 \end{array}\right)\left(\begin{array}{c}\omega_1 \\ \omega_2 \\ \omega_3 \end{array}\right) }
従って微小要素の運動エネルギーは
{\displaystyle dT_2 = \frac{1}{2}\rho\left||\vec{v}\right||^2 dxdydz = \frac{1}{2}\rho\vec{\omega}^T\left(\begin{array}{ccc}r_3^2+r_2^2 & -r_1r_2 & -r_1r_3 \\ -r_1r_2 & r_3^2+r_1^2 & -r_2r_3 \\ -r_1r_3 & -r_2r_3 & r_1^2+r_2^2 \end{array}\right)\vec{\omega} dxdydz }
{\displaystyle = \frac{1}{2}\rho\vec{\omega}^T\left\{(r_1^2+r_2^2+r_3^2)E - \left(\begin{array}{c}r_1 \\ r_2 \\ r_3\end{array}\right)\left(\begin{array}{ccc}r_1 & r_2 & r_3\end{array}\right)\right\}\vec{\omega} dxdydz }
（※ E は単位行列）
剛体全体の回転運動によるエネルギーは、上式を剛体全体にわたって積分すれば得られます。
{\displaystyle T_2 = \int dT_2 = \frac{1}{2}\vec{\omega}^T I \vec{\omega} }
ただし、慣性テンソル
{\displaystyle I = \int \rho\left\{(r_1^2+r_2^2+r_3^2)E - \left(\begin{array}{c}r_1 \\ r_2 \\ r_3\end{array}\right)\left(\begin{array}{ccc}r_1 & r_2 & r_3\end{array}\right)\right\} dxdydz }

感想

慣性テンソル懐かしい…。学生の時に工業力学で習ったなぁ…。

修正しました（2014-10-01 15:17）

@gm3d2 さんに幾つかミスをご指摘頂いて修正しました。ありがとうございました。

「池袋物理学勉強会(1)」の復習 #ikeph

池袋物理学勉強会（池袋物理学勉強会(1) - connpass）の復習です。第一回の開催からだいぶ日が経ったので今更感ありますが。この勉強会は講師の方が多大な労力を割いて下さっているので、心より感謝するとともに、ブログでも書いて微力ながら勉強会を盛り上げていこうなどと思った次第。

ラグランジュ方程式の \dot{q}の捉え方

下記にラグランジュ方程式を示します。
{\displaystyle \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{\vec{q}}}\right)-\frac{\partial L}{\partial \vec{q}}=\vec{0} }
この式の中に \dot{\vec{q}} による偏微分が現れています。この偏微分についてどう考えればよいのか、勉強会の最後に質問があったように記憶しています。私も疑問に思った箇所だったので質問者の方には共感しました。
疑問の内容はこうです。まず \dot{\vec{q}} というのは一般化座標 \vec{q} の時間微分です。つまり \dot{\vec{q}} は \vec{q} に従属する概念のはずです。しかしラグランジュ方程式では \dot{\vec{q}} があたかも独立変数であるかのように扱って偏微分を行っています。これはどのように捉えればよいのか、という疑問です。
この疑問については私なりの結論に至ったので、それが正しいかどうかを問う意味も込めてここで説明してみることにしましょう。
まず、一般化座標の自由度を f とすると一般化座標 \vec{q} は f 次元ベクトルです。
{\displaystyle \vec{q}=\left(q_1, q_2, \ldots, q_f\right)^T }
これを用いてラグランジアンは \vec{q} と \dot{\vec{q}} の関数として表記されています。
{\displaystyle L=L(\vec{q},\dot{\vec{q}}) }
この表記が既に変な感じなわけです。 \dot{\vec{q}}=d\vec{q}/dt という従属関係があるにも関わらず、明らかに \dot{\vec{q}} を独立変数として扱った表記になっています。独立なのか従属なのはハッキリしろと言いたくなるわけです。
私の結論としては、ここは独立変数として考えるべきです。ですから厳密には \dot{\vec{q}} と表記するのは正しくなくて、本来は独立した f 次元ベクトルの変数 \vec{r} を用いて
{\displaystyle L=L(\vec{q},{\vec{r}}) }
と表記するべきだと思うわけです。つまりラグランジアンLは、 \dot{\vec{q}} と \dot{\vec{r}} が張る 2f 次元空間に定義された実数値関数と捉えるべきではないでしょうか。そしてラグランジュ方程式は次のような2本の式の連立方程式として表記されるのが正しい姿ではないでしょうか。

  • {\displaystyle \frac{d}{dt}\left(\frac{\partial L}{\partial \vec{r}}\right)-\frac{\partial L}{\partial \vec{q}}=\vec{0}}
  • {\displaystyle \vec{r}=\frac{d\vec{q}}{dt}}

ラグランジュ方程式の考え方やメリット

これは私もまだ漠然としか分かりませんが、「こういうことかな？」と思ったことを書いておきます。
まず、ニュートンの運動方程式と比べて運動の捉え方が大きく変わります。例えば2粒子系の力学系を考えたとき、通常のニュートンの運動方程式では3次元空間上を飛び回る2つの粒子として運動を捉えることになります。しかしラグランジュ方程式では、2粒子系の自由度は6であるためラグランジアンLが定義される空間は12次元となり、この12次元空間中の一本の軌跡線として運動を捉えることになります。これはニュートンの運動方程式では決して得られなかった視点です。12次元空間ですから物理的イメージも幾何的イメージも湧きにくい抽象的な捉え方ですが、代わりに数学的には多粒子系でも一本の方程式で表せるわけですからメリットがありそうな気がします。
これは私の個人的な感覚ですが…、ニュートンの運動方程式は運動する粒子一つ一つに着目していますから、粒子の数が増えると着目する対象が沢山の点に分散してしまい系全体を俯瞰できないような気がします。木を見て森を見ず、という感じでしょうか。これに対してラグランジアンは多粒子系の系全体を捉えた物理量と言えますから、系全体としての特性を俯瞰して捉えようとする雰囲気を感じます。
教科書の先の方を読むと変分原理でラグランジュ方程式を捉え直すことが出来ます。つまり作用積分の最小化問題として運動方程式を捉え直すことが出来ます。コンピュータによる数値計算の手法として最小化問題を解くアプローチは色々ありますから、そういった手法を利用して運動方程式を解くことが出来るのかもしれません。つまり、運動方程式を別の形で定式化出来れば別の数値計算手法が利用可能になるメリットがある、といえるかもしれません。（この辺は完全に私の想像にすぎないので、大外ししているかもしれません）

ラグランジュ方程式の座標変換

自分で手を動かして計算してみたところ、ベクトルの表記を使うと教科書よりもシンプルに計算できるように思ったので、その計算をここ記しておきます。
教科書では次のような一般化座標の座標変換について書かれています。
{\displaystyle \vec{q}=\vec{q}(Q_1, Q_2, \ldots, Q_f)=\vec{q}(\vec{Q}) } ... 式(★)
しかしここでは、まずはより一般化した座標変換を考えて、後から上記の座標変換に戻ってくることにしましょう。つまり次のような一般的な座標変換を考えていくことにします。
{\displaystyle \vec{q}=\vec{q}(\vec{Q},\dot{\vec{Q}}) }
一般化座標の時間微分 \dot{\vec{q}} も同様に \vec{Q} と \dot{\vec{Q}} の関数となります。
{\displaystyle \dot{\vec{q}}=\dot{\vec{q}}(\vec{Q},\dot{\vec{Q}}) }
これら2つをまとめて微分すると次のようになります。
{\displaystyle \left(\begin{array}{c}d\vec{q} \\ d\dot{\vec{q}}\end{array}\right)= \left[\begin{array}{cc} \frac{\partial \vec{q}}{\partial \vec{Q}} & \frac{\partial \vec{q}}{\partial \dot{\vec{Q}}} \\ \frac{\partial \dot{\vec{q}}}{\partial \vec{Q}} & \frac{\partial \dot{\vec{q}}}{\partial \dot{\vec{Q}}} \\ \end{array}\right] \left(\begin{array}{c}d\vec{Q} \\ d\dot{\vec{Q}}\end{array}\right) =J\left(\begin{array}{c}d\vec{Q} \\ d\dot{\vec{Q}}\end{array}\right) }
ここで右辺に登場した座標変換行列を J とおきました。この行列は一般にヤコビ行列と呼ばれているものです。
このヤコビ行列を用いて、任意の実数値関数 f=f(\vec{q},\dot{\vec{q}})=f(\vec{Q},\dot{\vec{Q}}) の微分を考えていきましょう。次のような式変形が出来ます。
{\displaystyle df=\frac{\partial f}{\partial \vec{q}}\cdot d\vec{q}+\frac{\partial f}{\partial \dot{\vec{q}}}\cdot d\dot{\vec{q}} \\ =\left(\begin{array}{cc}d\vec{q}^T & d\dot{\vec{q}}^T\end{array}\right) \left(\begin{array}{c}\frac{\partial f}{\partial \vec{q}} \\ \frac{\partial f}{\partial \dot{\vec{q}}}\end{array}\right) =\left(\begin{array}{cc}d\vec{Q}^T & d\dot{\vec{Q}}^T\end{array}\right) J^T\left(\begin{array}{c}\frac{\partial f}{\partial \vec{q}} \\ \frac{\partial f}{\partial \dot{\vec{q}}}\end{array}\right) }
従って次の関係が得られます。
{\displaystyle \left(\begin{array}{c}\frac{\partial}{\partial \vec{Q}} \\ \frac{\partial}{\partial \dot{\vec{Q}}}\end{array}\right)=J^T\left(\begin{array}{c}\frac{\partial}{\partial \vec{q}} \\ \frac{\partial}{\partial \dot{\vec{q}}}\end{array}\right) }
さて、ここで最初の座標変換式(★)に戻ってみましょう。この変換式では\vec{q}は\vec{Q}のみの関数であり、\dot{\vec{Q}}の関数にはなっていません。従って次式が成り立ちます。
{\displaystyle \frac{\partial \vec{q}}{\partial \dot{\vec{Q}}}=0 }
また、この\vec{q}を時間微分して\dot{\vec{q}}を計算すると
{\displaystyle \dot{\vec{q}}=\left[\frac{\partial\vec{q}}{\partial\vec{Q}}\right]\dot{\vec{Q}} }
が得られ、この両辺を \dot{\vec{Q}}で微分することにより次式が得られます。
{\displaystyle \left[\frac{\partial\dot{\vec{q}}}{\partial\dot{\vec{Q}}}\right]=\left[\frac{\partial\vec{q}}{\partial\vec{Q}}\right] }
これらを用いるとヤコビ行列 J は次のように書き換えられます。
{\displaystyle J=\left[\begin{array}{cc} \frac{\partial \vec{q}}{\partial \vec{Q}} & 0 \\ \frac{\partial \dot{\vec{q}}}{\partial \vec{Q}} & \frac{\partial \vec{q}}{\partial \vec{Q}} \\ \end{array}\right] }
さあ、これらの結果を用いていよいよ、新しい座標系 \vec{Q}でもラグランジュ方程式が成立するかどうかを調べていくことにしましょう。準備として次のような作用素を定義しておきます。
{\displaystyle \Phi_q=\frac{d}{dt}\left(\frac{\partial}{\partial \dot{\vec{q}}}\right)-\frac{\partial}{\partial \vec{q}} }
これを用いるとラグランジュ方程式は次のように書けます。
{\displaystyle \Phi_qL=\vec{0} }
同様にして、新しい座標系 \vec{Q} では作用素は次のように定義されます。
{\displaystyle \Phi_Q=\frac{d}{dt}\left(\frac{\partial}{\partial \dot{\vec{Q}}}\right)-\frac{\partial}{\partial \vec{Q}} }
この作用素を、上で求めたヤコビ行列による変換式を用いて式変形していくことにしましょう。まずは第一項に注目すると、
{\displaystyle \frac{d}{dt}\left(\frac{\partial}{\partial \dot{\vec{Q}}}\right)= \frac{d}{dt}\left(\left[\frac{\partial \vec{q}}{\partial \vec{Q}}\right]^T\frac{\partial}{\partial \dot{\vec{q}}}\right)= \left[\frac{\partial \dot{\vec{q}}}{\partial \vec{Q}}\right]^T\frac{\partial}{\partial \dot{\vec{q}}}+ \left[\frac{\partial \vec{q}}{\partial \vec{Q}}\right]^T\frac{d}{dt}\left(\frac{\partial}{\partial \dot{\vec{q}}}\right) }
第二項は
{\displaystyle \frac{\partial}{\partial \vec{Q}}= \left[\frac{\partial\vec{q}}{\partial\vec{Q}}\right]^T\frac{\partial}{\partial\vec{q}}+ \left[\frac{\partial\dot{\vec{q}}}{\partial\vec{Q}}\right]^T\frac{\partial}{\partial\dot{\vec{q}}} }
従って作用素は次のようになります。
{\displaystyle \Phi_Q=\left[\frac{\partial\vec{q}}{\partial\vec{Q}}\right]^T\left(\frac{d}{dt}\left(\frac{\partial}{\partial\dot{\vec{q}}}\right)-\frac{\partial}{\partial\vec{q}}\right)= \left[\frac{\partial\vec{q}}{\partial\vec{Q}}\right]^T\Phi_q }
こうして \Phi_q と \Phi_Q の間にはとてもシンプルな変換式があることが分かりました。これによりラグランジュ方程式 \Phi_qL=\vec{0} が成り立てば新しい座標系でも方程式 \Phi_QL=\vec{0} が成り立つことが一目瞭然に分かります。

最後に

ブログで数式たくさん書くのツライ…(´・ω・｀)

続・Weak Event パターン

先日こんなエントリ（Weak Event パターン - かたちづくり）を書きましたが、Observable モジュールを拡張して弱参照でサブスクライブできるようにしておくと応用が効いて便利そうかな、などと思い作ってみました。

module Observable =

  type private WeakSubscription<'a> (callback : 'a -> unit, subscription : IDisposable) =
    member __.Callback = callback
    override __.Finalize () = subscription.Dispose ()
    interface IDisposable with member __.Dispose () = subscription.Dispose ()

  let subscribeWeakly (callback : 'a -> unit) source =
    let callbackRef = WeakReference callback
    let subscription = source |> Observable.subscribe (fun x ->
      match callbackRef.Target with :? ('a -> unit) as callback -> callback x | _ -> ())
    new WeakSubscription<'a> (callback, subscription) :> IDisposable

しかし subscribeWeakly という関数名はもうちっとなんとかならんかなー、と。これじゃあまるで週刊で購読してるみたいだ。（いちおうスペル違うけどさ）