2016-12-20

逐次代入法を用いた超シンプルな例と収束因子

数学

$x = 1-x$

を考える。もちろん答えは $x = 0.5$ である。
これを敢えて逐次代入法で解いてみよう。この方法は固体物理や第一原理計算で用いるSCF法と同じノリである。

最初に初期値を設定する。答えに近ければ近い程良いが、とりあえず $x_0 = 0.7$ とおく。これを右辺に代入したものが左辺を与えるとして、

$x_1 = 1 - x_0 = 0.3$
$x_2 = 1 - x_1 = 0.7$

と次々に項が求まる。が、よく見ると全然0.5に近づいておらず、このままだとずっと振動して $x$ が求まらない。

そこで収束因子 $\alpha$ を導入する。 $\alpha$ は新しい $x$ をどれくらい古い $x$ に混ぜるかというもので、例えば $x_1$ は次のようになる。

$x'_0 = 1 - x_0$
$x_1 = \alpha x'_0 + ( 1 - \alpha) x_0$

$\alpha$ が小さい程、古い $x$ の情報が残り、 $x$ の変化が緩やかになる。
具体的に値を入れてその振舞を見てみると、

$\alpha = 0.2$
$x_0 = 0.7$
$x'_0 = 0.3$
$x_1 = 0.2 \cdot 0.3 + 0.8 \cdot 0.7 = 0.06 + 0.56 = 0.62$
$x_2 = 0.2 \cdot 0.38 + 0.8 \cdot 0.62 = 0.072 + 0.496 = 0.568$

というように振動が抑えられ、0.5に近づいていくのが分かる。
では本当にこれを続けて0.5になるかどうかを解析的に解いてみよう。

$x_1 = \alpha ( 1 - x_0) + ( 1 - \alpha) x_0 = x_0 (1 - 2\alpha) + \alpha$
$x_2 = x_1 (1 - 2\alpha) + \alpha = x_0 ( 1 - 2\alpha)^2 + \alpha (1 - 2\alpha) + \alpha$
$\vdots$
$x_n = x_0 ( 1 - 2\alpha)^n + \alpha ( 1 + (1 - 2\alpha) + (1 - 2\alpha)^2 + \cdots + (1 - 2\alpha)^{n-1})$

よって、これを無限に繰り返すと、 $|1 - 2\alpha| < 1$ であれば、

$x_{\infty} = \alpha \frac{1}{1-(1-2\alpha)} = \frac{1}{2}$

となり、初期値にも収束因子にも依存せず答えがちゃんと求まることが示せた。
$|1 - 2\alpha| < 1$ 、つまり $0 < \alpha < 1$ という条件は、元々が新しい $x$ と古い $x$ の割合ということを思い出せば全く自然な範囲である。
ちなみに、最初に $x_0=0.5$ を代入すると一発で求まるし、 $\alpha=0.5$ でも一発で求まるので、問題によって適切な初期値と収束因子を選ぶことが重要であることがわかる。しかし、一般に収束因子に関しては物理的なものではないので類推は困難であり、経験的に適切だと思われる値を入れることが多いと個人的には思う。
最適な収束因子を求める問題とか数値計算の方でありそうな気がするが、不用意に飛び込むと全ての数値計算法に満足出来なくなりそうで怖いのでやめておこう。。。

2016-12-19

初等的な縮退の話

量子力学

サクライの章末問題に載っててフムフムとなった縮退の話。

演算子 $A_1, A_2$ が互いに交換せず( $[ A_1, A_2 ] \neq 0$ )、かつそれぞれがハミルトニアンと同時固有状態を作るとする $[ A_1, H ] = 0, [ A_2, H ] = 0$ 。
このときに、一般に縮退が存在することが証明出来る。

これを真っ向から挑む場合に、必要十分であることを言うのはなかなかしんどい気がする。ホワンとした証明しか思いつかん。
これは背理法を使うと結構スッキリ話を持って行けて気持ちいい。

背理法として、縮退が無いとまず仮定する。あるエネルギー固有状態 $|E_i>$ を $A_1$ の固有値 $a^{(1)}_i$ でラベル付けるか、 $A_2$ の固有値 $a^{(2)}_i$ でラベル付けるか好きな方を選択出来るが、 $[ A_1, A_2 ] \neq 0$ の前提条件により $a^{(1)}_i$ と $a^{(2)}_i$ を一つの状態に同時に指定することは出来ない。しかし、今の仮定では縮退が無いため、

