以下のサイトの下の方に、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" )
最初の固有値25個くらいは値が負なので、束縛状態であることがわかる。
最初の5個くらいの波動関数を確認すると、しっかり中心の有限井戸に束縛されている。
プロットしてみるとわかるが、energyの図の真ん中らへんは、無限井戸も含めた全体に広がる解で、さらに高いエネルギーの解は有限井戸を避けて片側に強く局在した状態になる。
片側だけに寄るのは対称性的におかしいが、そこはまぁ数値解の御愛嬌という感じだろうか。
とりあえず、デカイ箱を用意しておいて、その中に引力ポテンシャルを入れれば、束縛解をパッと出せることがわかった。