c======================================================= c se: Program to solve schrodinger equation: c c -u''(x) + (V(x) - E)u = 0 c c Boundary conditions????? c c c usage: se c c c output to std out is c c x_it u_it u'_it c c it = 1,2,... nout c where: c x_1 = xmin c x_nxout = xmax c nxout = 2**olevel + 1 c c======================================================= program se implicit none c include 'fcn.inc' integer iargc, i4arg real*8 r8arg c------------------------------------------------------- c Command-line arguments: c c xmin: initial integration time c xmax: final integration time c u0: initial value: u(0) c du0: initial value: u'(0) c tol: error tolerance (this program uses LSODA's c pure absolute error control) c olevel: output level: c dxout = (xmax-xmin)/2**olevel c------------------------------------------------------- real*8 xmin,xmax,u0,du0,tol integer olevel real*8 r8_never parameter ( r8_never = -1.0d-60 ) c======================================================= c start of LSODA declarations c------------------------------------------------------- c----------------------------------------------------------- c Note that 'fcn' and 'jac' are user supplied SUBROUTINES c (not functions) which evaluate the RHSs of the ODEs and c the Jacobian of the system. Under normal operation, c (as in this case), the Jacobian evaluator can be a c 'dummy' routine; if and when needed, LSODA will compute c a finite-difference approximation to the Jacobian. c----------------------------------------------------------- external fcn, jac c----------------------------------------------------------- c number of ODEs (when written in canonical first c order form). c----------------------------------------------------------- integer neq parameter ( neq = 2) c----------------------------------------------------------- c y(neq): Storage for approximate solution c x: Initial time for LSODA integration sub-interval c xout: Final time for LSODA integration sub-interval c----------------------------------------------------------- real*8 y(neq) real*8 x, xout c----------------------------------------------------------- c Tolerance parameters: c c The following comment block is extracted from the c LSODA documentation. c----------------------------------------------------------- c rtol = relative tolerance parameter (scalar). c atol = absolute tolerance parameter (scalar or array). c the estimated local error in y(i) will be controlled so c as to be less than c ewt(i) = rtol*abs(y(i)) + atol if itol = 1, or c ewt(i) = rtol*abs(y(i)) + atol(i) if itol = 2. c thus the local error test passes if, in each component, c either the absolute error is less than atol (or atol(i)), c or the relative error is less than rtol. c use rtol = 0.0 for pure absolute error control, and c use atol = 0.0 (or atol(i) = 0.0) for pure relative error c control. CAUTION.. actual (global) errors may exceed c these local tolerances, so choose them CONSERVATIVELY. c----------------------------------------------------------- real*8 rtol, atol integer itol c----------------------------------------------------------- c Control parameters and return code (see below). c----------------------------------------------------------- integer itask, istate, iopt c----------------------------------------------------------- c Work arrays. c----------------------------------------------------------- integer lrw parameter ( lrw = 22 + neq * 16 ) real*8 rwork(lrw) integer liw parameter ( liw = 20 + neq ) integer iwork(liw) c----------------------------------------------------------- c 'jt' defines which type of Jacobian is supplied or c computed; we use jt = 2 here which, as mentioned c above, instructs LSODA to compute a finite-difference c approximation to the Jacobian if and when needed. c----------------------------------------------------------- integer jt c=========================================================== c End of LSODA declarations c=========================================================== c----------------------------------------------------------- c miscellaneous variables c----------------------------------------------------------- real*8 dxout integer it, nxout c----------------------------------------------------------- c argument parsing c----------------------------------------------------------- if( iargc() .ne. 6 ) go to 900 xmin = r8arg(1, r8_never) xmax = r8arg(2, r8_never) u0 = r8arg(3, r8_never) du0 = r8arg(4, r8_never) tol = r8arg(5, r8_never) olevel = i4arg(6, -1) if (xmax .le. xmin .or. u0 .eq. r8_never .or. & du0 .eq. r8_never .or. tol .eq. r8_never & .or. tol .gt. 10d-2 .or. tol .lt. 10d-12) & go to 800 c----------------------------------------------------------- c set LSODA parameters ... see LSODA documentation c for more complete description c----------------------------------------------------------- itol = 1 ! indicates that 'atol' is scalar rtol = 0.0d0 ! use pure absolute tolerance atol = tol ! absolute tolerance itask = 1 ! normal computation iopt = 0 ! indicates no optional inputs jt = 2 ! jacobian type c----------------------------------------------------------- c compute number of outputs and output interval c and initialize sub-interval start time and solution c estimate. c----------------------------------------------------------- nxout = 2**olevel +1 dxout = (xmax - xmin) / (nxout -1) x = xmin y(1) = u0 y(2) = du0 c----------------------------------------------------------- c output initial solution c----------------------------------------------------------- write(*,*) x, y(1) c----------------------------------------------------------- c Loop over requested outputs ... c c Set istate to 1 to indicate initial call, istate c should be set to 2 for subsequent calls, but lsoda c will automatically do this so long as the initial c call is successful. c----------------------------------------------------------- istate = 1 do it = 2, nxout if (x .eq. 0.0d0) then x = x + dxout end if c set final integration position for current interval xout = x + dxout c call lsoda to integrate system on [x ... xout] call lsoda(fcn, neq, y, x, xout, & itol, rtol, atol, itask, & istate,iopt,rwork,lrw,iwork,liw,jac,jt) c----------------------------------------------------------- c Check return code and exit with error message if c there was trouble. c----------------------------------------------------------- if( istate .lt. 0 ) then write(0,1000) istate, it, nxout, x, x + dxout 1000 format(/' sode: Error return ',i2, & ' from integrator LSODA.'/ & ' sode: At output time ',i5,' of ',i5/ & ' sode: Interval ',1p,e11.3,0p, & ' .. ',1p,e11.3,0p/) go to 500 end if write(*,*) x, y(1) end do 500 continue stop stop 800 continue write(0,*) 'xmax > xmin' write(0,*) '10e-12 <= tol <= 10e-2' write(0,*) 'olevel > 0' write(0,*) 900 continue write(0,*) 'usage: se '// & ' ' stop end subroutine fcn(neq,x,y,yprime) implicit none real*8 calcv2 external calcv2 integer neq real*8 x, y(neq),yprime(neq) real*8 v integer nmax parameter ( nmax = 100 ) call calcv(v, x, nmax) yprime(1) = y(2) yprime(2) = (v - 1.0d0) * y(1) return end subroutine jac implicit none return end