2019-01-09

微分係数の逆数は逆の微分？（続き）

数学

前回、 $x \equiv x( u ), \, x_v=0$ の時に、 $u_x = x_u^{-1}$ であることを示した。
koideforest.hatenadiary.com

しかし、前回得られた式をよく見ると、 $x_v \neq 0, \, y_u = 0$ でも $u_x = x_u^{-1}$ になることがわかる。
$\displaystyle u_x = \frac{ y_v }{ x_u y_v - x_v y_u } = \frac{ y_v }{ x_u y_v } = x_u^{-1}$

$x_v \neq 0$ はつまり $x \equiv x( u, v )$ であるわけだが、 $x$ が $v$ に依存しているにも関わらず、 $y_u =0$ のために $v_x = 0$ となることがわかる。
$\displaystyle v_x = - \frac{ y_u }{ x_u y_v - x_v y_u } = - \frac{ 0 }{ x_u y_v } = 0 \neq x_v$
これは、結構直感に反するのではないだろうか？

そこで、 $y$ の微分についても調べることにする。
$u, v$ の微分は連鎖律を用いて、 $x, y$ の微分で表せる。
$\displaystyle \frac{ d }{ du } = x_u \frac{ d }{ dx } + y_u \frac{ d }{ dy } \\ \displaystyle \frac{ d }{ dv } = x_v \frac{ d }{ dx } + y_v \frac{ d }{ dy }$

ここから逆に $x, y$ の微分に分離すると、
$\displaystyle y_v \frac{ d }{ du } - y_u \frac{ d }{ dv } = \frac{ \partial ( x, y ) }{ \partial ( u, v ) } \frac{ d }{ dx } \\ \displaystyle x_v \frac{ d }{ du } - x_u \frac{ d }{ dv } = - \frac{ \partial ( x, y ) }{ \partial ( u, v ) } \frac{ d }{ dy } \\ \displaystyle \frac{ \partial ( x, y ) }{ \partial ( u, v ) } \equiv x_u y_v - x_v y_u$

整理すると、
$\displaystyle \frac{ d }{ dx } = \left( \frac{ \partial ( x, y ) }{ \partial ( u, v ) } \right)^{-1} y_v \frac{ d }{ du } - \left( \frac{ \partial ( x, y ) }{ \partial ( u, v ) } \right)^{-1} y_u \frac{ d }{ dv } \\ \displaystyle \frac{ d }{ dy } = - \left( \frac{ \partial ( x, y ) }{ \partial ( u, v ) } \right)^{-1} x_v \frac{ d }{ du } + \left( \frac{ \partial ( x, y ) }{ \partial ( u, v ) } \right)^{-1} x_u \frac{ d }{ dv }$

ここまでは一般的。ここで、 $y_u = 0$ の条件を課すと、
$\displaystyle \frac{ d }{ dx } = x_u^{-1} \frac{ d }{ du } \\ \displaystyle \frac{ d }{ dy } = - \frac{ x_v }{ x_u y_v } \frac{ d }{ du } + y_v^{-1} \frac{ d }{ dv }$
と表せる。
$x$ が $u, v$ 両方に依存する影響が、 $y$ の方に出ていることがわかる。

これが合っているのか、簡単な計算で確かめることにする。
$\displaystyle x = 2 u + v, \quad y = v \\ \displaystyle f = x + y = 2u + 2v \\ \displaystyle \frac{ d f }{ dx } = 1 = \frac{ 1 }{ 2 } \cdot 2 = x_u^{-1} \frac{ d f }{ du } \\ \displaystyle \frac{ d f }{ dy } = 1 = - 1 + 2 = - \frac{ 1 }{ 2 \cdot 1 } \cdot 2 + \frac{ 1 }{ 1 } \cdot 2 = - \frac{ x_v }{ x_u y_v } \frac{ df }{ du } + y_v^{-1} \frac{ df }{ dv }$
ちゃんと答えが一致した。
どちらかが一変数にしか依存しない場合、対角成分では逆数の関係が成り立つ（ $u_x = x_u^{-1}, \, v_y = y_v^{-1}$ ）が、非対角成分には注意が必要である。
もちろん、 $x,y$ の両方が $u,v$ に依存する場合には、逆数の関係は成り立たないので、逆数に期待し過ぎるとちょっと危険かもしれない。

