nano_exit

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

Path operator と Coherent Potential Approximation (CPA)

Path operator を次のように用意しておく。


\displaystyle
\tau = G + GtG + GtGtG + ...
\\
\displaystyle
= G \left( 1 - tG \right)^{-1} = \left( G^{-1} - t \right)^{-1}
\\
\displaystyle
= t^{-1} \left[ \left( 1 - tG \right)^{-1} - 1 \right]= \left( t^{-1} - G \right)^{-1} - t^{-1}

 Gはサイト間の移動、 tはサイト上での繰り込まれた相互作用(site T matrix)を表すとする。
後のために tを一応、site potential  vを使って定義しておく。

\displaystyle
t = v + v g v + v + v g v g v + ... = \left( v^{-1} - g \right)^{-1}
 gはサイト内の伝搬を表すと定義するが、後の議論には出てこないのでここでは割愛。

サイトと tが一対一対応すれば良いが、例えば合金のように一つのサイトに対して複数の原子がある確率を持って占有する場合、 tをsite potentialの平均から作る近似(仮想結晶近似:VCA)は粗過ぎる。
何故なら、本来は固定された tのセットを全ての可能な配置に対して求め、それを配置数で割るという平均の取り方をするべきであり、複数種類の tが同時にそのサイトに重みを持って存在している訳ではない。

したがって、配置に対して、すなわちpathの取り方に対して、ある種の平均操作を行うのが妥当だと思われる。
そのような平均操作によって得られたpath operatorを \bar{ \tau }と定義する。
仮に \bar{ \tau }が求めまっていると仮定する。
この \bar{ \tau }を与えるような \bar{ t }を上でpath operatorを用意したように、次の様に定義する。

\displaystyle
\bar{ \tau } = G + G \bar{ t } G + ... = \left( G^{-1} - \bar{ t } \right)^{-1}
\\
\displaystyle
\therefore \; \bar{ t } = G^{-1} - \bar{ \tau }^{-1}
\\
\displaystyle
or
\\
\displaystyle
\bar{ t }^{-1} = \left( G^{-1} - \bar{ \tau }^{-1} \right)^{-1}
= G G^{-1} \left( G^{-1} - \bar{ \tau }^{-1} \right)^{-1} \bar{ \tau }^{-1} \bar{ \tau }
= G \left( \bar{ \tau } - G \right)^{-1} \bar{ \tau }

仮に求まっているとした平均的な場(と呼んで良いかは微妙だが、ここでは有効場と呼ぶことにする)に、注目している相互作用 t_{\alpha}が埋まっているとする。
その場合のpath operator  \bar{ \tau }_{\alpha}は次の様に求められる*(厳密になんでそうなるかは不勉強)。

\displaystyle
\bar{ \tau }_{\alpha}
= \bar{\tau} \left[ 1 + \left( \bar{ t }^{-1} - t^{-1}_{\alpha} \right) \bar{\tau} \right]

これの意味は、次の変形を踏まえると分かり易くなる。

\displaystyle
\bar{ t }^{-1} - t^{-1}_{\alpha}
=  \left( \bar{ v }^{-1} - g \right) -  \left( v^{-1}_{\alpha} - g \right)  
= \bar{ v }^{-1} - v^{-1}_{\alpha}
= \frac{ v_{ \alpha } - \bar{ v } }{ v_{ \alpha } \bar{ v } }
= \frac{ \Delta v }{ v_{ \alpha } \bar{ v } }

つまり、以下の様に、有効場を与える様なsite potential(決してただの平均ではない(と思う))を考慮したいsite potentialで置き換えることで、「埋まる」ことを表現している(なんで v_{\alpha} \bar{ v }で割るのかは不明。しかもこれだとただの一次摂動に思える)。

\displaystyle
\bar{ \tau }_{\alpha}
= \bar{\tau} + \bar{\tau} \frac{ \Delta v }{ v_{ \alpha } \bar{ v } } \bar{\tau}

