nano_exit

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

多面体の内側に点があるかどうかを判定する。

多面体の頂点の座標が分かっているとする。
多面体の内側の点を原点に取り、対象とする点が多面体の内側に入っているかどうかを判断したい。

もし対象点が多面体の外側にあるとき、「多面体の頂点ベクトルを(原点から見た)対象点の方向に射影させると、射影成分は常に(原点からの)対象点の距離より小さい」が成り立つ。
逆に言えば、各頂点について調べていき、一つでも射影成分の方が大きければ、多面体の内側にあることが分かる。

訂正(2021/0812)
多面体の頂点との内積では上手く行かないことが分かった。
原点から各面に垂線を下し、面と垂線の交点と垂点と定義する。
原点から対象点へ向かうベクトルを各垂点方向へ射影したときに、射影成分が垂線の長さよりも短ければ、多面体の内側にいると判断できる。

以下、テストコード。

import numpy as np
import itertools as it

def make_vertical_points( vertices ):
    vps = []
    for v1, v2, v3 in it.combinations( vertices, 3 ):
        vv = np.cross( v2 - v1, v3 - v1 )
        ev = vv / np.linalg.norm( vv )
        pv = np.dot( v1, ev )
        nv = abs( pv )
        if np.all( nv >= np.array( [ np.dot( vertex, ev ) for vertex in vertices ] ) ) and \
           not np.any( [ np.all( pv * ev == vp ) for vp in vps ] ):
            vps.append( pv * ev )
    return vps

def check_inside_polyhedron( point, vps ):
    return np.all( [ np.dot( point, vp/np.linalg.norm(vp) ) <= np.linalg.norm( vp ) for vp in vps ] )

# example
vertices = np.array( list( it.product( [-1,1], repeat=3 ) ) )
vps = make_vertical_points( vertices )
check_inside_polyhedron( np.array( [0.9, 0, 0] ), vps )

pythonでLebedev quadratureを使う(更新)

以前にquadpyについて紹介した。
koideforest.hatenadiary.com
しかし、記事が古くなっていて、仕様がいろいろ変更されたようなので、更新する。
詳しい内容は公式のGitHubページを参照。
GitHub - nschloe/quadpy: Numerical integration (quadrature, cubature) in Python

テストコード。

import numpy as np
import quadpy

def f_cart( xyz ):
    # np.shape( xyz ) -> ( 3, number of points )
    n_points = np.shape( xyz )[1]
    return np.ones( n_points )
    # DO NOT WRITE "return 1"

def f_sphere( theta_phi ):
    # np.shape( theta_phi ) -> ( 2, number of points )
    # theta: polar : theta_phi[0]
    # phi: azimuth : theta_phi[1]
    n_points = np.size( theta_phi[0] )
    return np.ones( n_points )
    # DO NOT WRITE "return 1"

#scheme = quadpy.u3._lebedev.lebedev_019()
scheme = quadpy.u3.get_good_scheme(19)

# \int^\pi_0 d\theta \sin\theta \int^{2\pi}_0 d\phi = 4 \pi

xyz_center = ( 0., 0., 0. )
radius = 1.
val_cart = scheme.integrate( f_cart , xyz_center, radius )

val_sphere  = scheme.integrate_spherical( f_sphere )

print( val_cart, val_sphere, 4*np.pi )

エネルギー等分配則を並進・回転・振動に分けずに理解する。

エネルギー等分配則の本質は、

\displaystyle
H \equiv \sum_i h_i = \sum_i c_i \xi_i^2
\\
\displaystyle
\varepsilon_i = \left< h_i \right> = \frac{1}{2\beta}
である。
期待値の計算は、ここでは省略する。

もう少し具体的にハミルトニアンを表すと,

