c======================================================================= c History: newt1.f c c This program computes eigenvalues for the triangular well: c c H0 = p**2/(2*m) -V0*x c c The routine uses a Newton's method approach to solve for zeros c of the Airy function, and uses these to calculate eigenvalues c for n = 1,...,maxst. c======================================================================= program triwev implicit none c----------------------------------------------------------------------- c Argument parsing variables c----------------------------------------------------------------------- integer iargc real*8 r8arg, r8_never parameter (r8_never = -1.0d-60) real*8 v0, l, tol real*8 default_tol parameter (default_tol = 1.0d-12) c----------------------------------------------------------------------- c Filename where eigenvalues will be written c----------------------------------------------------------------------- character*(*) Efile parameter (Efile = 'trianglE') integer getu, ufrom c----------------------------------------------------------------------- c Hard-code parameters hbar and m c----------------------------------------------------------------------- real*8 hbar, m parameter (hbar = 1.0d0) parameter (m = 1.0d0) c----------------------------------------------------------------------- c maxst: maximum state for which E will be computed c E: energy eigenvalue c----------------------------------------------------------------------- integer ist, maxst parameter (maxst = 1 000) real*8 e, s, pi c----------------------------------------------------------------------- c Newton variables: c res: residual values c x: current trial value c dx: correction to trial value c iter: current iteration number c mxiter: maximum number of Newton iterations c nrm2_: 1-2 norm of quantity _ (absolute value for 1-d) c----------------------------------------------------------------------- real*8 fprime, res, x, dx integer mxiter parameter (mxiter = 10) integer iter real*8 nrm2res, nrm2dx, nrm2x real*8 drelabs c----------------------------------------------------------------------- c Airy function and derivative from imsl library c----------------------------------------------------------------------- real*8 dai, daid external dai, daid c----------------------------------------------------------------------- c Argument parsing c----------------------------------------------------------------------- if (iargc() .lt. 2) go to 900 v0 = r8arg(1, r8_never) l = r8arg(2, r8_never) if (v0 .eq. r8_never .or. l .eq. r8_never) go to 900 tol = r8arg(3, default_tol) if (tol .le. 0.0d0) tol = default_tol pi = 4.0d0*atan(1.0d0) c----------------------------------------------------------------------- c Get available unit number and open file fname c----------------------------------------------------------------------- ufrom = getu() open(ufrom,file=Efile,status='new') c----------------------------------------------------------------------- c Loop over states 1 to maxst. Each pass calculates the energy c eigenvalue for the state n = ist c----------------------------------------------------------------------- do ist = 1, maxst c----------------------------------------------------------------------- c Initialize value for zero of airy function. We use a power c series approximation to the zeros (see Handbook of Math. Func.) c----------------------------------------------------------------------- s = 3.0d0*pi*(4.0d0*dble(ist)-1.0d0)/8.0d0 x = -s**(2.0d0/3.0d0)*(1.0d0 + (5.0d0/(4.8d1*s**2)) - & (5.0d0/(3.6d1*s**4)) + (7.7125d4/(8.2944d4*s**6)) & - (1.08056875d8/(6.967296d6*s**8)) + & (1.62375596875d11/(3.34430208d8*s**10))) write(0,*) 'Newton iteration for state n = ',ist write(0,*) 'Iter x log10(dx) '// & ' log10(res)' c----------------------------------------------------------------------- c Begin Newton loop for state n = ist c----------------------------------------------------------------------- do iter = 1, mxiter c----------------------------------------------------------------------- c Evaluate residual and correction c----------------------------------------------------------------------- res = dai(x) nrm2res = abs(res) fprime = daid(x) dx = res/fprime nrm2dx = abs(dx) nrm2x = abs(x) c----------------------------------------------------------------------- c Update solution value c----------------------------------------------------------------------- x = x - dx c----------------------------------------------------------------------- c Tracing of iterations c----------------------------------------------------------------------- write(0,1000) iter, x, log10(max(nrm2dx,1.0d-60)), & log10(max(nrm2res,1.0d-60)) 1000 format(i2, 1p, e24.16, 2f8.2) c----------------------------------------------------------------------- c Check for convergence c----------------------------------------------------------------------- if (drelabs(nrm2dx,nrm2x,1.0d-6) .lt. tol) go to 100 end do c----------------------------------------------------------------------- c No convergence exit c----------------------------------------------------------------------- write(0,*) 'Computation terminated at n = ',ist write(0,*) 'No convergence after ',mxiter,' iterations.' write(0,*) tol close(ufrom) stop c----------------------------------------------------------------------- c Convergence exit c----------------------------------------------------------------------- 100 continue c----------------------------------------------------------------------- c Use solution of Ai(x)=0 to calculate energy eigenvalue c----------------------------------------------------------------------- e = -((v0**2*hbar**2/(2*m))**(1.0d0/3.0d0))*x - v0*l write(ufrom,*) ist, e end do c----------------------------------------------------------------------- c Close file fname c----------------------------------------------------------------------- close(ufrom) stop c----------------------------------------------------------------------- c Usage message c----------------------------------------------------------------------- 900 continue write(0,*) 'usage: triwev []' write(0,*) write(0,*) 'writes energies to ',Efile stop end c======================================================================= c drelabs: function for 'relativizing' quantity being c monitored for detection of convergence (by Matt Choptuik) c======================================================================= real*8 function drelabs(dx,x,xfloor) implicit none real*8 dx, x, xfloor if (abs(x) .lt. abs(xfloor)) then drelabs = abs(dx) else drelabs = abs(dx/x) end if return end