有効場に埋まっている各相互作用、ここでは二元系として A, B、に対して、有効path operator \tau_A, \tau_Bが求まったが、これらが有効場を構成しているわけだから、この平均がもともとの有効場に一致すると考えれば、

\displaystyle
x_A \bar{ \tau }_A + x_B \bar{ \tau }_B = \bar{ \tau }
 x_A, x_B A, Bの割合である。

これで変数の閉じた自己無撞着方程式が得られた。
この方法でpath operatorを得る方法をCoherent Potential Approximation (CPA)と呼ぶ。

やっぱり、 \bar{ \tau }_{\alpha}の導出を理解しないと、わかった気になれんなぁ。。。

*Ref: V. Popescu et al., J. Synchrotron Rad. 6 (1999) 711.

グランドカノニカル平均からのKeldysh形式に移る前の導入の不満

参考文献:"Nonequilibrium Green's Functions Approach to Inhomogeneous Systems", Karsten Blazer and Michael Bonitz (Springer 2013).
該当箇所:2.1.1 Keldysh Contour

この手のもので自分が最初に詰まってしまうのは、

  • 時間に依存する場合、時間の異なるハミルトニアンの交換を仮定しているのかどうか

が曖昧な点である。

最終的には、一般化された順序演算子を導入して好き勝手に順番を弄れるようにするから何でも良いのだろうが、導入の部分だけを真面目に考えると、適用範囲が狭過ぎて「これ意味あるの??」と思ってしまった。でも導入だからあまり真面目に書かれていない気がする。

時間 tにおけるグランドカノニカル平均を次の様に定義。

\displaystyle
< \hat{A} > (t)
= \frac{ 1 }{ Z_0 } {\rm Tr } \left\{ e^{ - \beta ( \hat{H}(t) - \mu \hat{N} ) } \hat{A}_S  \right\}
= \frac{ 1 }{ Z_0 } {\rm Tr } \left\{ e^{ - \beta \hat{K}( t ) } \hat{A}_S  \right\},
\\
\displaystyle
Z_0 =  {\rm Tr } \left\{ e^{ - \beta \hat{K}( t_0 ) }  \right\}
\\
\displaystyle
\hat{K}( t ) = \hat{H}( t ) - \mu \hat{N}
ハットは演算子であることを表す(ハットの位置がズレて見えるのが不満で仕方がない)。

一般的な時間発展演算子は時間順序演算子 \hat{T}を使って表すと、
{ \displaystyle
\hat{U}( t_1, t_0 ) = \hat{T} {\rm exp} \left( - \frac{ i }{ \hbar } \int^{t_1}_{t_0} dt \; \hat{H}( t ) \right)
}

ハミルトニアンが時間に依存するが、 \left[ \hat{H}(t_1), \hat{H}(t_2) \right] = 0 の時(時間毎にエネルギーがwell defined)、
{ \displaystyle
\hat{U}( t_1, t_0 ) = {\rm exp} \left( - \frac{ i }{ \hbar } \int^{t_1}_{t_0} dt \; \hat{H}( t ) \right)
}

最後に、ハミルトニアンが時間に依存しない場合、
{ \displaystyle
\hat{U}( t_1, t_0 ) = {\rm exp} \left( - \frac{ i ( t_1 - t_0 ) }{ \hbar } \hat{H}  \right)
}

ここで、 t_0における温度因子を次の様に定義。
{\displaystyle
\hat{U}( t_0 - i \hbar \beta, t_0 ) = {\rm exp} \left( - \beta \hat{K}( t_0 ) \right)
}
ここの取り扱いは \hat{K}( t ) = \hat{K}( t_0 ) として、ハミルトニアン(もしくは \hat{K})が時間に依存しないとした時の時間発展演算子に対応させたものになる。

 \hat{U}の中のハミルトニアン \hat{K}に読み替え、ハミルトニアン \hat{K})が時間に依存しなければ、
{\displaystyle
\hat{U}( t_0 - i \hbar \beta, t_0 ) = \hat{U}( t, t_0 ) \hat{U}( t_0, t ) \hat{U}( t_0 - i \hbar \beta, t_0 ) = \hat{U}( t, t_0 ) \hat{U}( t_0 - i \hbar \beta, t_0 ) \hat{U}( t_0, t )
}

