ゆるふわブログ

プログラミング,大学の勉強,日常生活で感じたことをゆるふわに書いていけたらなと思います.技術的に拙いところがあっても温かい目で見守っていただければ幸いです.

pV 図上の吸発熱の判定

こんにちは,Ysmr-Ry です.
今回は僕が高校のとき曖昧だったことをまとめてみたいと思います.
意外とコレがなかなか教えてくれなかったりする.

話題は高校の熱力学です.(高校の熱力学の基礎的な知識を仮定します)
pV 図というのを書くじゃないですか.囲む面積が仕事になるので…
大学に行くとエントロピーという状態量を習って,TS 図とかいうのもかけたりするんですね.
それは囲む面積が熱を表すんです.pV と TS どちらも次元がエネルギーになってるからなんですけど.

つまり高校ではエネルギーを供給する熱と仕事のうち仕事にフォーカスを当てて,そちらをわかりやすいようにしてるのです.物事というのは往々にしてトレードオフであり,そうすると熱がわかりにくくなってしまうんですね.

今回はそのわかりにくい熱のお話です.

熱力学とは

熱力学というのは系の熱平衡状態を仮定して,状態量というマクロな少数パラメータで熱力学的状態を記述し,その上での法則を考えようというものです.モル数を固定すると,状態量としては,圧力  p,体積  V絶対温度  T,内部エネルギー  U などがあります.(高校で出てくるのはたぶんこれだけ)
大学に行くともっとたくさん出てくるのですが… そのうちで最も重要なのがエントロピー  S です.まぁざっくり言うと熱を司り,エントロピーが変化しないのなら熱がなく断熱過程です.(ここではその程度の認識で十分)
で,状態量の自由度は 2 なのです.状態量としてはいっぱいあるのですが,そのうちの 2 つの量がわかれば,他の状態量は既知の状態量との関係式をもとに原理的に決定できるのです.(僕はこういう理解なのですが間違ってたら教えてください)
冒頭で pV 図や TS 図が出てきましたが,これらは 2 次元平面です.平面上の 1 点を与えれば,実数 2 つ分の情報量を与えることができ,それが熱力学的な状態と 1 対 1 に対応しているというわけです.
で,高校での主役は pV 図です.その上での熱の話をしたいのです.熱を司るのはエントロピーなのです.
あらわに目にする状態量はこの場合  p, V なわけですが,その裏で  S が働いているという構図です.
 p, V S の関係を熱力学的な法則から考えたいのです.
まずは,"あらわでない状態量が裏で働いている" 様を体感してもらいましょう.

等温過程

あらわでない状態量として,絶対温度  T を考えます.
例えば, T=298K の状態を集めてくると pV 図上ではどうなるのでしょうか.熱力学的状態は pV 図上の 1 点に対応しているので,pV 図上の領域になることは容易にわかります.
それに加えて,自由度が 2 であることも考えると,そのうちの 1 つである  T を固定しているため,実質的な自由度は 1 になります.この場合自由度を次元と置き換えることができ,1 次元,すなわち曲線になることが予想されます*1
高校物理の熱力学で考える作業物質は理想気体がほとんどですから,具体的な曲線は Boyle の法則が与えます.
Boyle の法則により,等温であるとき, pV が一定となります.
 xy 座標平面で言えば  xy=const. であり,これは反比例のグラフです.(直角双曲線) (曲線にちゃんとなった)
また,理想気体の状態方程式から, pV の値は  T に比例します.
よって, T をとびとびに増やしてやると,このように平面を覆っているのがわかります.(一つ一つが  T が一定な状態を集めた双曲線)

右上に行くほど一定である  T の値が大きくなります.
これが "あらわでない状態量が裏で働いている" 例です.

断熱過程

