グロッソ(中)に載ってた。
井戸型 vs クーロン
今更ながらIPythonデータサイエンスクックブックを購入した。もちろん私費である。
せっかくなので一次元のシュレーディンガー方程式でも解こうかと思った。
一階微分方程式に対しては、かくあきさんのサイトで紹介されているローレンツアトラクターのソースが簡潔で分かり易くて非常に参考になった。
Python NumPy SciPy : 1 階常微分方程式の解法 | org-技術
Hartree unitで、
井戸型ポテンシャルでは、
クーロンポテンシャルではゼロ割を防ぐためにを分母に足して、
初期条件として、を用いた。理由は特にない。それっぽいだけ。実はいろいろ理由があったが、微妙に長くなるので割愛。軽く言及すると、井戸型は純粋に一次元の問題だが、クーロンポテンシャルは三次元の問題を球対称として動径部分のみの方程式に落とし込んでいて、を求めていることになる。
エネルギーを振って計算。
結果:
波動関数の大きさは規格化していないので置いといて、やはりクーロンポテンシャルの波の引きずり込み具合は素晴らしい。
Thomas-Fermiの遮蔽ポテンシャルも一瞬考えたが、Fermiエネルギーとかの設定がめんどいので今日はこの辺で。
井戸型ポテンシャルの場合のソースも一応曝しておく。
import numpy as np import scipy.integrate as spi import matplotlib.pyplot as plt import prettyplotlib as ppl def well( x ): a = 1. if -a < x < a: v = -1. else: v = 0. return v """ Initialization """ q0 = np.zeros(2) q0[0] = 0. q0[1] = 1. def f( q, x, e ): p, pdot = q[0], q[1] v = well( x ) pdotdot = - 2 * ( e - v ) * p return np.r_[ pdot, pdotdot ] x = np.linspace( 0., 10., 100 ) for e in np.linspace( 0.2, 1., 5 ): q = spi.odeint( f, q0, x, args=(e,) ) ppl.plot( x[:], q[:,0], 'o-', mew=1, ms=8, mec='w', label='e={:4.2f}'.format(e) ) ppl.plot( x[:], map( well, x[:] ), 'o-', mew=1, ms=8, mec='w', label='v' ) ppl.legend() plt.ylim( -1.0, 1.0 ) plt.show()
ストークスの定理の最も簡単な例
xy平面内でのみ回転している、大きさ1の回転を考える。
つまり、
ここから元のベクトルを復元すると、
ただし、微分でしか定義されていないので、一意には決まらずに任意性が残る。
例えば以下のようなものも可能である。
今回は上のものを考える。これをgnuplotで図示すると、
グルグルとベクトルが渦を巻いているのが良くわかる。
原点から遠ざかるにつれて、そのベクトルの大きさは大きくなる。
実は、「最も簡単な回転しているベクトル」を想像したとき、これの全部のベクトルの大きさが1のver.を思いついたが、それはが1にならずになって原点で発散、外側で減衰していくので、個人的にはちょっとした発見があった。
これを用いてストークスの定理を考察する。
積分範囲を半径1の円に取る。そのため、線積分の方は、単位円の周りにそって一周する形になる。
面積積分の方は、何も考えずに
と求まる。
線積分の方は体積素片を系に合わせて考察する必要があり、ここが読んでるだけではピンと来ないところではある。
単位円一周なので、の極座標に対してをまで動かして積分すればよい。積分変数をに変換するのを忘れないように気を付けると、
という感じで、面積積分と同じになることが確認できた。
ちなみに、rotationを計算したpythonのコードは以下のものを使用。
import numpy as np def a( x, y ): x1 = 0.2 * x y1 = 0.2 * y ax = -0.5 * y1 ay = 0.5 * x1 aa = np.sqrt( ax**2 + ay**2 ) s = '{:f} {:f} {:f} {:f} {:f}\n'.format( x1, y1, ax, ay, aa ) return s f = open( 'simple_rot.dat', 'w') for x in range(6): for y in range(6): if np.sqrt( ( 0.2 * x )**2 + ( 0.2 * y )**2 ) > 1. : continue f.write( a( x, y ) ) f.write( a( -x, y ) ) f.write( a( x, -y ) ) f.write( a( -x, -y ) ) f.close()
gnuplotのコマンドは以下の通り。
set terminal png
set output 'simple_rot.png'
set xr[-1.2:1.2]
set yr[-1.2:1.2]
set palette rgbformulae 31,13,10
pl 'simple_rot.dat' u 1:2:3:4:5 w vector lc palette
unset output
条件反射シリーズ:ガンマ関数
とが一緒に入っている積分:
ガンマ関数に持って行けないか、疑ってみましょう。
ちなみにガンマ関数の定義は、
の指数がだったかだったかややこしいし、の肩がプラスだったかマイナスだったか忘れるけど、とにかく
との組み合わせ =>
と思っておけば見通しは良いはず。
ガンマ関数を使いそうと思ってからググっても遅くはないはず。
例:Sommerfeld展開で見かけるあいつ
ガンマ関数にまとめるために、をの形に持って行けるように変形します。
の形にするために、分数を級数展開します。
これによりとまとまりますが、前に掛かっていると指数の肩を合わせるために、の積分変数の変換を行うと、に注意して、
という感じで、ガンマ関数にまとめることが出来ました。めでたし。
それで、前に掛かってる奴は、見るからにゼータ関数を含んでいるので、この形に持って行くと、
となり、具体的なの値に対して関数は求まっているので、計算が出来るということになります。
追記:
から分かる通り、ゼータ関数とガンマ関数は繋がる。
積分の分母がじゃなくて、のとき、分数の級数展開にはマイナスは出て来ないから、もの凄くストレートに
となる。
周期的境界条件の位相の和の関係
前回、Tight-Bindingで使った位相の和の関係についてまとめておく。
koideforest.hatenadiary.com
これらを証明する。
簡単のため、格子定数の一次元格子を考える。
格子点個を含む範囲(今の場合は長さと表現できる。2,3次元では箱とよく呼ばれる)を設定し、このが周期的に繰り返しているとする。ここに割と自分は嵌まった。もともと格子定数で既に周期的であるが、わざわざそれよりも大きい入れ物で周期性を考えるのである。このの中で、格子点はが含まれるとする。とが同じ値を取るとするのが周期的境界条件である。
この一次元格子上で平面波が満たすべき条件は、
より、(ただし、は整数)となり、(結晶)運動量が定まる。としない限り、運動量は離散化する。
ここで、この平面波を内で和を取る。と出来るから、
この級数は、
ここでより、
よって、であり、のとき。
のときとなる。
これをまとめると、
となる。kについても同様に、和の範囲をに制限すると、
また、は整数であるため、
あとは、の和と同様の級数の扱いで求まる。
この導出はよく忘れるので、個人的にやっとまとめられたなぁと感慨深くあります。
第二量子化でTight-Binding
※2019/12/09 修正
位置表示での対角項は、後で行うフーリエ変換(級数)しても対角的なので、省略する。
電子の移動のみをハミルトニアンで扱うと、
これを、最近接のみに制限(近似)する。
を最近接のペアを表すとすれば、
ここまでは系に依らない(この項だけを扱うのが正当化されるかどうかはもちろん系によるが、それはまた次元の違う近似の話である)。
ここで一次元格子に制限すると、と具体化される。正直、ここまで具体化しないと全然良くわからん。
注意として、
であるため、逆方向への飛び移りはエルミート共役に対応することがわかる。
したがって、 は逆方向の飛び移りだから、
これを対角化したいというのが人間の性というものである。
位置基底が固有状態 => 局在(動かない)
運動量基底が固有状態 => 非局在(動く)
という発想の元、フーリエ変換(級数)すると、
変換の変換が元に戻るか確かめると、
ここで、
を使用。
iの和の場合は、
細かい話は省略。
この変換により、例えば上の一次元話でi+1からiへの移動の項を考えると、
は隣のサイトの距離である。今、一次元の格子だから格子定数と言って良い。
がサイト間の相対座標にしか依らないとすると、で和が取れるから、
となり、kに対して対角的になる。つまり、動くものは運動量の固有状態になり、隣のサイトに行った分、位相がずれる。
ただし、無限小変位に対して不変ではないため、運動量は離散的であり、結晶運動量と呼ばれる。
結局、
よって、対角化が完了。
が絶対座標に依存したら並進対称性がないわけで、それはバンドにならんことがわかります。
逐次代入法を用いた最もシンプルな例と収束因子:図のテスト
前回、簡単な例で逐次代入法を考察した。koideforest.hatenadiary.com
図の使い方のテストもかねて、その内容を図にしてみた。
収束因子を使わない場合、は以下のようになる。
ただの長方形をグルグルまわるだけで一向に収束しないのが分かる。
収束因子を使うと、
のようにとの交点にドンドン近づいていくのが分かる。
ちなみに、これらはgnuplotで出力しているが、直線の連続はデータ点を読み込ませていて、とはgnuplotをそのまま使っている。つまり、
pl 'iteration.dat' u 1:2 w l
rep x
rep 1-x
と入力している。
各ファイルの中身は、
iteration.dat:
-0.50000 -0.50000
-0.50000 1.50000
1.50000 1.50000
1.50000 -0.50000
-0.50000 -0.50000
-0.50000 1.50000
1.50000 1.50000
1.50000 -0.50000
-0.50000 -0.50000
-0.50000 1.50000
iteration_broyden.dat:
-0.50000 -0.50000
-0.50000 1.50000
0.90000 0.90000
0.90000 0.10000
0.34000 0.34000
0.34000 0.66000
0.56400 0.56400
0.56400 0.43600
0.47440 0.47440
0.47440 0.52560
である。
これらのデータはpythonで計算させて出力させていて、収束因子を使った方は、
#!/usr/bin/env python def fn( x ): return 1.0 - x f = open( 'iteration_broyden.dat', 'w' ) x = -0.5 alpha = 0.7 for i in map( float, range(5) ): f.write( '{:10.5f} {:10.5f}\n'.format( x, x ) ) f.write( '{:10.5f} {:10.5f}\n'.format( x, fn( x ) ) ) x = ( 1.0 - alpha ) * x + alpha * fn( x ) f.close()
というような感じになっている。