これによって、
 
\displaystyle
< \hat{A} > (t) = \frac{ 1 }{ Z_0 } {\rm Tr } \left\{ \hat{U}( t, t_0 ) \hat{U}( t_0 - i \hbar \beta, t_0 ) \hat{U}( t_0, t ) \hat{A}_S  \right\}
\\
\displaystyle
= \frac{ 1 }{ Z_0 } {\rm Tr } \left\{ \hat{U}( t_0 - i \hbar \beta, t_0 ) \hat{U}( t_0, t ) \hat{A}_S  \hat{U}( t, t_0 ) \right\}

これは t_0 \rightarrow t \rightarrow t_0 \rightarrow t_0 - i \hbar \beta という時間発展を考えましょうということを表している。
でもこれ自身は、ハミルトニアンが時間に依存していない時の話だから、こんな仰々しいことしないで松原先生にお願いすることになると思われる。

多分大事なのは、もっと一般のハミルトニアンの時には、直接こうには書けないけれども、「後で順番整理するからとりあえず好きなように並べていいよ!」として同じような時間経路を組んだらどう計算出来るか、という流れなのだと思う。
まあ時間に依存しなきゃ順番なんて何でも良いから、その時に順序演算子は要らなくて、厳密にその形でかけるようになるのは当たり前と言えば当たり前ですが。。。

Formulationは触ったことありますが、実際に運用したことはないので、参考文献を読み進める次第です。

調和振動子と平面波の比較

零点振動している調和振動子と同じエネルギーを持つ平面波とで波動関数を比較する。

\displaystyle
l_0
= \sqrt{ \frac{ \hbar }{ m \omega } }
\\
\displaystyle
\xi
= \frac{ x }{ l_0 }
\\
\displaystyle
\psi_n( x )
= A_n H_n\left( \frac{ x }{ l_0 } \right) e^{ - \frac{1}{2} \left( \frac{ x }{ l_0 } \right)^2 }
= A_n H_n( \xi ) e^{ - \frac{ \xi^2 }{ 2 } }
\\
\displaystyle
A_n
= \sqrt{ \frac{ 1 }{ 2^n \; n! } } \; \sqrt{ \frac{ l_0 }{ \sqrt{ \pi } } }
\\
\displaystyle
H_0
= 1
\\
\displaystyle
E_0
= \frac{ \hbar \omega }{ 2 }
最後ら辺、サボりました。

ここで、 E_0 = E_kとして自由平面波のエネルギーを代入して、波数と有効距離 l_0を結びつけると、

\displaystyle
l_0
= \sqrt{ \frac{ \hbar }{ m } \frac{ \hbar }{ 2 E_k } }
= \sqrt{ \frac{ \hbar^2 }{ m } \frac{ 2 m }{ 2 \hbar^2 k^2_0 } }
= \frac{ 1 }{ k_0 }
綺麗に有効距離の逆数になる。

これで求めた k_0を持つ平面波( cos( k_0 x ) )を n=0調和振動子波動関数と比較する(規格化定数は無視)と、全体で一周期分をちょっと超えたかなくらいの波が調和振動子の中に入っていることがわかる。
f:id:koideforest:20180110005333p:plain

もう少し厳密に見てみると、 k_0を逆に固定して、 k_0 x_0 = \piとなる x_0における調和振動子(係数は無視)の値を求めると、

\displaystyle
e^{ - \frac{1}{2} \left( \frac{ x }{ l_0 } \right)^2 }
= e^{ - \frac{1}{2} ( k_0 x )^2 }
\rightarrow e^{ - \frac{1}{2} ( k_0 x_0 )^2 }
= e^{ - \frac{\pi^2}{2} } = 0.00719... \sim 0.01 \ll 1
したがって、十分小さいと見なせるところまでの範囲で、平面波一周期分が入っていることがわかる。