では,裏で働いている状態量がエントロピー  S である場合を考えましょう.
熱の議論をするためには基本的に熱力学第一法則から出発しなければなりません.
熱力学第一法則は内部エネルギーの変化を  dU,気体がする仕事を  d'W,熱を  d'Q として, d'Q=dU+d'W となります.高校では  \Delta を使って書いていましたが,厳密には微小量です*2.ここで  d'W = pdV です.作業物質は理想気体なので, dU = nc_VdT です.( c_V は定積モル比熱)
すると,
\begin{align*}
\frac{d'Q}{T} &= nc_V\frac{dT}{T}+\frac{p}{T}dV \\
&= nc_V\frac{dT}{T}+nR\frac{dV}{V}
\end{align*}
1 行目から 2 行目へは状態方程式  pV=nRT を使いました.
両辺を積分すると,
 \begin{align*}
\int \frac{d'Q}{T} &= nc_V\log T+nR\log V \\
&= nc_V\log TV^\frac{R}{c_V} \\
&= nc_V\log \frac{pV^\frac{c_V+R}{c_V}}{nR}
\end{align*}
これが断熱で 0 ならば,Mayer の関係  c_P = c_V+R と比熱比  \gamma=\frac{c_P}{c_V} を使って, pV^\gamma = const. です.
Poisson の式として教科書にぽつりと載っていたりします.この式の表す pV 図上の曲線を断熱曲線といいますが,そのある接点での接線の傾きは, \frac{dp}{dV} = -\gamma\frac{p}{V} になります.接点の表す状態から接線に沿って進むことは断熱過程を表します.その点からの過程においてこのラインが吸発熱の境界線になります.
ちなみにエントロピーの定義が  \int \frac{d'Q}{T} だったりするので,いわば "裏でエントロピーが働いている" 曲線なわけです.

変形の仕方を少し変えると,
\begin{align*}
d'Q &= nc_VdT+pdV \\
&= nc_V\frac{Vdp+pdV}{nR}+pdV \\
&= \frac{c_V}{R}Vdp+\left(\frac{c_V}{R}+1\right)pdV
\end{align*}
となります.1 行目から 2 行目へは,状態方程式から, dT = \frac{d(pV)}{nR} = \frac{Vdp+pdV}{nR} を使いました.
 d'Q > 0 で吸熱のときは,
 \begin{align*}
\frac{c_V}{R}Vdp+\left(\frac{c_V}{R}+1\right)pdV &> 0 \\
\frac{dp}{dV} &> -\gamma\frac{p}{V}
\end{align*}
となります.
 d'Q < 0 で発熱のときは  \frac{dp}{dV} < -\gamma\frac{p}{V} です.
つまり,ある 1 点の状態からの過程は,そこを接点とする断熱曲線の接線よりも右上ならば吸熱,左下ならば発熱なのです.

図の紫の曲線が断熱曲線です.オレンジの点の状態から,紫の接線に沿って進めば断熱過程,青い線に沿って進めば吸熱,赤い線に沿って進めば発熱です.

つまり,

このように,青い吸熱の領域と赤い発熱の領域に分けられているイメージです.(この場合はオレンジの 1 点に関しての話です.)

おまけ

以下のような熱サイクルを考えてみましょう.

f:id:Ysmr_Ry:20171231030407p:plain

先ほどの話から,左,上に進む過程はそれぞれ発熱,吸熱と決まりますが,斜めの過程はそうとは限りません.
特にこの場合,途中で切り替わるようにできています*3.切り替わるのはこの斜めの線と断熱曲線が接している点です.

右下の方のオレンジの点が切り替わる点です.試しにそれよりも左上のもう 1 つのオレンジの点を見てみましょう.
この点を通る断熱曲線の接線がオレンジの線です.この線よりも熱サイクルの青の斜めの線の方が右上にあるので,この点では吸熱です.接点をすぎた先の点では発熱に切り替わります.
これは, \frac{d^2p}{dV^2} = \gamma\frac{p}{V^2} > 0 より断熱曲線が下に凸なことによります.

理想気体の場合  \gamma=\frac{5}{3} なので*4,この斜め線の方程式と Poisson の式を連立すれば切り替わる接点の座標がわかります.計算すると  (\frac{15}{8}V_0,\frac{9}{8}p_0) であり,この熱サイクルの熱効率は  \frac{48}{97} となります.(間違ってたら教えてください)

*1:ここで直線といいたくなるかもしれませんが,途中折れ曲がっていてもよく,折れ線を限りなく細かくしていくと曲線になるので一般には曲線です.

*2:僕の先生は熱力学第一法則を  \Delta で書くやつを信用するなとまで言っていました.

*3:過去に東工大がこのサイクルを出題したそうです.途中で切り替わるのは出題ミスではと僕の先生が言ってました.

*4:ここらへんは統計力学の結論です.

高速フーリエ変換 (FFT)

お久しぶりです,Ysmr-Ry です.
試験範囲なので高速フーリエ変換 (FFT) について僕の理解を記しておこうと思います.
高速フーリエ変換とは離散フーリエ変換を高速に計算するためのアルゴリズムで,コンピューターで (離散) フーリエ変換をするために重要な手法です.
フーリエ変換・離散フーリエ変換については,下記のようなサイトで調べればわかると思います.たぶん.僕読んでないけど.

www.ic.is.tohoku.ac.jp

成果物

見て楽しくないとあれかなぁと思ったので一応 Visualizer を作ってみました.
矩形波 (もどき) を Fourier 変換して表示してます.Enter を押すとマイクから音を拾ってきます.
もう一度 Enter を押すとその波形で固定して Fourier 変換してくれます.

CodePen - Fourier Visualizer

また,高速フーリエ変換の理解に欠かせないバタフライダイアグラムを自動で作ってくれるのも作ったので記事と合わせてみると良いかもしれません.↑↓ で数を増減できます.

CodePen - Butterfly Diagram Maker

4 つのポイント

以下の 4 点を押さえていればおそらく FFT を理解できるであろうと思います.

  • 分配法則
  • 約分
  • 二進数
  • 再帰的 (漸化式的) 定義

上の 3 つはすぐわかると思います.4 つ目が曲者なのでそのつもりで.
この 4 点を順に確認した後,離散フーリエ変換の高速化を考えます.

分配法則

ご存知だと思いますが一応.分配法則とは,
 2 \times (1+3) = 2\times1+2\times3 = 2+6 = 8 という風に,かっこの中の各項に掛け算を分配することを言います.

約分

これも大丈夫でしょう.
 \frac{2}{8} = \frac{1}{4} のようにすることです.

二進数

コンピューターと相性のいい数の表し方です.
我々が普段使っているのは十進数で,桁が上がるごとに桁の担っている数が 10 倍になっていきます.一の位,十の位,百の位,千の位…という風に.
二進数はこれが 2 倍になっていきます.一の位,ニの位,四の位,八の位…という風です.
例えば,11 を表したいなら,1+2+8 より,1011 となります.

再帰的 (漸化式的) 定義

ここまでは「何を当たり前のことを」という感じだったかもしれませんが,これはちょっと難しいです.
これは高校で習ったであろう漸化式を思い出してもらえるとわかりやすいです.
ご存知フィボナッチ数列が良さそうですね.
フィボナッチ数列とは,1, 1 から始めて直前の 2 つの数を足し合わせて次を作っていく数列です.
 1, 1, 2, 3, 5, 8, 13, \cdots
左から数えて n 番目のフィボナッチ数列の数を  f_n と書くと,「直前の 2 つの数を足し合わせて次を作っていく」は,
 f_n = f_{n-1} + f_{n-2}\quad(n\geq3)
と表せます.これが漸化式です.
フィボナッチ数列の作り方のように,最初 ( f_1, f_2) から順に作っていけば  f_n を得られますが,今は何もない状態からぽっと  f_n を聞かれた状況を考えてみましょう.最初から数列を作っていないので何もわかりません.
その時この漸化式を見ると,「n 番目を求めるには n-1 番目と n-2 番目を求めて足せばいいのさ!」と語りかけてきます.
そう,左から n 番目を求める問題を,n-1 番目と n-2 番目を求めると言い換えられるのです.
n-1 番目は n-2 番目と n-3 番目がわかればよく,n-2 番目は…と繰り返していけば,いつか 1 番目と 2 番目がわかればわかるというところまできます.1, 2 番目は 1 と決まっていますから,問題が解けたことになります.
このように,より小さな (この場合は左側) な問題の組み合わせに言い換えていくことができるのです.
よって,問題を解くのに必要なのは,より 小さな問題に言い換えるルール (漸化式)もっとも単純な問題の答え (1, 2 番目が 1) だけです.
これを覚えておいてください.

離散フーリエ変換 (DFT)

 W_N = e^{-i\frac{2\pi}{N}} とし,これを回転子と言います.(いきなり何)
難しそうですがなんのことはなく,名前の通り 360° を N 等分した分時計回りに回ることを表します.
例えば  W_4 によって円周を一周すると図のようになります.

f:id:Ysmr_Ry:20171108232907p:plain

はい.ちなみに円の右端を 0° として時計回りに回ってます.番号の順に回ってます.
 W_N^k W_N の回転を  k 回連続で行うことを表します.
ここでは簡単のため, W_N^k を円の中心から回転で行き着いた点を結ぶ矢印で表します.
例えば, W_4^2 ならば,円の右端から半周するので左端にいます.図の 3 の位置です.
よって,矢印 ← で表します.
すると,大きさ 8 の離散フーリエ変換は与えられた  x_i\,(i=1,2,\cdots,8) に対し,次の行列計算を行うことを指します.

 \begin{align*}
\begin{pmatrix}
X_0 \\ X_1 \\ X_2 \\ X_3 \\ X_4 \\ X_5 \\ X_6 \\ X_7
\end{pmatrix} = 
\begin{pmatrix}
\mbox{→} & \mbox{→} & \mbox{→} & \mbox{→} & \mbox{→} & \mbox{→} & \mbox{→} & \mbox{→} \\
\mbox{→} & \mbox{↘} & \mbox{↓} & \mbox{↙} & \mbox{←} & \mbox{↖} & \mbox{↑} & \mbox{↗} \\
\mbox{→} & \mbox{↓} & \mbox{←} & \mbox{↑} & \mbox{→} & \mbox{↓} & \mbox{←} & \mbox{↑} \\
\mbox{→} & \mbox{↙} & \mbox{↑} & \mbox{↘} & \mbox{←} & \mbox{↗} & \mbox{↓} & \mbox{↖} \\
\mbox{→} & \mbox{←} & \mbox{→} & \mbox{←} & \mbox{→} & \mbox{←} & \mbox{→} & \mbox{←} \\
\mbox{→} & \mbox{↖} & \mbox{↓} & \mbox{↗} & \mbox{←} & \mbox{↘} & \mbox{↑} & \mbox{↙} \\
\mbox{→} & \mbox{↑} & \mbox{←} & \mbox{↓} & \mbox{→} & \mbox{↑} & \mbox{←} & \mbox{↓} \\
\mbox{→} & \mbox{↗} & \mbox{↑} & \mbox{↖} & \mbox{←} & \mbox{↙} & \mbox{↓} & \mbox{↘} \\
\end{pmatrix}\begin{pmatrix}
x_0 \\ x_1 \\ x_2 \\ x_3 \\ x_4 \\ x_5 \\ x_6 \\ x_7
\end{pmatrix}
\end{align*}

行列の 2 行目は 8 等分を 1 個飛ばしで,k 行目は 8 等分を k-1 個飛ばしで回っていることがわかります.
これは通常  O(n^2) の計算量になるのですが,隠れた対称性を見出すことでこの計算量を  O(n\log n) まで落とすのが高速フーリエ変換です.
バブルソートクイックソートぐらい違いますね.

高速フーリエ変換 (FFT)

高速フーリエ変換を一言でいうと,「 \frac{2}{8} って約分して  \frac{1}{4} にできんじゃん!」です.たぶん.
つまり, W_8^2 で 8 等分を 2 個飛ばしで回るのと, W_4^1 で 4 等分を順に回るのとでは区別がつきません.

f:id:Ysmr_Ry:20171109002725p:plain

f:id:Ysmr_Ry:20171108232907p:plain

これです.
これによって,行を半分に分けた時,回転していく様子が "対等" となり分配法則によってくくれるのです.
サイズの小さい場合から具体的に見てみましょう.

n = 2 の場合

離散フーリエ変換は,

 \begin{align*}
\begin{pmatrix} X_0 \\ X_1 \end{pmatrix} = \begin{pmatrix}
\mbox{→} & \mbox{→} \\ \mbox{→} & \mbox{←}
\end{pmatrix}\begin{pmatrix}
  x_0 \\ x_1
\end{pmatrix}
\end{align*}

となります.
ここで,行を半分 (1 個ずつ) に分けてみましょう.
1 行目はどちらも同じです.
2 行目は右側は左側を半周 (←) 回したものです.
よって,
 \begin{align*}
\begin{pmatrix} X_0 \\ X_1 \end{pmatrix} &= \begin{pmatrix}
\mbox{→} x_0 + \mbox{→} x_1 \\
\mbox{→} x_0 + \mbox{←} x_1 \\
\end{pmatrix} \\
&= \begin{pmatrix}
\mbox{→}(x_0 +x_1) \\
\mbox{→}(x_0 + \mbox{←} x_1) \\
\end{pmatrix} \\
&= \mbox{→}
\begin{pmatrix}
x_0 +x_1 \\ x_0 + \mbox{←} x_1
\end{pmatrix} \\
&=
\begin{pmatrix}
x_0 +x_1 \\ x_0 + \mbox{←} x_1
\end{pmatrix}
\end{align*}

 \mbox{→} = W_2^0 = 1 ですからこうなります.
これを図にしてみるとこうなります.

f:id:Ysmr_Ry:20171109004714p:plain

左から線をたどっていきます.合流したら足します.
三角があったら流れてきた数に傍に書いてある回転子を書きます.
この図は, X_0 = x_0 + W_2^0 x_1 = x_0 + x_1 X_1 = x_0 + W_2^1 x_1 = x_0 + \mbox{←} x_1 を表しています.

この図を,蝶のようなのでバタフライダイアグラムといいます.

n = 4 の場合

離散フーリエ変換は,

 \begin{align*}
\begin{pmatrix}
X_0 \\ X_1 \\ X_2 \\ X_3
\end{pmatrix} = \begin{pmatrix}
\mbox{→} & \mbox{→} & \mbox{→} & \mbox{→} \\
\mbox{→} & \mbox{↓} & \mbox{←} & \mbox{↑} \\
\mbox{→} & \mbox{←} & \mbox{→} & \mbox{←} \\
\mbox{→} & \mbox{↑} & \mbox{←} & \mbox{↓} \\
\end{pmatrix} \begin{pmatrix}
x_0 \\ x_1 \\ x_2 \\ x_3
\end{pmatrix}
\end{align*}

今度は偶数行目と奇数行目で考えてみましょう.偶数行目は行を半分に分けた時に同じものの繰り返しですが,奇数行目は半周分ずれています.
よって.

 \begin{align*}
\begin{pmatrix}
X_0 \\ X_2
\end{pmatrix} &= \begin{pmatrix}
\mbox{→} & \mbox{→} & \mbox{→} & \mbox{→} \\
\mbox{→} & \mbox{←} & \mbox{→} & \mbox{←} \\
\end{pmatrix} \begin{pmatrix}
x_0 \\ x_1 \\ x_2 \\ x_3
\end{pmatrix} \\
&= \begin{pmatrix}
\mbox{→} x_0 + \mbox{→} x_1 + \mbox{→} x_2 + \mbox{→} x_3 \\
\mbox{→} x_0 + \mbox{←} x_1 + \mbox{→} x_2 + \mbox{←} x_3
\end{pmatrix} \\
&= \begin{pmatrix}
\mbox{→} (x_0+x_2) + \mbox{→} (x_1+x_3) \\
\mbox{→} (x_0+x_2) + \mbox{←} (x_1+x_3) \\
\end{pmatrix} \\
&= \begin{pmatrix}
\mbox{→} & \mbox{→} \\
\mbox{→} & \mbox{←} \\
\end{pmatrix} \begin{pmatrix}
x_0+x_2 \\ x_1+x_3
\end{pmatrix} \\
\begin{pmatrix}
X_1 \\ X_3
\end{pmatrix} &= \begin{pmatrix}
\mbox{→} & \mbox{↓} & \mbox{←} & \mbox{↑} \\
\mbox{→} & \mbox{↑} & \mbox{←} & \mbox{↓} \\
\end{pmatrix} \begin{pmatrix}
x_0 \\ x_1 \\ x_2 \\ x_3
\end{pmatrix} \\
&= \begin{pmatrix}
\mbox{→} x_0 + \mbox{↓} x_1 + \mbox{←} x_2 + \mbox{↑} x_3 \\
\mbox{→} x_0 + \mbox{↑} x_1 + \mbox{←} x_2 + \mbox{↓} x_3
\end{pmatrix} \\
&= \begin{pmatrix}
\mbox{→} (x_0+\mbox{←}x_2) + \mbox{↓} (x_1+\mbox{←}x_3) \\
\mbox{→} (x_0+\mbox{←}x_2) + \mbox{↑} (x_1+\mbox{←}x_3) \\
\end{pmatrix} \\
&= \begin{pmatrix}
\mbox{→} & \mbox{↓} \\
\mbox{→} & \mbox{↑} \\
\end{pmatrix} \begin{pmatrix}
x_0+\mbox{←}x_2 \\ x_1+\mbox{←}x_3
\end{pmatrix} \\
\end{align*}

ここまでが第 1 段階です.これは n = 2 の場合とやっていることが全く同じです.
n = 2 の場合の  x_0 +x_1 が偶数行目に  x_0 + \mbox{←} x_1 が奇数行目に対応しています.(n = 2 の場合はそれぞれ 1 個ずつだったのが,今は 2 個ずつ)
しかし,添字の番号が違っていますね. x_0 + x_1 に対応するのが  x_0+x_2, x_1+x_3 です.ここにどのような対応があるのでしょう?
答えは二進数にあります.添字を二進数で表してみると,(0,2) が (00,10), (1,3) が (01,11) です.二進数の一番上の桁に注目すると,どちらの組も (0,1) となっていることがわかります.これがルールです.
これが第 k 段階ならば上から k 桁だけに注目します.あとで詳しく考えます.
次に,
 \begin{align*}
\begin{pmatrix}
X_0 \\ X_2
\end{pmatrix} &= \begin{pmatrix}
x_0+x_2 + \mbox{→}(x_1+x_3) \\ x_0+x_2 + \mbox{←}(x_1+x_3)
\end{pmatrix} \\
\begin{pmatrix}
X_1 \\ X_3
\end{pmatrix} &= \begin{pmatrix}
x_0+\mbox{←}x_2 + \mbox{↓}(x_1+\mbox{←}x_3) \\ x_0+\mbox{←}x_2 + \mbox{↑}(x_1+\mbox{←}x_3)
\end{pmatrix}
\end{align*}

で完成です.(第 2 段階)
これをバタフライダイアグラムで表すと,

f:id:Ysmr_Ry:20171109093939p:plain

となります.左側の列の蝶の集まりが第 1 段階,右側が第 2 段階に対応しています.
第 1 段階は n = 2 の場合に対応すると言ったように,n = 2 の時のバタフライが 2 つ入っています.
サイズ n の FFT を実行することを, fft_n(x_i) と表すとすると, fft_4(x_0,x_2,x_1,x_3) = fft_2(x_0,x_2) + fft_2(x_1,x_3) + (\mbox{一番右のバタフライ}) です.これが漸化式です.

n = 8 の場合

これも同様です.バタフライダイアグラムだけ示すと,

f:id:Ysmr_Ry:20171109094547p:plain

となります.1 列目 (第 1 段階) には n = 2 の場合のバタフライが 4 つ,1, 2 列目を合わせると n = 4 の場合のバタフライが 2 つ隠れています.
この場合も, fft_8(x_0,x_4,x_2,x_6,x_1,x_5,x_3,x_7) = fft_4(x_0,x_4,x_2,x_6) + fft_4(x_1,x_5,x_3,x_7) + (\mbox{一番右のバタフライ}) が漸化式になります.

先ほど確認したように,問題を解くには漸化式の他にもっとも単純な問題の答えが必要です.
もっとも単純なのはこの場合 n = 1 の場合を考えると簡単です.n = 1 の場合は,フーリエ変換する前もなくそのままになります.
 fft_1(x_i) = x_i です.これで問題が解けます.

解けるんですが,あと問題なのは,バタフライダイアグラムの左側の  x_i の並びです.これにはある規則があります.

なぜビットリバースの順なのか

バタフライダイアグラムの左側の  x_i の右に書いてあるかっこの中の数は添字を二進数で表したものです.
これをそれぞれ逆から読んでみると,000, 001, 010, 011, 100, … で,0, 1, 2, 3, 4, … と順に増えていることがわかります.これをビットリバースの順と言います.
なぜこうなるのか考えてみましょう.漸化式に沿って考えてみると,右から 2 番目のバタフライが 2 つある段階 (第  \log_2 n-1 段階) において,添字を二進数で表した時の上から  \log_2 n-1 桁だけ見ています.言い換えると,一番右の桁だけ無視しています.なので,( \log_2 n-1 桁)0 と ( \log_2 n-1 桁)1 の 2 つのセットがあります.これが 2 つのバタフライに対応します.
漸化式から一番右のバタフライと 2 つの半分のサイズの  fft に言い換えられるため,半分のサイズのバタフライダイアグラムに今と同じことを順に適用していくと,n = 2 の場合,一番上の桁が 0, 1 なので,ビットリバースの順に増えていくことに対応します.
例えば,n = 4 の場合,第 1 段階は 0, 1 ですが,第 2 段階まで考えると,その後ろに 0 をつけるか,1 をつけるかの 2 セットあります.これを順に書くと,00, 10, 01, 11 でビットリバースの順になっています.
n = 8 ならば,00, 10, 01, 11 の後ろに 0 をつけるか 1 をつけるかで,000, 100, 010, 110, 001, 101, 011, 111 です.
ビットリバースの順に増えていくのは,そこまでの列の後ろに 0 を加えたものと 1 を加えたものをくっつけて拡張するということを繰り返していくと作れます.
ここにも再帰的定義が隠れているのです.
(伝わる?)

まとめ

個人的にネットに転がってる FFT の説明がわかりにくかったので自分なりにまとめてみました.
徐々に難しくしていったつもりなんですが,ビットリバースのところがうまく伝わったか心配です…
不明な点があったらコメントください.間違っていた場合もお願いします.

行列と複素数 (小ネタ)

とりあえずなんでもいいから書かなきゃなぁと思ったので簡単な小ネタ.


複素数って回転拡大を表すじゃないですか.
 r(\cos\theta+i\sin\theta) って  r 倍拡大, \theta 回転を表しますね*1
ちょっと書き方を変えると  a+bi = \sqrt{a^2+b^2}\left(\frac{a}{\sqrt{a^2+b^2}}+\frac{b}{\sqrt{a^2+b^2}}i\right) という複素数 \sqrt{a^2+b^2} 倍拡大, \phi = \arctan \frac{b}{a} 回転を表します.


回転で行列といえば回転行列というのもありますね.
 \theta 回転ならば  \begin{pmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{pmatrix} です.
拡大はこの行列の前に拡大率をかければよいので, a+bi の表す回転拡大に対応する行列は  \sqrt{a^2+b^2}\begin{pmatrix} \cos\phi & -\sin\phi \\ \sin\phi & \cos\phi \end{pmatrix} = \begin{pmatrix} a & -b \\ b & a \end{pmatrix} となります.(というか基底を直交する  (a,b) (-b,a) に張り替えると考えた方が早い.)


逆にこの行列を  a+bi とみなせば見通しがよくなります.

例えば,線形代数で †随伴行列† ってやるじゃないですか.中二病みたいなマーク使うやつ.
複素行列  A に対して  A^† A を転置して各成分の複素共軛をとります.
これは実数行列の転置に当たる操作なんですけど,複素数バージョンの場合はなんで共軛をとるんだろうって不思議じゃないですか.(初めて習ったときそう思った)
これが各成分の複素数を対応する行列に書き換えるとわかるんです.
例えば,
 \begin{align*}
\begin{pmatrix}
  1 & i \\ i & 1
\end{pmatrix}
\end{align*}
は,対応する行列が  1 \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} i \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} なので,
 \begin{align*}
\begin{pmatrix}
  1 & i \\ i & 1
\end{pmatrix} &= \begin{pmatrix}
  \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} & \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} \\
  \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} & \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}