$|a^{(1)}_i, E_i > = |a^{(2)}_i, E_i >$

である。これは次の式で見るように、 $[ A_1, A_2 ] \neq 0$ に抵触している。

$( AB - BA ) |a^{(1)}_i, E_i > = (AB-Ba^{(1)}_i) |a^{(1)}_i, E_i > \\ = (AB-Ba^{(1)}_i)|a^{(2)}_i, E_i > = (A a^{(2)}_i - a^{(2)}_i a^{(1)}_i) |a^{(2)}_i, E_i > \\ = (A a^{(2)}_i - a^{(2)}_i a^{(1)}_i) |a^{(1)}_i, E_i > = (a^{(1)}_i a^{(2)}_i - a^{(2)}_i a^{(1)}_i) |a^{(1)}_i, E_i > \\ = 0$

よって、縮退が無いという仮定が誤りであり、実際には縮退していることが言える。

これの具体例で良さそうだなぁと思ったのが点群で、例えば $C_{3v}$ を考えて、 $A_1 = C_3$ 、 $A_2 = \sigma_v$ と対応付けられて、図を描くとわかるが、 $C_3 \sigma_v = \sigma''_v$ 、 $\sigma_v C_3 = \sigma'_v$ （鏡映面の異なる鏡映操作）となり、交換しない。加えて、系が $C_{3v}$ に属していればハミルトニアンは $C_{3v}$ の対称操作に対して変わらない。よって、エネルギー固有状態は縮退している。
実際、 $C_{3v}$ は二次元の既約表現 $E$ を持ち、二重に縮退する。もちろん、一次元の全対称既約表現 $A_1$ （どの対称操作に対しても不変）しかないような部分空間を用意するならば、話は別で縮退は無い（この部分空間においては $[ C_3, \sigma_v ] = 0$ が成立）。
でもやっぱ点群とかは図が無いと説明が絶望的に伝わらんなぁ。。。でもサクライだと角運動量 $L_z$ と $L_x$ で説明してて、交換しないのがただ数式的な感じがして個人的には縮退のイメージがピンとこなかったんだよなぁ。。。悩ましい。。。

2016-12-16

二次元ハミルトニアンとパウリ行列

量子力学

$H = a ( |1><1| - |2><2| + |1><2| + |2><1| )$ の固有値を求める問題をパウリ行列使って解く方法が割と楽しい（サクライでは章末問題になっている）。
これを行列で書くと、

$\displaystyle H = a \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix} = a \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} + a \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \\ = a ( \sigma_z + \sigma_x ) = a \sqrt{2}( \sigma \cdot \hat{n} )$

ただし、 $\hat{n} = \frac{1}{\sqrt{2}}(1,0,1)$ の単位ベクトルである。
せっかくなので、パウリ行列をスピン演算子に直してあげれば、 $S_i = \frac{\hbar}{2}\sigma_i$ より

$\displaystyle H = \sqrt{2} a \frac{2}{\hbar} ( {\bf S} \cdot \hat{n} )$

で、変形して何なんだという話であるが、 ${\bf S} \cdot \hat{n}$ の固有値は $\pm \frac{\hbar}{2}$ で物理的に明らか（ $\hat{n}$ 方向を向いたスピンの固有値を表している。なぜそうなるのは自明というよりかは自然（実験事実）からの要請に近い。）であるため、このハミルトニアンの固有値は $E = \pm \sqrt{2}a$ であることが行列式を使わずにわかってしまう。

ちなみに固有ベクトルだが、違う章末問題で ${\bf S} \cdot \hat{n}$ の固有値を回転行列を使わないで真面目に固有値方程式を解く趣旨のものがあり、そっちも教育的で楽しい。固有値はもうわかっているから展開係数を直接求めるのだが、そういうプロセスはあんまりないなぁとボンヤリ思ったりもした。いずれそれもまとめたい。

2016-12-15

ハミルトニアンの行列表示

量子力学

J.J. Sakuraiのゼミをしていたとき、ハミルトニアンの行列表示について議論になった。

$\displaystyle H \doteq \begin{pmatrix} H_{11} & H_{12} & \cdots & \\ H_{21} & H_{22} & & \\ \vdots & & \ddots & \end{pmatrix}$

Diracのbraketを用いた表示で書けば、

$\displaystyle H \doteq \begin{pmatrix} <1|H|1> & <1|H|2> & \cdots & \\ <2|H|1> & <2|H|2> & & \\ \vdots & & \ddots & \end{pmatrix}$

