Я пишу код на питоне для развития зависящего от времени уравнения Шредингера с использованием схемы Крэнка-Николсона. Я не знал, как справиться с этим потенциалом, поэтому огляделся и нашел путь от -crank-nicolson-method">этот вопрос, который я проверил из нескольких других источников. Согласно им, для гармонического потенциала осциллятора схема C-N дает
AΨn+1=A∗Ψn
где элементы на главной диагонали A равны dj=1+[(i∆t)/(2m(∆x)^2)]+[(i∆t(xj)^2)/4], а элементы на верхней и нижней диагоналях a=−i∆t/[4m(∆x)^2]
Насколько я понимаю, я должен задать начальное условие (я выбрал когерентное состояние) в виде матрицы Ψn и мне нужно вычислить матрицу Ψn+ 1 , которая представляет собой волновую функцию после времени Δt. Чтобы получить Ψn+1 для данного шага, я инвертирую матрицу A и умножаю ее на матрицу A*, а затем умножаю результат с Ψn. Результирующая матрица становится Ψn для следующего шага.
Но когда я это делаю, я получаю неправильную анимацию. Предполагается, что волновой пакет колеблется между границами, но в моей анимации он едва отклоняется от своего начального среднего значения. Я просто не понимаю, что я делаю неправильно. Мое понимание проблемы неверно? Или это ошибка в моем коде? Пожалуйста, помогите! Я разместил свой код ниже и видео своей анимации здесь . Извините за длину кода и вопрос, но это сводит меня с ума, не зная, в чем моя ошибка.
import numpy as np
import matplotlib.pyplot as plt
L = 30.0
x0 = -5.0
sig = 0.5
dx = 0.5
dt = 0.02
k = 1.0
w=2
K=w**2
a=np.power(K,0.25)
xs = np.arange(-L,L,dx)
nn = len(xs)
mu = k*dt/(dx)**2
dd = 1.0+mu
ee = 1.0-mu
ti = 0.0
tf = 100.0
t = ti
V=np.zeros(len(xs))
u=np.zeros(nn,dtype="complex")
V=K*(xs)**2/2 #harmonic oscillator potential
u=(np.sqrt(a)/1.33)*np.exp(-(a*(xs - x0))**2)+0j #initial condition for wave function
u[0]=0.0 #boundary condition
u[-1] = 0.0 #boundary condition
A = np.zeros((nn-2,nn-2),dtype="complex") #define A
for i in range(nn-3):
A[i,i] = 1+1j*(mu/2+w*dt*xs[i]**2/4)
A[i,i+1] = -1j*mu/4.
A[i+1,i] = -1j*mu/4.
A[nn-3,nn-3] = 1+1j*mu/2+1j*dt*xs[nn-3]**2/4
B = np.zeros((nn-2,nn-2),dtype="complex") #define A*
for i in range(nn-3):
B[i,i] = 1-1j*mu/2-1j*w*dt*xs[i]**2/4
B[i,i+1] = 1j*mu/4.
B[i+1,i] = 1j*mu/4.
B[nn-3,nn-3] = 1-1j*(mu/2)-1j*dt*xs[nn-3]**2/4
X = np.linalg.inv(A) #take inverse of A
plt.ion()
l, = plt.plot(xs,np.abs(u),lw=2,color='blue') #plot initial wave function
T=np.matmul(X,B) #multiply A inverse with A*
while t<tf:
u[1:-1]=np.matmul(T,u[1:-1]) #updating u but leaving the boundary conditions unchanged
l.set_ydata((abs(u))) #update plot with new u
t += dt
plt.pause(0.00001)
