ポテンシャルがある場合の1Dシュレーディンガー方程式の解

一次元系のシュレーディンガー方程式は、

\[\left( -\frac{\hbar^2}{2m}\frac{d^2}{dx^2} + V(x) \right) \psi(x) = \epsilon \psi(x)\]

である。前回はポテンシャル$V(x)$がゼロの場合を解いたが、次はポテンシャルを入れて解いてみよう。

前回のおさらい

座標を離散化すると、シュレーディンガー方程式は、

\[\sum_j \left( -\frac{\hbar^2}{2m} d_j + V(x_j) \delta_{ij} \right) \psi(x_j) = \epsilon \psi(x_i)\]

と書くことができる。ここで、

\[\frac{d^2}{dx^2}\psi |_{x=x_i} \rightarrow \sum_j d_j \psi(x_j) \]

である。この方程式は、先ほどと同様に行列とベクトルで表現することができて、

\[\hat{H} {\bf \psi} = \epsilon {\bf \psi}\]

となる。これがシュレーディンガー方程式の行列表示である。ここで、${\bf \psi}_n = U {\bf c}_n$を使えば、

\[\hat{H} U {\bf c}_n = \epsilon U {\bf c}_n \]

\[U^+ \hat{H} U {\bf c}_n = \epsilon U^+ U {\bf c}_n \]

\[\hat{H}' {\bf c}_n = \epsilon {\bf c}_n\]

となる。ベクトルはユニタリー変換で色々変換できるので、kでもxでも好きな表示で計算できる。

みてわかるように、$a$の大きさを小さくすると左右の均衡が少し改善された。つまり、数値計算上のエラーであることがわかる。$a$が小さくなればさらに改善されるだろう。矩形ポテンシャルをきちんと中心においた場合$center = (L+1)/2$であれば、改善される。このような違いは、矩形のポテンシャルに不連続な変化があるからである。いま、差分化している微分演算子の精度が悪いのである。これを確認するために、ポテンシャルの形状をなだらかにしてみよう。

function calc_V(N,V0)
    vec_V = zeros(Float64,N)
    dx = N/6
    center = (N)/2
    for i in 1:N
        vec_V[i] = V0*exp(-(i-center)^2/(dx^2))
    end
    return vec_V
end

N = 1000
xaxis = Int64[]
for i in 1:N
    push!(xaxis,i)
end

V0=1.0
vec_V = calc_V(N,V0)
plot(xaxis[1:N],vec_V[1:N],label="Potential")

fig1

ポテンシャルを、上のようにガウス関数にした。このとき、最小固有値の固有関数は

N = 1000
a = 0.01
function gs2(N,a)
    groundstates = []
    labels = []
    for v in 1:10
        V0 = v*0.5
        mat_H = make_H1dv(N,a,V0)
        ε,ψ = eigen(mat_H)
        println("Potential = ",V0," Minimum eigenvalue = ",ε[1])
        push!(groundstates,ψ[:,1])
        if v == 1
            labels = [string(V0)]
        else
            labels = [labels string(V0)]
        end
    end    
    return groundstates,labels
end
groundstates,labels = gs2(N,a)

Potential = 0.5 Minimum eigenvalue = 0.3216661503241833
Potential = 1.0 Minimum eigenvalue = 0.47496666834409224
Potential = 1.5 Minimum eigenvalue = 0.5787666683372826
Potential = 2.0 Minimum eigenvalue = 0.6538526217544764
Potential = 2.5 Minimum eigenvalue = 0.7124623403946451
Potential = 3.0 Minimum eigenvalue = 0.7609213437478689
Potential = 3.5 Minimum eigenvalue = 0.8025976518969778
Potential = 4.0 Minimum eigenvalue = 0.83945543554293
Potential = 4.5 Minimum eigenvalue = 0.8727520445958912
Potential = 5.0 Minimum eigenvalue = 0.9033145515618101

plot(xaxis[1:N],groundstates,label=labels)

fig1

先ほどよりもましになったが、相変わらず$V0=5$では、左右対称が崩れてしまっている。これは、ポテンシャルの変化が急激すぎることに起因している。最後に、ポテンシャルの中心を厳密に系の中心とおくと

function calc_V(N,V0)
    vec_V = zeros(Float64,N)
    dx = N/6
    center = (N+1)/2
    for i in 1:N
        vec_V[i] = V0*exp(-(i-center)^2/(dx^2))
    end
    return vec_V
end

N = 1000
a = 0.01
function gs3(N,a)
    groundstates = []
    labels = []
    for v in 1:10
        V0 = v*0.5
        mat_H = make_H1dv(N,a,V0)
        ε,ψ = eigen(mat_H)
        println("Potential = ",V0," Minimum eigenvalue = ",ε[1])
        push!(groundstates,ψ[:,1])
        if v == 1
            labels = [string(V0)]
        else
            labels = [labels string(V0)]
        end
    end    
    return groundstates,labels
end
groundstates,labels = gs3(N,a)

Potential = 0.5 Minimum eigenvalue = 0.3216680437594529
Potential = 1.0 Minimum eigenvalue = 0.47497702625749055
Potential = 1.5 Minimum eigenvalue = 0.5788014430373154
Potential = 2.0 Minimum eigenvalue = 0.6539436690715744
Potential = 2.5 Minimum eigenvalue = 0.7126667576418749
Potential = 3.0 Minimum eigenvalue = 0.7613285627430311
Potential = 3.5 Minimum eigenvalue = 0.8033158575112843
Potential = 4.0 Minimum eigenvalue = 0.8405613523837149
Potential = 4.5 Minimum eigenvalue = 0.8742498997697548
Potential = 5.0 Minimum eigenvalue = 0.905155534659135

plot(xaxis[1:N],groundstates,label=labels)

fig1

きれいに左右対称になった。このことから、厳密に左右対称なものが欲しい場合には、数値計算でもきちんと左右対称にすべきであることがわかる。

ポテンシャルがある場合の1Dシュレーディンガー方程式の解

前回のおさらい

離散化座標表示での解

ポテンシャルがない場合

ポテンシャルがある場合