Google ColabのJulia言語で1次元シュレーディンガー方程式の井戸型ポテンシャルによる波束の散乱を計算してアニメーションにしてみる。空間方向は普通の差分、時間は自前の13段8次のルンゲクッタ法を使用。

今回はシュレーディンガー方程式をやってみる。

昔Excel VBAで計算して、Pythonでもやってみた1次元シュレーディンガー方程式のポテンシャルによる波束の散乱を計算しよう。
これはシッフの量子力学の例題で出ていたものを「計算物理」という本でfortranで書かれていたものの移植。

ただ、空間微分は普通に(φ(i+1) + φ(i-1) - 2 φ(i))/h²で差分化するが、時間方向は13段8次のルンゲクッタ法を使う。DifferentialEquations.jlを使えば一発だが、Google Colabで使うときは毎回Pkg.addとしないといけなくて、それが45分かかるのでもうやだ、ということで自前のもの（Excel VBA移植版）。

コードはこんな感じで。ルンゲクッタ部分の係数が長い…

using Plots

using Printf

function main()

    n = 2000 #空間分割数

    m = 5500 #時間分割数

    #波動関数の実部u、虚部v

    u = zeros(Float64, n + 2)

    v = zeros(Float64, n + 2)

    ud = zeros(Float64, n + 2)

    vd = zeros(Float64, n + 2)

    #ポテンシャル

    u_pot = zeros(Float64, n + 2)

    #物理パラメータ

    v_width = 0.064

    v0 = 0.6 * (70.7 * pi) ^ 2

    deltax = 0.035

    x0 = -0.3

    p0 = pi * 50.0

    t = 0.0

    dx = 0.001

    dt = dx * dx / 2.0

    xmax = dx * (n - 1) / 2.0

    xmin = -dx * (n - 1) / 2.0

    x = zeros(Float64, n + 2)

    #表示用

    results = []

    ndiv = 25

    #8次のルンゲクッタ法の係数

    a = zeros(Float64, 13)

    b = zeros(Float64, 13, 13)

    c = zeros(Float64, 13)

    ku = zeros(Float64, n + 1, 13)

    kv = zeros(Float64, n + 1, 13)

    a[1] = 0.0

    a[2] = 2.0 / 27.0

    a[3] = 1.0 / 9.0

    a[4] = 1.0 / 6.0

    a[5] = 5.0 / 12.0

    a[6] = 0.5

    a[7] = 5.0 / 6.0

    a[8] = 1.0 / 6.0

    a[9] = 2.0 / 3.0

    a[10] = 1.0 / 3.0

    a[11] = 1.0

    a[12] = 0.0

    a[13] = 1

    b[2, 1] = 2.0 / 27.0

    b[3, 1] = 1.0 / 36

    b[3, 2] = 1.0 / 12.0

    b[4, 1] = 1.0 / 24.0

    b[4, 3] = 1.0 / 8.0

    b[5, 1] = 5.0 / 12.0

    b[5, 3] = -25.0 / 16.0

    b[5, 4] = 25.0 / 16.0

    b[6, 1] = 1.0 / 20.0

    b[6, 4] = 1.0 / 4.0

    b[6, 5] = 1.0 / 5.0

    b[7, 1] = -25.0 / 108.0

    b[7, 4] = 125.0 / 108.0

    b[7, 5] = -65.0 / 27.0

    b[7, 6] = 125.0 / 54.0

    b[8, 1] = 31.0 / 300.0

    b[8, 5] = 61.0 / 225.0

    b[8, 6] = -2.0 / 9.0

    b[8, 7] = 13.0 / 900.0

    b[9, 1] = 2.0

    b[9, 4] = -53.0 / 6.0

    b[9, 5] = 704.0 / 45.0

    b[9, 6] = -107.0 / 9.0

    b[9, 7] = 67.0 / 90.0

    b[9, 8] = 3.0

    b[10, 1] = -91.0 / 108.0

    b[10, 4] = 23.0 / 108.0

    b[10, 5] = -976.0 / 135.0

    b[10, 6] = 311.0 / 54.0

    b[10, 7] = -19.0 / 60.0

    b[10, 8] = 17.0 / 6.0

    b[10, 9] = -1.0 / 12.0

    b[11, 1] = 2383.0 / 4100.0

    b[11, 4] = -341.0 / 164.0

    b[11, 5] = 4496.0 / 1025.0

    b[11, 6] = -301.0 / 82.0

    b[11, 7] = 2133.0 / 4100.0

    b[11, 8] = 45.0 / 82.0

    b[11, 9] = 45.0 / 164.0

    b[11, 10] = 18.0 / 41.0

    b[12, 1] = 3.0 / 205.0

    b[12, 6] = -6.0 / 41.0

    b[12, 7] = -3.0 / 205.0

    b[12, 8] = -3.0 / 41.0

    b[12, 9] = 3.0 / 41.0

    b[12, 10] = 6.0 / 41.0

    b[13, 1] = -1777.0 / 4100.0

    b[13, 4] = -341.0 / 164.0

    b[13, 5] = 4496.0 / 1025.0

    b[13, 6] = -289.0 / 82.0

    b[13, 7] = 2193.0 / 4100.0

    b[13, 8] = 51.0 / 82.0

    b[13, 9] = 33.0 / 164.0

    b[13, 10] = 12.0 / 41.0

    b[13, 12] = 1.0

    c[6] = 34.0 / 105.0

    c[7] = 9.0 / 35.0

    c[8] = 9.0 / 35.0

    c[9] = 9.0 / 280.0

    c[10] = 9.0 / 280.0

    c[12] = 41.0 / 840.0

    c[13] = 41.0 / 840.0

    #ポテンシャルと波動関数の初期値（ガウシアン波束）設定

    for i in 2:(n + 1)

        x[i] = xmin + (i - 1) * dx

        if x[i] >= -v_width / 2.0 && x[i] <= v_width / 2.0

            u_pot[i] = v0

        else

            u_pot[i] = 0.0

        end

        u[i] = exp(-((x[i] - x0) ^ 2) / (4 * deltax * deltax)) * cos(p0 * x[i]) / ((2 * pi * deltax * deltax) ^ 0.25)

        v[i] = exp(-((x[i] - x0) ^ 2) / (4 * deltax * deltax)) * sin(p0 * x[i]) / ((2 * pi * deltax * deltax) ^ 0.25)

    end

    u[1] = u[2]

    u[n + 2] = u[n + 1]

    v[1] = v[2]

    v[n + 2] = v[n + 1]

    x[1] = xmin - dx

    x[n + 2] = xmax + dx

    for i in 1:m

        # 8次のルンゲクッタ法計算

         @inbounds for ii in 2:(n + 1)

            ku[ii - 1, 1] = f(ii, t + dt, u, v, dx, u_pot, n) * dt

            kv[ii - 1, 1] = g(ii, t + dt, u, v, dx, u_pot, n) * dt

        end

        @inbounds for j in 2:13

            @inbounds @simd for ii in 2:(n + 1)

                ud[ii] = u[ii]

                vd[ii] = v[ii]

            end

            @inbounds for k in 1:(j - 1)

                @inbounds @simd for ii in 2:(n + 1)

                    ud[ii] = ud[ii] + b[j, k] * ku[ii - 1, k]

                    vd[ii] = vd[ii] + b[j, k] * kv[ii - 1, k]

                end

            end

            @inbounds for ii in 2:(n + 1)

                ku[ii - 1, j] = f(ii, t + a[j] * dt, ud, vd, dx, u_pot, n) * dt

                kv[ii - 1, j] = g(ii, t + a[j] * dt, ud, vd, dx, u_pot, n) * dt

            end

        end

        @inbounds @simd for j in 1:13

            for ii in 2:(n + 1)

                u[ii] = u[ii] + c[j] * ku[ii - 1, j]

                v[ii] = v[ii] + c[j] * kv[ii - 1, j]

            end

        end

        t = t + dt

        # ndivステップごとに結果(波動関数の大きさ)を配列に格納

        if i % ndiv == 0

            push!(results, sqrt.(u.^2 + v.^2))

        end

    end

    # 計算結果をアニメーションで表示

    anim = @animate for i in 1:length(results)

        plot(x, results[i], ylim = (0, 4), linewidth = 5,title="Tunnelling Effect", label="Wave function", xlabel="x", ylabel="phi", size=(900,500))

        plot!(x, u_pot, label = "Potential")

    end

    gif(anim, "Tunnel.gif", fps = 30)    

end

function f(i, t, u, v, dx, u_pot, n)

    v[1] = v[2]

    v[n + 2] = v[n + 1]

    d2 = (v[i + 1] + v[i - 1] - 2.0 * v[i]) / (dx * dx)

    return -d2 + u_pot[i] * v[i]

end

function g(i, t, u, v, dx, u_pot, n)

    u[1] = u[2]

    u[n + 2] = u[n + 1]

    d2 = (u[i + 1] + u[i - 1] - 2.0 * u[i]) / (dx * dx)

    return d2 - u_pot[i] * u[i]

end

実行結果。一部は通過して一部は反射している。

次はトンネル効果やってみよう。