\end{pmatrix} \\
&= \begin{pmatrix}
 1 & 0 & 0 & -1 \\
 0 & 1 & 1 & 0 \\
 0 & -1 & 1 & 0 \\
 1 & 0 & 0 & 1
\end{pmatrix}
\end{align*}

こいつの転置をとると,
 \begin{align*}
\begin{pmatrix}
 1 & 0 & 0 & 1 \\
 0 & 1 & -1 & 0 \\
 0 & 1 & 1 & 0 \\
 -1 & 0 & 0 & 1
\end{pmatrix} &=  \begin{pmatrix}
  \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} & \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \\
  \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} & \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}
\end{pmatrix} \\
&= \begin{pmatrix}
1 & -i \\ -i & 1
\end{pmatrix} 
\end{align*}

 \begin{pmatrix}
  1 & i \\ i & 1
\end{pmatrix}^† になってます.複素共軛は対応する行列の転置に当たるので,複素数という行列のブロックを転置して,各々の共軛を取れば全てを展開した実行列の単なる転置に当たるのです.†納得†.

あとは指数行列ですね.指数関数に行列をぶち込むっていう.
 e^x = 1+x+\frac{x^2}{2}+\cdots から類推して,
 e^A = E+A+\frac{A^2}{2}+\cdots と定義します.