思っていたよりも調和振動子と平面波の対応関係がしっかりしていてビックリした。
原点近傍でポテンシャルがゼロに近くなるから、その近傍で自由解に近づくから当然といえば当然であるが。
自由解だけ波動関数があるような領域は、運動エネルギー(波動関数)ではなくポテンシャル項によってエネルギーが補填されているんだなぁとしみじみ感じた。

PythonのLaguerre陪関数について

Wikipediaには日頃からよくお世話になっているが、自分で確認するのは大事だなと感じた事件。

水素原子の波動関数の可視化とかカッコイイなぁと思い、ネットサーフィン中に以下のサイト様を訪問。
scipyの特殊関数 - 篠突く雨の日記

自分も早速コピーしてやってみて、ふむふむと散布図の便利さに舌鼓を打っていた。
ちょっとアレンジして、波動関数の正負で色を変えるようにしたら気付いてしまった。
「ん?1s波動関数に節があんぞ??」
波動関数の部分はscipy.special.assoc_laguerreで実装されていて、原因はこいつであることがわかった。
何か引数がおかしいのかと思い、Wikipediaで水素原子の解析解におけるLaguerre陪関数を確認。
水素原子におけるシュレーディンガー方程式の解 - Wikipedia

\displaystyle
\rho \frac{ d^{ 2 } u( \rho ) }{ d\rho^{ 2 } } + ( 2l + 2 - \rho ) \frac{ d u( \rho ) }{ dx } + ( n - l - 1 ) u( \rho ) = 0 \\
\displaystyle
\rho \frac{ d^{ 2 } u( \rho ) }{ d\rho^{ 2 } } + ( ( 2l + 1 ) + 1 - \rho ) \frac{ d u( \rho ) }{ dx } + ( ( n + l ) - ( 2l + 1 ) ) u( \rho ) = 0 \\
\displaystyle
u( \rho ) = c L^{ 2 l + 1 }_{ n + l }( \rho ) \\
\displaystyle
L^m_k( \rho ) = \frac{ d^m }{ d\rho^m } e^{ \rho } \frac{ d^k }{ d\rho^k } ( e^{ -\rho } \rho^k ) \\
Legendreではlとmは直接打ち込めば良かったけど、Laguerreはそういえばややこしかったんだよなぁ、と遠い記憶が蘇ってきた。
しかし、可視化スクリプトの引数はこれと対応したものが入っている。
実は上付きと下付きは逆の表記になっているとかあるか??と思い、他も参照。

ラゲールの陪多項式 - Wikipedia

\displaystyle
\left( x \frac{ d^{ k + 2 } }{ dx^{ k + 2 } } + ( k + 1 - x ) \frac{ d^{ k + 1 } }{ dx^{ k + 1 } } + ( n - k ) \frac{ d^k }{ dx^k } \right) L^k_n( x ) = 0 \\
\displaystyle
L^k_n( x ) = \frac{ d^{ k } }{ dx^{ k } } L_n( x ) = \frac{ d^{ k } }{ dx^{ k } } \left( e^x \frac{ d^{ n } }{ dx^{ n } } ( x^n e^{ -x } ) \right)
上がkだったり下がkだったりしてめちゃくちゃややこしいが、とりあえずLaguerre多項式(下付き添字のみ)を上付き添字の分だけ微分しまくるのが陪関数であることがわかった。
その意味で、水素原子におけるシュレーディンガー方程式の解 - Wikipediaの上付き下付きの意味はラゲールの陪多項式 - Wikipediaと同じ。

scipy.special.assoc_laguerreのサイトに行ってみた。
scipy.special.assoc_laguerre — SciPy v1.0.0 Reference Guide
見る限り、スクリプトの引数の順番も間違っていなさそうだった。

しかし、

scipy.special.assoc_laguerre( x, n + l, 2 * l + 1 )

はn=1, l=0で -x + 1のプロットを示し、意味不明。

そして英語版のwikipediaに辿り着く。
Laguerre polynomials - Wikipedia

