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
最近のコメント