🍻

Juliaで線形代数：線形変分法

2022/01/29に公開

Julia

量子力学の支配方程式はSchrödinger方程式という微分方程式である. 一見, 線形代数とは関係無いように思われるが, 適当な基底関数で展開することによって, 行列の固有値問題に帰着できる. この記事では, その方法の一つである線形変分法について, Juliaによる具体的な計算例を挙げながら解説する.

パッケージ
線形変分法（Rayleigh-Ritzの変分法）
無限に深い井戸型ポテンシャル
試行波動関数と基底関数
行列要素の導出
変分計算
おわりに
参考文献

パッケージ

次の3つのパッケージを利用する. Printf.jlとLinearAlgebra.jlは標準パッケージであるため, 初めからインストールされている. Plots.jlはこちらの記事に従ってインストールしておこう.

パッケージの読み込み

# using Pkg
# Pkg.add("Plots")
using Plots
using Printf
using LinearAlgebra

線形変分法（Rayleigh-Ritzの変分法）

時間に依存しないSchrödinger方程式

\hat{H}\psi(\boldsymbol{r}) = E\psi(\boldsymbol{r})

を非直交な基底関数 $\phi_1, \phi_2, \phi_3, \cdots$ で展開すると, 行列の一般化固有値問題

\boldsymbol{H}\boldsymbol{C} = E\boldsymbol{S}\boldsymbol{C}

に帰着できることが知られている（証明は参考文献に譲る）. ただし, $\boldsymbol{H}$ はハミルトニアン行列といい, その要素は $H_{ij} = \langle\phi_i|\hat{H}|\phi_j\rangle$ である. $\boldsymbol{S}$ は重なり行列といい, その要素は $S_{ij} = \langle\phi_i|\phi_j\rangle$ である. $\boldsymbol{C}$ は固有ベクトル(縦ベクトル)を横に並べた行列であり, $E$ は固有値を並べた縦ベクトルである. Juliaのeigen()は $(\boldsymbol{H},\boldsymbol{S})\mapsto(\boldsymbol{C},E)$ , つまり入力として $\boldsymbol{H}$ と $\boldsymbol{S}$ を渡すと, 出力として $\boldsymbol{C}$ と $E$ が得られる仕様である. 一般固有値問題を標準固有値問題に変形する方法も, 標準固有値問題を解く方法も既に多く考案され実装されているから, 我々に残された仕事は $H_{ij}$ や $S_{ij}$ の具体的な値を求めることである.

無限に深い井戸型ポテンシャル