2019-01-09

微分係数の逆数は逆の微分？

数学

以下の式は一般に正しいだろうか？
$\displaystyle \frac{ d u }{ d x } = \left( \frac{ d x }{ d u } \right)^{-1} \quad (?)$

$x \equiv x( u )$ の時は正しい。
$\displaystyle \frac{ d }{ d u } = \frac{ d x }{ d u } \frac{ d }{ d x } \\ \displaystyle \frac{ d }{ d x } = \frac{ d u }{ d x } \frac{ d }{ d u } \\ \displaystyle \therefore \frac{ d u }{ d x } = \left( \frac{ d x }{ d u } \right)^{-1}$
$x$ が $u$ で定義されている場合、 $\frac{ d x }{ d u }$ の方が求め易いことが多い。
なので、このような関係式をなるべく使いたい訳である。

では、多変数の場合にはどうか？
$\displaystyle x \equiv x( u, v ), \quad y \equiv y( u, v ) \\ \displaystyle \frac{ d }{ d u } = x_u \frac{ d }{ d x } + y_u \frac{ d }{ d y } \\ \displaystyle \frac{ d }{ d v } = x_v \frac{ d }{ d x } + y_v \frac{ d }{ d y }$

これを $\frac{ d }{ dx }$ だけを抜き出すように式変形すると、
$\displaystyle y_v \frac{ d }{ d u } - y_u \frac{ d }{ d v } = ( x_u y_v - x_v y_u ) \frac{ d }{ d x } \equiv \frac{ \partial ( x, y ) }{ \partial ( u, v ) } \frac{ d }{ d x }$
$\frac{ \partial ( x, y ) }{ \partial ( u, v ) }$ がヤコビアンと呼ばれているものである。
したがって、
$\displaystyle \frac{ d }{ d x } = \left( \frac{ \partial ( x, y ) }{ \partial ( u, v ) } \right)^{-1} y_v \frac{ d }{ d u } - \left( \frac{ \partial ( x, y ) }{ \partial ( u, v ) } \right)^{-1} y_u \frac{ d }{ d v }$

一方で、 $\frac{d}{dx}$ を $u,v$ の微分で直接表せば、
$\displaystyle \frac{ d }{ d x } = u_x \frac{ d }{ d u } + v_x \frac{ d }{ d v }$

したがって、
$\displaystyle u_x = \left( \frac{ \partial ( x, y ) }{ \partial ( u, v ) } \right)^{-1} y_v \\ \displaystyle v_x = - \left( \frac{ \partial ( x, y ) }{ \partial ( u, v ) } \right)^{-1} y_u$

よって、一般には $u_x \neq x_u^{-1}$ である。

$x_v = 0$ 、すなわち、 $x \equiv x( u )$ のとき~~に限り~~を確認すると、
$\displaystyle u_x = \frac{ y_v }{ x_u y_v - x_v y_u } = \frac{ y_v }{ x_u y_v } = x_u^{-1}$
が成り立つことがわかる。

次回、もう少し複雑な場合を考える。

2019-01-09

微分の連鎖律：係数は前？後？

数学

微分の連鎖律
$\displaystyle x = x( u, v ), \quad y = y( u, v ) \\ \displaystyle \frac{d}{du} = \frac{d}{dx} \frac{dx}{du} + \frac{d}{dy} \frac{dy}{du} \\ \displaystyle \frac{d}{dv} = \frac{d}{dx} \frac{dx}{dv} + \frac{d}{dy} \frac{dy}{dv}$

例えば、
$\displaystyle x = u^2, \quad y = v^2 \\ \displaystyle \frac{dx}{du} = 2u, \quad \frac{dx}{dv} = 0 \\ \displaystyle \frac{dy}{du} = 0, \quad \frac{dy}{dv} = 2v$

この時、
$\displaystyle \frac{d}{du} = 2\sqrt{ x } \frac{d}{dx} \, or \, \frac{d}{dx} \, 2\sqrt{ x } \\ \displaystyle \frac{d}{dv} = 2\sqrt{ y } \frac{d}{dy} \, or \, \frac{d}{dy} \, 2\sqrt{ y }$
で若干迷ったりすることはないだろうか？
後者だと、もう一回微分しないといけなくなり、答えは同じにならない。