例えば  e^{\begin{pmatrix} a & 0 \\ 0 & a \end{pmatrix}} = \begin{pmatrix} e^a & 0 \\ 0 & e^a \end{pmatrix} ってなります.
で, e^{\begin{pmatrix} 0 & -\theta \\ \theta & 0 \end{pmatrix}} = \begin{pmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{pmatrix} になるんですね.定義に代入して確かめればいいんですけど,直感的に収束先を捉えたくて.そこでこの小ネタの出番です.


 1 \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} i \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} だったじゃないですか.それぞれ  E,\ I と表します.すると,左辺は  e^{I\theta},右辺は  E\cos\theta + I\sin\theta じゃないですか.どっかで見た事ありますよね?*2 そういう事です.

 e^{\begin{pmatrix} a & -b \\ b & a \end{pmatrix}} = e^a\begin{pmatrix} \cos b & -\sin b \\ \sin b & \cos b \end{pmatrix} もわかると思います.


ちなみに  \begin{pmatrix} 0 & \theta \\ \theta & 0 \end{pmatrix} みたいなのは複素数を拡張した分解型複素数なる "双曲線の" 複素数を考えれば納得できそうです.いつもの複素数は "円の" です. j^2=1 なる虚数単位チックな  j を導入して  a+bj とかするみたいです. j J=\begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} に対応付ければ 2 乗と 1 ( E) なのでつじつまが合います.

 e^{j\theta}=\cosh\theta+\sinh\theta なる関係が成り立つので同じように  e^{\begin{pmatrix} 0 & \theta \\ \theta & 0 \end{pmatrix}} = e^{J\theta} = E\cosh\theta+J\sinh\theta = \begin{pmatrix} \cosh\theta & \sinh\theta \\ \sinh\theta & \cosh\theta \end{pmatrix} です.

