井戸型 vs クーロン

今更ながらIPythonデータサイエンスクックブックを購入した。もちろん私費である。

せっかくなので一次元のシュレーディンガー方程式でも解こうかと思った。
一階微分方程式に対しては、かくあきさんのサイトで紹介されているローレンツアトラクターのソースが簡潔で分かり易くて非常に参考になった。
Python NumPy SciPy : 1 階常微分方程式の解法 | org-技術

Hartree unitで、

$\phi''(x) = 2(E-V(x))\phi(x)$

井戸型ポテンシャルでは、

$V(x) = - 1 ( 0 < x < 1 )$
$V(x) = 0 ( otherwise )$

クーロンポテンシャルではゼロ割を防ぐために $10^{-7}$ を分母に足して、

$V(x) = - 1 / ( x + 10^{-7} )$

初期条件として、 $\phi(0) = 0, \phi'(x) = 1$ を用いた。~~理由は特にない。それっぽいだけ。~~実はいろいろ理由があったが、微妙に長くなるので割愛。軽く言及すると、井戸型は純粋に一次元の問題だが、クーロンポテンシャルは三次元の問題を球対称として動径部分のみの方程式に落とし込んでいて、 $\phi(x)= x \psi(x)$ を求めていることになる。

エネルギーを振って計算。

結果：
f:id:koideforest:20170122004507p:plain
f:id:koideforest:20170122004031p:plain

波動関数の大きさは規格化していないので置いといて、やはりクーロンポテンシャルの波の引きずり込み具合は素晴らしい。
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()

nano_exit

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

井戸型 vs クーロン