1. Not finding help here? Sign up for a free 30min tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Fortran Help with Complex Variables - Values Blowing Up

  1. May 29, 2012 #1
    1. The problem statement, all variables and given/known data

    The Fortran 90 code shown below was written to solve a Crank-Nicholson algorithm that describes the motion of a quantum particle (hence the need for complex numbers). Probably not that crucial to know the details of the math/physics, just the code...

    The problem is that the numbers blow up and become infinite, and I can't figure out why.

    The fact that it uses a 'black box' solver for a tridiagonal matrix system from LAPACK doesn't make it any easier...

    2. Relevant equations

    This is the input data file that contains the values of the input variables:

    -1.0E01
    1.0E00
    1.0E00
    -2.0E01
    2.0E01
    399
    1.0E-02
    2000
    10
    as10output.dat

    3. The attempt at a solution

    ! Assignment 10 - Crank-Nicholson soln of Wavefunction
    !
    ! Takes input of matrices, boundary conditions and steps
    ! and solves a complex wavefunction using CGTSV and the
    ! Crank-Nicholson numerical method
    !
    ! Some code inspired by example code by M.Bromley, 2012
    !
    !23456
    PROGRAM Cranky
    IMPLICIT NONE
    !
    INTERFACE
    SUBROUTINE CGTSV(N, NRHS, DL, D, DU, B, LDB, INFO)
    INTEGER, INTENT(IN) :: N, LDB, NRHS
    COMPLEX, INTENT(INOUT) :: DL(N-1), D(N), DU(N-1), B(LDB,NRHS)
    INTEGER, INTENT(OUT) :: INFO
    END SUBROUTINE
    END INTERFACE
    !
    COMPLEX, ALLOCATABLE, DIMENSION(:) :: psivec,DL,D,DU,B
    REAL, ALLOCATABLE, DIMENSION(:) :: space
    !
    REAL :: &
    tfin & ! The final time
    ,x0 & ! Position at t=0
    ,xfin & ! Position at tfin
    ,alpha & ! Constant containing variables for tridiagonal
    ,sigmax & ! Initial width of Gaussian
    ,alphax & ! Initial position of Gaussian
    ,deltat & ! Width of time slices
    ,deltax & ! Width of space slices
    ,t & ! Current value of time
    ,k & ! Value of momentum kick
    ,gausstemp & ! Initial shape of the Gaussian
    ,psixave & ! Average value for psi in x
    ,psix2ave & ! Average value for psi^2 in x
    ,psinorm & ! Norm for the working estimate
    ,normfactor & ! For normalising Gaussian
    ,sigmatime & ! Average for time spread
    ,sigmaave ! Average for sigma
    !
    REAL, PARAMETER :: pi=2.0E00*ACOS(0.0E00) ! Accurate calc of pi
    !
    INTEGER :: &
    Nt & ! Number of time steps
    ,LDB & ! Some kinda LAPACK variable
    ,NRHS & ! Ditto
    ,Nx & ! Number of position slices
    ,INFO & ! CGTSV check variable
    ,tstep & ! Counting variable for t
    ,interv & ! Variable for periodic outputs
    ,jx ! Counting variable for x
    !
    CHARACTER (LEN=20) :: outfile
    !
    ! Set up file formats for output and input
    !
    100 FORMAT (7x,'tstep',4x,'time',10x,'psinorm',14x,'<x>',10x,'sigmaave',8x,'sigmatime')
    101 FORMAT (2x,i8,2x,f12.8,6(2x,g15.8))
    102 FORMAT (ES12.3/ES12.3/ES12.3/ES12.3/ES12.3/ &
    i8.3/ES12.3/i8.3/i8.3/A20)
    !
    OPEN (UNIT=20,FILE="as10input.dat",ACTION="READ") !Open inputfile
    READ (20,102) alphax,sigmax,k,x0,xfin,Nx,deltat,Nt,interv,outfile
    !
    OPEN (UNIT=21,FILE=outfile,ACTION='WRITE') ! Open output file
    !
    ! Allocate matrices
    !
    ALLOCATE (DL(1:Nx-1),D(1:Nx),DU(1:Nx-1))
    ALLOCATE (space(1:Nx),psivec(1:Nx),B(1:Nx))
    !
    ! Initialise things
    !
    deltax=(xfin-x0)/REAL(Nx+1) ! Calculate width of space slices
    alpha=(-2.5E-1*deltat)/(deltax**2) ! Calculate constant term
    normfactor=sqrt(1.0E00/(sqrt(pi)*sigmax)) ! Normalisation factor
    !
    ! Initialise space vector
    !
    DO jx=1,Nx
    space(jx)=x0+deltax*real(jx)
    END DO
    !
    ! Initialise psi vector
    !
    DO jx=1,Nx
    gausstemp=exp(-5.0E-1*((space(jx)-alphax)/sigmax)**2)
    psivec(jx)=cmplx(gausstemp*cos(k*space(jx)),gausstemp*sin(k*space(jx)))
    END DO
    psivec=cmplx(normfactor,0.0E00)*psivec
    !
    ! Calculate and write t=0 values
    !
    tstep=0
    psinorm=0.0E00
    psixave=0.0E00
    psix2ave=0.0E00
    DO jx=1,Nx
    psinorm=psinorm+cabs(psivec(jx))**2
    psixave=psixave+space(jx)*cabs(psivec(jx))**2
    psix2ave=psix2ave+space(jx)**2*cabs(psivec(jx))**2
    END DO
    psinorm=psinorm*deltax
    psixave=psixave*deltax
    psix2ave=psix2ave*deltax
    sigmaave=sqrt(psix2ave-psixave**2)
    sigmatime=sqrt(5.0E-1*(sigmax**2 + tstep**2/sigmax**2))
    !
    WRITE (21,100) ! Write column headers to file
    WRITE (21,101) 0,0.0E00,psinorm,psixave,sigmaave,sigmatime
    !
    DO tstep=1,Nt
    !
    DU = cmplx(0.0E00,alpha)
    D = cmplx(1.0E00,-2.0E00*alpha)
    DL = cmplx(0.0E00,alpha)
    LDB=Nx
    NRHS=1
    !
    ! More initialising
    !
    psinorm=0.0E00
    psixave=0.0E00
    psix2ave=0.0E00
    !
    ! Build B vector
    !
    jx=1
    B(jx)=cmplx(1.0E00,2.0E00*alpha)*psivec(jx) + cmplx(0.0E00,-alpha)*psivec(jx+1)
    DO jx=2,(Nx-1)
    B(jx)=cmplx(1.0E00,2.0E00*alpha)*psivec(jx) + cmplx(0.0E00,-alpha)*psivec(jx-1) + psivec(jx+1)
    END DO
    jx=Nx
    B(jx)=cmplx(1.0E00,2.0E00*alpha)*psivec(jx) + cmplx(0.0E00,-alpha)*psivec(jx-1)
    !
    CALL CGTSV(Nx,1,DL,D,DU,B,LDB,INFO) ! Call subroutine - solns
    !
    ! Use INFO variable to check that calculation worked:
    !
    IF (INFO .NE. 0) THEN
    WRITE (6,*) "There has been a solving error, code=", INFO
    STOP
    END IF
    !
    t=tstep*deltat
    !
    ! Copy values from B to psivec for next iteration
    !
    DO jx=1,Nx
    psivec(jx)=B(jx)
    END DO
    !
    IF (mod(tstep,interv) .EQ. 0) THEN
    !
    ! Calculate norms, averages and SDs
    !
    DO jx=1,Nx
    psinorm=psinorm+cabs(psivec(jx))**2
    psixave=psixave+space(jx)*cabs(psivec(jx))**2
    psix2ave=psix2ave + space(jx)**2*cabs(psivec(jx))**2
    END DO
    psinorm=psinorm*deltax
    psixave=psixave*deltax
    psix2ave=psix2ave*deltax
    sigmaave=sqrt(psix2ave-psixave**2)
    sigmatime=sqrt(5.0E-1*(sigmax**2 + tstep**2/sigmax**2))
    !
    ! Write results to output file
    !
    WRITE (21,101) tstep,t,psinorm,psixave,sigmaave,sigmatime
    !
    END IF
    !
    END DO
    !
    ! Deallocate all matrices
    !
    DEALLOCATE (DU,D,DL,psivec,B,space)
    !
    CLOSE(21)
    CLOSE(20)
    !
    END PROGRAM Cranky
     
  2. jcsd
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?



Similar Discussions: Fortran Help with Complex Variables - Values Blowing Up
  1. Fortran Project (Replies: 0)

  2. Do Loops in Fortran 95 (Replies: 0)

Loading...