となる。
Sakuraiだと、イコールではないことを強調して一応 $\doteq$ を用いてはいるものの、結構ポッと出な感もあるし、何よりも $<1|H|1>$ の中の $H$ は外の $H$ とどう違うの？とか、もう一回中に代入するってこと？みたいな誤解が起きることが分かった。確かに $H \psi(x) = E \psi(x)$ に直接行列表現ぶち込もうとすると訳分からんことになる。

ポイントは一本の二階微分方程式を解くモードか、微分を含まない連立方程式で解くモードかということだと思っている。

シュレーディンガー方程式から出発して、全然分からん固有状態 $\phi$ を素性が良くわかっている固有状態 $\phi_i$ で展開する（例えば何にも相互作用の無いときの固有状態である平面波とか）。

ここで、例えば $<\psi_1|$ を左から掛けることで、 $H|\psi>$ にどれだけ $\phi_1$ の情報が含まれているかを見る様なことをすると、

$\displaystyle <\phi_1| H |\psi> = <\phi_1| E |\psi> \\ C_1 <\phi_1|H|\phi_1> + C_2 <\phi_1|H|\phi_2> + \cdots = E C_1$

これを他の状態 $<\phi_i|$ でもやってあげれば、展開係数 $C_i$ に関する連立方程式が出来る。

$\displaystyle C_1 <\phi_1|H|\phi_1> + C_2 <\phi_1|H|\phi_2> + \cdots = E C_1 \\ C_1 <\phi_2|H|\phi_1> + C_2 <\phi_2|H|\phi_2> + \cdots = E C_2 \\ \vdots$

ここまで、 $H$ は微分方程式（シュレーディンガー方程式）のハミルトニアンのままであり、別に変なことはしていない。この連立方程式を行列でまとめて書くと、

$\displaystyle \begin{pmatrix} <\phi_1|H|\phi_1> & <\phi_1|H|\phi_2> & \cdots & \\ <\phi_2|H|\phi_1> & <\phi_2|H|\phi_2> & & \\ \vdots & & \ddots & \end{pmatrix} \begin{pmatrix} C_1 \\ C_2 \\ \vdots \end{pmatrix} = \begin{pmatrix} E & & & O \\ & E & & \\ O & & \ddots & \end{pmatrix} \begin{pmatrix} C_1 \\ C_2 \\ \vdots \end{pmatrix}$

という具合になる。
なので、 $E,\psi$ を求める問題が、 $E, C_i$ を求める問題にすり替わり、特に $E$ を求めたければ $C_i$ は必要なく、

$\displaystyle \begin{vmatrix} <\phi_1|H|\phi_1> - E & <\phi_1|H|\phi_2> & \cdots & \\ <\phi_2|H|\phi_1> & <\phi_2|H|\phi_2> - E & & \\ \vdots & & \ddots & \end{vmatrix} = 0$

を解きさえすれば良い（これがゼロでない場合には、上の行列に逆行列が存在して、 $C_i = 0$ （つまり $\psi = 0$ ）という全然興味のない答えしか得られないため、ゼロ出ないことを条件付ける訳である）。
そうすると、 $H$ の微分方程式における具体的な表式は、「行列要素の値が求まっていれば」どうでも良くなり、行列形式で問題を扱って行くときにいちいち違う記号を用いるのは煩雑になる（人によって違う記号を使いだすと多分ヤバい）ので、同じ $H$ を使っているのだと理解している。
最初の話に戻ると、行列の中の $H$ は実際のシュレーディンガー方程式の $H$ で、 $\doteq$ の左側にいる $H$ はただ $H$ という記号を使っている行列ということになる。

2016-12-15

fortranとpython（読み込み）

プログラミング

fortranユーザーからpythonに移って来ると、ただテキスト編集したいのにも関わらず、fortranの癖に慣れてしまったが故にpythonがなかなか使えるようにならんかったのでちょっと整理。

「ファイルを開いて中身を読んで閉じる」までをfortranとpythonで比較。

fortran

implicit none
integer, parameter :: INPUT = 10
real(kind(0d0))    :: InputData( 50, 50 ), ...
!----------
InputData
open( INPUT, 'test.inp', status='old' )
do i = 1, 50
  read(INPUT, *) InputData( i, : )
enddo
close( INPUT )

