目前我面临的问题是数值求解双阱谐波势的以下方程：

$iℏ\frac{\partial}{\partial t}ψ(x,t)=−\frac{ℏ}{2m}\frac{\partial ^2}{∂x^2}ψ(x,t)+V(x,t)ψ(x,t)$

在哪里

$V(x,t)=\frac{mω^2_0}{2}[a(t)x^4−b(t)x^2]$

其中和的形式为。我开始使用分步傅里叶方法算法，这是我的代码：$a(t)$$b(t)$$\frac{1}{2*\sqrt{const_1}}(const_1-const_2 \cos(\omega t))$

import numpy as np
import matplotlib.pyplot as plt

tfinal=10000#total simulation time
dt=0.001#time step
M=np.round(tfinal/dt)
time_steps=M.astype(int)# number of time steps

b=10.# final x
a=-10.# initial x
N=512# num of pieces

dx=(b-a)/N` # x-space step size
dk=2*np.pi/(b-a)
n=np.linspace(-N/2,N/2 -1,N) # this is gonna be my n's in x_n=a+dx*n

x=n*dx
k=n*dk
t=0

alpha=0.0005
beta=0.0001
acoef=(alpha-beta*np.cos(t))/(2*alpha**0.5)
bcoef=(alpha-beta*np.cos(t))**0.5 / (2*alpha**0.5)
V=0.5*(acoef*x**4 -bcoef*x**2)

u0=np.exp(-(x+4.72)**2 /2)*(1/np.pi)**0.25
u=u0

for i in xrange(1,time_steps):
    u=np.exp(-1j*dt*V/2)*u
    c=np.fft.fftshift(np.fft.fft(u))
    c=np.exp(-1j*dt*k**2 /2)*c
    u=np.fft.ifft(np.fft.fftshift(c))
    u=np.exp(-1j*dt*V/2)*u
    t=t+dt

plt.plot(x,u*np.conjugate(u),'r--')
plt.xlim(-10,10)
plt.ylim(0,0.6)
plt.show

然而出现的问题是在python中执行需要很长时间。我不确定我做错了什么，因为我是 python 新手，或者我是否正确实施了算法。我将不胜感激任何见解或建议。