\displaystyle
xy'' + ( \alpha + 1 - x ) y' + \beta y = 0 \\
\displaystyle
L^{ ( \alpha ) }_0( x ) = 1 \\
\displaystyle
L^{ ( \alpha ) }_1( x ) = 1 + \alpha - x
これや、、、これやったんや。。。

つまり、scipy.special.assoc_laguerreを水素様原子の波動関数に使うには、
 
\displaystyle
L^{ \alpha }_{ \beta } = L^{ 2 l + 1 }_{ n - l - 1 }
ということになる。
確かにこれで1s (n=1, l=0)の時 L^1_0 = 1となり、めでたく波動関数の節が無くなる。

お手軽に特殊関数が使える様になったとはいえ、油断は禁物ですな。

モンテカルロ積分で水素原子のエネルギー期待値を確認

「スレーター型軌道の数値積分は大変」と呪文のように聞かされて続けて心にハードルが出来てしまったので、ここらで打破することを試みた。
前回、多次元のモンテカルロ積分をまとめた。
koideforest.hatenadiary.com
今回はスクリプトを組んで、どれくらいちゃんと動くのか見てみることにした。
目標は、非相対論水素原子の解析解を使って、ハミルトニアンの期待値を出す。
以下、スクリプトのポイントをまとめる。

  • 運動エネルギーを出すところの微分は、中心差分法を使った。簡単便利。
  • 積分範囲は、水素原子波動関数の値が0.01を下回る距離をrmaxと定義し、その範囲で水素原子の8分の1が入る箱を考えた。対称性から、その箱の中の積分値を8倍すれば全体の積分値が求まる。
  • そんなことをしたのは、乱数の範囲が [ 0, 1 ] の0始まりなので、全体を平行移動したりして水素原子全体を箱の中に無理矢理入れるのは余り得策ではなく感じたから。今回は素直に水素原子の原点を0にした。
  • 箱の中の全体で積分値を求めるから、乱数点に関する制限は必要ない。下手にif分を入れると遅くなるらしいので、そこまで一般化してスクリプトを作るのは余り良くないかもしれない。今回はそこの一般化は行っていない。

以下ソース。

import numpy as np
import time

def MonteCalro3d_ws( func, scales ): # ws: whole space
    M = 10**6
    vec_rand = np.random.rand( M * 3 ).reshape( M, 3 )
    # Caution: the number of samplings is still M, not M * 3
    vec_rand = np.array( scales ) * vec_rand
    tot = np.sum( [ func( v[0], v[1], v[2] ) for v in vec_rand ] )
    return np.prod( scales ) * tot / float( M )

def h1s( x, y, z ):
    vec = np.array( [ x, y, z] )
    r   = np.linalg.norm( vec )
    return 2. * np.exp( - r ) / np.sqrt( 4. * np.pi )

def d2func_dx2( func, x ):
    h = 1.0e-7
    return ( func( x + 2.*h ) - 2.*func(x) + func( x - 2.*h ) ) / (2.*h)**2

def integrand( x, y, z ):
    v = np.array( [ x, y, z] )
    r = np.linalg.norm( v )
    deriv2x = d2func_dx2( lambda x: h1s( x,     v[1], v[2] ), v[0] )
    deriv2y = d2func_dx2( lambda y: h1s( v[0],     y, v[2] ), v[1] )
    deriv2z = d2func_dx2( lambda z: h1s( v[0],  v[1],    z ), v[2] )
    ke = - ( deriv2x + deriv2y + deriv2z ) / 2.
    pot_n = - h1s( v2[0], v2[1], v2[2] ) / r
    return h1s( v[0], v[1], v[2] ) * ( ke + pot_n )

rmax = max( [ x_ for x_ in np.linspace( 0, 5, 1000 ) \
              if h1s( x_, 0., 0. ) / h1s( 0., 0., 0. ) > 0.01 ] )
scales = [ rmax, rmax, rmax ]

start = time.time()
result = 8. * MonteCalro3d_ws( integrand, scales )
print( "Hydrogen atom energy = {}".format( result )  )
end = time.time()
print( end - start )

