Я решаю эволюцию нестационарного уравнения Шредингера под гармоническим потенциалом и начальной гауссовской волновой функцией.Обрабатывая hcut=1 и 2m=1 и разделяя волновую функцию в действительной и мнимой частях, получаются два связанных уравнения в терминах действительной и мнимой частей, обозначаемых как yr и yc соответственно.
xrange is [xi,xf] trange is [0,tf]
Метод, который я использовал:
сначала разделяет волновую функцию в действительной и мнимой части, а именно yr (x, t) и yc(х, т).затем обработав hcut = 2m = 1 и записав волновую функцию как yr[x,t]+i*yc[x,t], мы получим два связанных уравнения из TDSE.
1.D[yr[x,t],t]=-D[yc[x,t],x,x]+V[x]*yc[x,t]
2.D[yc[x,t],t]=D[yc[x,t],x,x]-V[x]*yc[x,t]
Затем я указал начальную волновую функцию как
yr[x,0]=exp[-x^2]
yc[x,0]=0
После этого, используя схему конечных разностей, я попытался найти y [x, t] из y [x, 0], то есть
yr[x,t+d]=yr[x,t]+d*D[yr[x,t],t]
=yr[x,t]+d*(D[yc[x,t],x,x]+V[x]*yc[x,t])
, проиндексированных как «a» в коде и «b» длясложная часть.
Значения y[x[i],t[j]] сохраняются как y[i,j], поскольку массив не может иметь действительный индекс.
Код, который я использовал, приведен ниже:
function v(x) result(s)
real::s,x
s=x**2
end function v
real::t(10000),x(10000),yc(10000,10000),yr(10000,10000),tf,xi,xf,d
integer::i,j,k,l,m,n,time
write(*,*) "time of plot divided by step size"
read(*,*) time
tf=50
xi=-10
xf=10
d=0.01
x(1)=xi
t(1)=0
i=1
do while(x(i).lt.xf) !input all values of x in x(i) array
x(i+1)=x(i)+d
i=i+1
end do
n=1
do while(t(n).lt.tf) !input all values of t in t(i) array
t(n+1)=t(n)+d
n=n+1
end do
do j=1,i
yr(j,1)=exp(-(x(j)-5)**2) !input of initial wavefunction
yc(j,1)=0
end do
!calculation of wavefunction at higher time using finite element technique[y[x,t+d]=y[x]+d*D[y[x,t],t] and then replacing the partial derivative with time
!using equation 1 and 2 .
l=1
do while(l.le.i-2)
k=1
do while(t(k).lt.tf)
yr(l,k+1)=yr(l,k)-(yc(l+2,k)-2*yc(l+1,k)+yc(l,k))/d&
+v(x(l))*yc(l,k)*d
yc(l,k+1)=yc(l,k)+(yr(l+2,k)-2*yr(l+1,k)+yr(l,k))/d&
-v(x(l))*yr(l,k)*d
k=k+1
end do
l=l+1
end do
open(1,file="q.dat")
do m=1,i-2
write(1,*) x(m),(yr(m,time))**2+(yc(m,time))**2
end do
close(1)
end
Ожидаемый результат: y (x, t) ^ 2 = yr (x, t) ^ 2 + yc (x, t) ^ 2
Ошибка: волновая функция не остается регулярной, только после t=0.05или 0.06, волновая функция становится огромной, а максимумы становятся порядка e30, несмотря на то, что гауссовская форма остается практически неизменной, как и ожидалось, поскольку прошло всего 0.05 секунд.