これを $f( x, y ) = x + y = u^2 + v^2$ の微分で確認してみると、
$\displaystyle \frac{df}{du} = 2u = 2u \frac{df}{dx} = \frac{dx}{du} \frac{df}{dx} \\ \displaystyle \frac{df}{dv} = 2v = 2v \frac{df}{dy} = \frac{dy}{dv} \frac{df}{dy}$

したがって、係数は微分の前に出しておくのが正解である。
公式を覚えてしまうのも良いが、自分で簡単な例を作って確認出来る方が安心感が強いと個人的には思う。

2019-01-05

Green関数の遅延条件（因果律）について

力学

以下のサイトで、一次元のGreen関数が分かり易くまとめられている。
slpr.sakura.ne.jp

$G(t,t')$ において $t \neq t'$ のところでは簡単に一般形が求まるし、更に遅延条件（因果律）を課すことで具体的な形が求まるというのは、なるほどと思った。
それによって、自分が以前求めたものと一致したものが得られていた。
koideforest.hatenadiary.com

他サイトの導出と自分の導出を比較して、気が付いたことがある。
自分の場合には、問題に合った定義域の下端が設定される（例えば $t \ge 0$ ）ことによって、Green関数に制限が加わり、結果的に遅延条件になった。
一方、他サイトの方法では、遅延条件は完全に人間の都合によるものであるとして導入し、定義域は基本的に無限空間内( $-\infty < t < \infty$ )を動くとしている。

もちろん、定義域の設定そのものが「人間の都合によるもの」でもあるのだが、Green関数に課される因果律よりも、むしろ定義域の方が問題設定として土台になっているのではないかと感じた。
そうでなければ、 $t_0$ で自由落下を始めるとするとき、 $t < t_0$ での物体が自由落下し始める前のことを考えなければならなくなってしまう。
その意味で、因果律そのものは、我々が普段自然に設定する初期条件と切っても切れない関係にあると思う。

2018-12-22

Crude adiabatic (Crude-Born-Oppenheimer) 近似

量子力学

以前にBorn-Oppenheimer近似（BO近似）及び断熱(adiabatic)近似について言及した。
koideforest.hatenadiary.com

ここでは、またちょっと微妙に違うCrude-BO近似やCrude adiabatic近似を紹介する。

前回で肝となっていたのは、パラメータに依存した演算子を定義したところである。
電子-原子核相互作用演算子を $\hat{V}_{en}$ とすると（ $\vec{r} = (\vec{r}_1, \vec{r}_2, \cdots)$ 、 $\vec{R} = (\vec{R}_1, \vec{R}_2, \cdots)$ ）、
$\displaystyle \hat{V}^{R_0}_{en} \neq \hat{V}_{en} \\ \left. \hat{V}^{R_0}_{en}( \vec{r} ) \right|_{R_0 = \vec{R} } = \hat{V}_{en}( \vec{r}, \vec{R} )$
つまり、パラメータ $R_0$ を上手く調節すると、元の演算子に一致する電子位置空間の演算子 $\hat{V}^{R_0}_{en}$ を定義した訳である。

この演算子を含む電子ハミルトニアンの固有状態 $\Phi^{R_0}_{e}( \vec{r} )$ で、任意の状態 $\Psi( \vec{r}, \vec{R} )$ を展開すると、
$\displaystyle \Psi( \vec{r}, \vec{R} ) = \sum_k C^{R_0}_{k}( \vec{R} ) \, \Phi^{R_0}_{e,k}( \vec{r} )$
この展開係数は、結局は原子核の波動関数そのもの $C^{R_0}_{k}( \vec{R} ) = \Phi^{R_0}_{n,k}( \vec{R} )$ であることがわかる。

前回（BO or adiabatic）では $R_0$ を $\vec{R}$ に一致するように選ぶことで、原子核に対する微分方程式を得た。
Crude-BO or Crude-adiabatic では、 $R_0$ を任意のベクトルに固定する（平衡位置に選ぶことがほとんどであろう）。
この場合、電子波動関数は $\vec{R}$ に依存していないので、原子核位置で微分してもゼロになる。
$\displaystyle \nabla_{R_i} \Phi^{R_0}_{e,k}( \vec{r} ) = 0$
そのため、BO近似やadiabatic近似で落とす項は元々含まれないことになる。

一方、電子-原子核相互作用は $\vec{R_0}$ でしか完全に取り込めないので、余剰分が出てくる。
$\displaystyle \Delta V^{R_0}( \vec{r}, \vec{R} ) \equiv V_{en}( \vec{r}, \vec{R} ) - V^{R_0}_{en}( \vec{r} ) = V_{en}( \vec{r}, \vec{R} ) - V_{en}( \vec{r}, \vec{R}_0 )$
この余剰分が、非対角項をもたらすことになる。

これらを元に、原子核に対する微分方程式を求めると、
$\displaystyle \sum_k \left( \left( -\sum_i \frac{\hbar^2}{2M_i} \nabla^2_{R_i} + E^{R_0}_{e,k} + V_{nn}(\vec{R}) \right) \delta_{k'k} + \Delta V^{R_0}_{k'k}( \vec{R} ) \right) \Phi^{R_0}_{n,k}( \vec{R} ) = E \Phi^{R_0}_{n,k'}( \vec{R} ) \\ \displaystyle \Delta V^{R_0}_{k'k}( \vec{R} ) \equiv \int d\vec{r} \, \left( \Phi^{R_0}_{e,k'}( \vec{r} ) \right)^* \, \Delta V^{R_0}( \vec{r}, \vec{R} ) \, \Phi^{R_0}_{e,k}( \vec{r} )$

そのため、Crude-BOでは $\Delta V^{R_0}_{k'k}( \vec{R} )$ を完全に無視したり、Crude adiabaticでは対角項のみを拾ったりする。

こうしてみると、この原子核に対する方程式は、透熱的な場合における一つの極限のような状態に対応していることがわかる。
koideforest.hatenadiary.com
通常のBOやadiabaticの意味で落とす項は、今の場合には全く含まれていないため、通常のBOやadiabaticによる影響はゼロと言える。
その代わりに、ポテンシャル残差による非対角項が生じていて、ここをどうするかで近似の度合いを調節出来る。

$\vec{R} \approx \vec{R}_0$ の時には良い近似となると容易に想像付くが、 $\vec{R}_0$ から離れると非対角項の影響が顕著になるため、平衡位置からの変位が大きい分子振動よりは、固体中のフォノンの扱いに適していると言える。

2018-12-21

透熱的電子基底とBorn-Oppenheimer近似

量子力学

前回は、断熱的電子基底（電子ハミルトニアンに対して対角）を用いたBorn-Oppenpheimer近似（BO近似）について解説した。
koideforest.hatenadiary.com
今回は、もう少し一般化した透熱的な場合を紹介する。

前回、ハミルトニアンを
$\displaystyle \hat{H} = \hat{H}_e + \hat{H}_n + \hat{V}_{en} \\ \hat{H}_e = \hat{T}_e + \hat{V}_{ee} \\ \hat{H}_n = \hat{T}_n + \hat{V}_{nn}$

と定義した。

ここで、 $\hat{H}_e$ の中のポテンシャル項を都合の良いように分ける（or ポテンシャルを付け足す）。
$\displaystyle \hat{H}_e = \hat{T}_e + \hat{V}^{(0)}_{e} + \hat{V}^{(1)}_{e} \equiv \hat{H}^{(0)}_e + \hat{V}^{(1)}_{e}$

$\vec{r} = ( \vec{r}_1, \vec{r}_2, \cdots )$ 、 $\vec{R} = ( \vec{R}_1, \vec{R}_2, \cdots )$ とし、 $R$ を原子核の位置をパラメータ化したものとすると、
$\displaystyle \hat{V}^R_{en} \neq \hat{V}_{en} \\ \displaystyle V^R_{en}( \vec{r} ) \equiv V_{en}( \vec{r}, \vec{R}) \\ \displaystyle \hat{H}^R_e \equiv \hat{H}_e + \hat{V}^R_{en} \\ \displaystyle \hat{H}^{(0),R}_e \equiv \hat{H}^{(0)}_e + \hat{V}^R_{en}$
と定義する。

前回の断熱的な場合(adiabatic)では、 $\hat{H}^R_e$ の固有状態で展開したが、透熱的な場合(diabatic)ではより恣意的な $\hat{H}^{(0),R}_e$ の基底で展開する。
$\displaystyle H^R_e( \vec{r} ) \Phi^{a,R}_{e,k}( \vec{r} ) = E^{a,R}_{e,k} \Phi^{a,R}_{e,k}( \vec{r} ) \\ \displaystyle H^{(0),R}_e( \vec{r} ) \Phi^{d,R}_{e,k}( \vec{r} ) = E^{d,R}_{e,k} \Phi^{d,R}_{e,k}( \vec{r} )$

よって、任意の状態 $\Psi$ は次のように展開出来る。
$\displaystyle \Psi^{a,R}_{e,k}( \vec{r}, \vec{R} ) = \sum_{k} C^{a, R}_k \Phi^{a,R}_{e,k}( \vec{r} ) = \sum_{k} C^{d, R}_k \Phi^{d,R}_{e,k}( \vec{r} )$

一見、ただ $a$ と $d$ が入れ替わっただけに見えるが、違いを強調しておくと、
$\displaystyle < \Phi^{a,R}_{e,k'} | \hat{H}^R_e | \Phi^{a,R}_{e,k} > = E^{a,R}_{e,k} \delta_{k'k} \\ \displaystyle < \Phi^{d,R}_{e,k'} | \hat{H}^R_e | \Phi^{d,R}_{e,k} > = E^{d,R}_{e,k} \delta_{k'k} + V^{(1),d,R}_{e,k'k}$
となり、透熱的な場合にはポテンシャルの非対角項が残る。

透熱的電子基底を用いた展開を、全ハミルトニアン $\hat{H}$ のシュレディンガー方程式に適用すると、
$\displaystyle \sum_k \left( T^d_{n,k'k}( \vec{R} ) + \left( E^{d}_{e,k}( \vec{R} ) + V_{nn}( \vec{R} ) \right) \delta_{k'k} + V^{(1),d}_{e,k'k}( \vec{R} ) \right) \Phi^{(k)}_n( \vec{R} ) = \Phi^{(k)}_n( \vec{R} ) \\ \displaystyle T^d_{n,k'k}( \vec{R} ) = -\sum_i \frac{\hbar^2}{2M_i} \nabla^2_{R_i} + T^{(1),d}_{n,k'k}( \vec{R} ) \nabla_{R_i} + T^{(2),d}_{n,k'k}( \vec{R} )$
$V^{(1),d}_{e,k'k}( \vec{R} )$ が前回には無かった項である。

新しい項が加わって、問題が複雑になっただけのように感じられるが、 $T^{(1)}$ や $T^{(2)}$ も基底の影響を受けている。
そのため、 $V^{(1),d}_{e,k'k}( \vec{R} )$ が導入される代わりに、 $|T^{(i),d}| < |T^{(i),a}|$ と出来るような基底を見つければ、BO近似（ $T^{(i)} = 0$ ）の影響を小さくする（精度を高める）ことが出来る。

どんな基底（もしくはポテンシャルの分け方）を選択すれば良いかは問題に依存するであろうが、基底を弄って近似を調節出来るのは面白いと思う。

2018-12-21

ポテンシャルのルジャンドル関数展開

量子力学プログラミング

「ポテンシャルの角運動量展開」の二次元版だと思って頂いて差し支えない。

$\displaystyle V( r, \theta ) = \sum_l V_l( r ) P_l( \cos\theta ) \\ \displaystyle V_l( r ) = \frac{ 2 l + 1 }{ 2 } \int^\pi_0 d\theta \, \sin\theta \, P_l( \cos\theta ) V( r, \theta )$

係数 $(2l+1)/2$ はルジャンドル関数の直交関係から出て来る。
koideforest.hatenadiary.com

上記の式を用いて、 $l_{max}=10$ まで展開してみる。

import numpy as np
from math import radians
from scipy.special import legendre
from scipy.integrate import simps
from scipy.interpolate import interp1d, interp2d
from matplotlib import pyplot as plt

def Plx( l, x ):
    # legendre( l ) >>> poly1d([ ... ])
    # poly1d([ ... ])( x ) >>> float
    return legendre( l )( x )

def pw_expand( func, l, r ):
    N = 300
    theta = np.linspace( 0., np.pi, N )
    coeff = ( 2 * l + 1 ) / 2.
    pw = []
    for rr in r:
        integrand = []
        for t in theta:
            x = rr * np.cos( t )
            y = rr * np.sin( t )
            integrand.append( coeff * np.sin( t ) * func( x, y ) * Plx( l, np.cos( t ) ) ) 
        print( np.array( integrand ) )
        pw.append( simps( np.array( integrand ), theta ) )
    return np.array( pw )

# For an interpolated function
def pw_expand_interp( func, l, r ):
    N = 300
    theta = np.linspace( 0., np.pi, N )
    coeff = ( 2 * l + 1 ) / 2.
    pw = []
    for rr in r:
        integrand = []
        for t in theta:
            x = rr * np.cos( t )
            y = rr * np.sin( t )
            integrand.append( coeff * np.sin( t ) * func( x, y )[0] * Plx( l, np.cos( t ) ) ) 
        pw.append( simps( np.array( integrand ), theta ) )
    return np.array( pw )

最初に、真四角の定数ポテンシャルを展開してみる。

N = 500
x = np.linspace( -2, 2, N )
y = np.linspace( -2, 2, N )
X, Y = np.meshgrid( x, y )

mat_v_const = np.zeros( ( N, N ) )
for ix, xx in enumerate( x ):
    for iy, yy in enumerate( y ):
        mat_v_const[ iy ][ ix ] = potential_const( xx, yy )

fig = plt.figure( figsize = ( 5, 4 ) )
ax = fig.add_subplot( 1, 1, 1 )
pcm = plt.pcolormesh( X, Y, mat_v_const, cmap = 'Greys' )
fig.colorbar( pcm )
plt.savefig( "v_const.png" )
plt.show()

f:id:koideforest:20181217080205p:plain — v_const

角運動量は、球対称の時に保存量となるため、球対称に近い時に展開の収束が良いと言える。
これは、フーリエ変換の時に正弦波に近いほど収束が良いのと同じ話である。
そのため、不連続な境界を円周上に持つ真四角ポテンシャルの収束は、最悪だと予想出来る。
実際、全然収束しない。

r = np.linspace( 0, 2 * np.sqrt(2), 200 )
pws = []
for l in range( 10 + 1 ):
    pws.append( pw_expand( potential_const, l, r ) )

r = np.linspace( 0, 2 * np.sqrt(2), 200 )
f_const = [ interp1d( r, pws[ l ], kind = "cubic" ) for l in range( 10 + 1 ) ]
mat_pw_const = np.zeros( ( N, N ) )
l_max = 10
for ix, xx in enumerate( x ):
    for iy, yy in enumerate( y ):
        rr = np.linalg.norm([ xx, yy ])
        tt = np.arctan2( yy, xx )
        for l in range( l_max + 1 ):
            if l % 2 == 0:
                mat_pw_const[ iy ][ ix ] += f_const[l]( rr ) * Plx( l, np.cos( tt ) )

fig = plt.figure( figsize = ( 5, 4 ) )
ax = fig.add_subplot( 1, 1, 1 )
pcm = ax.pcolormesh( X, Y, mat_pw_const, cmap = "bwr", vmin = -1, vmax = 1 )
fig.colorbar( pcm )
plt.savefig( "pw_const_plot.png" )
plt.show()

f:id:koideforest:20181217080304p:plain — pw_const_plot

元のポテンシャルとの差を取ると、展開が上手く行っていない部分がより分かり易くなる。

fig = plt.figure( figsize = ( 5, 4 ) )
ax = fig.add_subplot( 1, 1, 1 )
pcm = plt.pcolormesh( X, Y, mat_pw_const - mat_v_const, cmap = "bwr", vmin = -1, vmax = 1 )
fig.colorbar( pcm )
plt.savefig( "pw_const_plot_diff.png" )
plt.show()

f:id:koideforest:20181217080328p:plain — pw_const_diff_plot

各 $l$ の成分に分解した $V_l( r )$ を見てみると、 $l$ が大きくなっても寄与が消える気配がない。

for ipw, pw in enumerate( pws ):
    plt.plot( r, pw, label = "{:d}".format( ipw ) )
plt.legend()
plt.savefig( "pw_const.png" )
plt.show()

f:id:koideforest:20181217080356p:plain — pw_const

真四角ポテンシャルが隅のパリティ（空間反転に対して対称）を持つため、 $l$ が偶数の時だけ寄与する（ $l$ が奇数の時、パリティは奇）。

次に、境界は球対称だが、中身はランダムポテンシャルで対称性もクソもない場合。（後でポテンシャルを再利用するため、cubic splineで補間して関数として取っておく。）

import random

def potential_random( x, y ):
    if x**2 + y**2 <=1:
        return random.random()
    else:
        return 0.

mat_v_rand = np.zeros( ( N, N ) )
for ix, xx in enumerate( x ):
    for iy, yy in enumerate( y ):
        mat_v_rand[ iy ][ ix ] = potential_random( xx, yy )

f_mat_v = interp2d( x, y, mat_v_rand, kind = "cubic" )
mat_v_interp = np.zeros( ( N, N ) )
for ix, xx in enumerate( x ):
    for iy, yy in enumerate( y ):
        mat_v_interp[ iy ][ ix ] = f_mat_v( xx, yy )

fig = plt.figure( figsize = ( 5, 4 ) )
ax = fig.add_subplot( 1, 1, 1 )
pcm = plt.pcolormesh( X, Y, mat_v_interp, cmap = 'Greys' )
fig.colorbar( pcm )
plt.savefig( "v_rand.png" )
plt.show()

f:id:koideforest:20181217080418p:plain — v_rand

境界が丸いため、その部分は問題が無いが、やはり中身は苦労する。

pws_rand = []
for l in range( 10 + 1 ):
    pws_rand.append( pw_expand_interp( f_mat_v, l, r ) )

f_rand = [ interp1d( r, pws_rand[ l ], kind = "cubic" ) for l in range( 10 + 1 ) ]
mat_rand_pw = np.zeros( ( N, N ) )
l_max = 10
for ix, xx in enumerate( x ):
    for iy, yy in enumerate( y ):
        rr = np.linalg.norm([ xx, yy ])
        tt = np.arctan2( yy, xx )
        for l in range( l_max + 1 ):
            mat_rand_pw[ iy ][ ix ] += f_rand[l]( rr ) * Plx( l, np.cos( tt ) )

fig = plt.figure( figsize = ( 5, 4 ) )
ax = fig.add_subplot( 1, 1, 1 )
pcm = plt.pcolormesh( X, Y, mat_rand_pw, cmap = 'Greys', vmin = 0, vmax = 1 )
fig.colorbar( pcm )
plt.savefig( )
plt.show()

f:id:koideforest:20181217080443p:plain — pw_rand_plot

全然似ても似つかない。
差を取ると、何を再現出来ているのか不明なくらい差が大きい。

fig = plt.figure( figsize = ( 5, 4 ) )
ax = fig.add_subplot( 1, 1, 1 )
pcm = plt.pcolormesh( X, Y, mat_rand_pw - mat_v_interp, cmap = 'bwr', vmin = -1, vmax = 1 )
fig.colorbar( pcm )
plt.show()

f:id:koideforest:20181217080512p:plain — pw_rand_diff_plot

$V_l(r)$ もほぼ全ての $l$ が寄与を持つ。

for ipw, pw in enumerate( pws_rand ):
        plt.plot( r, pw, label = "{:d}".format( ipw ) )
plt.legend( loc = 1 )
plt.savefig( "pw_rand.png" )
plt.show()

f:id:koideforest:20181217080537p:plain — pw_rand

これらは、かなり非現実的なポテンシャルであるため、次回はもう少しマシなポテンシャルを展開してみたいと思う。

nano_exit

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

微分係数の逆数は逆の微分？（続き）

微分係数の逆数は逆の微分？

微分の連鎖律：係数は前？後？

Green関数の遅延条件（因果律）について

Crude adiabatic (Crude-Born-Oppenheimer) 近似

透熱的電子基底とBorn-Oppenheimer近似

ポテンシャルのルジャンドル関数展開