\displaystyle
H = \sum_n \frac{\vec p_n \cdot \vec p_n}{2m_n} + \sum_{n \neq n'} \frac{\kappa_{n,n'}}{2}( \vec r_n - \vec r_{n'} )^2
第二項はいわゆる「振動」を与える調和ポテンシャルである。
これがエネルギー等分配則を語る上で、本来最初に宣言するべきハミルトニアンである。
物理化学の教科書では、運動を並進・回転・振動に自由度を分けて、分子が複雑になって行くにつれて少しずつ自由度を足していくような論法が見受けられるが、これでは何が前提になっているのかサッパリ分からない。
エネルギー等分配則自体は、もうこのハミルトニアンから大体計算出来て、

  • 単原子分子: H = \frac{p_x^2 + p_y^2 + p_z^2 }{2m} \beta \varepsilon= 3/2
  • 二原子分子: H = \sum_n \frac{p_{x,n}^2 + p_{y,n}^2 + p_{z,n}^2 }{2m_n} + \frac{\kappa}{2}( \vec r_1- \vec r_2 )^2 \beta \varepsilon= (3/2 + 3/2) + 1/2 = 6/2 + 1/2
  • 直線三原子分子: H = \sum_n \frac{p_{x,n}^2 + p_{y,n}^2 + p_{z,n}^2 }{2m_n} + \sum_{n,n'} \frac{\kappa_{n,n'}}{2}( \vec r_n- \vec r_{n'} )^2 + V_{linear} \beta \varepsilon= 9/2 + 4/2
  • 非直線三原子分子: H = \sum_n \frac{p_{x,n}^2 + p_{y,n}^2 + p_{z,n}^2 }{2m_n} + \sum_{n,n'} \frac{\kappa_{n,n'}}{2}( \vec r_n- \vec r_{n'} )^2 \beta \varepsilon= 9/2 + 3/2

ここで、 V_{linear}は原子を直線状に拘束するための追加のポテンシャルである。その意味で、直線分子ではそもそものポテンシャルの作り方に気を付けなければならないかも知れない。
いずれにせよ、このように直線分子には拘束ポテンシャルが必要と考えると、物理化学で言われる並進・回転・振動に分解する方法と一致する。
これらを変数変換していき、適切な近似(例えば、剛体回転時の慣性半径に振動の効果を考慮しない等)を施せば、並進・回転・振動の表記にできる。

調和ポテンシャルが振動に寄与するのは明らかなので、ハミルトニアンから始める方が「何故振動の自由度だけ \beta \varepsilon_{vib} = 1 \neq 1/2なのか?」を説明し易い。

ちなみに、二原子分子のエネルギー等分配則で \beta \varepsilon = 5/2と説明されるが、それは実験に合わせるために振動の自由度を凍結させるからであり、「全ての自由度に等しいエネルギーが分配される」という真の意味のエネルギー等分配則からすると近似表現である。

さらに、この方法を一般化した場合、例えば非直線五原子分子だと、 \beta \varepsilon= 5 \times 3/2 + {}_5C_2 \times 1/2 = 25/2となるが、並進・回転・振動に分解する方法だと 3/2 + 3/2 + (15 - 6) = 24/2となる。
これは、並進・回転・振動に分解したときに、どこかのポテンシャル(おそらく一番遠い原子間のポテンシャル)を無視していることに対応しており、並進・回転・振動に分解する方法がやはり近似的な手法であることがわかる。

とはいえ、例えば「室温付近では振動の自由度は考えないことで実験を説明できる」等、並進・回転・振動の分解は有益な情報を齎すし、何より一般的なハミルトニアンが与える実際の運動を理解するために結局は並進・回転・振動に分解するので、最終的には必要である。

参考文献
http://phys.sci.hokudai.ac.jp/~kita/StatisticalMechanicsI/Stat6.pdf
エネルギー等配分の法則 - Wikipedia

カノニカル相関

ハミルトニアン Hが摂動部分 Vを持つとし、統計平均を以下で定める。

\displaystyle
H = H_0 + V
\\
\displaystyle
Z \equiv {\rm tr}[ e^{- \beta H } ] 
\\
\displaystyle
\left< A \right>
  \equiv \frac{ {\rm tr}[ e^{- \beta H } A ] }{ Z }
\\
\displaystyle
Z_0 \equiv {\rm tr}[ e^{- \beta H_0 } ] 
\\
\displaystyle
\left< A \right>_0
  \equiv \frac{ {\rm tr}[ e^{- \beta H_0 } A ] }{ Z_0 }

このとき、カノニカル相関は次のように定義される。

\displaystyle
\left< A ; B \right>
  \equiv \int_0^\beta \frac{d \lambda}{\beta} \left< e^{\lambda H_0} A e^{-\lambda H_0} B \right>_0
非摂動状態による統計平均で書かれているところがポイントである。

一次摂動における等温感受率は、カノニカル相関を用いて表現できる。
まず、ハミルトニアンがパラメータ xに依存する摂動を持つとすると、

\displaystyle
H( x )
  \equiv H_0( x_0 ) + V( x - x_0 ) = H_0 + V( x )
\\
\displaystyle
Z(x) \equiv {\rm tr}[ e^{- \beta H(x) } ] 
\\
\displaystyle
\left< A \right>_x
  \equiv \frac{ {\rm tr}[ e^{- \beta H(x) } A ] }{ Z(x) }
等温感受率は、統計平均によって以下の様に書ける。

\displaystyle
\chi_T( x ) = \frac{\partial}{\partial x} \left< A \right>_x

厳密に計算すると、

\displaystyle
\frac{\partial}{\partial x} Z^{-1}(x)
  = Z^{-2}(x) {\rm tr}\left[ - \beta e^{- \beta H(x)} \frac{\partial H(x) }{\partial x} \right]
\\
\displaystyle
\quad
  = Z^{-2}(x) {\rm tr}\left[ - \beta e^{- \beta H(x)} \frac{\partial V(x) }{\partial x} \right]
\\
\displaystyle
\quad
  \equiv Z^{-2}(x) {\rm tr}[ \beta e^{- \beta H(x)} R(x) ]
  = \beta Z^{-1}(x) \left< R(x) \right>
\\
\displaystyle
R(x) = - \frac{\partial V(x) }{\partial x}
\\
\displaystyle
\frac{\partial}{\partial x} {\rm tr}\left[ e^{- \beta H(x)} A \right]
  = \beta {\rm tr}\left[ e^{- \beta H(x)} R(x) A \right]
\\
\displaystyle
\therefore
\frac{\partial}{\partial x} \left< A \right>_x
  = \beta \left< R(x) A \right>_x - \beta \left< R(x) \right>_x \left< A \right>_x
\\
\displaystyle
\quad
  = \beta \left< (R(x) - \left< R(x) \right>_x) (A - \left< A \right>_x ) \right>_x
\\
\displaystyle
\quad
  \equiv \beta \left< \Delta_x R(x) \Delta_x A \right>_x
\\
\displaystyle
\because
\left< A B \right> - \left< A \right> \left< B \right> = ( A - \left< A \right>)( B - \left< B \right>)
\\
\displaystyle
\Delta_x A = A - \left< A \right>_x
しかし、fullに摂動を考慮した期待値を計算できないから困っているわけで、ここで行き詰ってしまう。

摂動展開を用いて、 Vの一次まで拾うようにしていく。
koideforest.hatenadiary.com

\displaystyle
e^{-\beta(H_0+V)}
  = e^{-\beta H_0} \left( 1 - \int_0^\beta d\lambda\, e^{\lambda H_0} V e^{- \lambda (H_0 + V)} \right)
\\
\displaystyle
\quad
  \approx e^{-\beta H_0} \left( 1 - \int_0^\beta d\lambda\, e^{\lambda H_0} V e^{- \lambda H_0} \right)
\\
\displaystyle
Z(x)
  \approx Z_0 - \int_0^\beta d\lambda\, Z_0 \left< V(x) \right>_0
  = Z_0 \left( 1 - \int_0^\beta d\lambda\, \left< V(x) \right>_0 \right)
\\
\displaystyle
\because
{\rm tr }[ e^{- \beta H_0} e^{\lambda H_0} V(x) e^{- \lambda H_0} ]
  = {\rm tr }[ e^{- \beta H_0} V(x) ] = Z_0 \left< V(x) \right>_0
\\
\displaystyle
Z^{-1}(x)
  \approx \left[ Z_0 \left( 1 - \int_0^\beta d\lambda\, \left< V(x) \right>_0 \right) \right]^{-1}
  \approx Z^{-1}_0 \left( 1 + \int_0^\beta d\lambda\, \left< V(x) \right>_0 \right)
\\
\displaystyle
\therefore
\frac{ e^{- \beta ( H_0 + V(x) ) } }{ Z(x) }
  \approx \frac{ e^{- \beta H_0 } }{ Z_0 } \left( 1 - \int_0^\beta d\lambda\, e^{\lambda H_0} V(x) e^{- \lambda H_0} \right) \left( 1 + \int_0^\beta d\lambda\, \left< V(x) \right>_0 \right)
\\
\displaystyle
\quad
  \approx \frac{ e^{- \beta H_0 } }{ Z_0 } \left( 1 - \int_0^\beta d\lambda\, \left[ e^{\lambda H_0} V(x) e^{- \lambda H_0} - \left< V(x) \right>_0 \right] \right)
\\
\displaystyle
\therefore
\left< A \right>_x
  \approx  \left< A \right>_0 - \int_0^\beta d\lambda\, \left[ \left< e^{\lambda H_0} V(x) e^{- \lambda H_0} A \right>_0 - \left< V(x) \right>_0 \left< A \right>_0 \right]
\\
\displaystyle
\quad
  =  \left< A \right>_0 - \int_0^\beta d\lambda\, \left< ( e^{\lambda H_0} V(x) e^{- \lambda H_0} - \left< V(x) \right>_0 ) ( A - \left< A \right>_0 ) \right>_0
\\
\displaystyle
\quad
  =  \left< A \right>_0 - \int_0^\beta d\lambda\, \left< e^{\lambda H_0} ( V(x) - \left< V(x) \right>_0 ) e^{- \lambda H_0} ( A - \left< A \right>_0 ) \right>_0
したがって、

\displaystyle
\chi_T (x) = \frac{\partial}{\partial x} \left< A \right>_x
  \approx \int_0^\beta d\lambda\, \left< e^{\lambda H_0} ( R(x) - \left< R(x) \right>_0 ) e^{- \lambda H_0} ( A - \left< A \right>_0 ) \right>_0
\\
\displaystyle
\quad
  = \beta \left< \Delta_0 R(x); \Delta_0 A \right>
\\
\displaystyle
\Delta_0 A = A - \left< A \right>_0

ブール領域を用いたエンタングルメントの説明。

ブール領域  \mathbb{B}:=\{0, 1\} の直積の部分集合として、以下のものを考える。

\mathcal{C} := \{ (0, 0), (1, 1) \} \subset \mathbb{B} \times \mathbb{B}

これは、次の様に書くことが出来ない。

\mathcal{C} \neq \{ (b, b') | b \in \mathbb{B}, b' \in \mathbb{B} \} = \mathbb{B} \times \mathbb{B}
これを、ここでは「分離不可能」と言うことにする。

一方で、次のようなものは分離可能である。

\mathcal{C}_1 := \{ (0, 0) \}
\\
\mathcal{C}_2 := \{ (1, 1) \}
\\
\mathcal{C}_3 := \{ (0, 1) \}
\\
\mathcal{C}_4 := \{ (1, 0) \}
\\
\mathcal{C}_5 := \{ ( 0, 0 ), ( 0, 1 ), ( 1, 0 ), ( 1, 1 ) \} = \mathbb{B} \times \mathbb{B}

よって、部分集合を取れば何でもよい訳ではなく、 \mathcal{C} は特殊な状態であることが分かる。

参考文献:圏論量子力学入門

追記(2021/03/28)
上の定義で言うと、分離不可能な部分集合は結構ある。

\mathcal{C}^c := \{ (0, 1), (1,0) \}
\\
\mathcal{C}_1^c := \{ (0, 1), (1,0), (1, 1) \}
\\
\mathcal{C}_2^c := \{ (0, 0), (0, 1), (1, 0) \}
\\
\mathcal{C}_3^c := \{ (0, 0), (1, 0), (1, 1) \}
\\
\mathcal{C}_4^c := \{ (0, 0), (0, 1), (1, 1) \}
水素分子的に言えば、 \mathcal{C}はイオン状態で、 \mathcal{C}^cは原子状態に対応させることが出来る。
なので、もつれている部分集合の中でも、性質の良いものだけを抜き出して何か数学的な構造が出来ると面白そう。
こうしてみると、量子論だけでなく、もっと普遍的にエンタングルメントが存在して何かしらの機能を果たしているような気もする。

上極限集合と下極限集合。

個人的には、数列を作った方が分かり易い。

  •  \limsup


\displaystyle
\lim_{m \to \infty} S_m \equiv \limsup_{n \to \infty} A_n = \bigcap^\infty_{k=1} \bigcup^\infty_{n=k} A_n
\\
\displaystyle
S_m \equiv \bigcap^{m}_{k=1} \bigcup^\infty_{n=k} A_n
\\
\displaystyle
S_1 = \bigcup^\infty_{n=1} A_n
\\
\displaystyle
S_{l+1} = S_l \cap \left( \bigcup^\infty_{n=l+1} A_n \right) \subset S_l
よって、数列 \{S_m\}は単調減少数列である。
また、 \bigcup^\infty_{n=l+1} A_n \supset A_{\infty}より、 \lim_{m \to \infty} S_m \supset A_{\infty}

  •  \liminf


\displaystyle
\lim_{m \to \infty} I_m \equiv  \liminf_{n \to \infty} A_n = \bigcup^\infty_{k=1} \bigcap^\infty_{n=k} A_n
\\
\displaystyle
I_m \equiv \bigcup^{m}_{k=1} \bigcap^\infty_{n=k} A_n
\\
\displaystyle
I_1 = \bigcap^\infty_{n=1} A_n
\\
\displaystyle
I_{l+1} = I_l \cup \left( \bigcup^\infty_{n=2} A_n \right) \supset I_l
よって、数列 \{I_m\}は単調増加数列である。
また、 \bigcap^\infty_{n=l+1} A_n \subset A_{\infty}より、 \lim_{m \to \infty} I_m \subset A_{\infty}


したがって、 \liminf_{n \to \infty} A_n \subset \limsup_{n \to \infty} A_n が自然に言える。


実際に計算するときにも、数列を考えた方が扱い易いと思う。

例: A_n = [ 1/n, 1 ]


\displaystyle
\forall k \ge 1, \bigcup^\infty_{n=k}  [ 1/n, 1 ] = ( 0, 1]
\\
\therefore
\displaystyle
\forall m \ge1, S_m = ( 0, 1]
\\
 \limsup_{n \to \infty} A_n = \lim_{m \to \infty} S_m= ( 0, 1 ]


\displaystyle
I_1 = \bigcap^\infty_{n=1}  [ 1/n, 1 ] = [1]
\\
\displaystyle
I_2 = I_1 \cup \left( \bigcap^\infty_{n=2}  [ 1/n, 1 ] \right) = [ 1 ] \cup [1/2, 1] = [1/2, 1]
\\
\displaystyle
\therefore
I_m =  [1/m, 1]
\\
 \liminf_{n \to \infty} A_n = \lim_{m \to \infty} I_m= ( 0, 1 ]


参考サイト:
極限集合[数学についてのwebノート]
極限集合の解釈(イプシロンデルタ風に) - 岡竜之介のブログ

時間反転演算子のエルミート共役について。

時間反転演算子 \Thetaは反ユニタリー演算子である。
反ユニタリー演算子は、ユニタリー演算子 Uと反線形演算子 Kの積で表される。
教科書によって、「時間反転演算子のエルミート共役は取ってはいけない」と書いてあったり、普通にエルミート共役 \Theta^\daggerが定義されていたりする。
ここでは、 \Theta^\daggerを何故考えてはいけないかを考察する。

エルミート共役をわざわざ取りたいのは、以下のような関係を利用したいからである。

\displaystyle
\langle \psi | O | \phi \rangle
  = \langle \psi | O \phi \rangle
  = \langle O^{\dagger} \psi | \phi \rangle

内積は、適当に基底を選べば具体的に計算できる。

\displaystyle
\langle \psi | \phi \rangle
  = \sum_a \langle \psi | a \rangle \langle a | \phi \rangle
  = \sum_a \int \psi^*(a) \phi(a)

見た目的には、 \psi \phiは対等な関係のように見える。
しかし、数学的には、ケットベクトルはブラベクトルを入力して内積を返す線形汎関数で定義されている。
多分 F_{\psi} \equiv \langle \psi |と書き換えると分かり易いかも。

\displaystyle
F_{\psi}[ | \phi \rangle ] \equiv \langle \psi | \phi \rangle \in \mathbb{C}

ここで、 \phiが線形結合で表されているとすると、

\displaystyle
F_{\psi}[ c_1 | \phi_1 \rangle + c_2 | \phi_2 \rangle ]
\\
\displaystyle
\quad
  = c_1 F_{\psi}[ | \phi_1 \rangle ]
  + c_2 F_{\psi}[ | \phi_2 \rangle ]
\\
\displaystyle
\quad
  = c_1 \langle \psi | \phi_1 \rangle
  + c_2 \langle \psi | \phi_2 \rangle
線形汎関数の感じが伝わると思う。

ここで、反ユニタリー演算子 \phiに掛けると、

\displaystyle
F_{ \psi }[ UK ( c_1 | \phi_1 \rangle + c_2 | \phi_2 \rangle ) ]
\\
\quad
\displaystyle
  =F_{ \psi }[ U ( c^*_1 | \phi_1 \rangle + c^*_2 | \phi_2 \rangle ) ]
\\
\displaystyle
\quad
  = c^{*}_1 \langle \psi | U | \phi_1 \rangle
  + c^{*}_2 \langle \psi | U | \phi_2 \rangle
上では \phi_{1,2}が基底ベクトルと仮定した( K| \phi_{1,2} \rangle = | \phi_{1,2} \rangle )。
基底ベクトルでなくても、ユニタリー変換でいつでも基底の線形結合で書けるので、一般性は失われていない。

それで、今度はブラベクトルに作用させたときに上記と同じものを与えるにはどうしたら良いを考える。

\displaystyle
F_{ (UK)^\dagger \psi }[ c_1 | \phi_1 \rangle + c_2 | \phi_2 \rangle ]
\\
\displaystyle
\quad
  = c_1 F_{ (UK)^\dagger \psi}[ | \phi_1 \rangle ]
  + c_2 F_{ (UK)^\dagger \psi}[ | \phi_2 \rangle ]
だがしかし、どう足掻いても \psi (UK)^\daggerを作用させるだけでは c^{*}_{1,2}は作れない。

実際には、「全体の複素共役を取る」などをすると対応出来るは出来るのだが、それは元々の線形汎関数の能力を超えているので、上手く定義出来ないのである。
なので、 \langle \psi | \Theta | \phi \rangleは常に \phiに掛けるようにして、 \psiには掛けないからエルミート共役 \Theta^\daggerはわざわざ定義しない、という形を取るのである。

参考文献:J. J. サクライ「現代の量子力学・下」