結果:
-0.500215806919
148.051234961

  • 原子単位系は原子単位系でも、Hartree unitで計算しているので、理論値は-0.5でちゃんと計算できている。
  • ちなみにHartreeかRydbergかは運動エネルギーが 1/2されているかどうかで判断できる。
  •  1/2の無いRydbergの場合には、核ポテンシャルの項も二倍して理論値は-1になる。
  • 細かい値を忘れたが、ざっくり 1 Hartree ~ 27 eV, 1 Rydberg ~ 13 eV, 1 Hartree = 2 Rydberg である。最後の等式は厳密。
  • 値の安定性はサンプリング数 Mに依存するが、 M = 10^6であれば有効数字一桁は安心できると思う。
  • cythonとか使えばもっと速くなる。

N次元のモンテカルロ積分

一般的なアイデアとして、積分を 体積 x 密度(割合) に焼き直すところがスタート。

\displaystyle
 I = \int_{ \Omega } d \textbf{ r } \; f( \textbf{ r } ) \
 = \left( \int_{ \Omega } d \textbf{ r } \right) \
     \left( \frac{ \int_{ \Omega } d \textbf{ r } \; f( \textbf{ r } ) }{ \int_{ \Omega } d \textbf{ r } } \right) \
 = V_{\Omega} < f >_{\Omega}

密度(割合)は単位量当たりにどれだけ対象となる事物があるかを表し、いわゆる期待値として焼き直せる。連続量の期待値なので、(確率)密度 or 重みをかけた積分に対応。それを離散量で近似しましょうというノリ。このあたりは本当は測度論とかが顔を出すんだろうなぁと妄想していますが、不勉強なので回避。
 \displaystyle
 < f >_{\Omega} = \int_{ \Omega } d \textbf{ r } \; \rho_{ \Omega }( \textbf{ r } ) f( \textbf{ r } ) \
 \sim \sum^{ m }_{ \mu_1 } \cdots \sum^{ m }_{ \mu_N } p_{\Omega}( r^{ (1) }_{ \mu_1 }, \cdots, r^{ (1) }_{ \mu_N } )  f( r^{ (1) }_{ \mu_1 }, \cdots, r^{ (N) }_{ \mu_N } )

 mはサンプリング数またはメッシュの細かさ。トータルのメッシュ数は M \equiv m^N で定義する。
一次元と三次元に関して先ほどの和を解して書いてみると、

\displaystyle
\sum^{ m }_{ \mu } p_{\Omega}( x_{ \mu } ) f( x_{ \mu } ) = p_{\Omega}( x_1 )  f( x_1 ) + p_{\Omega}( x_2 )  f( x_2 ) + \cdots + p_{\Omega}( x_m )  f( x_m ) \\
\displaystyle
\sum^{ m }_{ \mu_x } \sum^{ m }_{ \mu_y } \sum^{ m }_{ \mu_z } p_{\Omega}( x_{ \mu_x }, y_{ \mu_y }, z_{ \mu_z } ) f( x_{ \mu_x }, y_{ \mu_y }, z_{ \mu_z } ) \\
\displaystyle
\qquad = p_{\Omega}( x_1, y_1, z_1 )  f( x_1, y_1, z_1 ) + p_{\Omega}( x_2, y_1, z_1 )  f( x_2, y_1, z_1 ) \\
\displaystyle
\qquad + \cdots + p_{\Omega}( x_m, y_m, z_{ m-1 } )  f(  x_m, y_m, z_{ m-1 } ) + p_{\Omega} ( x_m, y_m, z_m )  f(  x_m, y_m, z_m )

一次元において、 0 \leq x \leq a における積分を、等間隔のメッシュを切った長方形で近似する場合で書いてみれば、もう少し期待値の描像が明確になる。
 
