« 松屋で甘唐辛子のトロたまごろチキコンボ牛めしに七味唐辛子をたっぷりかけていただく。 | トップページ | 「博士を殺した数式」を読んだ。カオス理論の大数学者が殺され、残された数式とは?というものだが、長男が超弦理論の研究者で、著名な物理学者の醜聞(浮気癖)も書かれてたり。そして謝辞がウィッテンやシュワルツなど。 »

2020年7月20日 (月)

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アニメが始まります。

Sh01

 

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

%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)

« 松屋で甘唐辛子のトロたまごろチキコンボ牛めしに七味唐辛子をたっぷりかけていただく。 | トップページ | 「博士を殺した数式」を読んだ。カオス理論の大数学者が殺され、残された数式とは?というものだが、長男が超弦理論の研究者で、著名な物理学者の醜聞(浮気癖)も書かれてたり。そして謝辞がウィッテンやシュワルツなど。 »

学問・資格」カテゴリの記事

コメント

コメントを書く

(ウェブ上には掲載しません)

« 松屋で甘唐辛子のトロたまごろチキコンボ牛めしに七味唐辛子をたっぷりかけていただく。 | トップページ | 「博士を殺した数式」を読んだ。カオス理論の大数学者が殺され、残された数式とは?というものだが、長男が超弦理論の研究者で、著名な物理学者の醜聞(浮気癖)も書かれてたり。そして謝辞がウィッテンやシュワルツなど。 »

最近の記事

2021年5月
            1
2 3 4 5 6 7 8
9 10 11 12 13 14 15
16 17 18 19 20 21 22
23 24 25 26 27 28 29
30 31          
フォト
無料ブログはココログ