- #1

Bravus

- 23

- 0

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