Python＋Scipyでルンゲクッタ8次のDOP853(Dormand＆Prince)を使う（その16) 時間依存のシュレーディンガー方程式のトンネル効果を計算してGIFアニメに。

時間依存のシュレーディンガー方程式も空間微分を差分化すると普通の常微分方程式になるので、ルンゲクッタ8次のDormand-Prince法で計算できる。
具体的には空間の拡散項は2次の差分
(φi+1 - 2*φi + φi-1) / (dx*dx)
で近似。中央に井戸型ポテンシャルを置いて、左から波束が入ってくるときの計算。

ただし、複素数になるのでライブラリがodeからcomplex_odeに変わる。

ずっと前にExcel VBAでやったもののPythonでの書き直し。これは”計算物理”にでてたもので、もともとはシッフの例題にあるものだそうです。

クリックするとGIFアニメが始まります。

 

ソースはこんな感じで（インデントがないけどココログにコピペすると消えるので、、、）

%matplotlib notebook
import numpy as np
from scipy.integrate import complex_ode
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation
from matplotlib.animation import PillowWriter
import numba

INTERVAL = 0.05 # [s] interval time
FRAME_INTERVAL = 1000 * INTERVAL # [msec] interval between frames
FPS = 1000/FRAME_INTERVAL # frames per second

Nt=6000
Nx=2000

t0=0
dx = 0.001
dt = dx * dx / 2
xmax = dx * float(Nx) / 2
xmin = -dx * float(Nx) / 2
tmax=Nt*dt

v_width = 0.064
v0 = 0.6*((70.7 * np.pi) **2)
deltax = 0.035
x0 = -0.3
p0 = np.pi * 50
t=np.linspace(0,tmax,Nt)
x=np.linspace(xmin,xmax,Nx)
v=np.empty(Nx,dtype=np.complex)
u0=np.empty(Nx,dtype=np.complex)

for i in range(Nx):
u0[i]=np.exp(-((x[i]-x0)**2) / (4 * deltax * deltax))*np.exp(p0*x[i]*1j)/((2 * np.pi * deltax * deltax) ** 0.25)
if x[i]>=-v_width/2 and x[i] <=v_width/2:
v[i]=v0
else:
v[i]=0

sol= np.empty((Nt, Nx),dtype=np.complex)
sol[0] = u0

def Schroedinger(t, u): #odeintのときとt,xの並びが逆
f=np.zeros(Nx, dtype=np.complex)

for i in range(Nx):
if i < Nx-1:
p1=i + 1
else:
p1=0
if i > 0:
m1=i - 1
else:
m1=Nx-1

diff2 = (u[m1] + u[p1] -2.0* u[i]) / ( dx ** 2)
f[i] = (-1j)*(-diff2 + v[i]*u[i])

return f

numba_f = numba.jit(Schroedinger,nopython=True)

solver=complex_ode(Schroedinger)
solver.set_integrator('dop853')
solver.set_initial_value(u0,t0) #なぜか関数と並びが逆

k=1
while solver.successful() and solver.t < tmax:
solver.integrate(t[k])
sol[k] = solver.y
k+= 1

def update(i):
plt.cla() #現在描写されているグラフを消去
plt.xlim(-1,1)
plt.ylim(0,4)
plt.plot(x,np.abs(sol[i]),c='Red',lw=3)
plt.plot(x,4*abs(v)/v0,c='blue',lw=3)
plt.xlabel("x")
plt.ylabel("Phi")
plt.title("Schroedinger equation by DOP853")
plt.grid()
# Plot
fig = plt.figure(figsize=(10,6))
ani = FuncAnimation(fig, update, frames=range(0,Nt,15), interval=FRAME_INTERVAL, blit=True)
plt.show()
ani.save("sh01.gif", writer = 'pillow',fps=FPS)