3.6一维薛定谔方程的定态解
program main
implicit none
integer(8)::m,n,i,j
real(8)::xmin,gamma,tol,e,de,x(3000),y(3000),yz(3000),c,psi(3000,3000),a,b,xm,xmax,k
real(8),external::matchpoint,v1
open(unit=100,file="xable.csv")
open(unit=101,file="j=1.csv")
open(unit=102,file="j=2.csv")
open(unit=103,file="j=3.csv")
open(unit=105,file="e.csv")
xmin=-5.0
xmax=5.0
n=1000
gamma=sqrt(2.0)
tol=10e-10
e=-1.0
de=0.1
m=3
do j=1,m,1
e=e+de
xm=matchpoint(xmin,e,tol)
call equelx(xm,xmax,xmin,n,x)
call equely(xm,n,e,xmin,xmax,tol,gamma,y)
do while(abs(de)>tol)
e=e+de
xm=matchpoint(xmin,e,tol)
call equelx(xm,xmax,xmin,n,x)
call equely(xm,n,e,xmin,xmax,tol,gamma,yz)
a=(y(n+1)-y(n))/(x(n+1)-x(n))-(y(n+3)-y(n+2))/(x(n+3)-x(n+2))
b=(yz(n+1)-yz(n))/(x(n+1)-x(n))-(yz(n+3)-yz(n+2))/(x(n+3)-x(n+2))
if(a*b>0)then
else
e=e-de
de=de/2.0
end if
end do
de=0.1
!!the nomalization
c=0.0
do i=1,2*n+1,1
c=c+(x(i+1)-x(i))/2.0*(yz(i)*yz(i)+yz(i+1)*yz(i+1))
end do
do i=1,2*n+2,1
psi(i,j)=yz(i)/sqrt(c)
end do
do i=1,2*n+2,1
if(j==1)then
write(101,*)psi(i,j)
else if(j==2)then
write(102,*)psi(i,j)
else
write(103,*)psi(i,j)
end if
end do
write(105,*)e
end do
do i=1,2*n+2,1
write(100,*)x(i)
end do
close(100)
close(101)
close(102)
close(103)
close(105)
end program main
function matchpoint(xmin,e,tol)
implicit none
real(8)::xmin,e,tol,dx,fo
real(8)::matchpoint
real(8),external::v1
matchpoint=xmin
dx=0.01
fo=v1(matchpoint)
do while(abs(dx)>tol)
matchpoint=matchpoint+dx
if((fo-e)*(v1(matchpoint)-e)>0)then
else
matchpoint=matchpoint-dx
dx=dx/2.0
end if
end do
return
end
function v1(x)
implicit none
real(8)::x
real(8)::v1
v1=-(1.0-1.0/2.0*x*x)
return
end
subroutine left(delta,xm,n,e,xmin,tol,gamma,y)
implicit none
real(8),intent(in)::delta,xm,e,xmin,tol,gamma
integer(8),intent(in)::n
real(8),intent(out)::y(3000)
real(8)::h,x(3000)
integer(8)::i
real(8),external::v1
h=(xm-xmin)/n
y(1)=0.0
y(2)=delta
do i=1,n+1,1
x(i)=xmin+h*(i-1)
end do
do i=2,n,1
y(i+1)=2.0*(1.0-5.0/12.0*h*h*gamma*gamma*(e-v1(x(i))))*y(i)-(1.0+h*h/12.0*gamma*gamma*(e-v1(x(i-1))))*y(i-1)
y(i+1)=y(i+1)/(1+h*h/12.0*gamma*gamma*(e-v1(x(i+1))))
end do
end subroutine left
subroutine right(xm,n,e,xmin,xmax,tol,gamma,y)
implicit none
real(8),intent(in)::xm,e,xmin,tol,gamma,xmax
integer(8),intent(in)::n
real(8),intent(out)::y(3000)
real(8)::h,x(3000)
integer(8)::i
real(8),external::v1
h=(xmax-xm)/n
y(n+1)=0.0
y(n)=0.1
do i=1,n+1,1
x(i)=xm+h*(i-1)
end do
do i=n,2,-1
y(i-1)=2.0*(1.0-5.0/12.0*h*h*gamma*gamma*(e-v1(x(i))))*y(i)-(1.0+h*h/12.0*gamma*gamma*(e-v1(x(i+1))))*y(i+1)
y(i-1)=y(i-1)/(1.0+h*h/12.0*gamma*gamma*(e-v1(x(i-1))))
end do
end subroutine right
subroutine equely(xm,n,e,xmin,xmax,tol,gamma,y)
implicit none
real(8),intent(in)::xm,e,xmin,xmax,tol,gamma
integer(8),intent(in)::n
real(8),intent(out)::y(3000)
real(8)::delta,ddelta,y1(3000),y2(3000),yy1(3000)
integer(8)::i
delta=0.4
ddelta=0.1
call left(delta,xm,n,e,xmin,tol,gamma,y1)
call right(xm,n,e,xmin,xmax,tol,gamma,y2)
do while(abs(ddelta)>tol)
delta=delta-ddelta
call left(delta,xm,n,e,xmin,tol,gamma,yy1)
if((yy1(n+1)-y2(1))*(y1(n+1)-y2(1))>0)then
else
delta=delta+ddelta
ddelta=ddelta/2.0
end if
end do
do i=1,2*n+2,1
if(i<(n+2))then
y(i)=yy1(i)
else
y(i)=y2(i-1-n)
end if
end do
end subroutine equely
subroutine equelx(xm,xmax,xmin,n,x)
implicit none
real(8),intent(in)::xm,xmax,xmin
real(8),intent(out)::x(3000)
integer(8),intent(in)::n
integer(8)::i
real(8)::p1,p2
p1=(xm-xmin)/n
p2=(xmax-xm)/n
do i=1,2*n+1,1
if(i<(n+2))then
x(i)=xmin+p1*(i-1)
else
x(i)=xm+p2*(i-n-2)
end if
end do
end subroutine equelx
绘制图像如下: