RSS

Schrödinger equation

シュレーディンガー方程式(Schrödinger Equation)の時間発展の計算。

時間依存のシュレーディンガー方程式。 $$ \displaystyle i \hbar \frac{\partial \Psi(\mathbf{x}, t) }{\partial t} = \left( -\frac{\hbar^2}{2 m}\nabla^2 + V(\mathbf{x}, t) \right) \Psi(\mathbf{x},t) $$

ポテンシャルが$V(x) = x^2$ の場合の挙動。

アニメーションはこちら。

Pythonスクリプト

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation, rc
from IPython.display import HTML

N = 50
numIter = 10000
dt = 0.001
dx = 10.0 / N
m = 1
hb = 1

# allocate zero arrays
x = np.zeros(N, dtype=np.float)
phi0 = np.zeros(N, dtype=np.complex)
phi1 = np.zeros(N, dtype=np.complex)
V = np.zeros(N, dtype=np.float)

# create fig, ax
fig, ax = plt.subplots()
plt.close()
ax.set_xlim((-5, 5))
ax.set_ylim((-1, 1))
line1, = ax.plot([], [], '-', label='real', lw=1.0)
line2, = ax.plot([], [], '-', label='imag', lw=1.0)
line3, = ax.plot([], [], '-', label='abs', lw=1.0)
ax.grid()
fig.legend(bbox_to_anchor=(0.9, 0.9), loc='upper right', borderaxespad=1, fontsize=10)

# initialization
def init():
  print("init called")
  dp = 0.0
  for i in range(0, N):
    global x, V, phi0, phi1
    x[i] = dx*i - 5.0
    V[i] = x[i]**2
    phi0[i] = np.exp(-((x[i]-1.0))**2)
    dp += np.abs(phi0[i])

  # normalize
  for i in range(0, N):
    phi0[i] /= dx*dp
  

  line1.set_data(x, phi0.real)
  line2.set_data(x, phi0.imag)
  line3.set_data(x, np.abs(phi0))
  return (line1, line2, line3)

# animation
def animate(k):
  if k % 100 == 0:
    print(k)
  global dt, hb, m, dx, phi0, phi1, V
  for i in range(1, N-1):
    phi1[i] = phi0[i] + (dt/(1.j*hb))*((-hb*hb/(2*m*dx*dx))*(phi0[i-1] - 2*phi0[i] + phi0[i+1]) + V[i]*phi0[i])
  
  phi0 = phi1

  line1.set_data(x, phi0.real)
  line2.set_data(x, phi0.imag)
  line3.set_data(x, np.abs(phi0))
  return (line1, line2, line3)

plt.rcParams['animation.embed_limit'] = 2**128  # avoid limitation of embedded animation size

anim = animation.FuncAnimation(fig, animate, init_func=init, frames=numIter, interval=20, blit=True)

plt.rcParams['animation.ffmpeg_path'] = '/usr/bin/ffmpeg' # For google colab
HTML(anim.to_html5_video())

#rc('animation', html='jshtml')
#anim