test.inpが 50 x 50 以上のデータを持っていると仮定して、InputDataの配列サイズを指定。
'test.inp'に対応する番号として10(=INPUT)を割り当ててあげてファイルを開き、i行目の各列のデータをreadでInputDataに流し込む。
それを50行目まで繰り返してファイルを閉じる。

python

f = open( 'test.inp', 'r' )
lines = f.read().split('\n')
f.close()
InputData = []
for i, line in enumrate( lines ):
    temps = line.split(' ')
    while '' in temps: temps.remove('')
    InputData += temps

読み込み形式('r')で'test.inp'を開き、中身をlinesにいきなりリストとして突っ込みファイルは閉じる。.split('\n')によって改行ごとに（つまり各行が）違う要素としてlinesの中に収納されている。lines[0]は一行目全部、lines[1]は二行目全部、といった具合。
それを値として全部ぶつ切りに区別したデータとするため、linesの中の要素line（つまりある一つの行）を空白で区切って（.split(' ')）、空白が連続している場合には空の要素''が生まれてしまうためwhile文を使って除いている。
綺麗に区切ったデータをInputDataに突っ込む作業を全ての行に対して行う。

おそらくpythonでfortranと同じように読ませようとするとこんな感じになるのではないかと思う。
fortranは「行列」を最初から想定しているような節を感じる。型が決まっているなら、その型にハマる様に最初に設定してしまえばあんまり何も考えなくて済む（その代わり型から外れると酷いが）。
一方でpythonは、一本の長いリストに収納してからグニグニ変形させる。かなり自由が利く分、自分がやりたいことを実践するために出来る方法も多く、結果どうしたら良いかわからなくなることも多かった（ファイルを読み込むだけなのにいくつか方法があるというのも新参者には割とストレスになった。。。）。
とりあえず今はこんな感じで落ち着いているが、もっとスマートな方法が多分あると思う。先は蛇のように長い。。。

2016-12-14

奇数次元の虚数について

数学

虚数 $i$ は二次元の正方行列を用いて、

$i_2 = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}$

と書ける。
では、三次元の正方行列で二乗したら単位行列にマイナスの掛かったものが得られるかをボンヤリ考えたところ、

$i_3 = \begin{pmatrix} 0 & 0 & -1 \\ 0 & i & 0 \\ 1 & 0 & 0 \end{pmatrix}$

しか思いつかんかった。この $i$ を二次元 $i_2$ と思えば、super matrixと解釈して各行列要素が二次元正方行列に拡張出来るから、

$\begin{pmatrix} 0 & 0 & 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 & 0 & -1 \\ 0 & 0 & 0 & -1 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \end{pmatrix}$

ちゃんと二乗するとマイナス単位行列になる。もちろん、基底を組み替えたものに対応した、

$\begin{pmatrix} 0 & 0 & 0 & 0 & 0 & -1 \\ 0 & 0 & 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & -1 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}$

もちゃんとなるし、どちらかというとこっちの方が帰納的でしっくり来るだろう。

でも、最初の三次元の $i$ を自分自身の $i_3$ で置き換えると無限に行列の次元が増えていく。
世界の自己相似性に奇数次元であることが絡んでいると面白そうだなぁと思った次第である。

2016-12-13

2GBの壁

プログラミング

*.f90に対して、64bit対応でcompileする分には問題無いが、それを32bitでcompileしたり、*.fの場合には静的メモリに2GBの制限が付くらしい。

relocation truncated to fit: R_X86_64_32S against *...
additional relocation overflows omitted from the output.

的なerrorが出る。
これを回避するには、compile optionに

-mcmodel=large

を付けると良いらしい。ちなみに-mcmodel=mediumにすると、データへの制限は無くなるが、実行ファイルへの2GBの制限は残るそうな。
しかし、これを付けたところで、大元のgfortranなりifortなりなんなりのライブラリーがこれらによってcompileされてなければそこで詰まる。
最近は数十GB超えのPC（workstation)はほぼ当たり前だし、数値計算用に購入したworkstation or serverなんかは100GBのメモリを超えてることはザラだろうから、この辺融通が効くといいなぁと思う次第である。

nano_exit

基礎的なことこそ、簡単な例が必要だと思うのです。

逐次代入法を用いた超シンプルな例と収束因子

初等的な縮退の話

二次元ハミルトニアンとパウリ行列

ハミルトニアンの行列表示

fortranとpython（読み込み）

奇数次元の虚数について

2GBの壁