\displaystyle
\int^{a}_{0} dx \; f(x) \sim \sum^{ m }_{ i = 1 } \frac{ a }{ m } f( ( i - 1 ) a / m ) \\ 
\displaystyle
\int^{a}_{0} dx \; f(x) = a < f >_a \\
\displaystyle
\therefore \; < f >_a \sim \sum^{ m }_{ i = 1 } \frac{ 1 }{ m } f( ( i - 1 ) a / m ) \
 = \sum^{ m }_{ i = 1 } p f( a ( i - 1 ) / m ) \\
\displaystyle
\qquad \sim \sum^{ m }_{ i = 1 } p f( a R_i )

全空間において一様に離散化・サンプリングすると、それは平均値を取りましょうということになる。
そして、これの近似として、 0 < R_i < 1で規格化された昇順の乱数列 \{ R \}を定義して置き換えたものが、積分領域に制約のない場合のモンテカルロ積分となる。
近似関係が分かりやすいようにあえて昇順と言ったが、結局は和なので順番は関係ない。
とりあえず多重積分したいという場合にはこっちの考え方になるだろう。いちいちメッシュの切り方を自分で考えなくて良いので、手軽で良い。
 pが単純にサンプリング数の逆数なので、モンテカルロ積分の例題で有名な \piを求める問題(円の中に入った点の数を数える)は、この考えでは出来ない。

これまでの話は領域 \Omegaの中で閉じていた。
なので、今までの p_{\Omega}はあまり確率の意味はなく、単純に重みの役割が大きかった。
この領域 \Omegaを扱い易い任意の箱 \Boxに拡張するが、その際に余分な領域での値がゼロになるような関数 \Thetaを定義する。
 \displaystyle
\int_{ \Omega } d \textbf{ r } \; f( \textbf{ r } ) \
 = \int_{ \Box } d \textbf{ r } \; \Theta( \textbf{ r } ) f( \textbf{ r } )
 = V_{ \Box } < \Theta f >_{ \Box } \\