V(x) = \left\{ \begin{array}{ll} \infty & |x|\gt a \\ 0 & |x|\leq a \end{array} \right.

の中で運動する質量 $m$ の粒子のハミルトニアンは

\hat{H} = \left[ - \frac{\hbar^2}{2m} \frac{\mathrm{d}^2}{\mathrm{d}x ^2} + V(x) \right]

であり, これに対する時間に依存しないSchrödinger方程式

\hat{H}\psi(x) = E \psi(x)

の解は

\psi_n(x) = \sqrt{\frac{1}{a}} \sin \frac{n\pi(x+a)}{2a} \\ E_n = \frac{\hbar^2 n^2 \pi^2}{8 m a^2}

である(固有関数であることの証明と固有値の導出, 規格化定数の導出). なお, 井戸の外 $x<-a,a<x$ ではポテンシャル $V(x)$ が発散しているから, Schrödinger方程式はそもそも意味を持たない. 定義域を井戸の中 $-a\leq x \leq a$ に限り, 自由粒子のSchrödinger方程式を考え, 境界条件 $\psi_n(-a) = \psi_n(a) = 0$ を付加するものだと考えればよい.

厳密解

# V(x; L=1) = -L<x ? (x<L ? 0 : Inf) : Inf
ψ(n, x; a=1) = sqrt(1/a) * sin(n*π*(x+a)/(2*a))
E(n; ℏ=1, m=0.5, a=1) = (ℏ^2*n^2*π^2) / (8*m*a^2)

厳密解の描写

# using Plots
a = 1
plot(xlabel="\$x\$", ylabel="\$V(x)\$", xlims=(-3,3))
plot!([-a,-a,a,a], [68,0,0,68], label="\$V(x)\$", lc="#000000", lw=2)
for i in 1:5
    plot!([-a, a], fill(E(i,a=a), 2), label="\$E_$i\$", lc=i, lw=2)
    plot!(-a:0.01:a, x->3*ψ(i,x,a=a)+E(i,a=a), label="", lc=i, lw=2, ls=:dash)
end
plot!()

基底関数

J.M.ティッセン著『計算物理学』(丸善出版, 2012) 3.2.1を参考に, 試行波動関数 $\psi(x)$ と基底関数 $\phi_i(x)$ をそれぞれ

\begin{align} \psi(x) &\simeq \sum_{i=1}^{N} c_i \phi_i(x) \\ \phi_i(x) &= - x^{i-1} (x-a) (x+a) \end{align}

とする(基底関数は変えているので注意). この基底関数は厳密な固有関数と概ね似た形をしている.

基底関数

φ(n, x; a=1) = -x^(n-1) * (x-a) * (x+a)

基底関数の描写

# using Plots

a = 1

plot(xlabel="x", ylabel="\$\\psi_n(x), V(x)\$", xlims=(-1.2,1.2), ylims=(-1.1,1.1))
plot!([-a,-a,a,a], [100,0,0,100], label="", lc="#000000", lw=2)
for i in 1:5
    plot!(-a:0.01:a, x->ψ(i,x,a=a), label="", lc=i)
end
plt1 = plot!()

plot(xlabel="x", ylabel="\$\\phi_n(x), V(x)\$", xlims=(-1.2,1.2), ylims=(-1.1,1.1), legend=:bottomright)
plot!([-a,-a,a,a], [100,0,0,100], label="Potential", lc="#000000", lw=2)
for i in 1:5
    plot!(-a:0.01:a, x->φ(i,x,a=a), label="n = $i", lc=i)
end
plt2 = plot!()

plot(plt1, plt2)

行列要素

前述の通り, 基底関数 $\phi_i(x)$ を

\phi_i(x) = - x^{i-1} (x-a) (x+a)

とすると, 後に示すように, 行列要素 $H_{ij}$ , $S_{ij}$ はそれぞれ

\begin{align} S_{ij} &= \langle\phi_i|\phi_j\rangle = a^3 \left[a^{i+j} + (-a)^{i+j}\right] \left[ \frac{1}{i+j+3} - \frac{2}{i+j+1} + \frac{1}{i+j-1} \right] \\ H_{ij} &= \langle\phi_i|\hat{H}|\phi_j\rangle = - \frac{\hbar^2}{2m} a \left[a^{i+j}+(-a)^{i+j}\right] \left[ \frac{j^2+j}{i+j+1} - 2 \frac{j^2-j+1}{i+j-1} + \frac{j^2-3j+2}{i+j-3} \right] \\ \end{align}

と求められる（J.M.ティッセン著『計算物理学』(丸善出版, 2012) 3.2.1を参考にしているものの, $i=0,1,2,\cdots$ ではなく $i=1,2,3,\cdots$ としており, 基底関数も少し変えており, さらに無次元化しないまま行列要素を計算しているため, 表式は異なる）. これらを用いて, Schrödinger方程式は一般化固有値問題

\boldsymbol{H}\boldsymbol{C} = E\boldsymbol{S}\boldsymbol{C}

に帰着できる.

重なり行列の導出

まず, こちらの積分

\int_{-a}^{a} x^n \mathrm{d}x = \left[ \frac{x^{n+1}}{n+1} \right]_{-a}^{a} = \frac{a[a^n+(-a)^n]}{n+1}

を用いて, 重なり積分を計算する.

\begin{align} S_{ij} &= \langle\phi_i|\phi_j\rangle \\ &= \int_{-a}^{a} \phi_i^\ast(x) \phi_j(x) \mathrm{d}x \\ &= \int_{-a}^{a} \left[ -x^{i-1} (x-a) (x+a) \right] \left[ -x^{j-1} (x-a) (x+a) \right] \mathrm{d}x \\ &= \int_{-a}^{a} x^{i-1} x^{j-1} [(x-a) (x+a)]^2 \mathrm{d}x \\ &= \int_{-a}^{a} x^{i+j-2} (x^2-a^2)^2 \mathrm{d}x \\ &= \int_{-a}^{a} x^{i+j-2} (x^4 -2 a^2 x^2 + a^4) \mathrm{d}x \\ &= \int_{-a}^{a} x^{i+j+2} -2 a^2 x^{i+j} + a^4 x^{i+j-2} \mathrm{d}x \\ &= \frac{\left[x^{i+j+3}\right]_{-a}^{a}}{i+j+2+1} -2a^2 \frac{\left[x^{i+j+1}\right]_{-a}^{a}}{i+j+1} +a^4 \frac{\left[x^{i+j-1}\right]_{-a}^{a}}{i+j-2+1} \\ &= \frac{a\left[a^{i+j+2} + (-a)^{i+j+2}\right]}{i+j+3} - 2a^2 \frac{a\left[a^{i+j+0} + (-a)^{i+j+0}\right]}{i+j+1} + a^4 \frac{a\left[a^{i+j-2} + (-a)^{i+j-2}\right]}{i+j-1} \\ &= \frac{a^3 [a^{i+j} + (-a)^{i+j}]}{i+j+3} - 2a^2 \frac{a^1 [a^{i+j} + (-a)^{i+j}]}{i+j+1} + a^4 \frac{a^{-1}[a^{i+j} + (-a)^{i+j}]}{i+j-1} \\ &= a^3 \left[a^{i+j} + (-a)^{i+j}\right] \left[ \frac{1}{i+j+3} - \frac{2}{i+j+1} + \frac{1}{i+j-1} \right] \end{align}

この結果を数値積分によって検証する. 概ね一致しているので問題ない.

重なり積分の検証

# using Printf

int(func, x_min, x_max, dx) = sum(func.(x_min:dx:x_max))*dx

a = 1.5
N = 5
S = zeros(N,N)

println("  i  j\tnumerical\tanalytical")
for j in 1:N
    for i in 1:N
        numerical = integral(x->φ(i,x,a=a)*φ(j,x,a=a), -a, a, 0.001)
        S[i,j] = isodd(i+j) ? 0 : a^3 * (a^(i+j)+(-a)^(i+j)) * (1/(i+j+3) - 2/(i+j+1) + 1/(i+j-1))
        @printf("%3d%3d\t%.9f\t%.9f\n", i, j, numerical, S[i,j])
    end
end

出力

i  j	numerical	analytical
1  1	8.100000000	8.100000000
2  1	0.000000000	0.000000000
3  1	2.603571429	2.603571429
4  1	0.000000000	0.000000000
5  1	1.952678571	1.952678571
1  2	0.000000000	0.000000000
2  2	2.603571429	2.603571429
3  2	0.000000000	0.000000000
4  2	1.952678571	1.952678571
5  2	0.000000000	0.000000000
1  3	2.603571429	2.603571429
2  3	0.000000000	0.000000000
3  3	1.952678571	1.952678571
4  3	0.000000000	0.000000000
5  3	1.997057630	1.997057630
1  4	0.000000000	0.000000000
2  4	1.952678571	1.952678571
3  4	0.000000000	0.000000000
4  4	1.997057630	1.997057630
5  4	-0.000000000	0.000000000
1  5	1.952678571	1.952678571
2  5	0.000000000	0.000000000
3  5	1.997057630	1.997057630
4  5	-0.000000000	0.000000000
5  5	2.419512128	2.419512128

以上が重なり行列の要素である.

ハミルトニアン行列の導出

同様に

\begin{align} H_{ij} &= \langle\phi_i|\hat{H}|\phi_j\rangle \\ &= \int_{-a}^{a} \phi_i^\ast(x) \hat{H} \phi_j(x) \mathrm{d}x \\ &= \int_{-a}^{a} \left[ -x^{i-1} (x-a) (x+a) \right] \left[ - \frac{\hbar^2}{2m} \frac{\mathrm{d}^2}{\mathrm{d}x ^2} \right] \left[ -x^{j-1} (x-a) (x+a) \right] \mathrm{d}x \\ &= - \frac{\hbar^2}{2m} \int_{-a}^{a} \left[ x^{i-1} (x^2-a^2) \right] \frac{\mathrm{d}^2}{\mathrm{d}x ^2} \left[ x^{j-1} (x^2-a^2) \right] \mathrm{d}x \\ &= - \frac{\hbar^2}{2m} \int_{-a}^{a} \left[ x^{i+1} - a^2 x^{i-1} \right] \frac{\mathrm{d}^2}{\mathrm{d}x ^2} \left[ x^{j+1} - a^2 x^{j-1} \right] \mathrm{d}x \\ &= - \frac{\hbar^2}{2m} \int_{-a}^{a} \left[ x^{i+1} - a^2 x^{i-1} \right] \left[(j+1) (j) x^{j-1} - a^2 (j-1) (j-2) x^{j-3} \right] \mathrm{d}x \\ &= - \frac{\hbar^2}{2m} \int_{-a}^{a} \left[ j(j+1) x^{i+j} - a^2 j(j+1) x^{i+j-2} - a^2 (j-1)(j-2) x^{i+j-2} + a^4 (j-1)(j-2) x^{i+j-4} \right] \mathrm{d}x \\ &= - \frac{\hbar^2}{2m} \int_{-a}^{a} \left[ j(j+1) x^{i+j} - a^2 (j^2 + j) x^{i+j-2} - a^2 (j^2 -3j +2) x^{i+j-2} + a^4 (j-1)(j-2) x^{i+j-4} \right] \mathrm{d}x \\ &= - \frac{\hbar^2}{2m} \int_{-a}^{a} \left[ j(j+1) x^{i+j} - a^2 (2j^2 -2j +2) x^{i+j-2} + a^4 (j-1)(j-2) x^{i+j-4} \right] \mathrm{d}x \\ &= - \frac{\hbar^2}{2m} \left[ j(j+1) \frac{a \left[a^{i+j}+(-a)^{i+j}\right]}{i+j+1} - a^2 (2j^2 -2j +2) \frac{a \left[a^{i+j-2}+(-a)^{i+j-2}\right]}{i+j-2+1} + a^4 (j-1)(j-2) \frac{a \left[a^{i+j-4}+(-a)^{i+j-4}\right]}{i+j-4+1} \right] \\ &= - \frac{\hbar^2}{2m} \left[ j(j+1) \frac{a \left[a^{i+j}+(-a)^{i+j}\right]}{i+j+1} - (2j^2 -2j +2) \frac{a \left[a^{i+j}+(-a)^{i+j}\right]}{i+j-1} + (j-1)(j-2) \frac{a \left[a^{i+j}+(-a)^{i+j}\right]}{i+j-3} \right] \\ &= - \frac{\hbar^2}{2m} a \left[a^{i+j}+(-a)^{i+j}\right] \left[ \frac{j^2+j}{i+j+1} - 2 \frac{j^2-j+1}{i+j-1} + \frac{j^2-3j+2}{i+j-3} \right] \\ \end{align}

と求められる. 同じく数値積分によって検証したが, ほぼ一致しているので問題ない.

ハミルトニアン行列の検証

# using Printf

int(func, x_min, x_max, dx) = sum(func.(x_min:dx:x_max))*dx
diff2(func, x, dx) = ((func(x+dx)-func(x))/dx-(func(x)-func(x-dx))/dx)/dx

a = 1.5
m = 1.0
ℏ = 1.0
N = 5
H = zeros(N,N)

println("  i  j\tnumerical\tanalytical")
for j in 1:N
    for i in 1:N
        numerical = -ℏ^2/(2*m)*int(x->φ(i,x,a=a)*diff2(x->φ(j,x,a=a),x,0.001), -a, a, 0.001)
        H[i,j] = isodd(i+j) ? 0 : -ℏ^2/(2*m) * a * (a^(i+j)+(-a)^(i+j)) * ((j^2+j)/(i+j+1) - 2*(j^2-j+1)/(i+j-1) + (j-1)*(j-2)/(i+j-3))
        @printf("%3d%3d\t%.9f\t%.9f\n", i, j, numerical, H[i,j])
    end
end

出力

i  j	numerical	analytical
1  1	4.499999500	4.500000000
2  1	-0.000000000	0.000000000
3  1	2.024998875	2.025000000
4  1	-0.000000000	0.000000000
5  1	1.952676040	1.952678571
1  2	0.000000000	0.000000000
2  2	6.074996625	6.075000000
3  2	0.000000000	0.000000000
4  2	5.858028121	5.858035714
5  2	0.000000000	0.000000000
1  3	2.024998875	2.025000000
2  3	0.000000000	0.000000000
3  3	7.159810797	7.159821429
4  3	0.000000000	0.000000000
5  3	10.251535976	10.251562500
1  4	0.000000000	0.000000000
2  4	5.858028121	5.858035714
3  4	0.000000000	0.000000000
4  4	11.227871682	11.227901786
5  4	0.000000000	0.000000000
1  5	1.952676040	1.952678571
2  5	0.000000000	0.000000000
3  5	10.251535976	10.251562500
4  5	0.000000000	0.000000000
5  5	19.471228780	19.471311891

変分計算

eigen()に行列 $\boldsymbol{H}, \boldsymbol{S}$ を入力として渡すと, 固有ベクトルを並べた行列 $\boldsymbol{C}$ と固有値を並べたベクトル $E$ が出力として得られる.

変分計算のルーチン

# using LinearAlgebra

function variational(N; ℏ=1, a=1, m=0.5)
    S = zeros(N,N)
    H = zeros(N,N)
    for j in 1:N
        for i in 1:N
            # S[i,j] = integral(x->φ(i,x,a=a)*φ(j,x,a=a), -a, a, 0.001)                         # 数値積分版
            # H[i,j] = -ℏ^2/(2*m)*int(x->φ(i,x,a=a)*diff2(x->φ(j,x,a=a),x,0.001), -a, a, 0.001) # 数値積分版
            S[i,j] = isodd(i+j) ? 0 : a^3 * (a^(i+j)+(-a)^(i+j)) * (1/(i+j+3) - 2/(i+j+1) + 1/(i+j-1))
            H[i,j] = isodd(i+j) ? 0 : -ℏ^2/(2*m) * a * (a^(i+j)+(-a)^(i+j)) * ((j^2+j)/(i+j+1) - 2*(j^2-j+1)/(i+j-1) + (j-1)*(j-2)/(i+j-3))
        end
    end
    return eigen(Symmetric(H),Symmetric(S))
end

基底関数を増やすごとに厳密解に近づく様子が見て取れる. これはJ.M.ティッセン著『計算物理学』(丸善出版, 2012) 3.2.1 とも一致する結果である.

変分計算

# using Printf
ℏ = 1.0
a = 1.0
m = 0.5

println("N \tE₁ \t\tE₂ \t\tE₃ \t\tE₄ ")

for N in 4:2:16
    @printf("%d\t%.6f\t%.6f\t%.6f\t%.6f\n", N, variational(N, ℏ=ℏ, a=a, m=m).values[1:4]...)
end

@printf("%s\t%.6f\t%.6f\t%.6f\t%.6f\n", "∞", [E(i, ℏ=ℏ, a=a, m=m) for i in 1:4]...)

出力

N 	E₁ 		E₂ 		E₃ 		E₄ 
4	2.467437	9.875388	25.532563	50.124612
6	2.467401	9.869618	22.293406	39.997979
8	2.467401	9.869604	22.207367	39.489241
10	2.467401	9.869604	22.206612	39.478505
12	2.467401	9.869604	22.206610	39.478418
14	2.467401	9.869604	22.206610	39.478418
16	2.467401	9.869604	22.206610	39.478418
∞	2.467401	9.869604	22.206610	39.478418

最後に固有ベクトル, つまり波動関数の収束を見よう. 試行波動関数は基底関数の線形結合であり, その展開係数は固有ベクトルの各成分である. ただし, LAPACKでは固有ベクトル $\psi_n(x)$ の正負が逆になることがある（数学的には間違っていないがグラフを比較しにくい）ので, 2乗した確率密度関数 $|\psi_n(x)|^2$ を描写している. 基底関数が8個の時点で第3励起状態( $n=4$ )までは厳密解とほぼ一致して見える.

試行波動関数

trial(n,x,C;a=1) = sum(C[:,n] .* [φ(i, x; a=a) for i in 1:size(C)[1]])

試行波動関数と厳密解の比較

# using Plots

function plots(N)
    plts=[]
    C = variational(N).vectors
    for n in 1:4
        plot(xlabel="\$x\$", ylabel="\$|\\psi_$n(x)|^2\$", legend=(n==1 ? :bottom : nothing), background_color_legend=nothing, foreground_color_legend=nothing)
        plot!(-a:0.01:a, x->ψ(n,x)^2, label="exact", lw=2, lc="#B61972")
        plot!(-a:0.01:a, x->trial(n,x,C)^2, label="trial", lw=2, lc="#007AB7", ls=:dash)
        plt = plot!()
        push!(plts, plt)
    end
    "基底関数の数：$N" |> display
    plot(plts...) |> display
end

plots(4)
plots(6)
plots(8)

基底関数の数：4

基底関数の数：6

基底関数の数：8

おわりに

「Juliaで線形代数」シリーズの最後はRayleigh-Ritzの変分法による井戸型ポテンシャルの数値計算であった. お楽しみ頂けただろうか. これは線形代数なのか？という疑問はさておき, Trotter分解やエルミート行列についての解説はなかなか教育的な内容であると思う. 井戸型ポテンシャルの変分計算も, 行列要素の計算が簡単なため, 変分法の入門にうってつけだろう.

2022年3月15日(火)～3月19日(土)に開催される日本物理学会第77回年次大会の領域12(化学物理)にて, Rayleigh-Ritzの変分法をクーロン三体系に適用した結果について報告する予定である. 興味のある方はぜひ聴きに来てほしい.

参考文献

https://ja.wikipedia.org/wiki/井戸型ポテンシャル
http://hooktail.sub.jp/quantum/squarewell-inf/
J.M.ティッセン著『計算物理学』(丸善出版, 2012) 3.1 変分計算, 3.2.1 無限に深い井戸型ポテンシャル
J.P.キリンベック著『パーソナルコンピュータによる量子力学』(サイエンス社,1985) 8.2 量子力学における行列
A.ザボ, Ｎ.Ｓ.オストランド著『新しい量子化学上電子構造の理論入門』(東京大学出版会, 1987) 1.3.2 線形変分問題

この記事の元になったノートブックはGistにアップロードしてある.

パッケージ

線形変分法（Rayleigh-Ritzの変分法）

無限に深い井戸型ポテンシャル

基底関数

行列要素

重なり行列の導出

ハミルトニアン行列の導出

変分計算

基底関数の数：4

基底関数の数：6

基底関数の数：8

おわりに

参考文献

Discussion