Lorentz 変換とかに関係ありそう.


あと全然関係ないんですが,なんで共軛とるかわからないといえばあれです.複素ベクトルの内積をとるときって片方の共軛をとるじゃないですか. \vec{a}=\begin{pmatrix} a_1 \\ a_2 \\ a_3 \end{pmatrix}, \vec{b}=\begin{pmatrix} b_1 \\ b_2 \\ b_3 \end{pmatrix} なら, \vec{a}\cdot\vec{b} = \bar{a_1}b_1+\bar{a_2}b_2+\bar{a_3}b_3 みたいに.あれがなんでかって話ですね.まぁ  \bar{z}z=|z|^2 なので, \vec{a}\cdot\vec{a} = \lvert a_1\rvert^2+\lvert a_2\rvert^2+\lvert a_3\rvert^2 = \lvert\vec{a}\rvert^2 となるようになってるといえばおしまいです.1 次元複素ベクトル (複素数) 
 \alpha = \lvert\alpha\rvert e^{\arg\alpha}, \beta = \lvert\beta\rvert e^{\arg\beta}内積を考えると, \bar{\alpha}\beta = \lvert\alpha\rvert\lvert\beta\rvert e^{\arg(\beta-\alpha)} となります.これの実部をとると 大きさ x 大きさ x cos (なす角) なので内積っぽいです.

計算幾何とかで複素数クラスをベクトルに流用する場合,このようにして内積外積をスマートに記述したりします.

以上.行列を複素数に対応付ければ見通しが良い例,無数にありそうですが僕が観測した範囲ではこんな感じです.他にあったら教えてください.

もっと本格的な記事も気が向けば書きます (書くとは言ってない)

*1:知らない人は複素平面で調べればわかるかも

*2: e^{i\theta}=\cos\theta+i\sin\theta,Euler の公式で調べてみてください

ボロノイ図を描く! ~Fortune's Algorithm~

こんにちは! Ysmr-Ry です.
Twitter の知り合いにボロノイ図の論文を見せてもらって興味が湧いたのでボロノイ図の描き方を勉強してみました!
ボロノイ図とはなんなのか,そしてその衝撃の描き方を見てみましょう.
ビーチに波が押し寄せるようにボロノイ図ができてゆきます.(面白そう)
長くなってしまいましたが,なぜ放物線なのかがわかれば十分な気がします (ごめんなさい)

ボロノイ (Voronoi) 図とは

ボロノイ図とは,与えられたいくつかの母点のうちどれに一番近いかで領域を分けた図のことです.
実際にボロノイ図を見てみましょう.

f:id:Ysmr_Ry:20170906172456p:plain

これがボロノイ図です.黒い点が母点 (ボロノイ図を描くために与えられた点) で色分けされてるのがそれぞれのボロノイ領域です.
よく見てると気づくとは思いますが,隣接するボロノイ領域の境界線は対応する母点を結んだ線分の垂直二等分線です.

球面上で作ったり,いろいろ拡張はあるようです.
次のリンク先は地球上の空港を母点として球面上にボロノイ図を描いています.

World Airports Voronoi

Fortune's Algorithm

