nano_exit

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

Juliaで一次元井戸型ポテンシャル

以下のサイトの下の方に、Juliaで一次元のシュレーディンガー方程式を解くPDFが紹介されている。
物理ノートby永井

Juliaの練習としてやってみた。
PDF内では、無限の井戸の中に斥力ポテンシャルを入れた場合をやっているが、ここでは引力ポテンシャルに対してやってみる。
解き方としては、中心差分を用いて二階微分を近似した有限要素法に対応する(と思っている)。
とりあえず、引力ポテンシャルは滅茶苦茶深くした。

using Plots
using LinearAlgebra

function calc_V( N, a, V0 )
    vec_V = zeros( Float64, N )
    center = N * a / 2
    L = N * a / 2
    for i in 1 : N
        x = i * a
        if center - L / 2  <= x <= center + L / 2
            vec_V[ i ] = - V0
        end
    end
    return vec_V
end

function make_H1dv( N, a, V0 )
    mat_H = zeros( Float64, N, N )
    vec_V = calc_V( N, a, V0 )
    for i in 1 : N
        for dx in -1 : 1
            j = i + dx
            v = -1 / ( 2 * a^2 )
            if dx == 0
                v = 1 / a^2 + vec_V[ i ]
            end
            if 1 <= j <= N
                mat_H[ i, j ] = v
            end
        end
    end
    return mat_H
end

function main( N, a, V0 )
    N = 100
    a = 1 / N
    V0 = 1e4
    mat_H = make_H1dv( N, a, V0 )
    energy, phi = eigen( mat_H )
    return energy, phi
end

energy, phi = main()

plot( energy )
savefig( "energy.png" )

plot( phi[ :, 1:5 ] )
savefig( "phi.png" )

f:id:koideforest:20181211191822p:plain
energy
最初の固有値25個くらいは値が負なので、束縛状態であることがわかる。

最初の5個くらいの波動関数を確認すると、

f:id:koideforest:20181211191803p:plain
phi
しっかり中心の有限井戸に束縛されている。

プロットしてみるとわかるが、energyの図の真ん中らへんは、無限井戸も含めた全体に広がる解で、さらに高いエネルギーの解は有限井戸を避けて片側に強く局在した状態になる。
片側だけに寄るのは対称性的におかしいが、そこはまぁ数値解の御愛嬌という感じだろうか。

とりあえず、デカイ箱を用意しておいて、その中に引力ポテンシャルを入れれば、束縛解をパッと出せることがわかった。