FORTRAN+计算物理学学习日记（11）

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

绘制图像如下：