ボロノイ図の描き方は当然いろいろあるとは思いますが,今回紹介するのは Fortune's Algorithm という方法です.
Fortune's Algorithm は平面走査法の一種で,円錐曲線としての放物線に着想を得ています.
簡単言うと,ボロノイ図を描きたい平面上を横線が上から下に動いていってスキャンし,それに放物線を重ね合わせた "波" が続いて,波の継ぎ目がボロノイ領域の境界をなぞっていきます.
実際見た方が早いですね.こんな感じです.( 点の配置によってはバグるのでそういう時はうまくいくまでページを更新してください((((( )

Voronoi Visualizer

水色の線は放物線の組み合わせでできている "波" で,それよりも上の領域が "海" と考えると良いです.
黒い横線が平面上を走査して行って,緑色のボロノイ図が描かれていくのがわかります.
赤い丸はなんなのか,右の木 (こういう形のデータ構造を木と言います) は何を表しているのか,はこの後説明します.

前提知識 1 (円錐曲線としての放物線)

何と言ってもこのアルゴリズムの面白いところは放物線が生きてくるということです.
放物線,2 次関数自体は中学生でもやると思いますが,ここで使う放物線の性質は現在の課程では数 III の 2 次曲線のところで出てきます.
それは,「ある定直線  d とその上にない定点  F があるとき, d F から等距離にある点の軌跡は放物線である」です.
定直線  d を準線 (directrix),定点  F を焦点 (focus) と言います.
アニメーションで見るとこんな感じです.

紫の定点と定直線がそれぞれ焦点と準線,赤いのがそれによってできる放物線,オレンジ色の円は焦点を通り定直線に常に接しているような円です.
ということはこの円の中心は焦点と準線から同じ距離 (半径) にあるので,放物線の上を動くというわけです.
ちなみにこれを一般化して焦点と準線からの距離の比が  \varepsilon となるような点の軌跡を離心率  \varepsilon の二次曲線といい,放物線・楕円・双曲線になります.

前提知識 2 (Binary Search Tree)

先ほどの Voronoi Visualizer の右側に表示されていた木を二分探索木 (Binary Search Tree; BST) といいます.
データ構造とか言ってますがそんなに難しい話じゃないので見ていきましょう.
木というのは先ほど見ていただいた通りの形をしたデータ構造なのですが,中でも二分木はそれぞれの頂点 (丸のこと) から伸びる枝が高々二股に分かれているものを言います.
それぞれの頂点には例えば数字が入っています.二分探索木は,「ある頂点に注目した時にその中に入っている値に対して,そこから左に生えている枝の先の部分木に入っている値は全てそれよりも小さく,右は全て大きいという性質を保つもの」です.
図で見ていきましょう.次のような木は二分探索木と言えます.

f:id:Ysmr_Ry:20170906184935p:plain

一番上の 5 の入った頂点 (根と言います) に注目します.すると,左の部分木は下の四角で囲った部分です.

f:id:Ysmr_Ry:20170906185322p:plain

見てみるとこの中のどの値も 5 より小さくなっていることがわかります.(1, 2, 3 < 5)
また右の部分技についてはどの値も 5 より大きくなっていることがわかります.(5 < 6, 8, 11)
今回は根に注目しましたがこれの性質はどの頂点についても成り立っています.(例えば 2 の頂点なら 1 < 2 で 2 < 3 です)
こういう風に数値のデータを持っていて何がうれしいかというと,探索木ですから,この中からある値を探す時に効率がいいのです.
例えばこの二分探索木から 6 を探したいとしましょう.
根から順にたどっていきます.根に入っている値は 5 です.5 < 6 なので,先ほど見た性質から,6 は少なくとも右の部分木にあることがわかるので,右の枝を辿ります.
これを繰り返せば,最悪の場合でも木の深さ分だけ比較をすれば見つかることがわかります.

f:id:Ysmr_Ry:20170906190834p:plain

難しいことを言うと木の深さ  d に対してその深さにある頂点はだいたい  2^d 個あるので計算量は (平衡状態であれば)  O(\log n) となります.
Fortune's Algorithm では二分探索木を放物線の波の最前線を管理するのに使います.

Fortune's Algorithm (なぜ放物線を考えるのか)

では本題の Fortune's Algorithm について考えていきましょう.
Fortune's Algorithm は平面走査法なので,走査線 (Sweep Line) が平面上を上から下に動きます.線が掃いていくから Sweep です.
ただこのアルゴリズムの場合は特殊で,もう一つ放物線の組み合わせの "波" である汀線 (Beach Line) を考えます.
これがなぜかということです.
今回はボロノイ図を構成するためのアルゴリズムなので,走査線よりも上の掃き終わった領域においてはボロノイ図の一部が組みあがっていて欲しいのです.
しかしながら,走査線よりも上ならば下にある母点たちの影響を全く受けないとは言えません.
そこで,それよりも上ならば下の母点たちから一切影響を受けないというもう一つの境界線,汀線を考えるのです.
先ほど言ったように汀線の上側の領域を "海" と呼称すると,ボロノイ図の出来上がった領域である "海" がじわじわと平面全体を覆っていくということです.
ではなぜその汀線が放物線となるのか.それは前提知識 1 で見た放物線の性質によるものです.
汀線は走査線が掃き終わった,上にある全ての母点をそれぞれ焦点とし,現在の走査線を準線とした放物線よりも上の領域のうち,もっとも下にある境界線をつぎはぎしてできた最前線として決めます.(日本語が不自由)
例えば,先ほどの Visualizer で放物線を全て表示すると

f:id:Ysmr_Ry:20170906193730p:plain

となるのですが,実際に考えるのは一番下にある,最前線の放物線の組み合わせです.
そうするとそれらの放物線は「走査線と対応する母点から等距離にある点の軌跡」なので,放物線の上の領域は「走査線よりも母点の方が近い点の集まり」になります.走査線よりも下の点は当然走査線よりも遠いですから,"海" は走査線よりも下の点から影響を受けないことがわかります.
また,隣接する放物線の交点を Breakpoint と呼ぶのですが,その点は準線と二つの母点からの距離が等しいのでボロノイ図の辺の上に乗ります.

Fortune's Algorithm (2 つのイベント)

平面走査法はアニメーションとしては連続的に走査線が動いていくのですが,実際のプログラムでは特別なイベントが起こる時点をとびとびにたどっていきます.
今回の場合そのイベントは 2 つあります.

Circle Event

f:id:Ysmr_Ry:20170906201010p:plain

図のように 3 つの母点が同一円周上にあり,その円の下側に走査線が接する時点を言います.(したがって上の図は Circle Event が起こる少し前)
このまま走査線が下に行くと,赤い円の中心に集まってきている 3 つの放物線の一部のうち,真ん中が完全に汀線からなくなります.
また,円の中心は Circle Event のまさにその時点では 3 つの放物線に同時に乗っていると言うことができ,それに対応する 3 つの母点 (と走査線) から等距離にあるため,ボロノイ図における点 (ボロノイ点) になります.
したがってこのイベントが起きた時は真ん中の放物線を汀線の放物線のリストの中から消去し,また,ボロノイ点で終わっている辺の終わりにボロノイ点をセットし,ボロノイ点から始まる辺の始点をセットします.

Site Event (NewPoint Event)

f:id:Ysmr_Ry:20170906210749p:plain

Site Event は母点の一つに走査線が触れた時点で起きます.上の図は走査線が 3 の母点を掃いた直後です.
母点 3 を焦点とする放物線が生じるため,それを汀線に組み込まねばなりません.

Fortune's Algorithm (実装の詳細)

先ほどのイベントは上から順番に並べて順に処理します.この時 y 座標が最小のものを効率よく取り出すデータ構造が必要で,Priority Queue というのを使うみたいです.ちなみに最初は母点から Site Event だけを作って並べて置いて,それを処理ながら Circle Event を発生させてイベントリストに追加します.
このアルゴリズムで使うもう一つのデータ構造は先ほども紹介したように二分探索木です.
頂点の持つ値として何を入れているかというと,Breakpoint の x 座標です.
Visualizer で (1,2) と頂点に表示されていたら,母点 1 と母点 2 を焦点とする放物線が左からこの順で交わっている時の交点の x 座標が入っています.(実質的に)
Site Event のところの図がわかりやすいと思います.(再掲)
f:id:Ysmr_Ry:20170906210749p:plain

この場合は,(0,2) (2,3) (3,1) (1,0) の順で左から Breakpoint が並んでいるのがわかります.(このうち画面内に見えているのは (2,3) (3,1) だけですが…)
ちなみにこの順序は二分探索木を深さ優先探索 (Depth-First Search) すると得られます.

Circle Event においては真ん中の放物線を消し (つまり真ん中と左,真ん中と右の作る Breakpoint を消す),新たな Breakpoint として円の中心を追加します.
Site Event では新たに放物線ができるので,既存の汀線とどこでぶつかるかを調べて,新たに Breakpoint を追加しなければなりません.
関係する汀線の部分を探索するのに効率が良いのです.

また実装の際には母点の座標と準線の位置さえあれば放物線も Breakpoint も計算できるのでそれさえ持っておけば OK です.
辺は始点だろうが終点だろうがどのみち Circle Event の時の円の中心になるのでその時決定すれば良いです.

詳しくは Haskell で実装してる人がいるのでそちらでどうぞ.

github.com

また,他の言語での実装もあります.(英語だけど)

Fortune’s algorithm and implementation

おわり

日本語の資料少ないので苦労しました…

http://www.cs.sfu.ca/~binay/813.2011/Fortune.pdf
http://cs.nyu.edu/yap/bks/egc/09/7Vor.pdf

ちなみに僕の Visualizer のソースコードは以下にあります.Haskell の Haste というライブラリを使って JavaScript のコードを生成しています.
したがって描画は HTML5Canvas です.
僕が確認しただけでも以下のようなバグがあるので,Pull-Request をいただけると助かります…(((((

・最初の母点 0, 1 だけの放物線の場合,片方が表示されずかくかくする (治った?)
・点の配置による原因不明のフリーズ
・放物線の継ぎ目がはみだしてる (CircleEvent が発生していない) ← CircleEvent がおきその先は画面によって切り取られているときにおかしくなる

これアルゴリズムの性質上 Circle Event が起きないと辺が終わらないので画面に切り取られているような辺は自分で計算しなきゃいけないんですよね…
半直線的な.それ僕が改造したんですが,なんかボロノイ辺がありえない感じになる場合があるみたいで… (僕のせい)

github.com

だいぶ下調べはしたはずなんですけど,間が空いちゃって今は Haskell のソースちらちら見ながら書いてるだけなので,アレなところがあったらご指摘お願いします(((((

プラネタリウム作ってみた

こんにちは!Ysmr-Ry です!

まあ知っての通り「放課後のプレアデス*1を愛してやまないので星に関わることをしたくって.

この度ヒッパルコス星表のデータを使い,JavaScriptcanvasプラネタリウムというか動く天球を作ってみました!
成果物で遊ぶだけでも楽しいので見ていってください〜

成果物

これです.

https://ysmr-ry.github.io/Hipplanetarium/

f:id:Ysmr_Ry:20170806104120p:plain

操作方法は左の方に書いてある通りなんですが,一応.

キー 操作内容
A 星座強調 On / Off
S 星座線 On / Off
X 星の名前表示 On / Off
C 星座名 非表示 / 和名 / 英名
Q モード切り替え 天球 / 天球距離反映 / 宇宙空間
Z 拡大
星をクリック 星にフォーカス / 詳細を表示
Space 星座ツアー開始 / 次の星座へ
Shift + Space 前の星座へ
矢印 (←↑↓→) キー 天球を回転

ちなみに星をクリックしてフォーカスしている時に,その星がどこかの星座に属していたならばそのまま Space を押すとその次の星座から星座ツアーが始まります.
クリックしていた星が属している星座にフォーカスしたい時は,Space → Shift + Space とするとできます(
星座線表示しながらモードを SphereDist にすると,すばるの「よりそうように輝く星も、本当は、ひとつひとつが、何光年も遠く遠く離れています。」というセリフも理解できるんじゃないでしょうか.
あ,完全に PC でやることを前提としてます.スマホ勢のみなさんごめんなさい(((

ヒッパルコス星表

ヒッパルコス星表とは,1988 年に打ち上げられた位置天文衛星ヒッパルコスが 4 年間かけて収集した恒星のデータです.
ちなみにヒッパルコスとは古代ギリシア天文学者の名前なのだそう.(英語の綴り違うみたいだけど)
今回の成果物はそのヒッパルコス星表の基本的な Visualizer です.
ちなみに Hipplanetarium という名前は Hipparcos + Planetarium ということです.
(ネーミングセンスが絶望的にない)

ちょっと舞台裏

ちょっと中身の話についてすると,canvas で描画してるので,遠近法の計算は自前でやっていて,3 次元回転には四元数 (Quaternion; クォータニオン) を使ってます.
複素平面とか知ってる人はその 3 次元 Ver. と思っていただければいいと思います.
まあその程度しか知らないんですが.

あとは基本データをうまく処理して描画するだけだったんですが,ちょっと苦労したのが星の色ですね.
ヒッパルコス星表には安直に星の色が RGB でいくら,みたいなことは載ってないので,そこはあるデータから自分で計算せねばならんのです.
今回は bv Color という値から下記の式で表面温度を推定し,そこから Plank の黒体輻射の分布式を XYZ 等色関数で数値積分して XYZ 表色系での値を求め,RGB に線形変換してます.

f:id:Ysmr_Ry:20170806102538p:plain

このサイトを参考にしました.

www.geocities.jp

ただ,これで計算すると RGB が負になったりしてよくわかんなかったところはあるんですね.
ちょっとだけ色のお勉強をして見た*2んですが,どうも RGB ですべての色が表せるという単純な話ではないようで,黒体の色温度がディスプレイで再現できる色の範囲を出てしまってるってことなんだと思いますたぶん.こんな風に黒体軌跡が色域の三角形を出てる感じ?

f:id:Ysmr_Ry:20170806103148p:plain

あとは色相を保って明るさ最大にしてるので,色の違いが目立っていいのかなぁと.
星の色ちゃんとあってますかね… さそり座のアンタレスが赤いから大丈夫だろうという謎の基準()
おかしかったら星に詳しい人教えてください〜

あとの詳細が知りたい人は僕にコメントとかで聞いてもらうか,Github にソースが上がってるので解読して貰えばいいと思います.(拙くて申し訳ないですが)

github.com

これの docs ってフォルダーの中です.

ちなみに

やっぱり同じようなこと考えた人はいて,もっとやばいことやってるのがありました.

HippLiner -- A 3D interstellar spaceship simulator with constellation writing function

でも星の知識とか絶望的に足りないし今回はとりあえずこんな感じで… orz

おわり

まあ,僕大昔に同じようなの作ったことあるんですが,その時は星の色とかかなり適当だったので作り直してみようかと思いまして.
それに canvas 使えばまあ Web ブラウザが動けばどこでも見れますし.
あと夏休みで暇だったので()

*1:いい加減しつこいのですがいいものはいいのです笑
なんか昨日はいつきちゃんの誕生日だったらしくて上映会をやってたみたい (?) ですね.
よかったら見てみて〜
sbr-gx.jp

*2:このサイトがわかりやすかったです.
qiita.com
qiita.com
あと,このサイトが参考にしてる本も特に 3 章の表色系の話が数学と絡んできて面白かったです.
( というかそれしか読んでない() )

色彩工学入門 定量的な色の理解と活用 | 森北出版株式会社

相対性理論の本を読んだ

こんにちは,Ysmr-Ry です.

大学が夏休みに突入したので春休みに買って挫折していた相対論の本に再チャレンジしました.
これです.

f:id:Ysmr_Ry:20170730211259j:plain

www.amazon.co.jp

これ,発売日が 2017 年 3 月 27 日だったんですね.
大学の書籍部に電話して予約してたのを思い出します…
1 年の時周りに相対論の授業をとってる人がいて,自分も勉強してみたいな,でも授業受けるだけじゃどうせわからないだろうしな,自分で歩み寄らなきゃなぁと思っていた僕にとっては春休みにこの本が発売されるのは天啓のように感じられましたw
「知る限りこの分野で一番分かり易い本」という評判がついてるくらいですしw
ちなみに再チャレンジしようと思ったきっかけは,「放課後のプレアデス」を見たからなんですね.なんか劇中で,「ブラックホールのシュワルツシルト半径がこんなにくっきり見えるなんて」「宇宙項 Λ が…」とか言ってて,似非物理屋としては「やらねば!」と思ったんですね.
物理屋は食い尽くし,もちろん内容も素晴らしいので,「放課後のプレアデス」一度見てみてはいかがでしょう.

sbr-gx.jp

しかも書いているのは石井俊全さん.先生の線形代数微積分の本には本当に助けられました.僕が大学に入ってからこれらの数学を楽しんで勉強できたのも先生のおかげです!

www.amazon.co.jp

www.amazon.co.jp

また,先生の統計学の本のおかげで統計検定 2 級を S ランクで合格できました…
 \chi^2 分布や  t 分布, F 分布の定義は何か,なぜこの統計検定量がこれらの分布に従うのかが式でわかってよかったです.

www.amazon.co.jp

www.amazon.co.jp

大学入ってから石井さんの本にばっかり助けられてる…

まあ,宣伝はともかく,中身の話を簡単にしましょう.
どこまで理解できてるか怪しいですが… (-_-;)
おかしかったらっょぃ人教えてください!

まず最初の難関はテンソルが理解できるかどうかだと思うんです.
僕の理解はこうです.
僕らの視界は首を左にひねれば右に動きますよね?
つまり,"反"対方向に"変"わるです.また,当たり前ですが視点はひねった方向と"共"に"変"わるわけです.
でも,僕らが見ている世界は実体としては 1 つですよね?
これを逆に考えて,首を左にひねったら視点は同じ方向に,視界は逆の方向に動く.このようにつじつまのあった動き方をする時に,それで捉えているもの (この場合世界) は普遍の実体であることがわかるのです.それを書きあらわす数学的な道具がテンソル
相対性理論では首をひねるのが座標系を変えること,普遍の実体が物理法則です.
普遍な実体を指しているので,テンソルは物理法則を表す数式に用いるのにもってこいというわけです.
この本のテンソルの説明はとてもわかりやすいです.(と言っても僕はこれで初めてテンソルというものを学んだわけですが)
通常物理の本では成分の書き換えの規則で定義するようなのですが,この本では数学寄りらしく,基底を出して説明しています.わかりやすい.

また,高校時代にやった曲率半径の話が再登場するのは感動しましたね.
高校数学はやっぱりすごい.
それを拡張したガウス曲率,それにつながるリーマンの曲率テンソル,接続係数,計量テンソルが時空における我々の空間の曲がり具合の情報を持っていて,それが重力の式や Einstein の重力場方程式に現れるということは,重力は空間の曲がり具合であると言って差し支えないということでしょうか.
曲面人の例えもすごくわかりやすかったです.
3 次元空間内の曲面に住んでいる人は,自分が 3 次元空間にはめ込まれているなんて思わないけど,入手可能な計量テンソルからガウス曲率を計算してみると 0 ではなく,自らの世界が曲がっているのを自覚できると.ガウスはこれを指して "驚きの定理" と言ったそうです.
我々の 4 次元時空もより高次元の空間では曲がっていて,その曲がり具合が重力を生み出すってことですかね.
高次元から見たら低次元の見えないものが見える,っていうのは「正解するカド」の40 (41?) 次元の異方存在を思い出しました.
最後急展開で賛否が分かれますが見てみてはどうでしょう.

seikaisuru-kado.com


また,この本ではあまり触れてませんでしたが,Lorentz 変換で普遍な Minkowski 距離はいわば双曲線的な距離の測り方で,それは Lorentz 変換とその逆変換の対等性から出てくるという話は実は挫折する前に書いていました.よかったらご覧ください.

ysmr-ry.hatenablog.com

この中で登場する「小人さんの妄想」ってブログはものすごくわかりやすいです.
こんなのもあった.

d.hatena.ne.jp

空間の曲がり具合につながる計量テンソルが場所によって違うということは重力が場所によって違うってことですかね.
でも 1 点においては計量テンソルが Minkowski 計量な直線座標である局所ローレンツ系を取れる?
テンソルが座標系の境界,変換表ってどういうことなんですかね.っょぃ人教えて〜

生存数の分布とカノニカル分布の類似性

このような記事を拝見しました*1

www.procrasist.com

この「何ヶ月ブログが続いているか」のヒストグラムが,最近統計力学の試験勉強をしていたこともあり,カノニカル分布*2 (要するに指数分布) に似ていることに気づいたので,考えたことを綴ろうと思います.

統計力学の理解ゆるゆるなので,おかしかったら教えてください.

カノニカル分布とは

統計力学ではミクロな世界の仕組みでマクロな世界の熱力学的物理量を説明するわけですが,そこに確率モデルを用います.
例えば気体分子の数はアボガドロ数 10^{23} オーダーですから,ここに運動方程式を立式して解くというのは現実的ではありません.数が莫大すぎるのです.
そこで,その数の莫大さを逆手に取り,莫大がゆえに確定的な議論ができるように確率モデルを導入するのです.
個々の量子状態に対応する確率を重み付けし,その総和を取れば期待値 (平均) が求められます.ですが,確率を用いている以上,それには不確かさ (ゆらぎ) が存在します.しかしながら,数が莫大がゆえにそのゆらぎは限りなく小さくなり,確率といういわば "不確かな" ものを用いて計算しているにもかかわらず,実質的には確定的な物理量を議論できるのです.そのような期待値として,ミクロの理論からマクロな物理量をひねり出します.
さて,ここで,各量子状態に対して,どのように確率を重み付けするかということが問題になります.
熱平衡状態として実現するマクロな状態を司るミクロな量子状態は,全ての量子状態の中で圧倒的大多数を占めており,また,その中でミクロな量子状態が異なってもマクロには区別がつきません.
そこで一番愚直には,全ての量子状態に一様に同じ確率を重み付けするというのが考えられます.すなわち,全ての状態数を  W とすると, \frac{1}{W} で確率が一定として重み付けするのです.
これをミクロカノニカル分布といいます.
もう少し凝ると,熱平衡状態を注目している系に対して圧倒的に巨大な熱浴にくっつけて実現していると考えて,頑張ると,量子状態  i のエネルギーが  E_i として,次のような確率モデルが考えられます.

f:id:Ysmr_Ry:20170720234055p:plain

 Z(\beta) は分配関数といい,単に確率の総和を 1 にするためのつじつま合わせの定数です.
要するに,エネルギー  E を持つ確率は  e^{-\beta E} (指数関数) に比例するということです.
これをカノニカル分布といいます.

カノニカル分布と貧富の格差

このように,とある分布が非自明に指数分布に従うと言ったことは,他の場面でも起きていて,(ブログのヒストグラムもその一種ではないかと考えているのですが) そのなかで最たるものが "貧富の格差の分布" です.

www.procrasist.com

d.hatena.ne.jp

カノニカル分布は,大雑把には,全エネルギーが一定の条件下の元でランダムに分配されると,エネルギー  E を持っている確率が  e^{-\beta E} に比例するという分布 (指数分布) に落ち着くということです.
エネルギーをお金だと思うと,世の中に流通するお金の総和は一定として,ランダムに分配がなされると,その分布は指数分布となり,貧乏 ~ 平凡がたくさんいる一方,持っているお金が上がれば上がるほど,それこそ指数関数的に急激に数が減っていくという,この世の理不尽がわかると思います.
つまり,お金の例とカノニカル分布の裏には共通の理屈が走っていることになります.
これがカノニカル分布を "直感的" に捉える上で重要だと個人的には思っています.
では,その理屈とはなんなのか.
それは,「総和一定という条件付きで,一様分布にもっとも "近い" 分布が指数分布である」ということです.

カノニカル分布の情報理論的導出

まあ,全てこの資料を見ていただければ解決なのですが,

http://genkuroki.web.fc2.com/20160616KullbackLeibler.pdf

簡単に説明します.
母集団分布  q に従っている母集団から,標本を抽出し,その経験分布を  p とします.
大数の法則から,抽出する標本の数を無限に飛ばせば  q p は一致します.
しかしながら,今考えたいのは,条件付きで抽出する標本の数を無限に飛ばせば母集団分布  q に対して経験分布  p はどうなるのかということです.
そこで,2 つの分布  p,\ q の "距離" を測る,カルバック・ライブラー (Kullback-Leibler) 情報量  D(p||q) を導入します.

f:id:Ysmr_Ry:20170721002235p:plain

これがなぜ分布の距離を測るのかは,上記資料の 1.1 ~ 1.3 節を読んでいただければ,多項分布と Stirling の公式から鮮やかに導出される様子が見れると思うのですが,
要するに,経験分布がほぼ  p になる確率  r というのが, \log r = -nD(p||q) という関係で結ばれていて, D(p||q) が最小になれば  r が最大となり,その  p が経験分布として実現することを意味することになるのです.
直感的には  D(p||q) が分布  p q の距離を表しており,実際に実現する経験分布  p はこの量を最小化する (もっとも近い) 分布に落ち着くということです.
また,情報理論的な相対エントロピー S(p||q) = -D(p||q) と定義でき,これは最大化します.このエントロピーは熱力学・統計力学に登場するエントロピーと定数倍の差しかなく,エントロピー増大の法則がわかると思います.

今,考えたいのは,条件付き大数の法則ですから,つまり,条件付きで  D(p||q) を最小化するということです.
まあ,お馴染みの Lagrange の未定乗数法を使うわけです.
 p の総和が 1 であること, p によるエネルギーの期待値が内部エネルギー  U であること (期待値 (平均) ではなく,総和でも OK) を等号条件として課し, D(p||q) を最小化すると,

f:id:Ysmr_Ry:20170721004457p:plain

と,元の分布  q e^{-\beta E} をつけただけになります*3
このあたりの詳細が知りたい人は,資料の 2.1 ~ 2.2 節らへんを見ると幸せになれます.
つまり,統計力学で考えているカノニカル分布は,母集団分布  q が全て一定の一様分布になっていることがわかります.
ミクロカノニカル分布が一様分布ですから,ちょうどエネルギーの期待値が内部エネルギーと一致するようにした上で,ミクロカノニカルにもっとも近い分布を自らの拡張として選んだ,ということでしょうか.
よく統計力学の本とか授業とかで Stirling の公式と Lagrange の未定乗数法を組み合わせて導出しますが,その本質は,"分布の距離の最小化" なのです.

つまり何が言いたいの

つまり,ブログのヒストグラムは,ブログの生存数の総和がほぼ一定なのではないかということです.その上で完全に一様に分配されるのなら,上記で語ってきたメカニズムで持って  m ヶ月生き残ったブログがいる確率は  e^{-m} に比例しているのではないかと.
ブログの生存数は,気体分子の中の速度分布や貧富の格差の分布と,本質的に同じ形をしているかもしれないという話です.

*1:リンク貼らせてもらいましたが,問題があればお申し付けください (他のリンクも同様)

*2:統計力学でもっとも有用と言っても良い確率モデル

*3:これから,カノニカル分布にとって当重率の原理 (母集団分布を一様分布にするということ) は本質的ではないことがわかります.