Fortran Help with Complex Variables - Values Blowing Up

[Bravus]

Homework Statement

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

Homework 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

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

 

[KataKoniK]

it is important to approach any problem with a systematic and analytical mindset. In this case, the first step would be to understand the problem at hand and the code provided. It is important to have a good understanding of the underlying physics and mathematics involved in order to effectively troubleshoot the issue.

The first thing to check would be the input data and make sure it is correct. Are the values reasonable? Are there any typos or mistakes? It would also be helpful to have some knowledge of the expected output in order to compare it to the actual output.

Next, it would be important to understand the code and how it is solving the problem. It is mentioned that a black box solver for a tridiagonal matrix system from LAPACK is being used. It would be helpful to understand how this solver works and what kind of results it should produce.

One possible reason for the numbers blowing up could be due to numerical instability. This can happen when the algorithm is not well-conditioned or when there are errors in the computations. It would be helpful to check for any potential sources of error in the code, such as division by zero or round-off errors.

Another possible reason could be an issue with the initial conditions or boundary conditions. It would be helpful to check these and make sure they are correctly implemented in the code.

It may also be useful to try running the code with different input values to see if the issue persists. This could help narrow down the potential sources of error.

If the problem still cannot be resolved, it may be helpful to seek assistance from other experts in the field or consult additional resources for guidance. It is important to approach the problem with a methodical and patient mindset, as solving complex scientific problems can often require trial and error and multiple iterations.