シュレーディンガー方程式

シュレーディンガー方程式

シュレーディンガー方程式は量子力学における基本的な方程式の一つであり、粒子の波動関数の時間変化を記述します。

しかし、シュレーディンガー方程式を正確に解くことは通常、解析的には不可能であるため、数値的な近似解を求めることが一般的です。

以下に、一次元の無限井戸ポテンシャル下でのシュレーディンガー方程式の時間依存性を無視した定常状態の解を求めるPythonのコードを示します。

このコードは、一次元の無限井戸の中を動く粒子の波動関数を計算し、それをプロットします。

無限井戸は物理学における理想化されたモデルで、井戸の中を動く粒子は井戸の壁を超えて脱出することはありません。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import diags
from scipy.sparse.linalg import eigs

# 定数
N = 1000  # 格子点の数
mass = 1.0  # 粒子の質量
hbar = 1.0  # 換算プランク定数
L = 1.0  # 井戸の幅
x = np.linspace(0, L, N)  # 位置

# 差分化
dx = x[1] - x[0]
T = hbar**2 / (2.0 * mass * dx**2) * (diags([1, -2, 1], [-1, 0, 1], shape=(N, N)).toarray())

# エネルギー固有値と波動関数の計算
values, vecs = eigs(T, k=30, which='SM')
idx = values.real.argsort()
values = values.real[idx]
vecs = vecs.real[:, idx]

# 波動関数の正規化とプロット
for n in range(3):  # 最初の3つのエネルギー状態をプロット
    psi = vecs[:, n] / np.sqrt(dx)
    plt.plot(x, psi, label='E = %.2f' % values[n])

plt.title('Wavefunctions in a 1D Infinite Square Well')
plt.xlabel('x')
plt.ylabel('psi')
plt.legend(loc='upper right')
plt.grid(True)

plt.show()

このコードでは最初の3つのエネルギー状態の波動関数をプロットします。

エネルギーが増えるにつれて波のノード（波の振幅がゼロとなる点）の数が増えることが確認できます。

これは量子力学における特從性の一つであり、粒子のエネルギー状態を視覚化することができます。

[実行結果]

ソースコード解説

このコードは、1次元の無限長方形井戸内の粒子の波動関数を計算し、その結果をグラフ化するものです。

1. ライブラリのインポート

1
2
3
4
import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import diags
from scipy.sparse.linalg import eigs
  • numpyは数値計算のためのライブラリで、配列操作や数学関数を提供します。
  • matplotlib.pyplotはグラフ描画のためのライブラリで、グラフの作成やカスタマイズを行います。
  • scipy.sparse.diagsは疎行列を扱うためのライブラリで、特定のパターンの疎行列を生成します。
  • scipy.sparse.linalg.eigsは疎行列の固有値問題を解くためのライブラリで、固有値と固有ベクトルを計算します。

2. 定数の設定

1
2
3
4
5
N = 1000  # 格子点の数
mass = 1.0  # 粒子の質量
hbar = 1.0  # 換算プランク定数
L = 1.0  # 井戸の幅
x = np.linspace(0, L, N)  # 位置
  • Nは格子点の数、massは粒子の質量、hbarは換算プランク定数、Lは井戸の幅を表します。
  • np.linspaceは$0$から L までの範囲を N 分割した値の配列を生成します。
    これにより、位置 x の配列が作成されます。

3. 差分化

1
2
dx = x[1] - x[0]
T = hbar**2 / (2.0 * mass * dx**2) * (diags([1, -2, 1], [-1, 0, 1], shape=(N, N)).toarray())
  • dxは隣接する位置の差を表します。
  • diags関数は特定のパターンの疎行列を生成します。
    ここでは、2階微分の差分演算子を表す疎行列$ (T) $を生成しています。

4. エネルギー固有値と波動関数の計算

1
2
3
4
values, vecs = eigs(T, k=30, which='SM')
idx = values.real.argsort()
values = values.real[idx]
vecs = vecs.real[:, idx]
  • eigs関数を使用して疎行列$ (T) $の固有値と対応する固有ベクトル（波動関数）を計算します。
  • k=30は最初の$30個$の固有値と固有ベクトルを取得することを意味します。
  • which='SM'は最小の固有値から計算を行うことを指定します。

5. 波動関数の正規化とプロット

1
2
3
4
5
6
7
8
9
10
11
for n in range(3):  # 最初の3つのエネルギー状態をプロット
    psi = vecs[:, n] / np.sqrt(dx)
    plt.plot(x, psi, label='E = %.2f' % values[n])

plt.title('Wavefunctions in a 1D Infinite Square Well')
plt.xlabel('x')
plt.ylabel('psi')
plt.legend(loc='upper right')
plt.grid(True)

plt.show()
  • 取得した固有ベクトル（波動関数）を正規化し、最初の3つのエネルギー状態に対応する波動関数をグラフ化しています。
  • plt.plotで各波動関数を位置$ (x) $に対する波動関数の値$ (\psi) $としてプロットしています。
  • plt.title、plt.xlabel、plt.ylabelはそれぞれグラフのタイトルと軸ラベルを設定し、plt.legendで凡例を表示しています。
  • plt.grid(True)はグリッド線を表示し、plt.show()でグラフを表示しています。

グラフ解説

[実行結果]

このグラフは、1次元の無限長方形井戸内の粒子の波動関数を計算し、その結果を表示したものです。

まず、N個の格子点で井戸の幅が L である領域を考えています。
その後、この領域を離散化するために、差分化を行っています。
具体的には、2階微分を表す行列$ (T) $を作成しています。
この行列は、離散的な位置$ (x)$ における運動エネルギーの演算子（運動量の二乗に比例する部分）を表します。

次に、eigs関数を使用して$ (T) $の固有値と対応する固有ベクトル（波動関数）を計算しています。
ここで k=30 は、最初の$30$個の固有値と固有ベクトルを取得することを意味します。
そして which='SM' は、最小の固有値から計算を行うことを指定しています。

最後に、最初の3つのエネルギー状態に対応する波動関数を取り出し、それぞれを正規化してグラフ化しています。
波動関数は各エネルギー状態における粒子の振る舞いを示し、横軸が位置$ (x) $、縦軸が波動関数$ (\psi) $の値を表しています。

タイトルは「Wavefunctions in a 1D Infinite Square Well」であり、各波動関数のラベルには対応するエネルギー値 $ (E) $が表示されており、それぞれのエネルギー状態が異なるラインで示されています。