# Fortran Help with Complex Variables - Values Blowing Up

1. May 29, 2012

### Bravus

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=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