\Theta( \textbf{ r } ) = \left\{ \
\begin{array}{ ll }
 1 & ( \textbf{ r } \in \Omega ) \\
 0 & (otherwise)
\end{array}
\right.

これを先ほどと同様に密度関数で表したのちに離散化すると、

 \displaystyle
< \Theta f >_{\Box} = \int_{ \Box } d \textbf{ r } \; \rho_{ \Box }( \textbf{ r } ) \Theta( \textbf{ r } ) f( \textbf{ r } ) \
 = \int_{\Box} d \textbf{ r } \; \rho( \textbf{ r } ) f( \textbf{ r } ) \\
 \displaystyle
 \sim \sum^{ m }_{ \mu_1 } \cdots \sum^{ m }_{ \mu_N } p( r^{ (1) }_{ \mu_1 }, \cdots, r^{ (1) }_{ \mu_N } )  f( r^{ (1) }_{ \mu_1 }, \cdots, r^{ (N) }_{ \mu_N } )

上と同様に一次元の例で考えれば、
 
\displaystyle
\int^{a}_{0} dx \; f(x) = \int^{b}_{0} dx \; \theta( a - x ) f(x) \sim b \sum^{ m }_{ i = 1 } \frac{ \sigma( ( i - 1 ) b / m ) }{ m } f( ( i - 1 ) b / m ) \\ 
\displaystyle
\int^{a}_{0} dx \; f(x) = b < \theta_a f >_b \\
\displaystyle
\therefore \; < \theta_a f >_b \sim \sum^{ m }_{ i = 1 } \frac{ \sigma( ( i - 1 ) b / m ) }{ m } f( ( i - 1 ) b / m ) \
 = \sum^{ m }_{ i = 1 } p( b ( i - 1 ) / m ) f( b ( i - 1 ) / m ) \\
\displaystyle
\qquad \sim \sum^{ m }_{ i = 1 } p( b R_i ) \; f( b R_i ) \\

ここで、

\displaystyle
\sigma( x ) = \left\{
\begin{array}{ll}
 1 & ( x \in [ 0, \; a ] ) \\
 0 & (otherwise)
\end{array}
\right.
と定義した。

 \sigma及び pで、積分領域の調整を乱数点が領域に入った入らないかで行うことが出来る。積分領域が複雑な時には、この方法が威力を発揮することになる。
乱数自体は [ 0, 1 ]の間でしか値が出ないでの、計算したい系に合わせてスケール因子を掛ける必要がある。
一つの次元を c倍する度に結果を c倍すれば正しい積分結果が得られる。

ラグランジュの未定乗数法と等高線プロット

ラグランジュの未定乗数法を直感的に非常に分かりやすく描いてあるサイト様を発見。
ラグランジュの未定乗数法の解説と直感的な証明

せっかくなので、実際に等高線プロット使って図示してみた。
f:id:koideforest:20171116003805p:plain

図で考えているのは、

 f(x,y) = ( x - 3 )^2 + y^2
 g(x,y) = x^2 + y^2 - 4

のような、何かポテンシャル的な場が f(x,y)で与えられていて、軌道 g(x,y)=0の上で f(x,y)を最大もしくは最小にする (x,y)は何でしょう?という問題。
数学の問題っぽく言えば、束縛条件 g(x,y) = 0の下で f(x,y)を最大(最小)にする (x,y)を求める問題。
文章を読むと、最大最小化する f(x,y)が主役な感じがするが、図を見るとわかる様に、(当たり前だが) g(x,y)=0の軌道の中の点で良さそうなのを探す方がイメージに近い。
ちなみに、二次元の今の問題において得られる最大(最小)の座標群がなぜ零次元(つまり点)なのかというと g(x,y)=0の時点で x yは独立ではなく、既に一次元の問題になっているところが味噌である。そこに f(x,y)を介して条件式が加わると、零次元になるというわけである。三次元でやると得られる座標群は線(つまり軌道・一次元)になる。
単純に変数 x, yそれぞれに対して微分するから、方程式が未定乗数を含めた変数の数だけ出てきて、最終的に値がカッチリ決まるというもの。未定乗数の微分では g(x,y) = 0がそのまま出てくる。
この意味で、主役はなんとなく g(x,y)=0のような気がしてならない。

図を見れば明らかに
最大: ( -2, 0 )
最小: ( 2, 0 )
であるが、ラグランジュの未定乗数法を使うと、

 L = ( x - 3 )^2 + y^2 - \lambda [ x^2 + y^2 - 4 ]
 \frac{ \partial L }{ \partial x } = 2 ( 1 - \lambda ) x - 6 = 0
 \frac{ \partial L }{ \partial y } = 2 ( 1 - \lambda ) y = 0
 \frac{ \partial L }{ \partial \lambda } = - ( x^2 + y^2 - 4 )= 0

ここから、

 x  = - 2 / \lambda
 ( y, \lambda ) = ( 0, a ) \; or \; ( a, 1 ) \; for \; \forall a \in \mathbb{R}
 x^2 + y^2 = 4

三つ目のは元々の束縛条件まんまなので、当たり前は当たり前。
 y = 0を入れれば、 x = \pm 2が求まって、図から類推した答えと合っている。
 \lambda = 1 x = -2と等価であり、その時も y=0となるから、全ての答えが得られている。

以下、ソースの中身。

import numpy as np
import matplotlib.pyplot as pet

N = 100
x = y = np.linspace( -5, 5, N )
X, Y = np.meshgrid( x, y )
F = ( X - 3 )**2 + Y**2
G = X**2 + Y**2 - 4

plt.text( 0, 0, 'g(x,y)=0' )
plt.text( -2.5, 3.5, 'f(x,y)', color='#0068b7')

plt.contour( X, Y, G, [0] )
plt.gray()
plt.pcolor( X, Y, F )
plt.cool()
plt.gca().set_aspect( 'equal', adjustable='box' )
plt.show()

以下のサイト様にお世話になりました。

等高線プロット
matplotlibカラープロット - 週末はいつも晴れ

陰関数
Pythonのmatplotlibでグラフを簡単に表示する方法part2 | 東大卒ニートと普通のサラリーマンのお金の稼ぎ方